© The Authors 2024

1 Introduction

Gas-phase molecular emission and absorption is ubiquitous in the spectra of the disks and outflows associated with young pro-tostars. Analysis of their spectra provides a valuable probe into the physical conditions where molecules are formed, and the mechanism for their excitation. Mid-infrared (IR) spectroscopy, in particular, offers a unique window for observing rovibrational molecular features (e.g. van Dishoeck 2004; van Dishoeck et al. 2023; Hollenbach & McKee 1989). Key star formation tracers such as C₂H₂, CO₂, and CH₄ are among the most abundant carbon and oxygen carrying species, but they do not have observable rotational spectra due to their lack of a permanent electric dipole moment, and thus they can only be observed in the mid-IR through their rovibrational transitions. Furthermore, observations from space telescopes are essential for observing molecules such as CO₂ and H₂O, whose presence in our atmosphere makes their detection from the ground difficult.

With the launch of the James Webb Space Telescope (JWST), space-based mid-IR spectroscopy is once again possible with the Mid-Infrared Instrument (MIRI; Rieke et al. 2015; Wright et al. 2015, 2023) Medium Resolution Spectrometer (MRS; Wells et al. 2015). JWST provides several orders of magnitude in improvements to sensitivity and spatial resolution over previous observatories such as Spitzer and the Infrared Space Observatory (ISO). The improved spectral resolution of MIRI/MRS of R = Δλ/λ ~ 3400–1600 (Argyriou et al. 2023) versus 50–600 for Spitzer significantly boosts the line to continuum ratio for gas-phase lines, although MIRI/MRS does not reach as high a resolution as the R ~ 30 000 of ISO-SWS. The unmatched sensitivity of JWST also enables spectroscopy of far fainter and more distant targets than possible before. This is particularly useful for the study of massive protostars, which are rarer than their low-mass siblings and have short evolutionary timescales of a few 100 000 yr (e.g. Mottram et al. 2011; Motte et al. 2018). Furthermore, the sub-arcsecond spatial resolution of JWST is valuable for densely clustered star-forming regions where high-mass stars preferentially form.

Past mid-IR studies of high-mass protostars have largely focussed on those in the more evolved ‘hot core’ phase, where significant heating of the envelope by a massive protostar produces diverse complex organic molecules (COMs: carbon-bearing molecules with six or more atoms) and the formation of compact HII regions (Gerner et al. 2014; Cesaroni et al. 1997). Prior to the hot core phase, the gas-phase COM abundance is lower. Towards these protostars, transitions of molecules such as CO₂, C₂H₂, HCN, and H₂O are typically seen in absorption (e.g. Evans et al. 1991; Helmich et al. 1996; Boonman et al. 2003a,c; Boonman & van Dishoeck 2003). High-resolution spectroscopy at up to R = 100 000 with the Stratospheric Observatory for Infrared Astronomy Echelon-Cross-Echelle Spectrograph (SOFIA/EXES) and Very Large Telescope Cryogenic High-resolution Infrared Echelle Spectrograph (VLT/CRIRES) suggests that the absorption features often arise from the surface layers against the accretion-heated midplane of a massive protostellar disk (Indriolo et al. 2020; Barr et al. 2020, 2022). However, mid-IR molecular emission has also occasionally been detected in high-mass protostars, where it is associated with warm gas in the outflow excited by mid-IR pumping (e.g. the outflows of Cepheus A, Sonnentrucker et al. 2006, 2007 and Orion-IRc2/KL, Boonman et al. 2003c).

Mid-IR molecular emission from embedded disks of both low- and high-mass protostars has rarely been detected, likely due to extinction from the envelope and lower sensitivities of past observatories (Lahuis et al. 2010). Recent MIRI/MRS observations of class 0 protostars have identified gas-phase features of CO and H₂O possibly associated with the disk (Yang et al. 2022), and a stronger line forest from a variety of species in the eruptive protostar EX Lup (Kóspál et al. 2023). In the more evolved class II stage, molecular emission is much more frequently detected and better characterized, and previous ground-based and Spitzer surveys have shown it to originate from the hot inner disk and warm disk surface (e.g. Blake & Boogert 2004; Carr & Najita 2008, 2011; Salyk et al. 2011; Pontoppidan et al. 2010, 2014; Banzatti et al. 2022, 2023). JWST observations of class II disks have confirmed this, and also suggested links between the inner disk molecular inventory and dust transport within the disk (Grant et al. 2023; Tabone et al. 2023; Gasman et al. 2023; Perotti et al. 2023; Banzatti et al. 2023).

In this paper, we present an analysis of gas phase molecular features in JWST/MIRI observations of the high mass star-forming region IRAS 23385+6053. This paper is part of a series on the first results from the JWST Observations of Young Protostars (JOYS) collaboration: an overview of the IRAS 23385+6053 MIRI/MRS observations and the detection of the accretion-tracing Humpreys α line has been presented by Beuther et al. (2023), an analysis of the outflows and warm gas traced in the mid-IR and sub-mm by Gieser et al. (2023), and the ice absorption features by Rocha et al. (2024). Some of the data presented in this paper have also been discussed by van Dishoeck et al. (2023) in the context of the astrochemistry of planet-forming disks.

IRAS 23385+6053 is located in the outer Galaxy at a distance of 4.9 kpc (Molinari et al. 1998) and a Galactocentric distance of 11 kpc, with a luminosity of ~3000 L_⊙ and total envelope mass of 510 M_⊙ (Cesaroni et al. 2019). Previous observations suggest that IRAS 23385+6053 is extremely young, with a low bolometric temperature of ~40 K (Fontani et al. 2004), and no observed HII regions or free-free emission. Furthermore, IRAS 23385+6053 is not as rich in COMs as typical hot core regions (Cesaroni et al. 2019), consistent with a young age. The star formation in IRAS 23385+6053 is highly clustered, with up to 6 different cores detected in the sub-mm (Cesaroni et al. 2019) and a complex system of at least three high velocity outflows (Beuther et al. 2023). The detection of large-scale gas motions suggests a large disk may be present around a ~9 M_⊙ star within one of the brightest cores (Cesaroni et al. 2019).

The remainder of this paper is organized as follows: in Sect. 2, we discuss the data reduction, extraction of spectra in apertures containing the various protostars and outflow knots, and the LTE slab models used to identify and fit the observed molecular features. In Sect. 3, we then list the identified species and the results of the LTE model fits. We continue with a discussion of the origin of each detected species and the components of the protostellar system they most likely trace in Sect. 4. Finally, we summarize our overall model for the molecular features in IRAS 23385+6053, and conclude with some remarks on the future outlook for JWST studies of high mass star formation in Sect. 5.

2 Observations and methods

2.1 Data reduction

The full details of the data reduction are described in Beuther et al. (2023). In summary, we have run the JWST calibration pipeline (Bushouse et al. 2022) using the reference context jwst_l0l7.pmap of the JWST Calibration Reference Data System (CRDS; Greenfield & Miller 2016). Additionally, we have applied an astrometric correction using Gaia DR3 stars in the simultaneous MIRI imaging field of view. A dedicated background field was taken off-source. However, due to significant astronomical emission in the dedicated background observation, subtracting this background from the science data resulted in negative fluxes. The telescope background was therefore estimated by extracting a spectrum from within the primary science field of view, but off-source from the main infrared continuum sources at the position within the IFU where the background flux was the lowest (RA (J2000) 23^h40^m54.15^s, Dec (J2000) 61^d 10^m26.96^s). The aperture size was the same as was used for extracting the science data (Sect. 2.2). The background subtraction also resulted in the subtraction of the strong PAH emission features around 8.6 µm and 11.3 µm since their emission was about equally strong in the offset position, although some residuals remain. However, the molecular emission discussed in this paper does not overlap with any of the PAH bands and is therefore not affected by them.

2.2 Source overview and spectra extraction

An overview of our MIRI/MRS observations of IRAS 23385+6053 has been presented by Beuther et al. (2023), which we briefly summarize here. In Fig. 1, the emission from the H₂ S(7) line in channel 1 is shown in the first panel, while the continuum from the second and third MIRI/MRS channels is shown by the colour scale. For reference, we indicate the putative directions of the major outflows in the left panel. At mid-IR wavelengths, two mid-IR sources with a separation of ~0.67" / 3280 au are detected; it is not known if these sources are gravitationally bound. Following the naming of Beuther et al. (2023), we refer to the north-west and south-east sources as A and B respectively. source A and B are clearly resolved in channels 1 and 2 (4.9–11.7 µm) but begin to become blended in channel 3 (11.55–17.98 µm) and are completely blended in channel 4 (17.7–27.9 µm). Spatially extended H₂ emission tracing the outflows is detected in the S(1) to S(8) transitions over the majority of the MIRI field of view (Beuther et al. 2023). On the basis of their large projected separation and the presence of ice and silicate absorption features towards both sources, Beuther et al. (2023) argue that source A and B are two distinct proto-stars, which are likely responsible for driving at least two of the observed outflows in this region.

IRAS 23385+6053 has also been observed at mm wavelengths as part of the CORE programme (Beuther et al. 2018; Cesaroni et al. 2019; Gieser et al. 2021), where extended continuum emission likely tracing the positions of several deeply embedded protostars is detected (grey contours, middle panel of Fig. 1). source A is seen to be closely aligned with an existing mm peak, but source B is clearly offset, although still located within the cold envelope. Furthermore, although mid-IR continuum emission is only detected from source A and B, up to 6 cores have been identified from sub-mm molecular emission by Cesaroni et al. (2019). The lack of significant mid-IR emission towards the other cores may reflect a difference in mass or evolutionary state of protostars in the cluster, or simply a lack of sensitivity, as suggested by Beuther et al. (2023).

As mid-IR molecular features have previously been detected at both the locations of protostars and in the outflows (Beuther et al. 2023; Gieser et al. 2023), we extract spectra towards sources A and B, as well as the three brightest H₂ knots (see Fig. 1), which are named S1, S2, and S3 following Gieser et al. (2023). We note that the brightest knot (S1) is close to a bright mm peak possibly tracing a prestellar core or embedded pro-tostar, and that the H₂ emission here is seen over an extensive velocity range of ± 140 km s⁻¹, suggesting strong outflow activity (Beuther et al. 2023). Our extraction apertures for each source are shown by the circles in Fig. 1. Since the spatial resolution of MIRI/MRS varies significantly from ~0.27" in channel 1 to ~1.0" in channel 4 (Law et al. 2023), we sum the flux in a ‘conical’ aperture whose radius scales with wavelength as r = 1.22(λ/6.5m). Because of the blending of sources A and B in channels 3 and 4, we also extract a spectrum in a large aperture containing both sources with a fixed diameter of 2.5", identical to that used by Beuther et al. (2023). To all of our spectra, we apply a residual fringe correction procedure to remove fringing behaviour not corrected for by the pipeline (Kavanagh et al., in prep.).

Fig. 1 First three channels of MIRI/MRS observations of IRAS 23385+6053. The left panel shows the H2 S(7) transition in channel 1 at 5.511 µm, while the centre and right panels show the averaged continuum in channels 2 and 3. The colour scale in all panels is normalized to the peak intensity. The arrows in the left panel indicate the directions and velocity shift of the major outflows proposed by Beuther et al. (2023). The extraction apertures for our spectra are shown in each panel by the coloured circles. Aperture A+B has a constant radius of 2.5", while all other apertures have aradius that scales with wavelength as r = 2(1.22λ/6.5 m). The grey contours overlaid on channel 2 are 1.3 mm continuum emission (Gieser et al. 2021). The white circle in the bottom left indicates the mean beam size of JWST in the given channel, while the grey ellipse in the second panel indicates the FWHM of the mm beam.

2.3 Baseline subtraction and ice absorption correction

To characterize the molecular features, we first perform a baseline subtraction using univariate spline fits to line-free regions, in order to remove the thermal continuum, PAH emission, and ice absorption. We also model the thermal continuum using a second order polynomial fit to areas assumed to be free of ice-absorption. Two examples of local baselines and thermal continuum fits to the ~6.5 µm and ~15 µm regions are shown in Fig. A.1. More global continuum fits to study ices are presented in Rocha et al. (2024).

2.4 LTE slab models

We fitted the continuum-subtracted molecular emission using the local thermodynamic equilibrium (LTE) slab models described in Grant et al. (2023), Tabone et al. (2023), which are similar in functionality to the existing slabspec code (Salyk 2020). The molecular emission is assumed to originate in a plane-parallel uniform cylindrical slab of radius R and column density N, with a gas at a single temperature T_ex. In reality, proto-stellar environments should contain gas over a wide range of densities and temperatures in the disk, outflow, and envelope which such simple models cannot capture, and non-LTE effects may be particularly important for regions with low density and strong radiative pumping (e.g. outflow shocks). In particular, the effect of mid-IR pumping in particular can result in model excitation temperatures and column densities which do not reflect the kinetic temperature and column density of that gas, which we discuss further in Sect. 4.1. Nevertheless, LTE models are still valuable for identifying and characterizing the molecular emission; we discuss caveats of their use inSect. 4.

The molecular spectroscopy data for the models (transition wavelengths, Einstein A coefficients, statistical weights, and partition functions) are taken from the HITRAN database (Gordon et al. 2022). This information is used in combination with an assumed T_ex and N to solve for the level populations. The opacity at line centre is determined assuming a turbulent broadening by a Gaussian with full width at half maximum of ΔV = 4.7 km s⁻¹, which has typically been assumed for LTE models of molecular emission in T Tauri disks (e.g. Salyk et al. 2011). The typical mm line widths in the CORE survey are between 2 and 5 km s⁻¹, so this assumption is reasonable for IRAS 23385+6035 (see Fig. 7 of Cesaroni et al. 2019). The wavelength dependent line optical depth is computed on a high resolution (λ/Δλ) = 10⁶ grid assuming a Gaussian line profile, as there may be significant overlap between closely spaced transitions from molecules with a Q-branch (Tabone et al. 2023). The model flux F(λ) is then computed as $$F(\lambda )=\pi {\left(\frac{R}{d}\right)}^2{B}_v(T)\left(1-{\text{e}}^{-\tau (\lambda )}\right),$$(1)

where R is the radius of the emitting slab, d = 4.9 kpc is the distance to the source, B_ν is the Planck function, and τ(λ) is the line optical depth. The model is then convolved to the average MIRI instrumental resolution in each sub-band (Labiano et al. 2021; Jones et al. 2023). For comparisons with the observations, we velocity shift the convolved synthetic spectra by v_LSR = −50.2 km s⁻¹ for sources A and B (Beuther et al. 2018) and a by-eye-estimated −100 km s⁻¹ for source S1, then finally resample at the observed wavelengths of the MIRI spectra.

Correction for extinction from the local molecular cloud and protostellar envelope is important to obtain accurate column densities from our molecular emission models. We follow the approach of Gieser et al. (2023), who corrected the IRAS 23385+6053 spectra using the extinction curves of McClure (2009) with a K-band extinction of A_K = 7.0. However, we make one important modification to more accurately for the effect of ice absorption. The extinction curves of McClure (2009) include features from cloud absorption of H₂O and CO₂ ice, which may vary significantly from cloud to cloud in shape and depth, and extincts the molecular emission with a strong wavelength dependence. We thus remove the ice features from these extinction curves using a local polynomial fit to ice-free regions. To then include the effect of extinction by ices unique to IRAS 23385+6053 (located in either the protostellar envelope or surrounding molecular cloud), we calculate the optical depth of the ice feature τ = exp (−F_cont/F_baseline) as a function of wavelength, and apply it as an additional correction factor to our molecular emission models. An example of this correction is shown for CO₂ in Fig. A.2, where the ice absorption is seen to significantly weaken the Q and P-branch emission at λ > 14.9 µm. Similarly, from 5.5 to 7.8 µm, there is also strong absorption from H₂O, ${\text{NH}}_4^{+}$, and CH₄ ice (lower panel of Fig. A.1) which overlaps with the H₂O and CH₄ emission features, and thus an analogous correction is applied to the models. An important caveat of this approach is the implicit assumption that all of the absorbing ice lies in front of all of the emitting gas. This assumption may not be completely correct if a large fraction of the emitting gas is in an outflow outside the protostellar envelope, and may be further complicated by the presence of multiple sources experiencing varying degrees of extinction within a single aperture in the case of IRAS 23385+6053.

Molecular absorption features are detected in one case (H₂O in source B), and modelled in a similar way using LTE slab models based on the work of Helmich (1996). We follow the same procedure outlined above, but only calculating the line optical depth as a function of wavelength, followed immediately by convolution to the MIRI resolution, velocity shifting, and re-sampling. The models are compared to the observed optical depth, which is calculated as τ_obs = − log(F_obs/F_baseline). Unlike the case of emission, corrections for extinction by the cloud and envelope ices are not applied, as these do not influence the observed optical depth of the transitions.

2.5 Model-fitting procedure

To identify best-fit emission models, for each molecule detected in source A, B, and S1 (see Table 1), we compute a grid of models with varying N and T_ex, the details of which are provided in Appendix C. For each model, we calculate $${\chi }^2={\displaystyle \sum _{i=1}^{N_{\text{obs}}}{\left(\left(F_{\text{obs},i}-F_{\text{model },i}\right)/\sigma \right)}^2},$$(2)

where F_model,i· is determined from Eq. (1), σ is estimated from the emission in nearby line-free regions, and a minimization routine is used to determine a value of R that minimizes χ². The value of σ determined empirically from line-free spectral regions is 0.1 mJy. We only fit regions of the spectrum containing obvious molecular features, thus avoiding contaminating emission from bright atomic lines or other molecules. In some cases, emission from several molecules is blended together (e.g. HCN, OH, C₂H₂ and CO₂, see Grant et al. 2023). Here, we iteratively fit each molecule in order of decreasing flux, using the residuals of the previous best-fit as the input to the following one.

A similar approach is used for identifying best-fit absorption models, but in units of model and observed optical depth, with χ² calculated as $${\chi }^2={\displaystyle \sum _{i=1}^{N_{\text{obs}}}{\left(\left({\tau }_{\text{obs},i}-{\tau }_{\text{model },i}\right)/{\tau }_{\sigma }\right)}^2},$$(3)

where τ_σ = σ/F_obs,i·. We note that in the case of absorption, there is no emitting area radius R to constrain for the ‘pencil-beam’ of the absorbing column.

Uncertainties on the best-fit N and T_ex are estimated following Carr & Najita (2011) and Salyk et al. (2011) by producing χ² maps for the best-fit column density and temperature with contours of the best R value overlaid. Confidence intervals are computed assuming K = 2 degrees of freedom for the column density and temperature following the approach in Avni (1976). Specifically, for the map of ${\chi }_{\text{red}}^2={\chi }^2/K$, we overlay contours relative to the best-fit minimum at levels $\Delta {\chi }_{\text{red }}^2={\chi }_{\text{red }}^2-{\chi }_{\text{red,min }}^2$. The contour levels corresponding to 1,2, and 3σ are thus $\text{Δ}{\chi }_{\text{red }}^2\approx$ 2.3, 6.2, and 11.8.

Fig. 2 Spectra of sources A+B (top panel) and the brightest H2 knots (lower panel) extracted from the apertures shown in Fig. 1. Flux offsets of −2.5 and −5.0 mJy are applied to the spectra of S2 and S3 respectively. In contrast with the A+B spectrum, no molecular lines other than H2, CO, and OH in S1 are detected towards the H2 knots.

3 Results

3.1 Detections

An overview of the spectra in the first three MIRI channels is shown in Fig. 2 for the combined A+B aperture (upper panel) and the H₂ knots (lower panel). Mid-IR continuum emission is detected towards A + B at a level of up to ~ 10 mJy, orders of magnitude lower than the brighter and more evolved high mass sources studied with ISO (e.g. Boonman et al. 2003a), which cannot be observed by JWST due to rapid saturation of the detectors. No continuum emission is seen at the shock positions S1 to S3. Overlaid on the continuum from A + B are strong absorption features from various silicates and ices (Rocha et al. 2024). The aforementioned H₂ S(1) to S(8) transitions as well several shock-tracing ionized atomic lines are detected towards all sources (see Beuther et al. 2023).

A variety of molecular species are detected towards the different sources, which we summarize in Table 1. As previously noted, H₂ is detected over the entire MIRI/MRS field of view, and is clearly associated with the outflows (see Gieser et al. 2023). The (0–0) R(6), R(5), and R(4) transitions of HD are detected, which in combination with the H₂ lines, will be used as a probe of the D/H abundance in a future publication (Francis et al., in prep.). Strong molecular emission from CO₂, C₂H₂, and HCN is detected towards source A + B (inset, top panel of Figs. 2, 3), as well as fainter ~6 µm and ~16 µm H₂O features in the individual sources (Fig. 4). Towards the H₂ knots no such features are present (insets, lower panel of Fig. 2), however, emission from CO and OH is detected in source S1 (see also Fig. 5).

The spectra towards source A and B are particularly rich in comparison to the H₂ knot positions. In Fig. 3, we show a zoom-in on the A+B spectra in the 13.6–15.6µm range where emission from the CO₂, C₂H₂, and HCN is detected. The CO₂ emission is co-located with a strong absorption feature of CO₂ ice which extincts the Q- and P-branches of the emission (Sect. 2.4). In both sources, weak emission from OH is detected from 14 to 17 µm, but most transitions except those at ~ 16.04 and ~ 16.84 µm are blended with CO₂ and HCN emission. The iso-topologues ¹³CO₂ and ¹³CCH₂ are also detected, but are blended with emission from the P-branch of CO₂ and the Q-branch of C₂H₂.

Several differences are apparent between the spectra of source A and B as well. In channels 1 and 2 (4.9–11.7 µm), emission from sources A and B is spatially resolved, but partial blending and overlap of our extraction apertures begins to occur in channel 3 (11.55–17.98 µm). We therefore analyse source A and B independently in channels 1 and 2 only. In source A, rovi-brational emission from the CO P-branch is detected from 4.9 to 5.2 µm with υ = 1 − 0, J → J + 1 levels up to J = 45 (E_up = 9481 K), indicating extremely energetic excitation conditions. CO emission is entirely absent in source B, however (see Fig. 4).

H₂O is detected in both sources from 4.9 to 8 µm in its υ₂(1 – 0) vibrational bending mode, but in emission in source A and in absorption in source B. In source A, a handful of ‘cold’ transitions of pure rotational water emission at ~16– 18 µm are also detected, though blending with source B makes this association ambiguous.

In the H₂ hotspots, fewer molecules are detected. Emission from CO and OH is only seen towards the brightest source S1 (see Fig. 5), likely tracing the presence of strong, dense shocks.

Fig. 3 Zoom in on the molecular emission from CO2, C2H2, HCN, and OH emission seen in the A+B aperture (see Fig. 1). No continuum subtraction has been applied.

3.2 LTE slab model fits

As the emission from CO₂ C₂H₂, HCN, and OH lies in channel 3 where blending between sources A and B is significant, we fitted these species using spectra from the combined A+B aperture.

For all species in channels 1 and 2, we fitted the spectra of source A and B independently. We show an example of the combined best-fit model to source A+B for each molecule in Fig. 6; the fits to the remaining molecules and sources are shown in Figs. B.1–B.3.

For each model, we obtain a best-fit N, T_ex, and R, and also calculate the total number of molecules $\mathcal{N}=\pi *R_{\text{best }}^2*N_{\text{best }}$ and a column density ratio with respect to H₂. The column densities of H₂ for each source have been obtained with a rotational diagram analysis by Gieser et al. (2023). Two temperature components are used in the rotational diagram fitting, and we use the column density of the cooler T_Tot ~ 300 K H₂ for our abundance calculations (see Appendix of Gieser et al. 2023). We note that this is warmer than the best-fit temperatures for CO₂, C₂H₂, and HCN, and cooler than the temperatures for the H₂O absorption, and thus the warm H₂ may not trace the same reservoir of gas. The H₂ emission from IRAS 23385+6053 clearly traces the outflows as well as the protostar positions (see Fig. 1 left panel), and thus our inferred abundances may be higher in reality if only a fraction of the H₂ emission originates from the quiescent disk and envelope. However, Gieser et al. (2023) also determine H₂ column densities using the mm continuum emission and an assumed gas-to-dust ratio, and find column densities consistent with the warm H₂ within a factor of a few. The best-fit parameters of our models and estimated column density ratios are summarized in Table 2.

Confidence interval contours in the N-T_ex plane for each fit are overlaid on maps of normalized χ² in Figs. C.1–C.3. The confidence contours can be complex but typically show one of the following behaviours. 1) The emission is optically thick and the confidence contours are banana-shaped, indicating that R is well-constrained, while N and T_ex have a non-linear degeneracy. 2) The emission is optically thin and the confidence contours are tongue-shaped, indicating T_ex is well-constrained, but R and N are fully degenerate, and only lower and upper limits respectively can be determined. However, in this case the total number of molecules N is still well-constrained. 3) The S/N of the molecular features is low, resulting in confidence intervals showing a combination of behaviours 1 and 2. Here, solutions with either cooler and more optically thick emission or hotter and more optically thin emission are degenerate.

Fig. 4 As Fig. 3, but for the CO emission and 4.9–6.7 µm H2O absorption in the spectra of source A (green) and B (blue).

3.2.1 Source A and B

The high S/N of the CO₂ and C₂H₂ emission allows the parameters of the LTE model fits to be particularly well constrained. For the CO₂ emission we fitted only the Q and R-branch. The width of the Q-branch and distribution of the R-branch transitions is particularly sensitive to the excitation temperature in the model fits. We do not fit the P-branch as its S/N is lower and its flux is thus more sensitive to the continuum determination and correction from ice absorption. For C₂H₂, only the Q-branch is clearly detected and thus used for fitting. The slab model fits for both species indicate optically thick emission, with well-constrained temperatures of 120 ± 10 K and 180 ± 30 K respectively. The best-fit emitting area radii of CO₂ (${180}_{-30}^{+110}\text{au}$) is larger than that of C₂H₂ (${14}_{-2}^{+14}\text{ au}$) by a factor of ~13. The emission for HCN has a similar best-fit temperature of ~60 K, but is optically thin and has a poorly constrained emitting area radius of >400 au.

Optically thin emission from the isotopologues ¹³CO₂ and ¹³CCH₂ is tentatively detected, but due to the low S/N of these features, the exact column densities and temperatures cannot be well determined when fit without any constraints. We thus perform the fitting for ¹³CO₂ and ¹³CCH₂ with the model T_ex and R fixed to the best-fit values from CO₂ and C₂H₂ respectively, and report the resulting column densities as upper limits. The ¹³C/¹²C column density ratio is thus >14 using the CO₂ and ¹³CO₂ fits and >9 using the C₂H₂ and ¹³CCH₂ fits. This ratio is consistent with what is expected from the ISM abundance at a galactocentric distance of 11 kpc of 87 ± 15 from the models of Milam et al. (2005), but higher S/N data is needed to place better constraints on the ¹³C/¹²C ratio.

Faint OH emission is detected from both sources and is blended with other molecules except for the groups of transitions at 16.0 and 16.8 micron. The OH emission appears optically thin and has a high T_ex of ${1800}_{-500}^{+600} \text{K}$, as expected from its highly non-LTE formation process (see Sect. 4.5).

For molecules at shorter wavelengths where source A and B can be spatially resolved, there are stark differences between the two spectra. Emission from the P-branch of CO is only detected towards source A. The best-fit solutions suggests optically thin emission with N < 5.6 × 10¹⁸ cm⁻² and T_ex of ${2400}_{-1000}^{+700} \text{K}$, but solutions with optically thick emission and relatively cooler temperatures of 500–1000 K are also plausible at a 2σ level (see Sect. 4.4). H₂O is detected in emission from source A but in absorption from source B. The rovibrational emission towards source A is detected at low S/N, and the fitting is further complicated by the difficulty of determining the continuum for such faint features. In Fig. B.1, we thus only show a representative model with a by-eye fit to illustrate the detection. Here, the emitting area radius is set to the best-fit for CO₂ of 180 au, a column density of 10¹⁴ cm⁻², and excitation temperature of 500 K. A handful of pure rotational transitions of H₂O in source A are also present at 16–17 µm, though at this wavelength blending with source B is strong. Here, these transitions are consistent with optically thin emission at a temperature of ~350 K but and column densities of < 10¹⁸ cm⁻². In source B, the absorption model has a best-fit temperature of ${550}_{-300}^{+500} \text{K}$ and column density of ${8.2}_{-5.7}^{+77.3}\times {10}^{18} {\text{cm}}^{-2}$, which is not well constrained due to the low S/N and uncertainty in the continuum determination for H₂O transitions deep within ice absorption features.

Fig. 5 As Fig. 3, but for the CO and OH emission spectra in the S1 knot.

Fig. 6 Combined best-fit LTE model for 12CO2, 13CO2, C2H2, 13CCH2, HCN, and OH emission in source A+B (red shaded region) overlaid on the continuum subtracted data with an arbitrary flux offset (black). Prominent emission lines are labelled. Best-fit results for individual molecules are shown as the coloured lines below, with the horizontal bars indicating the regions of spectrum used for fitting. Extinction by CO2 ice has been included in the models, which suppresses the flux from the CO2 gas P-branch (see Sect. 2.4).

3.2.2 H₂ hotspots

In S1, we only detect CO P-branch, OH, and H₂ emission. Towards the other H₂ hotspots, no other molecular emission is detected. Notably, no CO₂ emission is detected, as was seen with ISO Short Wavelength Spectrometer (SWS) towards the Orion-IRC2/KL shock positions (Boonman et al. 2003c). The CO emission appears to be optically thin and consistent with a well-constrained temperature of ${1900}_{-200}^{+300} \text{K}$. The best-fit emitting radius is larger than a few au, but degenerate with column density. The OH emission has a best-fit solution with ~1300 K optically thick emission, but higher temperatures and an optically thin solution are also possible, as is the case for source A+B.

4 Discussion

Here, we briefly discuss some of the options for the origin of the molecular emission, before discussing each molecule in turn, as they are likely tracing distinct components in IRAS 23385+6053.

4.1 Origin of the emission

Identification of the origin of the emission is complicated by several factors. The spatial resolution of MIRI/MRS at the distance to IRAS 23385+6053 ranges from 1300 to 5000 au, and thus the extraction apertures for our spectra (Fig. 1) may contain emission from a diverse range of environments, including the hot inner disk, warm and colder envelope, irradiated outflow cavity, and gas in outflows with both C and J type shocks. Furthermore, the emission from sources A and B begin to overlap in channel 3, and are entirely blended in channel 4. The fits to source A+B thus indicate the properties of molecular emission in a large region, probably containing at least two disks.

The assumption of LTE at the gas temperature T_kin is likely not appropriate for several of the detected species, as this requires that collisions are the dominant excitation mechanism, and that the density of molecular hydrogen exceeds the critical density for the observed transitions, which are of order 10¹² cm⁻³ for rovibrational transitions (Bruderer et al. 2015). While the typical densities in disks may be large enough for colli-sional excitation, rovibrational molecular emission has also been observed in the much more tenuous gas of protostellar outflows (e.g. Sonnentrucker et al. 2006, 2007; Tappe et al. 2012). Here, radiative pumping from the mid-IR thermal dust continuum or UV emission produced in shocks may instead be responsible for the excitation, and a non-LTE model including these effects would be more appropriate. If mid-IR pumping is the dominant excitation mechanism, the best-fit LTE model temperature may reflect the temperature of the radiation field, which likely originates from the thermal dust continuum. Furthermore, pumping can cause the column density to be overestimated by several orders of magnitude. This effect has recently been identified in observations of SO₂ emission from protostars at both mid-IR and sub-mm wavelengths, which confirm the presence of strong mid-IR pumping (van Gelder et al. 2024). However, in general the strength and shape of the local radiation fields are highly uncertain, and depend on the local extinction and source geometry. Even if molecular emission does originate in the disk, large gradients in temperature and density should be present, and the assumption of a single temperature may not be valid. Nonetheless, the LTE slab models provide a good fit to the observations, except for the case of H₂O. We discuss in the following sections the case of each molecule in detail, and note where non-LTE effects may be important.

4.2 Simple organics – CO₂, C₂H₂, HCN

Emission from the simple organic molecules CO₂, C₂H₂, and HCN has often been detected in absorption towards high mass protostars in the hot core phase (e.g Lahuis & van Dishoeck 2000; Boonman et al. 2003c; Barr et al. 2020). However, IRAS 23385+6053 is likely in an earlier stage of evolution, as indicated by the lack of a compact HII region or free-free emission (Molinari et al. 1998), and the low number of detected COMs associated with high temperature chemistry and ice sublimation (Cesaroni et al. 2019; Gieser et al. 2021). The general lack of absorption, which requires colder absorbing gas in front of a warm continuum along the line of sight, is thus not surprising. The similar temperatures (100–200 K) and emitting areas (10s to 100s of au) of the LTE model fits suggest a common origin for the simple organics. For source A+B, the CO₂ and C₂H₂ emission is optically thick while the HCN emission is optically thin, which is consistent with models that predict the strongest HCN transitions become optically thick at about an order of magnitude higher column density than C₂H₂ (Lahuis & van Dishoeck 2000).

The simple organics may be tracing a variety of components of a young protostellar system. Previous Spitzer observations of protostars have found emission from these species to be associated with the outflows. In Spitzer Infrared Spectrograph observations of the Cepheus A star-forming region (d = 690 pc), Sonnentrucker et al. (2006, 2007) found spatially extended (~21 000 au) CO₂ and C₂H₂ emission with temperatures of 50–200 K and column densities of 10¹⁷ to 10¹⁹ cm⁻². In the Orion-IRc2/KL (d ~ 450 pc) region, Boonman et al. (2003c) similarly detected CO₂, C₂H₂, and HCN emission from H₂ bright shocked regions. In both cases, the low densities of the outflow regions relative to the critical densities of the simple organics supported mid-IR pumping by the local dust as the dominant excitation mechanism for the observed emission. However, while Sonnentrucker et al. (2006, 2007) find a strong correlation between the spatial extent of the simple organics and extended H₂ emission, in IRAS 23385+6053 these molecules are only detected towards source A and B, have emitting area radii < 200 au, and are not seen at the outflow knot positions (cf. Fig. 2 and Table 1). This would suggest that the simple organics do not primarily originate in the shocks, and are instead more likely to be associated with a protostellar disk and/or dense inner envelope. The prior detection towards source A of large scale gas motions resembling combined keplerian rotation and infall onto a nearly face on disk supports the disk scenario (Cesaroni et al. 2019).

The large emitting area radii (> 10s of au), temperatures of 100–200 K, and lack of strong water emission argues against the emission arising from the hot inner disk, as has typically been found in the more evolved T Tauris objects (Carr & Najita 2011). Instead, the emission may come from the extended warm disk surface and dense inner envelope that would be expected around young high mass protostars. The emitting radius for CO₂ of ~ 180 au is much larger than a typical T Tauri star disk, but is consistent with interferometric observations of high mass protostars that suggest disks of hundreds to thousands of au in radius (e.g. Johnston et al. 2015; Maud et al. 2019). Models of high-mass embedded disks by Nazari et al. (2022, see also van Dishoeck et al. 2023) indeed show typical temperatures of 150–200 K at radii of ~100–200 au. The dense environment of the disk should exceed the critical densities of the simple organic molecules, and in this case the assumption of LTE may be reasonable. However, we cannot rule out a contribution from mid-IR pumping from the dust in the disk or outflows. The lack of mid-IR continuum emission towards the H₂ knot positions suggests that pumping would only be important close to source A and B, but more modelling of the location radiation field is required to evaluate this scenario fully and is beyond the scope of this paper.

4.3 H₂O

Similar to the case of the simple organics, H₂O is likely not associated with large scale molecular outflows, as the emission and absorption is only detected towards sources A and B, and none of the H₂ knot positions (cf. Fig. 2 and Table 1). This is in stark contrast to the ubiquitous detection of H₂O far-IR and sub-mm lines in shocks and outflows (van Dishoeck et al. 2021). Mid-IR H₂O emission in protostars has not been well characterized previously owing to the relative lack of sensitivity of previous space missions compared to JWST. In the low mass protostars of the Spitzer c2d survey, H₂O emission was rarely detected (12/43 sources), and for all but one source, only in the more evolved class I stage (Lahuis et al. 2010). For these more evolved sources, the H₂O emission appeared to have similar properties to class II disks, with hot temperatures and small emitting radii. In the sole mid-IR detection towards a class 0 source, the H₂O emission appeared cold (T_ex ~ 170 K) and was suggested to originate from an accretion shock onto the disk (Watson et al. 2007), although this was later disputed by an analysis of far-IR H₂O emission with Herschel Herczeg et al. (2012). In source A, the parameters of the H₂O emission are not well determined, but the temperature of our representative model of ~500 K would be consistent with an origin in the disk surface.

In previous observations of high mass protostars, H₂O has typically been detected in absorption (Boonman et al. 2003a; Boonman & van Dishoeck 2003; Indriolo et al. 2020; Barr et al. 2022). Two non-mutually-exclusive possible scenarios involve a disk origin of the absorption and have been proposed for previous observations of protostars. In the first, absorption occurs in the disk atmosphere against the continuum from the midplane, which is viscously heated by accretion (e.g Barr et al. 2022). In the second, the disk is viewed nearly edge on relative to the observer, and the absorption occurs against the continuum from the hot inner disk on a line of sight intersecting the disk (e.g., Lahuis et al. 2006; Knez et al. 2009). A difference in inclination is an attractive possibility for explaining why we see H₂O in emission in source A but in absorption in source B. However, a nearly edge-on disk would also be expected to produce strong absorption features from CO₂, C₂H₂, and HCN which we do not detect (Lahuis et al. 2006; Knez et al. 2009). A third possibility is that the H₂O absorption originates in a dense inner disk wind against the continuum from the disk, but this may require a very specific inclination angle of the disk which matches the opening angle of the wind. We propose that the most likely scenario is that the disks in source A and B are both viewed at a low to moderate inclination, but that the viscous heating of the mid-plane is stronger in source B due to a much higher accretion rate, resulting in absorption. In source A, the weak H₂O emission then traces the warm inner disk, similar to what has been observed for class II disks (Carr & Najita 2011).

4.4 CO

In MIRI/MRS observations, the detection of any CO emission indicates highly energetic excitation conditions, as the shortest wavelength cutoff of 4.9 µm limits us to the detection of P-branch transitions with high upper energy levels (E_up > 4600 K, J_up > 24). The best-fit LTE CO models for source A and source S1 appear to indicate optically thin emission with high temperatures T > 1900 K. However, even with the simplifying assumption of LTE, it may not be possible to robustly determine the gas temperature due to the combination of optical depth effects and the lack of lower energy transitions (Herczeg et al. 2011). To illustrate this, we show in Fig. 7 synthetic rotational diagrams of the CO υ = 1 − 0 P-branch produced using a 1000 K LTE slab model at a variety of column densities, with the optical depth of each line indicated by the colour scale. At E_up ≳ 4600 K (dotted line), the υ = 1 − 0 P-branch transitions lie within the MIRI/MRS range at λ > 4.9 µm. A linear fit to the MIRI/MRS range transitions recovers the true temperature of 1000 K for a very optically thin model (N = 10¹⁷ cm⁻²), but is increasingly biased towards higher temperatures as the emission becomes optically thick. While our LTE models can account for optical depth effects, the curvature in the rotational diagram is only easily recognized in the lower upper-energy transitions located at < 4.9 µm. This effect similarly holds for the υ = 2 − 1 P-branch transitions in the MIRI range, which have even higher upper energy levels of > 7000 K. Observations with instruments covering shorter wavelengths (e.g. JWST/NIRSpec or VLT/CRIRES+) are therefore required for accurate modelling.

With this in mind, we can only conclude reliably that the observed CO emission in sources A and S1 originates from gas with a high excitation temperature, T ≳ 1000 K. What could explain this CO emission (and the non-detection of CO emission in source B)? CO rovibrational emission from protostars has been studied in detail from ground based high spectral resolution observations. Analysis of the emission line velocity profiles typically indicates that the emission originates from a broad and hot component (T ~ 1000 K) in the inner edge of the disk – similar to what is commonly observed for Class II disks (Salyk et al. 2011) – or within a narrow and cooler (T ~ 300 K), component associated with a disk wind or outflow shocks (Pontoppidan et al. 2002; Herczeg et al. 2011). Assuming that LTE holds, an inner disk-origin of the CO emission seems most reasonable for source A, and would be consistent with the other molecule detections. However, CO can be significantly pumped by strong mid-IR radiation (Blake & Boogert 2004; Thi et al. 2001), and an outflow origin cannot be ruled out. An outflow origin is more likely the case for the H₂ knot in source S1. Here, strong outflow activity is suggested by the high H₂ line velocities, but we do not detect mid-IR continuum emission, nor any of the H₂O or simple organics seen towards source A B that suggest a disk. For source B, the reason for the non-detection of CO is unclear, as a disk is likely present given the significant H₂O absorption, and there is significant outflow activity within this aperture as well. High spectral resolution studies of CO have shown that the line profiles for some sources can show a complex velocity structure with a combination of emission and absorption from the disk and disk winds (Herczeg et al. 2011; Thi et al. 2010). At low spectral resolution, blending of these features may smooth out the line profiles to a very low flux, making the CO difficult to detect. This scenario would be more likely if the H₂O absorption indeed originates from a disk wind. Alternatively, the emitting region for CO may be much smaller in source B, and thus the emission would simply be below our detection threshold.

Fig. 7 Rotational diagrams of CO υ = 1 − 0 P-branch transitions produced using our LTE slab model for four different column densities at a temperature of T = 1000 K; the colour of each symbol indicates the line optical depth. The dotted line at Eup = 4600 K is the energy above which transitions fall shortwards of the MIRI/MRS blue wavelength limit of 4.9 µm. For each model, the rotational temperature derived from a linear fit to transitions in the MIRI range (i.e. above Eup = 4600 K) is shown by the black lines.

4.5 OH

The OH transitions detected towards sources A+B and S1 are best fit by optically thin LTE models with high excitation temperatures of T_ex > 1300 K. However, the detected transitions from 13 to 17 µm have very high upper levels (E_uv ~ 8000–10 000 K, N = 21–17) as well as large Einstein A coefficients (A_ij ~ 100 s⁻¹), indicating that collisional excitation is not plausible as an excitation mechanism and that strong non-LTE processes are at play. Mid-IR OH emission has previously been detected in environments with strong UV irradiation, such as shocks from protostellar outflows (Tappe et al. 2008, 2012) and the surface layers of class II disks (Carr & Najita 2014). Here, OH transitions in extremely high rotational states (E_up ~ 40 000 K) but low lying vibrational states are observed. These transitions can be explained by the photodissociation of H₂O by UV photons in the $\tilde{B}=114-143 \text{nm}$ absorption band, which produces OH in extremely hot rotational states, followed by ‘prompt’ emission of photons associated with a cascade of pure rotational transitions to the ground state (Tabone et al. 2021; Zannese et al. 2023). However, these prompt transitions lie primarily at 9–11 µm and are not detected in our observations. Furthermore, the rotational lines of OH are split into a quadruplet by spin-orbit coupling and Λ-doubling, and for prompt OH emission, a strong asymmetry is expected between the components, but we do not observe this (Tabone et al., in prep.). A different non-LTE mechanism is the formation pumping of OH to highly excited rotational and vibrational states as the result of the O + H₂ → OH + H reaction (Liu et al. 2000), which can produce the emission long-wards of ~13 µm we observe. Further non-LTE modelling is required to understand the excitation mechanism, but is beyond the scope of this paper. The lack of detection of OH towards the H₂ knots S2 and S3 could be explained by weaker shock activity at these locations, resulting in OH emission below our detection limits.

4.6 Chemistry

Here, we comment on the chemistry for the simple organics CO₂, C₂H₂, and HCN, which have relatively precise LTE model fits. We caution that if significant mid-IR pumping is present, this would bias our LTE model fits towards higher column densities and abundances than are actually observed, and the excitation temperature would reflect the temperature of the radiation field rather than the kinetic temperature of the gas. The H₂ column densities used for our abundance determination also likely traces gas in the outflows as well as the disk and envelope, which would bias out abundances towards lower values.

We proceed here assuming that emission from these molecules originates in the disk and is indeed in LTE. What then is the origin of the molecules themselves? The presence of deep CO₂ ice absorption features at ~14.9 µm clearly indicates the presence of abundant CO₂ on the grain mantles along the line of sight in the larger scale envelopes. The dust temperature towards sources A and B estimated from a black-body fit of ~150 K exceeds the sublimation temperature of CO₂ of ~55 K, indicating thermal desorption closer to the protostars. The detection of [FeII], [SI] towards source A and B (Gieser et al. 2023) implies strong shocks capable of sputtering ices from the grain mantle are also present. Either or both mechanisms could release CO₂ into the gas phase, and a similar scenario could apply for the C₂H₂ and HCN emission, although both of these molecules have yet to be detected in ices.

Gas-phase chemistry may also contribute to the formation of the observed molecules. This may particularly be likely for CO₂, which has enhanced production by the reaction OH + CO → CO₂ + H at temperatures of 100–250 K. Above this temperature range, Oxygen is instead driven into H₂O by OH + H₂ → H₂O + H, and CO₂ enhancement is not expected (van Dishoeck et al. 2023). For HCN and C₂H₂, gas phase production can occur efficiently at temperatures >200 K via the CN + H₂ → HCN + H reaction for HCN, and through multiple pathways for C₂H₂ (Doty et al. 2002). However, our LTE model fits find excitation temperatures < 200 K for both molecules, suggesting that gas phase formation is not important.

In the context of hot core observations of gas-phase absorption, the abundance of C₂H₂ and HCN has been found to be correlated with T_ex (Lahuis & van Dishoeck 2000), while CO₂ shows no such trends (Boonman et al. 2003b). The abundance of HCN is only constrained by our modelling to <10^–5, but C₂H₂ and CO₂ abundances of ~10⁻⁷ are consistent with hot cores with T_ex ~ 200 K. The hot core chemistry models of Lahuis et al. (2007) provide predictions for the CO₂, C₂H₂, and HCN abundance over a range of excitation temperatures and X-ray fluxes. Comparing our CO₂ model fit results with Fig. 3 of Lahuis et al. (2007), at a T_ex of ~120 K the CO₂ abundance is very sensitive to X-ray flux. However, our best-fit abundance of ~2 × 10⁻⁷ is possible for some models with lower X-ray flux (<1 erg s⁻¹ cm⁻²). For C₂H₂ at ~180 K, very low abundances of <10⁻⁹ are expected, but we have ~3 × 10⁻⁷. Deeper JWST searches for both HCN and C₂H₂ ices are warranted, which could mitigate the C₂H₂ discrepancy.

5 Summary and future outlook

In Fig. 8, we summarize the proposed origins of the different molecular emission components observed towards IRAS 23385+6053. Taken together, the detected molecular species are largely consistent with an origin in a young disk and inner envelope, with a possible contribution from a dense inner disk wind. The major caveats to this picture are the limited spatial resolution at the large distance to the source, and the possibility of non-LTE effects, which we discuss in Sect. 4.1. With this in mind, the interpretation of our models is as follows:

• In the blended spectrum from source A+B, warm 120–180 K emission of CO₂, C₂H₂, and HCN is seen. The origin is most likely predominantly in the warm disk surface or dense inner envelope, but some contribution from the base of an outflow within the JWST beam cannot be ruled out. These species are not detected towards the shock positions in source S1–S3;

• H₂O is only detected towards source A and B, and not towards any shock positions, suggesting a disk rather than a large scale outflow origin. In source B, strong accretion heating in the disk midplane can produce a temperature inversion, resulting in hot 250–1050 K H₂O absorption. Alternatively, the absorption may occur in a dense inner disk wind, though this requires a very specific geometry. In source A, the accretion heating is weaker, and some faint H₂O emission is detected, which may originate in the inner disk.

• CO tracing an inner disk and/or outflow shocks is detected towards source A and S1. Only high-J P-branch states are detectable in the MIRI/MRS range of 4.9 – 28.3 µm, which suggests high temperatures T > 1000 K. Observations at shorter wavelengths covering the lower J lines are needed to determine the temperature if the emission is optically thick. CO is not detected towards source B, possibly due to absorption in a disk wind being filled in by emission from the disk;

• Highly non-LTE OH emission is detected towards source A, B, and S1, but not S2 or S3. Photo-dissociation or chemical formation pumping is likely responsible for the emission in high rotational states observed, although we do not detect the prompt OH transitions at 9–11 µm or Λ-doublet asymmetry associated with photo-dissociation.

There are many opportunities for complementary observations and modelling to improve the exploitation of MIRI/MRS data. For complex and distant high mass star-forming regions, spatial resolution is still a limiting factor. Long-baseline observations with ALMA can approach the scale of the disk, though these are not possible for IRAS 23385+6053 given its high latitude relative to the ALMA observatory. At shorter wavelengths, ground based observations of CO and H₂ O with high spectral resolution can provide valuable kinematic information for distinguishing between a disk and outflow origin, although sources like IRAS 23385+6053 are too faint for 8m class telescopes. Observations with JWST/NIRSpec would capture the full CO rovibrational band and a variety of accretion tracing HI lines accessible in the near-IR. On the modelling front, non-LTE analysis will be important for molecular emission detected in the high energy and low density environments of protostellar outflows. In particular, analysis of prompt OH emission produced by photo-dissociation of H₂O can provide valuable insights into the local UV fields and density of H₂O.

The unprecedented sensitivity of JWST is opening a new window into high mass star formation allowing extremely faint and distant sources to be detected and characterized for the first time. With the imminent arrival of additional MIRI/MRS observations of high mass protostars, our understanding of how a very young object such as IRAS 23385+6053 fits into the broader evolution will greatly improve.

Fig. 8 Schematic of our proposed origins for the molecular emission and absorption observed towards source A (right half), B (left half), and S1 (outflow shocks, centre) in IRAS 23385+6053.

Acknowledgements

We would like to thank the anonymous referee for their very constructive discussions on this paper. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. 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. These observations are associated with programme 1290. Astrochemistry in Leiden is supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program-meme (grant agreement no. 101019751 MOLDISK), by the Netherlands Research School for Astronomy (NOVA), and by grant TOP-1 614.001.751 from the Dutch Research Council (NWO).

Appendix A Thermal continuum and baseline modelling

Here, we show examples of the local continuum modelling procedure applied to our data and additional extinction corrections applied to our gas-phase emission models. In Figure A.1, examples of cubic spline fits to the local continuum at shown by the black points and solid lines overlaid on the data for source A in green. Additionally, the thermal dust continuum estimated with a cubic spline fit (dashed line) to the red data points is indicated. These estimates of the thermal dust continuum are used to calculate the wavelength dependent extinction of the ice features τ_ice = − log(F_dust/F_continuum). The wavelength dependent extinction is applied as a correction factor exp(−τ_ice) to our LTE models for emission from molecular species overlapping with the ice features, as demonstrated in Figure A.2 for the gas of gas-phase CO₂. The effect in this case is to suppress the emission from the Q– and P– branch.

Fig. A.1 Example of local continuum spline fitting and modelling of the CO2 (top panel), H2O (middle panel), and ${\text{NH}}_4^{+}$ (bottom panel) ice absorption feature’s optical depth. The dashed line shows the model for the thermal continuum fit to red points on either side of the ice feature, while the solid black line indicates the cubic spline fit. The data is shown in green for source A. The majority of unlabelled emission features in the lower panel are H2O emission.

Fig. A.2 Correction for differential extinction by CO2 ice of the CO2 gas emission features. The blue line shows the continuum subtracted data with the CO2 ice feature removed. The best-fit LTE model of CO2 with extinction correction, but without any correction for ice absorption (orange) overproduces the Q and P-branch emission. The estimated optical depth of the CO2 ice feature (dashed grey line) is applied as a scaling factor to the same best-fit model (green curve), which significantly improves the fit to the Q-branch, though the P-branch features are now somewhat under-produced.

Appendix B Model fits

In this section, we show the best-fit LTE models for each molecule in the same manner as Figure 6. In Figures, B.1, B.2, and B.3, the LTE model for each fit is shown by the coloured lines, and compared to the continuum-subtracted data (black) in the overlaid red-shaded region.

Fig. B.1 As Figure 6, but for the best-fit LTE models of CO, H2O emission at 5.5–7.7 µm and 16–18 µm in source A.

Fig. B.2 As Figure 6, but for the best-fit LTE model of CO and OH emission in source S1.

Fig. B.3 As Figure 6, but the observed data for source B is shown in units of optical depth with an offset of −0.5 relative to the best-fit LTE model for H2O absorption.

Appendix C LTE model χ² maps

The resolution of the LTE model grids used to fit our observations are the same in column density, but vary in temperature depending on the S/N of the emission and excitation conditions. For all species, the column density range is N = 10¹⁴ –10²⁴ cm⁻² in steps of log(N) = 0.166. For the relatively high S/N emission from cold CO₂, HCN, and C₂H₂, the temperature range is T = 10 – 500 K in steps of 10 K. For the highly excited emission from CO and OH, a range of T = 10 – 500 K in steps of 67 K is used. For all other species, the the temperature in the grid ranges from T = 25 – 1500 in steps of 25K. Our LTE model fits are shown for source A+B in Figure C.1, for source A and B separately in Figure C.2, and for source S1 in Figure C.3.

Fig. C.1 Normalized χ2 maps (colour scale) of the LTE model temperature and column density grids for CO2, C2H2, HCN, and OH model fits in source A+B. The best-fit emitting radius for each model is overlaid as white contours, while the 1, 2, and 3 σ confidence interval contours are overlaid in black. The best-fit model at the minimum χ2 is shown by a blue cross.

Fig. C.2 As Figure C.1, but for CO emission in source A and H2O absorption in source B

Fig. C.3 As Figure C.1, but for CO and OH emission in source S1.

References

All Tables

All Figures

Fig. 1 First three channels of MIRI/MRS observations of IRAS 23385+6053. The left panel shows the H2 S(7) transition in channel 1 at 5.511 µm, while the centre and right panels show the averaged continuum in channels 2 and 3. The colour scale in all panels is normalized to the peak intensity. The arrows in the left panel indicate the directions and velocity shift of the major outflows proposed by Beuther et al. (2023). The extraction apertures for our spectra are shown in each panel by the coloured circles. Aperture A+B has a constant radius of 2.5", while all other apertures have aradius that scales with wavelength as r = 2(1.22λ/6.5 m). The grey contours overlaid on channel 2 are 1.3 mm continuum emission (Gieser et al. 2021). The white circle in the bottom left indicates the mean beam size of JWST in the given channel, while the grey ellipse in the second panel indicates the FWHM of the mm beam.

In the text

	Fig. 2 Spectra of sources A+B (top panel) and the brightest H2 knots (lower panel) extracted from the apertures shown in Fig. 1. Flux offsets of −2.5 and −5.0 mJy are applied to the spectra of S2 and S3 respectively. In contrast with the A+B spectrum, no molecular lines other than H2, CO, and OH in S1 are detected towards the H2 knots.
In the text

	Fig. 3 Zoom in on the molecular emission from CO2, C2H2, HCN, and OH emission seen in the A+B aperture (see Fig. 1). No continuum subtraction has been applied.
In the text

	Fig. 4 As Fig. 3, but for the CO emission and 4.9–6.7 µm H2O absorption in the spectra of source A (green) and B (blue).
In the text

	Fig. 5 As Fig. 3, but for the CO and OH emission spectra in the S1 knot.
In the text

Fig. 6 Combined best-fit LTE model for 12CO2, 13CO2, C2H2, 13CCH2, HCN, and OH emission in source A+B (red shaded region) overlaid on the continuum subtracted data with an arbitrary flux offset (black). Prominent emission lines are labelled. Best-fit results for individual molecules are shown as the coloured lines below, with the horizontal bars indicating the regions of spectrum used for fitting. Extinction by CO2 ice has been included in the models, which suppresses the flux from the CO2 gas P-branch (see Sect. 2.4).

In the text

Fig. 7 Rotational diagrams of CO υ = 1 − 0 P-branch transitions produced using our LTE slab model for four different column densities at a temperature of T = 1000 K; the colour of each symbol indicates the line optical depth. The dotted line at Eup = 4600 K is the energy above which transitions fall shortwards of the MIRI/MRS blue wavelength limit of 4.9 µm. For each model, the rotational temperature derived from a linear fit to transitions in the MIRI range (i.e. above Eup = 4600 K) is shown by the black lines.

In the text

	Fig. 8 Schematic of our proposed origins for the molecular emission and absorption observed towards source A (right half), B (left half), and S1 (outflow shocks, centre) in IRAS 23385+6053.
In the text

Fig. A.1 Example of local continuum spline fitting and modelling of the CO2 (top panel), H2O (middle panel), and ${\text{NH}}_4^{+}$ (bottom panel) ice absorption feature’s optical depth. The dashed line shows the model for the thermal continuum fit to red points on either side of the ice feature, while the solid black line indicates the cubic spline fit. The data is shown in green for source A. The majority of unlabelled emission features in the lower panel are H2O emission.

In the text

Fig. A.2 Correction for differential extinction by CO2 ice of the CO2 gas emission features. The blue line shows the continuum subtracted data with the CO2 ice feature removed. The best-fit LTE model of CO2 with extinction correction, but without any correction for ice absorption (orange) overproduces the Q and P-branch emission. The estimated optical depth of the CO2 ice feature (dashed grey line) is applied as a scaling factor to the same best-fit model (green curve), which significantly improves the fit to the Q-branch, though the P-branch features are now somewhat under-produced.

In the text

	Fig. B.1 As Figure 6, but for the best-fit LTE models of CO, H2O emission at 5.5–7.7 µm and 16–18 µm in source A.
In the text

	Fig. B.2 As Figure 6, but for the best-fit LTE model of CO and OH emission in source S1.
In the text

	Fig. B.3 As Figure 6, but the observed data for source B is shown in units of optical depth with an offset of −0.5 relative to the best-fit LTE model for H2O absorption.
In the text

	Fig. C.1 Normalized χ2 maps (colour scale) of the LTE model temperature and column density grids for CO2, C2H2, HCN, and OH model fits in source A+B. The best-fit emitting radius for each model is overlaid as white contours, while the 1, 2, and 3 σ confidence interval contours are overlaid in black. The best-fit model at the minimum χ2 is shown by a blue cross.
In the text

	Fig. C.2 As Figure C.1, but for CO emission in source A and H2O absorption in source B
In the text

	Fig. C.3 As Figure C.1, but for CO and OH emission in source S1.
In the text