JOYS: Disentangling the warm and cold material in the high-mass IRAS 23385+6053 cluster

,


Introduction
The formation of high-mass stars (stellar masses M ⋆ > 8 M ⊙ ) still remains a puzzle, despite being studied for many decades, we refer to Motte et al. (2018) and Rosen et al. (2020) for recent reviews.High-mass star-forming regions (HMSFRs) are typically located at distances of several kiloparsecs and in order to achieve a spatial resolution of a few thousand au, sub-arcsecond resolution is necessary.
With the recently launched James Webb Space Telescope (JWST), infrared (IR) bright protostars can be characterized at sub-arcsecond resolution and high sensitivity in the near-infrared (NIR) and mid-infrared (MIR) regime.The observations analyzed in this work are part of the European Consortium guaranteed time program JWST Observations of Young protoStars (JOYS, PI: E. F. van Dishoeck, id. 1290).The project targets low-mass and high-mass protostars with the Mid-InfraRed Instrument (MIRI) using the Medium Resolution Spectrometer (MRS) mode.The first target observed as part of the JOYS project is the HMSFR IRAS 23385+6053 (IRAS 23385 hereafter), also referred to as "Mol 160" (Molinari et al. 1996).
IRAS 23385 is associated with H 2 O maser activity (Molinari et al. 1996), but remains undetected at cm wavelengths (Molinari et al. 1998a).This suggests that IRAS 23385 is at an evolutionary stage prior to the ultra-compact (UC) Hii region phase (Molinari et al. 1998b).The mass of the central object and the envelope are estimated to be ≈9 M ⊙ (Cesaroni et al. 2019) and 510 M ⊙ (Beuther et al. 2018), respectively.Molinari et al. (1998b) derived a mean visual extinction of the core of A V = 200 mag with a H 2 column density N(H 2 ) of 2 × 10 24 cm −2 .
The Two Micron All Sky Survey (2MASS) data toward the region reveal that the dense core is surrounded by NIR sources and is located in between two UCHii regions with bright cm emission (Molinari et al. 2002;Thompson & Macdonald 2003).A large-scale overview toward the region is shown in Fig. 1 highlighting the 250 µm, 70 µm, and 24 µm images in the region taken with the Herschel and Spitzer space telescopes.The central region is labeled "mm", while the two nearby UCHii regions are labeled "VLA1" and "VLA2" (Molinari et al. 2002).Largescale filamentary cold dust emission is revealed by the Herschel 250 µm data.The Herschel 70 µm and Spitzer 24 µm data reveal an emission arc between the UCHii region VLA2 and the mm core.Molinari et al. (2008) detect faint 24 µm emission (with a flux density of F 24µm = 0.14 Jy) toward the position of the central dense core (their source "A"), but no emission is detected 40 m 00 s In the center and right panels, the green star, labeled as "mm", marks the 1 mm continuum peak position (Fig. 2).The blue stars, labeled as "VLA1" and "VLA2", show the position of the two nearby UCHii regions (Molinari et al. 2002).The black rectangles in the right panel indicate the JWST/MIRI MRS field-of-view of the 4-pointing mosaic in ch1A (solid) and ch4C (dashed).
at wavelengths below confirming the pre-ZAMS (zero-age main sequence) nature of the source.By modeling the spectral energy distribution (SED) the authors derive a bolometric luminosity of L = 3.2 × 10 3 L ⊙ for the core, while the luminosity of the entire region, i.e. including the two nearby bright UCHii regions (Fig. 1), is L = 1.6 × 10 4 L ⊙ (Molinari et al. 1998b).Bipolar outflows are found in broad line wings in CO emission (southwest (redshifted)-northeast (blueshifted), Wu et al. 2005) and in SiO and HCO + emission (north (redshifted)-south (blueshifted), Molinari et al. 1998b).Since the emission is not extended and the outflow lobes overlap, Molinari et al. (1998b) conclude that the direction of the north-south outflow is nearly pole-on.
IRAS 23385 is part of the NOrthern Extended Millimeter Array (NOEMA) large program "CORE" (PI: H. Beuther) which aims at imaging the central dense core at an angular resolution of 0 ′′ .5 at 1 mm (Beuther et al. 2018).Cesaroni et al. (2019) carried out a case-study toward IRAS 23385 using the CORE data and suggested the presence of a northwest (redshifted)-southeast (blueshifted) outflow.These authors, based on the observed IR, mm, and cm emission, concluded that a cluster of massive stars is probably being formed in the central core of IRAS 23385.
While interferometric mm observations trace the cold material, the presence of outflows and MIR emission suggests an additional warm component that is heated by the protostars and/or shocks.As a homonuclear diatomic, the most abundant molecule in the interstellar medium, H 2 , has no permanent dipole and therefore lacks a pure rotational dipole transitions.However, rovibrational quadrupole transitions of H 2 can be excited at high temperatures, with upper energy levels E u /k B > 500 K.Such transitions can be accessed with JWST/MIRI for H 2 that covers the ∆J = 0 rotational lines in the vibrational ground state 0-0 from S(1) to S(7) (Table 1).In addition, JWST/MIRI covers emission lines from neutral and ionized atoms and molecules.
The first JOYS results of IRAS 23385 using JWST/MIRI MRS data are presented in Beuther et al. (2023).Toward the central dense core, two MIR sources are detected at short wavelengths, that become unresolved at longer wavelengths (λ ≳ 15 µm).Aside from continuum emission, the MIRI spectrum from 5 to 28 µm is dominated by emission of H 2 lines and forbidden transitions of atomic lines tracing the MIR sources and multiple outflows.In addition, the spectrum shows deep absorption features by ice species.The continuum emission at 24 µm rises to ≈0.1-0.2Jy, a similar value derived from past Spitzer data (0.14 Jy, Molinari et al. 2008).Using a faint Humphreys α line, an accretion rate of 0.9×10 −4 M ⊙ yr −1 is estimated (Beuther et al. 2023).Compact molecular line emission toward the MIR sources is being studied by Francis et al. (in prep.) and the ice properties will be presented in Rocha et al. (in prep.), see also van Dishoeck et al. (2023).
In this work, we characterize and compare in detail, for the first time, the warm and cold material, taking advantage of new JWST/MIRI MRS observations that can be compared with NOEMA data at mm wavelengths and to study the nature of the continuum sources.With NOEMA the cold molecular gas and dust properties can be studied, while with JWST/MIRI observations we probe the warmer environment.In Sect. 2 the observational parameters and calibration of the JWST and NOEMA data are explained.The analysis of the continuum and spectral line data is presented in Sect. 3 and in Sect. 4 the results are discussed.Our conclusions are summarized in Sect. 5.
The kinematic distance is estimated to be 4.9 kpc (Molinari et al. 1998b) with a Local Standard of Rest (LSR) velocity of LSR = −50.2km s −1 (Beuther et al. 2018).Typically, kinematic distance estimates suffer from high uncertainties.Molinari et al. (2008) modeled the SED with a zero-age main-sequence (ZAMS) star embedded in an envelope but were unsuccessful in explaining the observed SED properly.The authors argued that the region may actually be located at 8 kpc or that the source is at a younger evolutionary stage.On the other hand, using the Galactic rotation curve and spiral arm model by Reid et al. (2016Reid et al. ( , 2019)), a distance of 2.7±0.2kpc could be estimated 1 .
However, this distance estimate includes as a prior its association with a spiral arm which is unknown.Maser parallax measurements are not available for this source.Choi et al. (2014) measured maser parallaxes toward high-mass star-forming regions and compared the corresponding distances with kinematic distance estimates finding significant discrepancies.The authors concluded that for G111.23−1.23, the source in their sample that is closest to IRAS 23385 and has a similar LSR , the kinematic distance is 4.8 kpc ( LSR = −53 km s −1 ), while the maser parallax measurement suggests a distance of 3.33 kpc.However, since IRAS23385 is not covered by Choi et al. (2014), and given the uncertainties, we assume in this work a distance of 4.9 kpc, consistent with previous studies in the literature and in line with the distance used by the JOYS team (Beuther et al. 2023).The high uncertainty in distance does not affect the analysis in this work.
The observations were calibrated using the JWST Science Calibration Pipeline (version 1.11.2) and Calibration Reference Data System (CRDS) context "jwst_1100.pmap".In addition, astrometric corrections were applied by comparing detected sources in the MIRI imaging field, that was observed simultaneously with MIRI MRS, with Gaia data (Gaia Collaboration et al. 2021).We refer to Beuther et al. (2023) for a detailed description of the calibration of this data set.
The resulting spectral cube data product consists of continuum and line emission, as well as ice absorption features.In this work, we focus on the extended H 2 emission and forbidden atomic fine structure transitions.In Table 1 we give an overview of all lines that are analyzed in this work.With a spectral resolution of a few tens of km s −1 , the high-velocity outflow components can be resolved (Beuther et al. 2023), but a detailed kinematic analysis of the gas properties close to the protostars is difficult.Thus, in this work we focus on the spatial distribution of the integrated emission.
A 5.2 µm continuum map is created using the ch1A band at a wavelength range between 5.2 µm and 5.3 µm in which bright line emission or absorption features are absent.The angular resolution is ≈0 ′′ .2. Faint extended background emission with a median of 48 MJy sr −1 is subtracted from the 5.2 µm intensity map.The noise measured as the standard deviation of the 5.2 µm continuum map is 18 MJy sr −1 after background subtraction.
In individual line cubes (Table 1), the local continuum was subtracted by fitting a first order polynomial to the spectra masking out channels with line emission.The line noise, σ line , is estimated in line-free channels and decreases from 38 MJy sr −1 to 5 MJy sr −1 from ch1A to ch3C in which the noise increases again up to ch4C to 32 MJy sr −1 .Notes.The wavelength and upper energy level of the H 2 lines are taken from Jennings et al. (1987) and for the remaining transitions from the atomic line list by Van Hoof (2018).

NOEMA
Located in the Northern hemisphere, with a high declination of ≈61  2 and 3 in Ahmadi et al. 2018).A continuous spectrum along the full bandwidth is available with a channel width of 3.0 km s −1 .For the spectral line data, complementary short spacing observations with the IRAM 30m telescope were obtained recovering missing flux filtered out by the interferometer.The combination of the interferometric and single-dish ("merged") data is explained in Mottram et al. (2020).
The NOEMA data were calibrated using the CLIC package in GILDAS2 .The interferometric line and continuum data were self-calibrated in order to increase the signal-to-noise (S /N) ratio.Phase self-calibration was performed on the continuum data and then applied to the spectral line data.A detailed description of the self-calibration procedure for the NOEMA data can be found in Gieser et al. (2021).
The NOEMA continuum and continuum-subtracted merged spectral line data were deconvolved using the GILDAS/MAPPING package and the Clark algorithm.The weighting was set to a robust parameter of 0.1 to achieve the highest possible angular resolution.For the 1 mm continuum data, the resulting synthesized beam (with major and minor axis θ maj × θ min ) is 0 ′′ .48×0 ′′ .43 with a position angle (PA) of 58 • .The continuum noise is estimated to be σ cont,1mm = 0.15 mJy beam −1 .The properties of the 1 mm molecular lines presented in this work are summarized in Table 2 with the label "CORE" in the NOEMA project column.The line sensitivity σ line , estimated in channels without line emission, is ≈1.2K and ≈0.4 K in the highand low spectral resolution data, respectively, at an angular resolution of ≈0 ′′ . 5.

3 mm data
Using the upgraded PolyFiX correlator at NOEMA, CORE+ (PI: C. Gieser, project code W20AV) is a follow-up project of CORE targeting the 3 mm wavelength range.In a first case-study, three regions of the CORE sample, including IRAS 23385, were selected (Gieser et al., in prep.).The NOEMA 3 mm observations were obtained in the A and D configuration and short-spacing information was provided by complementary IRAM 30 m observations.The observations were carried out from February to April 2021 with NOEMA and in August 2021 with the IRAM 30m telescope.Using the high-resolution spectral units, a channel width δ of 0.8 km s −1 is achieved in the merged line data.
Phase self-calibration was performed on the continuum data using the same procedure as in the 1 mm CORE data (Gieser et al. 2021) and the gain solutions were applied to the NOEMA line data.The NOEMA continuum and continuumsubtracted merged spectral line data were deconvolved using the GILDAS/MAPPING package and the Clark algorithm.The weighting was set to a robust parameter of 1 to achieve a compromise between high angular resolution and good sensitivity.For the 3 mm continuum data, the resulting synthesized beam is 1 ′′ .1×0 ′′ .85 (PA 47 • ) and the noise is estimated to be σ cont,3mm = 0.020 mJy beam −1 .
The properties of the 3 mm molecular lines analyzed in this study are summarized in Table 2 with the label "CORE+" in the last column.The line sensitivity σ line at a channel width of 0.8 km s −1 is ≈ 0.2 K.
For all 1 mm and 3 mm NOEMA continuum and line data products we apply a primary beam correction taking into account that the noise levels increase toward the edge of the primary beam.The primary beam is ≈22 ′′ and ≈60 ′′ at 1 mm and 3 mm, respectively.Thus, a larger FoV is covered in the NOEMA data compared to the 4-pointing mosaic JWST/MIRI observations (≈ 15 ′′ at the longest wavelengths).Since we aim to compare the MIR and mm data, we only focus here on the central area (≈ 12 ′′ × 12 ′′ ) that is also covered by the JWST/MIRI FoV.The upper energy levels E u /k B of the mm lines vary between 4 K and 100 K, that is significantly lower than for the JWST/MIRI transitions (> 500 K).Thus, with the NOEMA continuum and line data, we are able to trace in the same region also the cold molecular gas.

Results
A comparison of the detected JWST/MIRI MRS and NOEMA mm lines allows us to study the properties of the warm and cold material, respectively, at sub-arcsecond resolution corresponding to spatial scales below 5 000 au for IRAS 23385.In addition, the continuum allows us to map the different sources present in the FoV.In Sect.3.1 the continuum emission, detected sources, and outflows are presented.The spatial morphology of atomic and molecular line emission is compared in Sect.3.2.The temperature of the warmer gas (>100 K) is derived from the JWST/MIRI H 2 0-0 lines in Sect.3.3, while the temperature of the cold gas (<100 K) is estimated by modeling CH 3 CN emission lines in Sect.3.4.

Continuum sources and outflows
A comparison of the JWST/MIRI 5.2 µm and NOEMA 1 mm continuum data is presented in Fig. 2 in the left panel with detected sources labeled.For the mm and MIR sources, we follow the nomenclature of Cesaroni et al. (2019) and Beuther et al. (2023), respectively.
The 1 mm continuum data clearly resolve two cores (mmA1 and mmB) that are separated by ≈1 ′′ .6, corresponding to a projected separation of ≈7 700 au.The mmA continuum peak is elongated, can be separated into two sources "mmA1" and "mmA2".
The 5.2 µm continuum data reveal two sources, A and B, that are barely resolved at an angular resolution of 0 ′′ .2. The MIR source A is co-spatial with the elongated mm core "mmA2", hereafter referred to as "A/mmA2".While in the shorter wave- The dash-dotted white circles show the aperture in which the 5.2 µm flux density F 5.2µm was derived (Table 3).In the bottom left corner, the JWST/MIRI 5.2 µm and NOEMA 1 mm angular resolution is indicated in black and grey, respectively.The mm (mmA1 and mmB, marked by green squares) and MIR (A/mmA2 and B) continuum sources are labeled in green.In the right panel, the line integrated intensity of the H 2 0-0 S( 7) line with S /N > 5 is presented in color.The JWST/MIRI 5.2 µm continuum as shown in the left panel is presented as blue contours with contour levels at 5, 10, 15, 20, 25×σ cont,5µm .The angular resolution of the 5.2 µm continuum and H 2 0-0 S( 7) line data is indicated in the bottom left and right, respectively.Several shock spots evident in the MIRI MRS line emission are marked by green crosses and labeled in green (Sect. 3.2).The red and blue arrows indicate three bipolar outflows, labeled I, II, and III (as presented in Beuther et al. 2023).In both panels, the black bar indicates a spatial scale of 10 000 au at the assumed source distance of 4.9 kpc.
length regime of MIRI the two sources can be resolved, beyond ≈15 µm this is not possible due to the lower angular resolution (Beuther et al. 2023).Notably, no additional source is detected in the MIRI data at longer wavelengths.The positions of the four continuum sources are summarized in Table 3.In addition, the flux density at 5.2 µm F 5.2µm , calculated within a 0 ′′ .5 circular aperture (Fig. 2), of the MIR sources A/mmA2 and B are listed.
For the non-detections of sources mmA1 and mmB, 3σ upper limits are given.Notes.The flux density F 5.2µm is computed within a 0 ′′ .5 aperture (Fig. 2).In case of non-detections, 3σ upper limits are listed.
The 1 mm continuum morphology surrounding the cores reveals complex substructure.The mm continuum emission is elongated toward the north from position mmA1, but also elongated toward the southeast from position mmB.The 1 mm (and also the 3 mm) continuum data trace cold dust emission that is optically thin: Gieser et al. (2021) estimate that toward mmA1 and mmB the values of the continuum optical depth at 1 mm, τ cont,1mm , are 9.3×10 −3 and 7.9×10 −3 , respectively.In the 3 mm continuum data, we also find τ cont,3mm ≪ 1 (Sect.3.4).The core masses of mmA1 and mmB are estimated to be 21.9 M ⊙ and 9.4 M ⊙ (with T kin = 73 K, Beuther et al. 2018).Thus, mmA1 contains about twice the mass compared to mmB.These mass estimates are lower limits due to spatial filtering of the extended emission in the interferometric observations.In summary, the NOEMA observations trace the optically thin dust emission of a cold envelope in which the four sources are embedded.The two MIR sources reveal two more evolved protostars with a hot thermal component.While source A has a clear counterpart at mm wavelengths (mmA2), source B does not coincide with a fragmented mm core, which may be due to a poorer angular resolution (0 ′′ .5) of the NOEMA data compared to the JWST/MIRI 5.2 µm data (0 ′′ .2).Still, source B is embedded within the cold dust envelope.The fact that multiple sources are detected toward the central region of IRAS 23385 reinforces the idea that a cluster is forming (Cesaroni et al. 2019).The sources mmA1 and mmB with no MIR-bright counterpart could be either younger compared to A/mmA2 and B or their MIR-extinction is high enough to completely block the protostellar MIR emission.The nature of the continuum sources is further discussed in Sect.

4.
The right panel in Fig. 2 shows the integrated intensity of the H 2 0-0 S(7) line at 5.511 µm (Table 1), showing H 2 emission at high angular resolution, revealing the presence of multiple bipolar outflows.Following Beuther et al. (2023), the outflows I, II, and III could be associated with the continuum sources B, mmA1, and A/mmA2, respectively.We note that due to the clustered nature of IRAS 23385, the outflow driving sources cannot be identified unambiguously.We refer here to the most likely driving source, as discussed in Beuther et al. (2023).Complementary NIR observations with JWST that have even higher angular resolution could provide a more reliable assignment.The outflow directions are indicated by arrows in Fig. 2. In addition, several positions where strong shocks occur are marked by green crosses and are labeled as S1, S2, S3, and S4.In addition, bright H 2 emission is detected toward the south of source B and the southwest of source A/mmA2 tracing shocks and heated material very close to the protostars.While these shocked locations are already evident in the H 2 0-0 S( 7) map, we find that additional forbidden transitions of atomic lines can be present there, as discussed in the following section.No outflow could be linked to source mmB, since the close-by bright H 2 knot (S1) toward the northeast of mmB is most likely connected to outflow III (Beuther et al. 2023).

Atomic and molecular line emission
In this work, we focus on the spatial distribution of H 2 and the neutral and ionized atomic forbidden fine-structure lines, while a detailed analysis of the MIR molecular emission lines other than H 2 covered by JWST/MIRI, for example CO, CO 2 , and HCN, will be presented in Francis et al. (in prep.).The H 2 transitions covered by JWST/MIRI only trace the warmer/hotter gas (Table 1, Sect.3.3), typically the heated and/or shocked gas along the jet or gas entrained by the jet.Therefore, we use CH 3 CN and the mm dust emission to trace the cold H 2 component (Sect.3.4).In this section, we first compare the spatial morphology of the emission of atoms and molecules detected by JWST and NOEMA.
All studied transitions covered by JWST/MIRI are listed in Table 1, but in the following, for H 2 we only present the spatial morphology of the 0-0 S(1) transition, while for the S( 7) transition it is shown in the right panel in Fig. 2, since the spatial morphology of the remaining H 2 lines is similar.The studied molecular transitions covered by NOEMA are summarized in Table 2.In order to obtain the spatial distribution of each emission line, we compute the line integrated intensity of the continuum-subtracted data cubes,

Molecular hydrogen H 2
The 0-0 S(1) transition of H 2 is detected everywhere within the field-of-view (FoV) of MIRI revealing at least three outflows (Beuther et al. 2023).Bright emission peaks are found slightly toward the south and southwest of the MIR-bright sources B and A/mmA2, respectively.Additional bright H 2 emission peaks are found toward the east and north of mmB and we refer to these positions as "S1" and "S2" in the following and they are also marked in Fig. 3 (top right panel).Another interesting feature is a bright emission arc toward the east at the edge of the MIRI FoV.This arc-like structure is not seen in the S( 7) integrated intensity map (Fig. 2), most likely due to a smaller FoV in the MIRI MRS data.This is not an artifact, but is most likely caused by the irradiation of the nearby UCHii region heating up the environment and is also seen in the Spitzer 24 µm image (Fig. 1).Otherwise, both H 2 transitions trace the same spatial structures in the dense core: heated and shocked gas along the outflow direction.

Atomic sulfur [S i]
Atomic sulfur, [S i], is typically associated with high-density material of 10 5 cm −3 due to a high critical density of the line in non-dissociative C-type shocks (Hollenbach & McKee 1989;Neufeld et al. 2007;Anderson et al. 2013).The atomic sulfur [S i] emission line follows a similar morphology as H 2 .The emission arc toward the east is not detected in [S i], however, another bright emission peak, "S3", toward the north of mmA1 and close to S2 has enhanced [S i] emission.This knot is also detected faintly in H 2 0-0 S(1), but more clearly in H 2 0-0 S(7) (Fig. 2).Both the H 2 and [S i] emission reveal shocked gas toward the center of the IRAS 23385 region caused by protostellar outflows.The north-south outflow (II) most likely originates from source mmA1 (Beuther et al. 2023) with strong H 2 and [S i] emission close to the source toward the south and the shocked gas at the positions S2 and S3 toward the north.

Iron [Fe ii] and Nickel [Ni ii]
Fine-structure lines from ionized atoms, such as iron (Fe) and nickel (Ni), primarily arise in dissociative shock fronts in protostellar outflows (Neufeld et al. 2009).The [Fe ii] and [Ni ii] lines have a similar spatial distribution to each other, both peak toward source A/mmA2, with elongated emission from northeast to southwest.Along this elongation another emission peak, S4, appears in these lines that is located right between the two mm continuum peak positions.This suggests that source A/mmA2 drives an outflow from the northeast to southwest (III) and that S4 and S1 reveal shocked gas along the outflow.Interestingly, the emission is not continuous, but shows emission knots suggesting an outflow with strong episodic outbursts.Outflow III is the only one that shows emission by these ionized atoms.

Neon ([Ne ii] and [Ne iii]) and Argon [Ar ii]
The spatial morphology of [Ne ii], [Ne iii], and [Ar ii] are similar, peaking toward the position of source A/mmA2 and with extended emission toward the west.Their spatial distribution is very different from that of [Fe ii] and [Ni ii].This indicates that toward this source, the internal UV and X-ray radiation is high enough to ionize these species and to excite such transitions with upper energy levels in the range of 1 000 K -3 500 K (Table 1).There is a secondary peak, 0 ′′ .9 toward the west of source A/mmA2, though only faintly seen in [Ar ii].This feature is further discussed in Sect. 4.
The [Ne ii] emission also shows an extended component with bright emission toward the edges of the FoV.While [Ne ii] can arise from J-type shocks (Hollenbach & McKee 1989) toward the protostars, the extended emission arises from energetic radiation of the nearby UCHii regions (Fig. 1).

Millimeter lines
The integrated intensity maps of the mm lines are shown in Fig. 4. The bipolar outflows can partially also be revealed in the mm line data.A common shock tracer in the mm regime is SiO, with silicon (Si) being sputtered off the grains and forming SiO in the gas-phase (e.g., Schilke et al. 1997).The SiO emission peak is toward the northern part close to the shocked positions S2 and S3 also bright in H 2 and [S i]. 13 CS also shows an elongated structure toward the northern part.SO, H 13 CO + , and H 2 CO show bright emission toward mmA1, but is also enhanced toward the northern direction toward S2 and S3.The north-south outflow (II) thus has a rich molecular component.
In SiO, a large-scale outflow (I) is seen from the northwest to the southeast, which is also revealed by the CORE 1 mm In the [Ne ii] and [Ne iii] panels, the dash-dotted grey circles show the aperture (0 ′′ .9) in which the flux density was derived (Sect.4).
OCS shows very compact emission toward mmA1.HC 3 N shows extended emission with many filamentary structures, hinting at not only outflowing but most likely also inflowing gas, and peaks toward mmA1.CH 3 CN emission is strong around mmA1 and A/mmA2.Interestingly, CH 3 OH emission significantly drops toward the locations of the MIR-bright sources A/mmA2 and B.
In none of the NOEMA emission maps, mmB shows significant emission peaks, except for CH 3 OH emission.Since all remaining detected continuum sources (mmA1, A/mmA2, and B) show protostellar activity through outflows and ionized gas, mmB is likely to be less evolved or less massive, although it cannot be completely ruled out that the H 2 and [S i] emission toward S1 could be an outflow originating from mmB that is directed along the line-of-sight.
Shocks near the protostars are caused by the emerging outflows and are seen in SiO emission, but also in sulfur-bearing species, such as atomic sulfur, SO, and 13 CS.However, OCS only emits strongly toward the source mmA1 suggesting that it may be related to ice sublimation, instead of shock chemistry.For SO, both high-temperature gas-phase chemistry close to the protostar and shock chemistry in the outflow (II) seem to be rel- evant to explain the spatial distribution.The entire region is embedded in a filamentary structure, seen in HC 3 N, where most of these converge toward the central dense core.This suggests that even though there exist at least three outflows seen in H 2 emission by MIRI (I, II, and III, Fig. 2), they have so far only little impact on disrupting the envelope in which the protostars are embedded in.Thus, further star formation activity could be still possible in the central dense core of IRAS 23385 toward the dense millimeter source mmB that currently does not show any signatures of active star formation.

Molecular hydrogen excitation diagram analysis
The JWST/MIRI MRS spectral range covers multiple transitions of the pure rotational 0-0 state of H 2 (Table 1) and with the excitation diagram analysis, the H 2 column density N(H 2 ) and ro-tational temperature T rot can be estimated (e.g., Neufeld et al. 2009).The lowest H 2 transitions are excited through collisions (e.g., Black & van Dishoeck 1987), while higher rotational levels and vibrational states can be excited through UV pumping.Thus the warm component linked to the collisions is a realistic tracer of the underlying gas temperature, while for any hotter component, the estimated rotation temperature may not correspond to a physical property.Since the angular resolution varies from 0 ′′ . 2 to 0 ′′ .8 for 5 µm to 28 µm along the full MIRI spectral range, we smooth the S(7) to S(1) transitions to a common resolution of 0 ′′ .7 that corresponds to the angular resolution of the S(1) transition and regrid the remaining lines to the same spatial grid as the H 2 S(1) data.Spectral binning is not performed.This allows us to carry out a pixel-by-pixel excitation diagram analysis and to obtain a complete N(H 2 ) and T rot map, since the H 2 transitions are detected within the entire FoV.
To provide reliable results, it is crucial to correct for extinction along the line-of-sight.We use an extinction curve derived by McClure (2009) that is valid between a K-band (2.2 µm) magnitude A K of 1 mag and 7 mag.Since the extinction toward the central dense core is expected to be high (Molinari et al. 1998b), we assume in the entire field-of-view a K-band extinction of A K = 7 mag that roughly corresponds to a visual extinction of ≈50 mag.We also tested lower values of A K , but found that the integrated intensity of the S(3) transition that is located within the broad 10 µm silicon-absorption feature (Beuther et al. 2023), is significantly decreased compared to the transitions outside of the absorption feature for these lower values of A K .
As a first step, in each spatial pixel the seven H 2 transitions using the extinction-corrected data are fitted by a 1-dimensional (1D) Gaussian and from the integral of the Gaussian fit the integrated intensity of each transition is computed.Example spectra and Gaussian fits of the seven transitions taken from the position of source mmA1 and source B are presented in Fig. A.1.The integrated line intensities toward all four continuum sources (Table 3) and four shock positions (S1, S2, S3, and S4) are summarized in Table A.1.
In general, the lines are barely resolved at the spectral resolution of MIRI.The uncertainty in the absolute flux calibration of the MIRI MRS instrument is only on the order of a few percent (Argyriou et al. 2023).A much higher uncertainty in the calculation of the line integrated intensities arises from the extinctioncorrection: We assume A K = 7 mag (A V = 54 mag) in the entire FoV, while there may be spatial differences.For comparison, we also present in Table A.1 the extinction-corrected integrated line intensities calculated with A K = 5 mag (A V = 39 mag) and A K = 3 mag (A V = 23 mag).From the 9.7 µm silicate absorption feature A V = 30 − 40 mag is estimated (Rocha et al., in prep.), whereas previous studies suggested even higher extinction values of A V = 200 mag (Molinari et al. 1998b).Given the high column densities toward the region, A K = 5 − 7 mag should be a reliable estimate of the extinction.The line integrated line intensities vary between a factor of 2−3 within this K-band extinction range (Table A.1).We thus assume in the following that the observed line integrated intensities are uncertain by a factor of two.
Moreover, only transitions with a peak intensity of 50 MJy sr −1 (in the non-extinction corrected data) and above are further considered which is above the noise (Table 1).With this threshold, we exclude areas with faint H 2 emission in the higher excited lines where the excitation diagram analysis would provide unreliable results.As an additional constraint, a fit to the excitation diagram is only performed when 5 or more transitions are detected in a spatial pixel above the threshold, as otherwise a two component fit would not provide reliable results.
The excitation diagram analysis is performed using the pdrtpy 3 python package.Typically it is found that the observed H 2 transitions are best described by two rotational temperature components, which we refer to as T warm and T hot in the following.The observed line integrated intensities are converted into upper state column densities normalized by its statistical weight, N u /g u .The H 2 rotation temperature and column density can be estimated by plotting the logarithm of N u /g u , as a function of upper state energy level, E u /k B and fitting the observed data points with a straight line for each component.The column density, log   N, is proportional to the y-intercept and the slope is proportional to −1/T .We assume that all H 2 transitions are optically thin (which is valid even for N(H 2 )>10 23 cm −2 , Bitner et al. 2008).
Examples of the excitation diagram toward the sources mmA1 and B are shown in Fig. 5.The fit results obtained from pdrtpy for the warm and hot component toward all four continuum sources (Table 3) and four shock positions (S1, S2, S3, and S4) are summarized in Table A.2. Considering all positions in the map, the uncertainties for the estimated temperature are between 10% − 30%, while for the column densities the uncertainties are on the order of a factor of 0.5 − 3.In absolute numbers, the median uncertainties are log ∆N warm = 0.3 cm −2 , log ∆N hot = 0.6 cm −2 , ∆T warm = 43 K, ∆T hot = 270 K.
The spectral line extraction and excitation diagram analysis toward the mmA1 and B sources shown in Figs.A.1 and 5, respectively, reveal bright line emission in all transitions of H 2 .The excitation diagram shows the necessity to fit the data with two temperature components.The warm and hot temperature component toward mmA1 are ≈400 K and ≈1 400 K, respectively.The total column density, considering the contribution from both temperature components is N warm+hot ≈ 5 × 10 23 cm −2 .Toward source B we find a higher column density, but lower temperature.This could be attributed to the fact that the lower excited H 2 lines may become optically thick toward the MIR sources, further discussed in Sect. 4.
In Table A.2, we also present excitation diagram fitting results with A K = 5 mag and A K = 3 mag assumed in the calculation of the extinction-corrected line integrated intensities (Table A.1).As expected, the estimated column densities are higher when the extinction is higher.However, given the high uncertainties, the differences are not significant.There is an effect on the derived temperatures.While the temperature of the hot component does not significantly change, the temperature of the warm component does show some differences, however not that significant given the high uncertainties.This is due to the fact that the S(3) line of H 2 lies in the silicate absorption feature and thus suffers from a higher extinction A λ compared to the remaining lines.Since this line lies at the border of the two components in the excitation diagram (Fig. 5), this data point does influence the fitted slope and thus the estimated temperature.
The complete temperature and column density maps of the warm and hot components are presented in Fig. 6 and in Sect. 4 we compare the results with rotation temperature estimates based on the cold gas traced by the NOEMA data.One caveat mentioned before is that the excitation diagram analysis assumes that the H 2 lines are optically thin, which may especially not be the case toward the positions MIR sources A/mmA2 and B.

Methyl cyanide line modeling
Molecular emission at mm wavelengths can be used to probe the temperature of the cold molecular gas.In Gieser et al. (2021), the temperature of the cold gas component in the CORE regions, including IRAS 23385, is estimated using H 2 CO and CH 3 CN emission lines.Unfortunately, the H 2 CO lines are optically thick toward the center of IRAS 23385 causing XCLASS to converge to unrealistic high temperatures > 200 K.However, this highlights again that toward the central core, the densities and extinction must be high.For CH 3 CN, the temperature toward mmA1 is estimated to be 100 K (Cesaroni et al. 2019;Gieser et al. 2021), however, the transitions of the J = 12 − 11 K-ladder at 1 mm have high upper state energies and thus, the emission is compact around mmA1 and barely resolved.With the complementary 3 mm CORE+ observations of the CH 3 CN J = 5 − 4 Kladder, the temperature of more extended colder gas can be investigated (Table 2).
We model the rotational transitions of the CH 3 CN J = 5 − 4 K-ladder with XCLASS.With XCLASS, the 1D radiative transfer equation is solved assuming local thermal equilibrium (LTE) conditions.The fit parameters are source size θ source , rotation temperature T rot , column density N(CH 3 CN), line width ∆ , and velocity offset off .In total, we model four transitions of CH 3 CN (J = 5 − 4 and K = 0, 1, 2, 3, Table 2) and from the best-fit, the rotation temperature can be estimated.Since these CH 3 CN transitions have upper state energy levels <100 K compared to the JWST/MIRI H 2 transitions with >1 000 K, we refer to this temperature component as the cold component, T cold .In each spatial pixel, the CH 3 CN spectrum is modeled with XCLASS if the peak brightness temperature is higher than a threshold of 0.7 K (corresponding to ≈3-4σ line , Table 2).
As an example, the observed spectrum and best-fit model derived with XCLASS is shown in Fig . A.2 toward the position of mmA1.The derived rotation temperature, column density, and line width of CH 3 CN is 46 K, 6.4 × 10 13 cm −2 , and 3.7 km s −1 , respectively.The full temperature map is presented in Fig. 6 and is discussed and compared to the temperature estimated from H 2 in Sect. 4.
Using the CH 3 CN temperature and the 3 mm continuum data, we can also estimate the H 2 column density of the cold component, if the 3 mm continuum stems from optically thin dust emission and if the gas and dust temperatures are coupled.The 3 mm continuum optical depth can be calculated according to . (1) I pix,3mm is the 3 mm intensity in a pixel, converted from Jy beam −1 to Jy pixel −1 considering the beam area Ω = π/(4ln2)θ maj θ min and the area covered by one pixel (0.25 ′′ ) 2 , and B ν (T cold ) is the Planck function.We find that τ cont,3mm < 0.02 in the entire map.Thus, the 3 mm continuum emission stems from optically thin dust emission.
Using the 3 mm continuum data, that has a similar angular resolution, and the CH 3 CN temperature map, the H 2 column density of the cold gas can be derived (Hildebrand 1983): We assume a gas-to-dust mass ratio of γ = 150 (see Beuther et al. 2018), a dust opacity of κ ν = 0.17 cm 2 g −1 (Table 1 in Ossenkopf & Henning 1994, extrapolated from 1.3 mm, and assuming thin icy mantles at a gas density of 10 6 cm −3 ), a mean molecular weight of µ = 2.8 and m H is the mass of a hydrogen atom.The flux calibration of NOEMA is expected to be accurate to a 10% level at 3 mm4 .The temperatures estimated with CH 3 CN in XCLASS have uncertainties on the order of 30% (Gieser et al. 2021).The biggest uncertainty when estimating N cold lies in the gas-to-dust mass ratio γ and the dust opacity κ ν .It can thus be estimated that the calculated N cold values can be uncertain up to a factor of 2 − 4 (Beuther et al. 2018).The results are for T cold , M cold , and N cold (H 2 ) are presented in Fig. 6 and are discussed in Sect. 4.

Discussion
The JWST/MIRI MRS observations of IRAS 23385 reveal exciting new insights into energetic processes during the formation of stars.Multiple outflows (I, II, and III, Fig. 2) can be identified, one of which shows in addition to H 2 and [S i] emission, enhanced emission by [Fe ii], [Ni ii], and [Ne ii].The outflows are also associated with cold molecular gas traced by mm lines.Here, we focus on a few main points based on the analysis (Sect.3) and discuss in detail the different gas components, the outflow properties, and the spatial emission and intensity ratio of [Ne ii] and [Ne iii].

Gas components toward the central dense core
A comparison of the warm and hot gas traced by H 2 in the JWST/MIRI observations (Sect.3.3) and the cold gas traced by CH 3 CN and the 3 mm continuum in the NOEMA observations (Sect.3.4) is presented in Fig. 6.The top and bottom panels show the H 2 rotation temperature and column density, respectively, of the cold (left), warm (center), and hot component (right).To our knowledge, this is the first time that spatially resolved maps of the cold, warm and hot components have been obtained in a high-mass cluster-forming region.The cold component consists of dense (10 23 -10 24 cm −2 ) and cold (10-60 K) gas.Interestingly, for the warm component we derive similar column densities, but spatially different, since the H 2 lines trace shocked and heated gas due to the outflows, while the mm line and continuum traces the cold envelope.The temperature of the warm component ranges between 200-400 K.
Comparing the temperature and column density, we find that these properties are anti-correlated, such that toward the high column densities the temperature is lower compared to the lower column density regions.This hints that toward the dense shocks, the lower H 2 lines could become optically thick and we only trace outer layers toward these shock positions.In the hot component, the column densities are about two orders of magnitude lower, N hot ≈ 10 22 cm −2 , compared to the cold and warm component and the temperatures are 1 000-1 800 K.
How do these temperatures and column densities compare with those found in nearby low-mass regions for which more extensive studies are available in the literature?Taking as an example the L1157 Class 0 protostar (L = 8.4 L ⊙ , Froebrich 2005), Di Francesco et al. (2020) studied the cold material using Herschel observations and find dust temperatures of 10-18 K and column densities up to 2 × 10 22 cm −2 .Nisini et al. (2010) created temperature and H 2 column density maps of the L1157 outflow and also find that the H 2 line fluxes are best described by two temperature components, with 300 − 400 K and 1 100 − 1 400 K, respectively.These results are similar to the temperatures derived in the IRAS 23385 region (Fig. 6), however, the H 2 column density map derived by Nisini et al. (2010) is about four orders of magnitude lower as expected in less massive outflows driven by low-mass protostars.Compared to previous IR studies toward HMSFRs, with JWST we are for the first time able to probe highcolumn density regions (> 10 23 cm −2 ) due to the higher angular resolution.The temperature in IRAS 23385 is also similar to low-mass protostars in NGC 1333 (Dionatos et al. 2020).In the regime of high-mass protostars, van den Ancker et al. (2000) find temperatures of the warm component of 500 K and 730 K for the HMSFR S106 and Cep A, respectively.In the Orion OMC-1 outflow, Rosenthal et al. (2000) find that bulk (≈ 72%) of warm gas has a temperature of ≈630 K.
A system that is comparable to IRAS 23385 is IRAS 20126+4104 (d = 1.7 kpc, Cesaroni et al. 1997) with a similar mass of the central protostar (7 M ⊙ ) and luminosity (10 4 L ⊙ ).Jet precession in this source also hints at underlying multiplicity (e.g., Cesaroni et al. 2005;Caratti o Garatti et al. 2008).Using higher excited transitions of H 2 in the NIR a high temperature component at 2 000 K is derived by these authors, while also a warm component is found at 500 K.Additional observations in the NIR, for example with the Near-Infrared Spectrograph (NIRSpec) instrument of JWST, are required in order to investigate the presence of potentially even higher temperature components in IRAS 23385.
Since toward IRAS 23385 the highest mass protostar is estimated to be 9 M ⊙ (Cesaroni et al. 2019), at the lower end of what is considered as a "high-mass" star, it is expected that more massive protostars can heat up the environment to even higher temperatures through powerful outflows and expanding UCHii regions.The fact that a cold component as traced by the NOEMA mm observations still exists toward the central dense core, stresses that star formation can efficiently occur in shielded dense gas.While outflows seem to impact and surrounding envelope through shocks and heating of the material, the feedback is not (yet) strong enough to completely dissipate the central dense core as outflows are collimated.

Composition of the outflows and evolutionary stage of the continuum sources
The northeast-southwest outflow III (Fig. 2) is not only seen in H 2 and [S i] emission knots, but these knots are also bright in [Fe ii] and [Ni ii] emission (Fig. 3).Dionatos et al. (2014) find [Fe ii] and [Ne ii] to trace the FUV irradiation on the disk surface and on the walls of the outflow cavity and Caratti o Garatti et al. ( 2015) also find that [Fe ii] and H 2 knots can be co-spatial.Toward the position of the source A/mmA2 we detect emission from ionized atoms from and [Ar ii] (Fig. 3).This suggests that this protostar is the most evolved with strong UV and X-ray radiation.The ionization energies are 7.6 eV (Ni), 7.9 eV (Fe), 15.8 eV (Ar), and 21.6 eV (Ne).In Fig is detected (the emission of this line is further discussed in the following).This suggests that the protostars and their outflows in IRAS 23385 are not energetic enough to produce and excite the double or higher ionized transitions.
We do not detect MIR continuum emission toward the mm continuum peak source mmA1.Either the extinction is too high toward the location of this source, for example, due to the geometry or high density, or the protostar itself is at an earlier evolutionary phase compared to the MIR-bright sources A/mmA2 and B. However, we could expect this source to become bright at the longer wavelengths within the MIRI coverage, which is not the case.The bipolar outflow (II) in the north-south direction only consists of molecules and neutral atomic lines (see Fig. 7 for the atomic line emission at the shock position S1).The fact that CH 3 CN and OCS emission is detected toward the location of mmA1, suggests that the central protostar is currently heating up the surrounding material creating a hot core region where ices sublimate from the grains.The high density may currently still prevent energetic UV and X-ray radiation and also the MIR radiation of the protostar itself to escape far from the protostar and thus remains undetected in the MIR.
The source mmB does not show any signatures of starforming activity.In fact, the CH 3 OH emission detected toward the location of mmB may be shocked gas caused by the northeast-southwest outflow (III).Though it may not be a coincidence that mmB is along this outflow direction as it may be cold material that was swept up by the outflow (III).Tychoniec et al. (2021) created an overview of molecular lines that trace different physical entities surrounding low-mass Class 0 and I protostars, for example, the hot core region, outflows, UV-irradiated cavity walls, jets, disks, and the surrounding envelope (see their Fig.11 and Table 2).All molecular lines at mm wavelengths presented in this work (Fig. 4), except for HC 3 N, are covered by their overview.We find a similar spatial distribution in IRAS 23385, though due to the large distance we are not able to resolve a potential disk or differentiate between the outflow and cavity walls.Interestingly, we do not find strong emission by complex organic molecules (COMs), except for CH 3 OH and CH 3 CN.The observed CH 3 OH emission traces the outflow (I), it is non-thermally desorbed and likely sputtered as observed in Perotti et al. (2021).Given that the temperature of the cold component is ≈60 K, the COMs could still mostly reside in ices on the dust grains.The composition of the MIRI MRS ice absorption spectra could reveal the presence of COMs in the ices and will be presented in Rocha et al. (in prep.).Early JWST observations toward low-mass protostars and their envelopes have already revealed a variety of different ice species (Yang et al. 2022;McClure et al. 2023).
For both outflows III and II (Fig. 2), the atomic sulfur is cospatial with H 2 knots caused by shocks due to the outflow (Fig. 3).This suggests that sulfur is being released by C-type shocks (similar to outflows of Class 0 protostars, Anderson et al. 2013).Since the observed [S i] line at 25.2 µm has a large critical density (10 5 cm −3 , Hollenbach & McKee 1989), it will be readily detected in high-density shocked gas.Its abundance can be enhanced with sputtering of the grains in shocks.
Shocks in star-forming regions can occur through various processes, e.g., cloud-cloud collisions, jets and molecular outflows, stellar winds, or expanding Hii regions.In C-type shocks, the maximum temperature is 2 000−3 000 K and shock velocities are <40 km s −1 , that is too low to destroy molecules, while in J-type shocks temperatures can reach 10 5 K with increased fine structure emission from atoms such as Fe, Mg, and Si that sputtered off the grains and shock velocities that are typically higher than 40 km s −1 (Hollenbach & McKee 1989;van Dishoeck 2004).With the relatively low spectral resolution of JWST/MIRI MRS (R ≈ 3 700 − 1 400) it is difficult to estimate the shock velocities, however, the presence of [Fe ii] (velocity-resolved maps from −185 km s −1 to 120 km s −1 are presented in Beuther et al. 2023, their Fig. 6) and [Ni ii] in the northeast-southwest outflow (III) suggests that there could be a contribution of J-type shocks, while for the remaining two outflows (II and I) only enhanced emission of H 2 and [S i] is detected, and ionized atomic fine structure lines are absent toward these outflow directions.With the low angular resolution of these older data often multiple physical entities are covered, for example, Hii regions, photon-dominated/photon-dissociation regions (PDRs), and outflows (e.g. in Orion IRc2 van Dishoeck et al. 1998).In addition, due to low filling factors toward distant HMSFRs such lines were often not detected.With JWST/MIRI MRS we are able to resolve individual components.In fact, both shocks and PDRs can contribute to many of the lines observed in this work.
The spatial morphology of the [Ne ii] line (Fig. 3) reveals extended emission within the full MIRI FoV, and is also detected toward the east at the arc-like structure that is a PDR illuminated by the nearby UCHii region.In PDRs the radiation field is dominated by FUV photons (6-13.6 eV) (e.g.Hollenbach & Tielens 1999).PDRs show emission by MIR H 2 transi- Fig. 7. JWST/MIRI MRS spectra of atomic line emission.In each panel, the spectra are presented extracted from the positions at the MIR source A/mmA2 (black, Table 3) and two shock locations (orange: S1 and green: S3, see Fig. 3).The bottom and top x-axis show the observed wavelength and velocity, respectively.The red vertical line indicates the wavelength of the line at a source velocity of −50.2 km s −1 .
tions, as well as neutral and single ionized atomic fine structure lines (van Dishoeck 2004, and references within).In the case of the IRAS 23385 PDR, we detect emission of H 2 and [Ne ii] in the arc highlighting the impact of the clustered mode of highmass star-formation including photons with energies >13.6 eV that can ionize Neon.
The [Ne ii] transition also traces both outflows and the inner disk toward low-mass protostars where these regions can be sufficiently resolved (Sacco et al. 2012).These authors find that the [Ne ii] emission is slightly blueshifted (2−12 km s −1 ) and argue that it may be caused by a disk wind.Glassgold et al. (2007) suggest that the emission could stem from the disk that is irradiated by stellar X-rays.[Ne ii] is commonly detected in Class II disks (Pascucci et al. 2022).While most atomic fine structure lines are indeed blueshifted by 50−100 km s −1 with respect to the source velocity, the [Ne ii] and [Ne iii] lines toward the position of source A/mmA2 are not blueshifted but more centered on the source velocity (Fig. 7), so the emission at source A/mmA2 could either stem from the protostar itself or a disk wind.
Considering the [Ne ii] and [Ne iii] peak toward the west of source A/mmA2 (Fig. 3), a disk wind is unlikely to explain the emission since it is ≈ 5 000 au offset from the MIR source.Another explanation could be the that we are tracing the UV and Xray irradiated cavity of the large scale northwest-southeast outflow (I) launched from source B. However, the emission is not perfectly aligned with that direction.Since most emission lines toward source A/mmA2 shown in Fig. 7  ] transitions, we perform aperture photometry with an aperture of 0 ′′ .9 toward both [Ne iii] peak positions (highlighted by grey dash-dotted circles in Fig. 3).Toward source A/mmA2 the line integrated flux density is 3.1 mJy for [Ne ii] and 0.23 mJy for [Ne iii], while to the western position it is 2.6 mJy and 0.19 mJy, respectively.While the flux density is higher at the position of source A/mmA2, the [Ne ii]/[Ne iii] flux density ratio at both positions is similar (13 and 14) suggesting that even at a projected distance of 5 000 au, the X-ray and UV-radiation from the protostar can efficiently funnel along the outflow.Simpson et al. (2012) find [Ne ii]/[Ne iii] ratios between 17 and 125 in a sample of massive protostars in the G333.2−0.4 giant molecular cloud and Rosenthal et al. (2000) find a ratio of 0.95 toward the brightest outflow knot in the Orion OMC-1 outflow.Toward the Orion IRc2 source that causes this outflow, a ratio of 1.4 is derived (van Dishoeck et al. 1998).In the Orion PDR a ratio of 1.3-1.9 is derived.Güdel et al. (2010) find a higher [Ne ii] luminosity in sources associated with jets compared to those without which would be consistent with our high ratio.These authors find a weak correlation with X-ray luminosity.When the full sample of the JOYS targets is observed, a statistical analysis of such properties can be conducted from low-to high-mass protostars and contributions for jets versus UV and/or X-rays can be disentangled.

Conclusions
In this study we compare the MIR properties of the warm gas in IRAS 23385 traced by JWST/MIRI MRS observations that are part of the JOYS project with the cold material traced by NOEMA observations at mm wavelengths.In total, four continuum sources are resolved that are all embedded within the dense cold envelope.Two sources are MIR bright, while the two bright mm cores remain undetected in the MIRI range between 5 and 28 µm.The spatial morphology of the atomic and molecular lines is investigated using integrated intensity maps.The gas temperature and column density of different components are estimated using H 2 MIR and CH 3 CN mm line emission and 3 mm continuum emission.Our conclusions are summarized as follows: 1. Combining the JWST/MIRI MRS and NOEMA observations allows us to disentangle several gas components, consisting of cold (<100 K), warm (200-400 K) and hot (≥1 000 K) material, traced by CH 3 CN and H 2 lines.While the warm and hot component, mostly caused by heated and/or shocked gas in the outflows, is extended, there is still a significant amount of cold material that allows for star formation to proceed.This is also highlighted by the H 2 column densities of the cold and warm components that are both on the order of 10 24 cm −2 .The column density of the hot component is in contrast significantly lower, on the order of 10 22 cm −2 .2. The H 2 0-0 S(1) to S( 7) transitions reveal C-type shocked and extended warm gas due to multiple bipolar outflows (I, II, and III).These are also evident in the atomic sulfur line [S i], stressing the importance of sulfur being released from the grains due to shocks.The outflows and their cavities can also be seen in many molecular lines at mm wavelengths (SiO, SO, 13 CS, H 13 CO + , H 2 CO, CH 3 OH).and [Ne iii] peaks at the source velocity.This suggests that in these cases, not only the outflow may cause the ionization and excitation of the line emission, but, for example, a disk wind or the disk atmosphere itself or an X-ray irradiated cavity, is the cause.However, the angular and spectral resolution are not sufficient to differentiate between these cases.5.The [Ne ii] emission is extended within the full MIRI FoV and also traces bright emission toward the edges of the MIRI FoV, a PDR region caused by the nearby UCHii regions (Fig. 1).
To fully understand the clustered properties of high-mass star formation it is important to conduct a multi-wavelength study tracing both the warm and the cold gas properties.With JWST we now have a new unprecedented view on the IR properties and diagnostics of protostars and their environments at high angular resolution and sensitivity.The line integrated intensities of the Gaussian fits are converted to upper state column densities divided by the statistical weight and when plotted against the upper energy level, the H 2 rotation temperature and column density can be estimated from the slope and intercept, respectively.Examples of the H 2 excitation diagrams are shown in Fig. 5 for sources mmA1 and B. Typically for star-forming regions, the profile is best described by two temperature components, which we refer to as the "warm" and "hot" component in contrast to the "cold" component traced by the NOEMA data.The fit results for the warm and hot component are summarized in Table A.2 for all continuum sources and shock positions.For comparison, we also present the results using extinction values of A K = 5 mag and A K = 3 mag.There are no significant differences for the derived column densities.A small effect is seen for the temperatures given that the S(3) transition that lies in the silicate absorption feature influences the two fitted slopes.
The temperature of the cold component is estimated using CH 3 CN line emission covered by the CORE+ project.In each spatial pixel, the CH 3 CN J = 5 − 4 K-ladder is modeled using XCLASS (Möller et al. 2017).An example toward the location of the mmA1 source is shown in Fig. A.2.In these examples, we trace for source mmA1 temperature components at 50 K, 400 K, and 1 400 K using the H 2 and CH 3 CN lines as thermometers.

Fig. 1 .
Fig. 1.Multi-wavelength overview of IRAS 23385.In color, the Herschel 250 µm (left), Herschel 70 µm (center), and Spitzer 24 µm (right) emission are shown in log-scale.The grey dashed squares show the field-of-view of the following panel.In the center and right panels, the green star, labeled as "mm", marks the 1 mm continuum peak position (Fig.2).The blue stars, labeled as "VLA1" and "VLA2", show the position of the two nearby UCHii regions(Molinari et al. 2002).The black rectangles in the right panel indicate the JWST/MIRI MRS field-of-view of the 4-pointing mosaic in ch1A (solid) and ch4C (dashed).

Fig. 2 .
Fig.2.MIR Continuum (left) and H 2 0-0 S(7) line emission (right) for IRAS 23385.In the left panel, the JWST/MIRI 5.2 µm continuum with emission > 5σ cont,5µm is presented in color.Grey contours show the NOEMA 1 mm continuum with contour levels at 5, 10, 20, 40, 80×σ cont,1mm .The dash-dotted white circles show the aperture in which the 5.2 µm flux density F 5.2µm was derived (Table3).In the bottom left corner, the JWST/MIRI 5.2 µm and NOEMA 1 mm angular resolution is indicated in black and grey, respectively.The mm (mmA1 and mmB, marked by green squares) and MIR (A/mmA2 and B) continuum sources are labeled in green.In the right panel, the line integrated intensity of the H 2 0-0 S(7) line with S /N > 5 is presented in color.The JWST/MIRI 5.2 µm continuum as shown in the left panel is presented as blue contours with contour levels at 5, 10, 15, 20, 25×σ cont,5µm .The angular resolution of the 5.2 µm continuum and H 2 0-0 S(7) line data is indicated in the bottom left and right, respectively.Several shock spots evident in the MIRI MRS line emission are marked by green crosses and labeled in green (Sect.3.2).The red and blue arrows indicate three bipolar outflows, labeled I, II, and III (as presented inBeuther et al. 2023).In both panels, the black bar indicates a spatial scale of 10 000 au at the assumed source distance of 4.9 kpc.
the JWST and NOEMA observations, respectively.Integration ranges were selected based on visual inspection of the emission.The integrated intensity maps of the emission lines are presented in Figs.3 (MIR lines) and 4 (mm lines) with a comparison to the JWST/MIRI 5.2 µm continuum, shock positions, and outflow directions.

[−[−Fig. 3 .
Fig.3.Integrated intensity maps of atomic and molecular lines detected with JWST/MIRI MRS.In color, the line integrated intensity is shown in a log-scale.The JWST/MIRI 5.2 µm continuum is presented as blue contours with contour levels at 5, 10, 15, 20, 25×σ cont,5µm .The two mm sources are indicated by green squares and several shock positions are highlighted by green crosses.In the bottom left and right corners, the angular resolution of the JWST/MIRI continuum and line data, respectively, is shown.In the top left panel, the bipolar outflows are indicated by red and blue dashed arrows.The shock locations (Sect.3.2) and continuum sources are labeled in the top right and center left panel, respectively.In the [Ne ii] and [Ne iii] panels, the dash-dotted grey circles show the aperture (0 ′′ .9) in which the flux density was derived (Sect.4).

Fig. 4 .
Fig.4.Integrated intensity maps molecular lines detected with NOEMA plotted over the same area as the JWST data in Figs.2 and 3.In color, the line integrated intensity is shown.The JWST/MIRI 5.2 µm continuum is presented as blue contours with contour levels at 5, 10, 15, 20, 25×σ cont,5µm .The two mm sources are indicated by green squares and several positions are highlighted by green crosses.In the bottom left and right corners, the synthesized beam of the NOEMA continuum and line data, respectively, is shown.In the top left panel, the bipolar outflows are indicated by red and blue dashed arrows(Beuther et al. 2023).The shock locations (Sect.3.2) and continuum sources are labeled in the top central and center left panel, respectively.

Fig. 5 .
Fig. 5. Example of the H 2 diagram analysis with pdrtpy of source mmA1 (top) and source B (bottom).The observed data is shown in black and the two component fit is shown by red and blue dots, corresponding to the warm and hot component, respectively and the total fit is indicated by a green line.

Fig. 6 .
Fig. 6.Temperature and H 2 column density of the gas components in IRAS 23385 derived using CH 3 CN and H 2 as a diagnostic tool (Sects.3.3 and 3.4).In the top and bottom panels, the rotation temperature and H 2 column density map, respectively, of the cold (left), warm (center), and hot (right) components are shown in color.The angular resolution of the line data is indicated by a purple ellipse in the bottom right corner.The JWST/MIRI 5.2 µm continuum is presented by blue contours with contour levels at 5, 10, 15, 20, 25×σ cont,5µm and the angular resolution is highlighted by a black ellipse in the bottom left corner.The NOEMA 3 mm continuum (top and bottom left panels) is highlighted by grey contours with levels at 5, 10, 20, 40, 80×σ cont,3mm and the synthesized beam is highlighted by a grey ellipse in the bottom left corner.All continuum sources are labeled in green in the bottom left panel and the mm continuum sources are marked by green squares.Several shock positions are indicated by green crosses (Sect.3.2).
. 7, example spectra of the atomic transitions extracted from the positions of the continuum source A/mmA2 and the shock positions S1 and S3 are shown.The [S i] line peaks at the source velocity toward source A/mmA2 and shock position S1.The shock position S3 caused by the redshifted part of outflow II is as expected redshifted.If detected, [Ni ii], [Ar ii], [Fe ii] lines are blueshifted compared to the source velocity tracing the blueshifted outflow components.The [Ne ii] and [Ne iii] lines are further discussed in Sect.4.3.Rosenthal et al. (2000) studied the MIR spectrum of one of the brightest H 2 shock position of the Orion OMC-1 outflow.In comparison with the IRAS 23385 outflows, these authors find emission from double or triple ionized fine structure lines by [Fe iii], [S iii], [P iii], [S iv], and [Ar iii], while in IRAS 23385 only [Ne iii]
are blueshifted compared to the source velocity, while [Ne ii] and [Ne iii] peak at the source velocity, this suggests that [Ne ii] and [Ne iii] are tracing a region closer to the protostar.It remains inconclusive what causes the [Ne ii], [Ne iii], and faint [Ar ii] emission to the west of source A/mmA2 (Fig. 3).High-mass (proto)stars have enough UV-intensity to ionize [Ne ii] and [Ne iii].To estimate the line integrated flux density of both [Ne ii] and [Ne iii

3 .
An energetic outflow (III) emerges probably from source A/mmA2 where the blueshifted outflow reveals emission knots by [Fe ii] and [Ni ii].Thus, the outflow seems to have undergone multiple outbursts in the past causing dissociative J-type shocks.4. The MIR source A/mmA2 also shows enhanced emission toward the central location of the protostar by [Fe ii], [Ni ii], as well as [Ne ii], [Ne iii], and [Ar ii].While the emission is blueshifted by 50−100 km s −1 compared to the source velocity for most of these species, the emission from [Ne ii], Fig. A.1.Examples of the observed and extinction-corrected H 2 lines (black) and Gaussian fit (red line) of source mmA1 (top) and source B (bottom).The vertical grey line indicates the wavelength of each transition.
Fig. A.2. Example of CH 3 CN line modeling with XCLASS.In black, the observed spectrum toward mmA1 is shown and in red the best-fit model derived by XCLASS.The modeled transitions of CH 3 CN (Table 2) are indicated by blue vertical bars.

Table 2 .
Molecular lines covered by NOEMA at 1 mm and 3 mm analyzed in this work.
3pdrtpy is developed by Marc Pound and Mark Wolfire.This project is supported by NASA Astrophysics Data Analysis Program grant 80NSSC19K0573; from JWST-ERS program ID 1288 provided through grants from the STScI under NASA contract NAS5-03127; and from the SOFIA C+ Legacy Project through a grant from NASA through award #208378 issued by USRA.
Notes.Transition properties are listed inTable 1. Article number, page 16 of 17 C. Gieser et al.: JOYS: Disentangling the warm and cold material in the high-mass IRAS 23385+6053 cluster Table A.2. Fit results from the H 2 excitation diagram analysis with pdrtpy (Sect.3.3) with a warm and hot component, adopting values of A K = 7, 5, and 3 mag.