# JOYS+: Mid-infrared detection of gas-phase SO 2 emission in a low-mass protostar

M. van Gelder, M. Ressler, E. van Dishoeck, P. Nazari, Benoît Tabone, J.

Black, Ł. Tychoniec, L. Francis, M. Barsony, H. Beuther, et al.

► To cite this version:

M. van Gelder, M. Ressler, E. van Dishoeck, P. Nazari, Benoît Tabone, et al. JOYS+: Mid-infrared detection of gas-phase SO 2 emission in a low-mass protostar. Astronomy & Astrophysics - A&A, 2024, 682, pp.A78. 10.1051/0004-6361/202348118. hal-05412365

HAL Id: hal-05412365

https://hal.science/hal-05412365v1

Submitted on 12 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# JOYS+: Mid-infrared detection of gas-phase SO₂ emission in a low-mass protostar

## The case of NGC 1333 IRAS 2A: Hot core or accretion shock?

M. L. van Gelder¹, M. E. Ressler², E. F. van Dishoeck¹,³, P. Nazari¹, B. Tabone⁴, J. H. Black⁵, Ł. Tychoniec⁶, L. Francis¹, M. Barsony⁷, H. Beuther⁸, A. Caratti o Garatti⁹, Y. Chen¹, C. Gieser³, V. J. M. le Gouellec¹⁰,*, P. J. Kavanagh¹¹, P. D. Klaassen¹², B. W. P. Lew¹³, H. Linnartz¹⁴, L. Majumdar¹⁵,¹⁶, G. Perotti⁸, and W. R. M. Rocha¹,¹⁴

*(Affiliations can be found after the references)*

Received 29 September 2023 / Accepted 24 November 2023

## ABSTRACT

**Context.** Thanks to the Mid-InfraRed Instrument (MIRI) on the James Webb Space Telescope (JWST), our ability to observe the star formation process in the infrared has greatly improved. Due to its unprecedented spatial and spectral resolution and sensitivity in the mid-infrared, JWST/MIRI can see through highly extincted protostellar envelopes and probe the warm inner regions. An abundant molecule in these warm inner regions is SO₂, which is a common tracer of both outflow and accretion shocks as well as hot core chemistry.

**Aims.** This paper presents the first mid-infrared detection of gaseous SO₂ emission in an embedded low-mass protostellar system rich in complex molecules and aims to determine the physical origin of the SO₂ emission.

**Methods.** JWST/MIRI observations taken with the Medium Resolution Spectrometer (MRS) of the low-mass protostellar binary NGC 1333 IRAS 2A in the JWST Observations of Young protoStars (JOYS+) program are presented. The observations reveal emission from the SO₂ v₃ asymmetric stretching mode at 7.35 µm. Using simple slab models and assuming local thermodynamic equilibrium (LTE), we derived the rotational temperature and total number of SO₂ molecules. We then compared the results to those derived from high-angular-resolution SO₂ data on the same scales (~50–100 au) obtained with the Atacama Large Millimeter/submillimeter Array (ALMA).

**Results.** The SO₂ emission from the v₃ band is predominantly located on ~50–100 au scales around the mid-infrared continuum peak of the main component of the binary, IRAS 2A1. A rotational temperature of 92 ± 8 K is derived from the v₃ lines. This is in good agreement with the rotational temperature derived from pure rotational lines in the vibrational ground state (i.e., v = 0) with ALMA (104 ± 5 K), which are extended over similar scales. However, the emission of the v₃ lines in the MIRI-MRS spectrum is not in LTE given that the total number of molecules predicted by a LTE model is found to be a factor of 2 × 10⁴ higher than what is derived for the v = 0 state from the ALMA data. This difference can be explained by a vibrational temperature that is ~100 K higher than the derived rotational temperature of the v = 0 state: T_vib ~ 200 K versus T_rot = 104 ± 5 K. The brightness temperature derived from the continuum around the v₃ band (~7.35 µm) of SO₂ is ~180 K, which confirms that the v₃ = 1 level is not collisionally populated but rather infrared-pumped by scattered radiation. This is also consistent with the non-detection of the v₂ bending mode at 18–20 µm. The similar rotational temperature derived from the MIRI-MRS and ALMA data implies that they are in fact tracing the same molecular gas. The inferred abundance of SO₂, determined using the LTE fit to the lines of the vibrational ground state in the ALMA data, is 1.0 ± 0.3 × 10⁻⁸ with respect to H₂, which is on the lower side compared to interstellar and cometary ices (10⁻⁸–10⁻⁷).

**Conclusions.** Given the rotational temperature, the extent of the emission (~100 au in radius), and the narrow line widths in the ALMA data (~3.5 km s⁻¹), the SO₂ in IRAS 2A likely originates from ice sublimation in the central hot core around the protostar rather than from an accretion shock at the disk-envelope boundary. Furthermore, this paper shows the importance of radiative pumping and of combining JWST observations with those from millimeter interferometers such as ALMA to probe the physics on disk scales and to infer molecular abundances.

**Key words.** astrochemistry – stars: formation – stars: low-mass – stars: protostars – ISM: molecules – ISM: individual objects: NGC 1333 IRAS 2A

## 1. Introduction

The embedded protostellar phase of star formation is very rich in terms of chemistry. The James Webb Space Telescope (JWST) provides unique new opportunities for probing these deeply embedded protostellar sources (Yang et al. 2022; Harsono et al. 2023; van Dishoeck et al. 2023; Beuther et al. 2023; Gieser et al. 2023). An interesting element to study is sulfur, as the

total volatile sulfur budget in protostars appears to be depleted by more than two orders of magnitude with respect to the diffuse clouds (Ruffle et al. 1999). The sulfur likely resides in unobservable refractory reservoirs, sulfur allotropes, or FeS inclusions (Woods et al. 2015; Kama et al. 2019). Even among the volatile species, the sulfur budget remains poorly constrained (e.g., Drozdovskaya et al. 2018; Navarro-Almaida et al. 2020; Kushwahaa et al. 2023). It is therefore important to constrain the physical conditions in which the volatile sulfur-bearing species reside in order to understand their chemistry and constrain the

<footnote>* NASA Postdoctoral Program Fellow.</footnote>
main sulfur reservoirs. JWST has proven that it can detect sulfur-
bearing species both in interstellar ices (e.g., SO₂ and OCS;
Yang et al. 2022; McClure et al. 2023; Rocha et al. 2024) and
in exoplanetary atmospheres (e.g., SO₂; Tsai et al. 2023). Here,
we present one of the first detected medium-resolution mid-
infrared (MIR) spectra of a low-mass Class 0 protostellar system,
NGC 1333 IRAS 2A, containing the first detection of gaseous
SO₂ in emission at MIR wavelengths.

One of the most frequently detected sulfur-bearing species toward low-mass protostellar systems in millimeter observations is SO₂ (e.g., Artur de la Villarmois et al. 2023). It is a shock tracer that is often present in outflow and jet shocks (e.g., Blake et al. 1987; Codella et al. 2014; Taquet et al. 2020; Tychoniec et al. 2021), where it is either sputtered from icy dust grains or formed through gas-phase chemistry (e.g., Pineau des Forêts et al. 1993). SO₂ is also often observed on smaller scales in accretion shocks at the disk-envelope boundary (e.g., Sakai et al. 2014; Oya et al. 2019; Artur de la Villarmois et al. 2019, 2022; Garufi et al. 2022), where it is either formed through gas-phase chemistry or thermally sublimated from the ices (Miura et al. 2017; van Gelder et al. 2021). Furthermore, it is also a good tracer of disk winds (Tabone et al. 2017) and hot core regions, where the bulk of the ices are sublimating (e.g., Drozdovskaya et al. 2018; Codella et al. 2021). Most of these studies are based on submillimeter observations with interferometers such as the Atacama Large Millimeter/submillimeter Array (ALMA) or the NOrthern Extended Millimeter Array (NOEMA), which can trace the pure rotational transitions of the vibrational ground state of SO₂ but are not able to detect ro-vibrational transitions.

These ro-vibrational lines are best traced at MIR wave-
lengths (i.e., ~5–30 µm). SO₂ has three fundamental vibrational
modes: the v₁ symmetrical stretching mode around 8.5–9 µm,
the v₂ bending mode around 18–20 µm, and the v₃ asymmetri-
cal stretching mode around 7.2–7.4 µm (Briggs 1970). At MIR
wavelengths, gaseous SO₂ has thus far only been observed in
absorption toward high-mass protostellar systems (Keane et al.
2001; Dungee et al. 2018; Nickerson et al. 2023). These stud-
ies have shown that the SO₂ in these high-mass systems resides
at typical temperatures of ~100–300 K, although higher tem-
peratures of up to 700 K have also been reported (Keane et al.
2001). The average abundances with respect to H₂ are >10⁻⁷,
which is consistent with the SO₂ abundance of interstellar ices
(Boogert et al. 1997, 2015; Zasowski et al. 2009; McClure et al.
2023; Rocha et al. 2024) and cometary ices (Altwegg et al. 2019;
Rubin et al. 2019), suggesting that gaseous SO₂ may originate
from ice sublimation in the inner hot cores of these high-mass
protostellar systems. A MIR detection of gaseous SO₂, either in
absorption or emission, toward low-mass protostellar systems is,
to the best of our knowledge, still lacking.

Most studies at submillimeter wavelengths assume local thermodynamic equilibrium (LTE) when deriving the physical conditions of SO₂ (i.e., column density and excitation temperature), which is a good approximation given that most pure rotational levels can be collisionally populated at typical inner envelope densities (10⁶–10⁸ cm⁻³) and can therefore be characterized by a single excitation temperature. However, when studying ro-vibrational lines at MIR wavelengths, it is important to be aware of non-LTE effects. The critical densities of these ro-vibrational transitions are typically >10¹⁰ cm⁻³ (e.g., for HCN and CO₂; Bruderer et al. 2015; Bosman et al. 2017), meaning that vibrationally excited levels are only collisionally populated in the inner ≤1 au around low-mass protostars. Furthermore, molecules can be pumped into a vibrationally excited state by a strong infrared radiation field which boosts the line

fluxes of MIR transitions far above what is expected from col-
lisional excitation (e.g., Boonman et al. 2003; Sonnentrucker
et al. 2007). For SO₂, no collisional rate coefficients are available
for its ro-vibrational transitions in the MIR, therefore prevent-
ing a full non-LTE analysis, but its MIR transitions likely have
high critical densities similar to those of HCN and CO₂. How-
ever, the importance of non-LTE effects can still be constrained
through the comparison of ro-vibrational transitions detected at
MIR wavelengths to pure rotational transitions in the vibrational
ground state measured at submillimeter wavelengths.

In this paper we present JWST/Mid-InfraRed Instrument (MIRI; Rieke et al. 2015; Wright et al. 2015, 2023) observations from the JWST Observations of Young protoStars (JOYS+) program, providing the first MIR detection of SO₂ in emission toward the low-mass protostellar system NGC 1333 IRAS 2A (hereafter IRAS 2A), and we compare it to the results of high-angular-resolution ALMA data of the same region. IRAS 2A is a binary Class 0 system with a separation between IRAS 2A1 (VLA1) and IRAS 2A2 (VLA2) of ~0.6" (i.e., ~180 au; Tobin et al. 2015, 2016; Jørgensen et al. 2022) located in the Perseus star-forming region at a distance of about 293 pc (Ortiz-León et al. 2018). It hosts one of the most famous hot corinos (e.g., Jørgensen et al. 2005; Bottinelli et al. 2007; Maury et al. 2014; Taquet et al. 2015) and drives two powerful almost perpendicu- lar large-scale outflows (e.g., Arce et al. 2010; Tobin et al. 2015; Taquet et al. 2020). Recently, Jørgensen et al. (2022) showed that the binary interaction results in a misalignment in the out- flow and infalling streamers around IRAS 2A1. Furthermore, IRAS 2A was part of the *Spitzer* Cores to Disks (c2d) survey, which revealed the main components of the ices (Boogert et al. 2008; Öberg et al. 2011), but the Infrared Spectrograph (IRS) of *Spitzer* only had a spectral resolving power of *R* = 60–120 at the critical wavelength range of 5–10 µm and a large aperture (>5"). Here, we present MIRI observations taken with the Medium Res- olution Spectrometer (MRS; Wells et al. 2015; Labiano et al. 2021; Argyriou et al. 2023; Jones et al. 2023) with a spectral resolving power of *R* ~ 3500 at 5–10 µm and subarcsecond res- olution (~0.2–0.4"). This is one of the first JWST/MIRI-MRS spectra of a low-mass Class 0 protostellar system to be detected (see also, e.g., IRAS 15398-3359; Yang et al. 2022).

This paper is structured as follows. In Sect. 2 the JWST/MIRI and ALMA observations are described together with the method of fitting the SO₂ emission. The results are presented in Sect. 3 and discussed in Sect. 4. Our main conclusions are summarized in Sect. 5.

## 2. Observations and analysis

2.1. *Observations*

2.1.1. MIRI-MRS

The MIRI-MRS observations were carried out as part of the Cycle 1 guaranteed time observations (GTO) program 1236 (PI: M. E. Ressler) on January 9, 2023, with a two-point dither pattern optimized for extended sources. A dedicated background observation was also taken with two-point dither pattern to allow for a proper subtraction of the telescope background and detector artifacts. In both cases, the FASTR1 readout mode was used with all three gratings (A, B, and C) in all four channels, providing the full wavelength coverage of MIRI-MRS (4.9–28.6 µm). The integration time in each grating was 111 s, resulting in a total integration time of 333 s. All observations included simultaneous off-source imaging using the F1500W filter.
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Spectrum (top panel) and continuum images at various wavelengths (bottom row) of IRAS 2A. In the top panel, four insets are presented: around the SO₂ v₃ asymmetric stretching mode (~7.35 µm), SO₂ v₁ symmetric stretching mode (~8.7 µm), SO₂ v₂ bending mode (~19 µm), and CO₂ v₂ bending mode (~15 µm). In the latter, the emission originates mostly from the outflow (see Fig. B.1). The 12 sub-bands of the MIRI-MRS spectrum are not stitched and show minor offsets in overlapping wavelengths (see, e.g., the inset on the SO₂ v₁ symmetric stretching mode). The bottom row shows, from left to right, the dust continuum around 5.3 µm, 8.1 µm, 12.5 µm, and 19.3 µm; a sqrt stretch has been used to enhance fainter features without over-saturating bright emission. The open white circle indicates the size of the aperture from which the spectrum was extracted (i.e., 1.1", 1.5", 2.1", and 3.0" in diameter at 5.3 µm, 8.1 µm, 12.5 µm, and 19.3 µm, respectively), which increases as a function of wavelength with the increase in the size of the PSF. A scale bar is displayed in the bottom left of each panel, and the size of the PSF is presented as the filled white circle in the bottom right.</figcaption></figure>

The observations were processed through all three stages of the JWST calibration pipeline version 1.11.0 (Bushouse et al. 2023) using the reference context *jwst_1097.pmap* of the JWST Calibration Reference Data System (CRDS; Greenfield & Miller 2016). The raw data were first processed through the *Detector1Pipeline* using the default settings. Following this, the *Spec2Pipeline* was performed, including the correction for fringes with the fringe flat for extended sources (Mueller et al., in prep.) and applying the detector level residual fringe correction (Kavanagh et al., in prep.). Furthermore, the background was also subtracted at this step using the rate files of the dedicated background. A bad-pixel routine was applied to the *Spec2Pipeline* products outside of the default MIRI-MRS pipeline using the Vortex Image Processing (VIP) package (Christiaens et al. 2023). The data were further processed using the *Spec3Pipeline* with both the master background and outlier rejection routines switched off. The background was already subtracted in the *Spec2Pipeline* and the outlier rejection routine was skipped because it did not significantly improve the quality of the data. Moreover, the data cubes created with the *Spec3Pipeline* were created for each band of each channel separately.

The main component of the binary, IRAS 2A1 (VLA1; Tobin et al. 2015), is clearly detected across all wavelengths whereas the companion protostar, IRAS 2A2 (VLA2), is not detected

(see Fig. 1). The reason for the non-detection could be because it is about an order of magnitude less massive and less luminous than IRAS 2A1 (Tobin et al. 2016, 2018), or because it is more embedded. Therefore, only one spectrum was manually extracted from the peak of the MIR continuum (RA (J2000) 03h28m55.57s, Dec (J2000) 31d14m36.76s) at 5.5 µm (channel 1A; see Fig. 1) and is assumed to only contain emission related to IRAS 2A1 (hereafter IRAS 2A). The diameter of the aperture was set to 4 × FWHM<sub>PSF</sub>, with FWHM<sub>PSF</sub> the empirically derived full width at half maximum of the MIRI-MRS point spread function (FWHM<sub>PSF</sub> = 0.033(λ/µm) + 0.106", i.e., 0.35" at 7.35 µm; Law et al. 2023), in order to capture as much emission as possible without adding too much noise. An additional spectrum level residual fringe correction was applied primarily to remove the high-frequency dichroic fringes in channel 3 and 4 (Kavanagh et al., in prep.). No additional spectral stitching was applied between the bands.

### 2.1.2. ALMA data

The ALMA data analyzed in this paper are taken from program 2021.1.01578.S (PI: B. Tabone), which targeted several Class 0 protostars in Perseus in Band 7, covering several transitions of SO₂ and ³⁴SO₂. These data were taken in an extended configuration (C-6, θ<sub>beam</sub> = 0.12" × 0.09") and in a more compact
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Subtracting the local continuum around the SO₂ v₃ band. Top: spectrum of IRAS 2A around 7.35 µm showing the gas-phase SO₂ emission superimposed at various ice absorption bands. The red line shows the estimated local continuum based on a fourth-order polynomial fit through the red dots. Middle: spectrum of IRAS 2A with the local continuum subtracted. The 7.25 µm complex organics ice absorption feature is clearly present. A simple Gaussian fit to this ice absorption feature based on fitting the spectrum in the shaded region is shown in green. Bottom: final continuum- and ice-feature-subtracted spectrum.</figcaption></figure>

configuration (C-3, $\theta_{\text{beam}} = 0.58'' \times 0.34''$) to include larger-scale emission. All spectral windows have a velocity resolution of 0.22 km s⁻<fnref n="1" /> except for two windows that have a resolution of 0.44 km s⁻<fnref n="1" /> and one continuum window at 0.87 km s⁻<fnref n="1" />. The data were pipeline calibrated and imaged with the Common Astronomy Software Applications<fnref n="1" /> (CASA; McMullin et al. 2007) version 6.4.1.12. The continuum was subtracted with the uvcontsub task using carefully selected line-free channels. Following this, the two configurations were combined using the tclean task for a Briggs weighting of 0.5 with a circular mask with a radius of 2'' centered on the main continuum peak. The synthesized beam of the final data cubes is $\theta_{\text{beam}} = 0.13'' \times 0.10''.$ In order for direct comparison to the MIRI data, spectra were extracted from a 1.4'' diameter aperture centered on the continuum position of IRAS 2A1. The noise level of the extracted spectrum is 0.015 Jy, and a flux calibration uncertainty of 5% was assumed.

## 2.2. Continuum subtraction around the SO₂ v₃ band

SO₂ is an asymmetric rotor and has three fundamental vibrational modes: a symmetrical stretching mode ($v_1 = 1151$ cm⁻<fnref n="1" />, $\lambda \sim 8.5–9.0$ µm), a bending mode around ($v_2 = 518$ cm⁻<fnref n="1" />, $\lambda \sim 18–20$ µm), and an asymmetrical stretching mode ($v_3 = 1362$ cm⁻<fnref n="1" />, $\lambda \sim 7.2–7.4$ µm; Briggs 1970; Person & Zerbi 1982). Around 7.35 µm, clear molecular emission originating

from the v₃ band is present (see the inset in Fig. 1), but both the v₁ and v₂ bands are not detected. The analysis of the MIRI-MRS data will therefore be focused on the v₃ band and the absence of the v₁ and v₂ bands will be further discussed in Sect. 4.1.1.

In order to fit the SO₂ emission in the v₃ band, the local continuum had to be subtracted. However, the spectral region surrounding the emission features is dominated by absorption of ices such as the 7.25 and 7.4 µm ice bands that are typically ascribed to complex organics including ethanol (CH₃CH₂OH) and acetaldehyde (CH₃CHO) or the formate ion (HCOO⁻; e.g., Schutte et al. 1999; Öberg et al. 2011; Boogert et al. 2015; Terwisscha van Scheltinga et al. 2018; Rocha et al. 2024). A fourth-order polynomial was fitted through obvious line-free channels (i.e., 7.015, 7.113, 7.162, 7.442, 7.489, 7.545, and 7.635 µm; see the top panel of Fig. 2) selected outside of the SO₂ emission range and the 7.25 and 7.4 µm ice bands to estimate the local continuum. Following the subtraction of this local continuum, a clear ice absorption feature remains around 7.25 µm. The generally equally strong 7.4 µm ice feature appears to be absent, but this likely originates from the superposition with the P-branch of the SO₂ emission (see Sect. 3). To estimate the strength of the 7.25 µm ice absorption band, a simple Gaussian model was fitted to the part of the spectrum where no clear molecular emission features are present (7.1–7.26 µm; green shaded area in Fig. 2). Using this Gaussian fit, the contribution of the 7.25 µm ice feature was subtracted, providing a spectrum with only the contribution of the molecular emission (Fig. 2 bottom).

<footnote n="1">¹ https://casa.nrao.edu/</footnote>
## 2.3. LTE slab model fitting

The molecular emission features in the continuum-subtracted spectra were fitted with simple slab models assuming LTE. This is a similar approach to what is applied to other low- and high-mass sources (Carr & Najita 2008; Salyk et al. 2011; Tabone et al. 2023; Grant et al. 2023; Perotti et al. 2023; Francis et al. 2024). It is important to note here that the assumption of LTE is likely not valid for the MIR ro-vibrational transitions of SO₂ detected in IRAS 2A. Nevertheless, LTE models can still constrain physical quantities such as the rotational temperature within a vibrationally excited state. The implications of the assumption of LTE will be further discussed in Sect. 4.1.

In the slab models, the emission is assumed to arise from a slab of gas with an excitation temperature $T_{ex}$ (which in the case of LTE is approximately equal to the kinetic temperature), a column density $N$, and an emitting area equal to $\pi R^2$. The latter was parameterized as a circular emitting area with a radius $R$, but could in reality originate from any shaped region with an area equal to $\pi R^2$. The intrinsic velocity broadening $\Delta V$ was fixed to $3.5 \text{ km s}^{-1}$ based on ALMA data covering SO₂ (see Sect. 3.2). It is important to note that in case of optically thin emission, the derived column densities are independent of $\Delta V$, whereas in case of optically thick emission the column density scales with $1/\Delta V$ (Tabone et al. 2023). Given that SO₂ has many ro-vibrational lines close to each other in wavelength in the MIR, it was important to take line overlap into account in order to derive an accurate column density and excitation temperature. Moreover, the spectral resolving power of MIRI-MRS around $7.35 \text{ µm}$ is about $\lambda/\Delta\lambda \sim 3500$ (Labiano et al. 2021; Jones et al. 2023), corresponding to a velocity resolution of $\sim 85 \text{ km s}^{-1}$, which means that the lines are not spectrally resolved and could be blended with other lines (of, e.g., SO₂ itself, H₂O) in the MIRI data.

The spectroscopic information (i.e., line wavelengths, Einstein $A_{ij}$, upper energy level $E_{up}$, degeneracy $g_{up}$) of the MIR ro-vibrational transitions of SO₂ was taken from the HITRAN database<fnref n="2" /> (Gordon et al. 2022) and converted to the Leiden Atomic and Molecular Database (LAMDA) format (van der Tak et al. 2020) in order to make it compatible with the slab model code. The partition function of SO₂ was calculated with the TIPS_2021_PYTHON code provided by the HITRAN database.

The best-fitting $N$ and $T_{ex}$ were computed by creating a large grid covering $10^{14}-10^{19} \text{ cm}^{-2}$ in steps of 0.05 in $\log_{10}$ and 50–250 K in steps of 1 K, respectively. Any higher column densities or temperatures were excluded by visual inspection of the spectrum. For each grid point, the LTE spectrum of SO₂ was calculated at a spectral resolving power of $R = 3500$ and the best-fitting emitting area was computed by minimizing the $\chi^2$ (see Appendix C of Grant et al. 2023, for more details). Only selected channels (i.e.,7.23–7.255, 7.258–7.268, and 7.2755–7.3535 µm) that do not include obvious emission or absorption features related to other molecules (i.e., H₂O) or artifacts are taken into account. Specifically, all wavelengths longer than $>7.3535 \text{ µm}$ were excluded because of the contribution of the 7.4 µm ice feature to the flux. Similarly, the noise level was estimated from line-free channels (7.06–7.13, 7.195–7.225, 7.477.50 µm) to be about 0.59 mJy, and a flux calibration uncertainty of 5% was assumed (Argyriou et al. 2023). The LTE models were corrected for an absolute extinction of $A_V = 55 \text{ mag}$; using the modified version of the extinction law of McClure (2009) introduced in Appendix C. The effect of

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Integrated intensity map of SO₂ observed with MIRI-MRS in color. The image is integrated over the *Q*-branch of the $v_3$ band: [7.34,7.35] µm and shown using a sqrt stretch to enhance fainter emission. The continuum around 7.35 µm is overlaid in white contours and peaks on the primary component, IRAS 2A1, extending somewhat toward the south in the direction of the scattered light. The positions of IRAS 2A1 and IRAS 2A2 are depicted with the white stars. A white scale bar is displayed in the bottom left, and the size of the PSF is presented as the filled white circle in the bottom right. The directions of the two outflows originating from IRAS 2A1 and IRAS 2A2 are indicated with the colored arrows (Tobin et al. 2015).</figcaption></figure>

differential extinction caused by the 7.25 µm and 7.4 µm ice absorption on the molecular emission is assumed to be negligible. The best-fit parameters ($N$, $T_{ex}$, and $R$) and their uncertainties were derived from the grid by minimizing the $\chi^2$.

# 3. Results

## 3.1. MIRI-MRS

The full spectrum of IRAS 2A is presented in the top panel of Fig. 1. In the bottom panels of Fig. 1, the continuum maps are presented at four increasing wavelengths, each obtained at a different MIRI-MRS channel. At the shortest wavelengths (~5.3 µm; left panel), the continuum clearly shows scattered light from the blue-shifted outflow cavity extending toward the south. This scattered light is visible up to wavelengths of ~8 µm (middle-left panel), but is no longer present at longer wavelengths (>8.5 µm; two panels on the right). Only the main component of the binary, IRAS 2A1 (VLA1; Tobin et al. 2015), is detected whereas the companion protostar, IRAS 2A2 (VLA2), is not detected. The non-detection of IRAS 2A could originate from it being an order of magnitude less massive and luminous (Tobin et al. 2016, 2018), or because it is more embedded.

The integrated intensity map over the SO₂ $v_3$ *Q*-branch in the (continuum-subtracted) MIRI-MRS data is presented in Fig. 3.

<footnote n="2">² https://hitran.org/</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Continuum-subtracted spectrum (blue) overlaid with the best-fit SO₂ model (orange). The orange bars in the bottom indicate the wavelength ranges that were included in the fit, and the black bars show the wavelength ranges over which σ is calculated. The *P*-branch lines around 7.4 µm appear to be overfitted since the 7.4 µm ice absorption could not be disentangled from the gas-phase emission.</figcaption></figure>

The SO₂ emission is mostly located around the continuum peak of IRAS 2A1. After deconvolution with the size of the PSF (Law et al. 2023), the extent of SO₂ emission is 0.5" × 0.3" in diameter (i.e., ~150 × 90 au). This very similar to the extent of H₂O emission lines that are likely tracing the inner hot corino (see Fig. B.1). Some very weak extended SO₂ emission (on ~500 au scales) is also present in the direction of the blue-shifted lobe of the north-south outflow (Tobin et al. 2015; Jørgensen et al. 2022). However, the extended SO₂ emission could also be related to IRAS 2A2, which is located about ~0.6" (i.e., ~180 au) from IRAS 2A1 to the southwest (Tobin et al. 2015, 2016, 2018; Jørgensen et al. 2022).

The SO₂ emission is clearly tracing a more compact component than the CO₂ emission around 15 µm, which peaks more toward the blue-shifted part of the north-south outflow than toward the central protostellar component (see Fig. B.1). Similarly, the lower *S*-branch lines of H₂ (e.g., S(1)) are also peaking in the blue-shifted outflow, but higher *S*-branch lines (e.g., S(7)) originating from higher *E*<sub>up</sub> levels are mostly peaking on the continuum position (see Fig. B.1). Furthermore, several ice absorption features, including those of SO₂ ice, are detected toward the bright continuum source and have been analyzed by Rocha et al. (2024). Any gas emission and absorption features other than SO₂ in this source, including the CO₂ emission in the outflow, will also be presented in future papers.

The best-fit LTE slab model to the $v_3$ band is presented in Fig. 4 overlaid on the continuum-subtracted spectrum. The fit slightly overestimates the *Q*-branch around 7.35 µm but provides a very reasonable fit to the *R*-branch around 7.3 µm. The *P*-branch lines around 7.4 µm were not included in the fit because the contribution of the 7.4 µm ice feature, which decreases the measured fluxes of the *P*-branch lines, could not be disentangled from the gas-phase emission of SO₂ a priori. Armed with a good model for the *Q*- and *R*-branches, Fig. B.2 clearly shows that the local continuum was estimated too high in the 7.4 µm region, proving that in fact the 7.4 µm ice absorption feature is present. The emission of the SO₂ *P*-branch is almost equally strong as the absorption caused by the ices.

The $\chi^2$ map is shown in Fig. 5. Based on the $\chi^2$ map, the best-fitting excitation temperature can be accurately constrained (with a $1\sigma$ error) at $T_{ex} = 92 \pm 8$ K. Any lower excitation temperature

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** $\chi^2$ map of the LTE slab model fit of SO₂. The $\chi^2$ values are inversely normalized by the minimum $\chi^2$ and shown in color. The 1, 2, and $3\sigma$ confidence intervals are presented as the black contours. The best-fitting emitting radius, *R* (in au), is represented by the white contours and is calculated for each grid point by minimizing the $\chi^2$ assuming a circular emitting area of $\pi R^2$. The $\chi^2$ map clearly reveals that the SO₂ emission is optically thin with an excitation temperature of $92 \pm 8$ K.</figcaption></figure>

results in a too strong and too narrow *Q*-branch compared to the *R*-branch, whereas higher excitation temperatures result in a too broad *Q*-branch and also start overshooting the *R*-branch. Furthermore, the SO₂ emission is optically thin since optically thick emission would give a lower ratio between the measured *R*- and *Q*-branch fluxes. The column density and emitting area are therefore completely degenerate with each other (i.e., the vertical profile of the contours in Fig. 5). Optically thick emission would appear as a more banana-shaped profile in Fig. 5 and allow both the column density and emitting area to be constrained (see, e.g., Pontoppidan et al. 2002; Salyk et al. 2011; Grant et al. 2023; Tabone et al. 2023). Given the degeneracy, no accurate column density of SO₂ can be derived, but the total number of molecules, $N = N\pi R^2$, can be determined: $N = 4.5 \pm 4.0 \times 10^{51}$ molecules.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Integrated intensity maps of the SO₂ 82,6–71,7 (left, *E*<sub>up</sub> = 43 K), SO₂ 164,12–163,13 (middle, *E*<sub>up</sub> = 164 K), and ³⁴SO₂ 174,14–173,15 (right, *E*<sub>up</sub> = 178 K) transitions in color. The images are integrated over [-2,2] km s⁻¹ with respect to the V<sub>1sr</sub> of 6.7 km s⁻¹. The white vertical bar in the colorbar on top of each image indicates the 3σ threshold. The 0.875 mm continuum is overlaid as black contours. The main continuum peak is associated with IRAS 2A1 and the secondary peak toward the south with IRAS 2A2. The directions of the two outflows originating from IRAS 2A1 and IRAS 2A2 are indicated with the colored arrows in the middle panel (Tobin et al. 2015). The size of the beam is shown in the bottom right, and in the bottom left a scale bar is displayed.</figcaption></figure>

However, the data can only be accurately fitted for a very large emitting area of ≥5000 au, which is in strong contrast with the extent of the emission in Fig. 3 (<150 au).

This is a clear indication that the assumption of LTE is not valid and thus that derived physical quantities, most notably the number of molecules, should be analyzed in more detail. Since the rotational distribution within the $v_3$ band is well-fitted by a single temperature, this temperature accurately describes the rotational temperature within this vibrationally excited state, but it does not necessarily represent the kinetic temperature. On the other hand, the values derived for the total number of molecules and the size of the emitting area are not valid, both because of the conflict with the measured size of the emitting area as well as the large discrepancy with the values derived for the vibrational ground state in the ALMA data (see Sect. 3.2). This will be further discussed in Sect. 4.1.

## 3.2. Comparison to ALMA data

Integrated intensity maps of the SO₂ 82,6–71,7 (*E*<sub>up</sub> = 43 K), SO₂ 164,12–163,13 (*E*<sub>up</sub> = 164 K), and ³⁴SO₂ 174,14–173,15 (*E*<sub>up</sub> = 178 K) pure rotational transitions in their vibrational ground states (i.e., *v* = 0) are presented in Fig. 6. The emission is spatially resolved and is extended in the northwest to southeast direction, following the direction of the outflow of IRAS 2A1 (e.g., Tobin et al. 2015; Jørgensen et al. 2022). However, most of the emission is concentrated in the inner ~1", similar to what is seen in the MIRI-MRS data (Fig. 3). After deconvolution with the beam, the size of the emission originating from the SO₂ 164,12–163,13 transition is ~0.5" × 0.3" in diameter, corresponding to ~150 × 90 au. The SO₂ 82,6–71,7 transition has a slightly larger deconvolved emitting size of ~0.6" × 0.4" (~180 × 120 au) likely because its lower *E*<sub>up</sub> makes it more sensitive to colder and more extended material in the envelope and outflow. The size of the ³⁴SO₂ emission is slightly smaller with a deconvolved size of ~0.4" × 0.2" (~120 × 60 au), but this likely originates from the lower signal-to-noise of this transition. These emitting areas are consistent with that of HDO 33,1–42,2 (*E*<sub>up</sub> = 335 K) and

complex organics such as methanol (CH₃OH) and methyl formate (CH₃OCHO; see Fig. A.3) and with the compact emission in the MIRI-MRS data (i.e., 0.5" × 0.3" in diameter; see Fig. 3).

The pure rotational lines of SO₂ and ³⁴SO₂ detected by ALMA were fitted using the same LTE slab model fitting procedure as for the MIRI-MRS data (Sect. 2.3). The submillimeter line lists of both SO₂ and ³⁴SO₂ were taken from the Cologne Database for Molecular Spectroscopy<fnref n="3" /> (CDMS; Müller et al. 2001, 2005; Endres et al. 2016), where the entry of SO₂ was mostly taken from Müller & Brünken (2005) and the entry of ³⁴SO₂ is based on several spectroscopic works (e.g., Lovas 1985; Belov et al. 1998, for transitions at ALMA Band 7 frequencies). The ALMA data cover in total five pure rotational (i.e., *v* = 0) transitions of SO₂ and also five transitions of the ³⁴SO₂ isotopologue (see Table A.1). However, several of these transitions are blended with strong emission of complex organics such as CH₃OCHO and CH₃CHO. Fortunately, two transitions of SO₂ and one transition of ³⁴SO₂ are (relatively) unblended. This is important to determine the column density and excitation temperature. All lines of SO₂ and ³⁴SO₂ are spectrally resolved with a FWHM of 3.5 km s⁻¹, which is similar to the FWHM of the lines of HDO 33,1–42,2 and the complex organics. No vibrational correction is necessary for the column densities derived from the pure rotational lines since the first excited vibrational states are included in the partition function.

Interestingly, the excitation temperature *T*<sub>ex</sub> (i.e., rotational temperature *T*<sub>rot</sub> since it is derived from pure rotational lines) derived for SO₂ from the ALMA data is very similar to that derived for MIRI: 104 ± 5 K (see Figs. A.1 and A.2). This suggests that ALMA is probing the same gas as that being probed with MIRI with the temperature reflecting the distribution of rotational levels within the *v* = 0 state. The column density can be accurately constrained at *N*<sub>SO₂,ALMA</sub> = 3.4 ± 1.1 × 10<sup>16</sup> cm<sup>-2</sup> for an emitting area with a radius of *R* = 85 au. In contrast to MIRI, the column density and emitting area can be fairly accurately constrained since the lines are marginally optically thick

<footnote n="3">³ https://cdms.astro.uni-koeln.de/</footnote>
(τ ~ 0.1). Moreover, the derived column density does not suffer severely from optical depth effects given that the column density derived from the optically thin ³⁴SO₂ isotopologue is 2.7 ± 0.5 × 10¹⁵ cm⁻², giving a ³²SO₂/³⁴SO₂ ratio of 13 ± 5, which is within a factor of 2 of the average ³²S/³⁴S derived for the local ISM (³²S/³⁴S = 22; Wilson 1999). However, the total number of SO₂ molecules measured by ALMA is N = 1.6 ± 0.2 × 10⁴⁷ molecules, which is ~4 orders of magnitude lower than what is detected by MIRI. This discrepancy cannot be explained by a difference in dust opacity since such an effect would only worsen the discrepancy in the detected number of molecules between ALMA and MIRI. The most logical explanation is the importance of non-LTE effects for the ro-vibrational transitions detected by MIRI.

# 4. Discussion

## 4.1. Importance of non-LTE effects

### 4.1.1. Absence of the v₁ and v₂ bands

If the emission originating from the v₃ asymmetric stretching mode were in LTE, both the v₁ symmetric stretching mode (λ ~ 8.5–9 µm) and v₂ bending mode (λ ~ 18–20 µm) should also have been detected. However, no SO₂ emission is detected in the v₁ and v₂ bands (see Fig. 1). The v₂ bending mode in particular has a significantly lower vibrational energy (518 cm⁻¹, 745 K) than the v₃ asymmetrical stretching mode (1362 cm⁻¹, 1960 K; Briggs 1970; Person & Zerbi 1982) and is therefore more easily collisionally excited. The expected SO₂ flux in the 18–20 µm region by employing the best-fit LTE models of the v₃ band around 7.35 µm (Sect. 3.1) and the vibrational ground state (i.e., (v₁, v₂, v₃) = (0, 0, 0), denoted as v = 0) in the ALMA data (Sect. 3.2) are presented in Fig. B.3. Similar to the analysis of the v₃ band, both models are corrected for an extinction of A<sub>V</sub> = 55 mag (Rocha et al. 2024) using a modified version of the McClure (2009) extinction law (see Appendix C). No obvious SO₂ emission features are present in the data whereas the best-fit v₃ LTE model clearly predicts that we should have seen the SO₂ at these wavelengths. On the other hand, the LTE model derived from the pure rotational lines in the ALMA data agrees with the non-detection of the v₂ band, predominantly due to the larger extinction at MIR wavelengths compared to millimeter wavelengths (e.g., McClure 2009; Chapman et al. 2009).

One solution for the lack of emission in the v₂ band could be that the line-to-continuum ratio is too low in this region, since the continuum flux level is on the order of ~1 Jy at 18–20 µm compared to ~20 mJy around 7.35 µm. The predicted line emission peaks at ~1 Jy for the strongest peaks of the v₃ model, corresponding to a line-to-continuum ratio of 1. This suggests that this cannot explain the absence of the v₂ band since a line-to-continuum ratio down to 0.01 should still be detectable on a strong continuum. Moreover, this would also not explain the discrepancy between the number of molecules needed in the LTE slab models to explain the emission in v = 0 state with ALMA and the v₃ band with MIRI. Scattered continuum radiation is present out to ~8 µm but not to longer wavelengths (see Fig. 1), suggesting that the v₃ band could be infrared-pumped on ~100 au scales whereas this does not occur for the v₂ band at 19 µm. The infrared radiation around 19 µm originates from thermal dust emission that is not extended (see Fig. 1).

Similarly to the v₂ bending mode, the v₁ symmetrical stretching mode between 8.5 and 9 µm is overproduced by the best-fit v₃ model (see Fig. B.4). The absence of the emission from the

v₁ band can be simply explained by the large extinction due to strong silicate absorption feature around these wavelengths

### 4.1.2. Infrared pumping

The critical densities of ro-vibrational transitions are typically >10¹⁰ cm⁻³ (e.g., for HCN and CO₂; Bruderer et al. 2015; Bosman et al. 2017, for SO₂ no collisional rates are available for the MIR transitions but they are expected to be similar), suggesting that LTE conditions are only valid for ro-vibrational transitions in the inner ≤1 au. However, SO₂ is present in IRAS 2A on ~100 au scales (see Figs. 3 and 6) where the densities are typically 10⁶–10⁸ cm⁻³ (e.g., Jørgensen et al. 2004, 2022; Kristensen et al. 2012), far below the critical densities of the MIR transitions, further suggesting that the vibrational level populations are not collisionally excited. The critical densities for the v = 0 transitions probed by ALMA are much lower, on the order of ~10⁶–10⁷ cm⁻³, meaning that the LTE assumption is valid for these pure rotational lines.

Despite the fact that the assumption of LTE is not valid for the MIR transitions, physical information can still be derived. The rotational temperature derived from the LTE models of the MIR v₃ transitions is well constrained given the accuracy of the fit to the R- and Q-branches. It is very similar to the rotational temperature derived for the v = 0 state, for which the LTE assumption is valid. Combined with the similar extent of the emission in the MIRI and ALMA images (Figs. 3 and 6, respectively), this strongly suggests that they are probing SO₂ gas on similar scales (i.e., 50–100 au). Because the rotational levels within the v = 0 state can be characterized by a single temperature, the column density of SO₂ derived from the ALMA data is also robust. However, the total number of SO₂ molecules derived from the v₃ MIR transitions is found to be >4 orders of magnitude higher than what is measured for the v = 0 state. This suggests that the v₃ vibrational state is more highly populated than just through collisional excitation by itself.

A likely explanation for populating the v₃ vibrational level could be infrared pumping. Vibrational levels can become more highly populated due to the absorption of infrared photons in the presence of a strong infrared radiation field, boosting the line fluxes far above that expected from collisional excitation alone (e.g., Bruderer et al. 2015; Bosman et al. 2017). During the infrared radiative pumping process, the distribution of rotational levels within the vibrational states is largely maintained since only ΔJ = 0, ±1 transitions are allowed. Infrared pumping thus leads to a more highly populated vibrational level without significantly changing the rotational distribution (i.e., rotational temperature), supporting the similarity in rotational temperatures measured in the v₃ = 1 and v = 0 states.

The importance of infrared pumping can be further quantified by comparing the results of the LTE models to the v₃ band and the v = 0 state. Assuming the v₃ = 1 vibrational level is radiatively pumped from the v = 0 state, the vibrational level population is set following the Boltzmann distribution by a vibrational temperature T<sub>vib</sub> that is different from the rotational temperature, T<sub>vib</sub> ≠ T<sub>rot</sub>. The difference between the total number of molecules predicted in LTE models of the MIRI and ALMA data can then be approximated as

$$ \frac{N_{SO_2,MIRI}}{N_{SO_2,ALMA}} \propto \frac{e^{-\hbar v/(k_B T_{vib})}}{e^{-\hbar v/(k_B T_{rot})}}, \quad (1) $$

where N<sub>ALMA</sub> and N<sub>MIRI</sub> are the total number of molecules needed in the LTE models of the ALMA and MIRI data, respectively, T<sub>rot</sub> the derived rotational temperature in the v = 0 state
with ALMA (i.e., 104 ± 5 K; Sect. 3.2), $T_{vib}$ the vibrational temperature, $v$ the frequency of the MIR transitions, and $h$ and $k_B$ Planck's and Boltzmann's constants, respectively. In Eq. (1), $g$ factors have been neglected for simplicity. In the case where both $v_3 = 1$ and $v = 0$ are collisionally excited, $T_{vib} = T_{rot}$ and the number of SO₂ molecules predicted by the LTE models of the MIRI and ALMA data should have been equal, $N_{ALMA} = N_{MIRI}$. However, for the derived $N_{SO_2,MIRI}/N_{SO_2,ALMA} \sim 2 \times 10^4$ and $T_{rot} = 104$ K, Eq. (1) results in $T_{vib} \sim 200$ K at a wavelength of 7.35 µm (i.e., $v = 40.79$ THz), which is ~80 K higher than the derived $T_{rot}$ of the vibrational ground state. We estimate an uncertainty on $T_{vib}$ of 50 K based on the assumption of similar rotational temperatures in the $v_3 = 1$ and $v = 0$ states and ignoring $g$ factors in Eq. (1).

If infrared pumping is indeed the cause for the elevated vibrational temperature, the brightness temperature of the infrared radiation around 7.35 µm, $T_{IR}$, has to be similar to $T_{vib}$,

$$T_{IR} = \frac{h\nu}{k_B} \ln^{-1} \left( 1 + \frac{2h\nu^3}{I_v c^2} \right), \quad (2)$$

where $I_v$ is the observed intensity (in Jy sr⁻¹) of the continuum at a frequency $v$ and $c$ is the speed of light. Given that $I_v \sim 6 \times 10^4$ MJy sr⁻¹ (i.e., 1.4 Jy arcsec⁻²) around 7.35 µm (i.e., 40.79 THz) and assuming an extinction of $A_V = 55$ mag (Rocha et al. 2024), this results in $T_{IR} \sim 180$ K. This is in very good agreement with $T_{vib} \sim 200$ K and therefore strongly suggests that the vibrational temperature of the $v_3$ band of SO₂ is indeed set by the infrared radiation field of the protostar rather than by collisions.

In turn, this means that both MIRI and ALMA are in fact tracing the same molecular gas with the same rotational temperature, but that the transitions in the $v_3$ band are just pumped by the strong infrared radiation field of the central star and originate from a region much smaller than 5000 au. ALMA thus measures the true number of SO₂ molecules (1.8 ± 0.2 × 10⁴⁷), and in order to directly derive the physical properties (i.e., excitation temperature, column density, emitting area) from the ro-vibrational transition of SO₂ at MIR wavelengths it is important to take into account non-LTE processes such as infrared radiative pumping.

## 4.2. Physical origin of the SO₂

The SO₂ emission detected with both ALMA and MIRI is clearly tracing warm inner regions of the protostellar system on disk scales of ≲100 au in radius. This implies that the emission either originates from the central hot core where many complex organs are detected (e.g., Jørgensen et al. 2005; Bottinelli et al. 2007; Maury et al. 2014; Taquet et al. 2015), or that the emission originates from an accretion shock at the disk-envelope interface (e.g., Artur de la Villarmois et al. 2019, 2022; van Gelder et al. 2021). Alternatively, the SO₂ could be located in a jet or outflow close to the source (e.g., Codella et al. 2014; Taquet et al. 2020; Tychoniec et al. 2021), or even in a disk wind (e.g., Tabone et al. 2017). However, since the emission in the MIRI data is clearly centrally peaked with only little extended emission (≈100 au in radius; see Fig. 3) and the line widths of the pure rotational lines detected by ALMA are ~3.5 km s⁻¹ (compared to >10 km s⁻¹ typically observed in outflows; e.g., Taquet et al. 2020; Tychoniec et al. 2019, 2021), an extended outflow origin is less likely. On the contrary, CO₂ is clearly detected in the outflow and not on the continuum peak (see Fig. B.1), which will be further discussed in a separate paper.

One way to investigate further the possible origin of the SO₂ emission is to estimate its abundance with respect to H₂. The most direct way to determine the number of H₂ molecules is using the detected H₂ MIR lines (see Appendix D). However, H₂ does not only trace the warm inner regions but is also present in outflows and disk winds (Tychoniec et al., in prep.). Indeed, in IRAS 2A H₂ is mostly located in the outflow toward the southwest (see Fig. B.1). Moreover, the rotational temperature of the warm component is >300 K (see Fig. D.2), which further suggests that H₂ may not be tracing the same gas as SO₂ ($T_{rot} = 104 \pm 5$ K).

Another method for determining the amount of ~100 K H₂ gas is to use the ALMA Band 7 continuum. The total gas mass, $M_{gas}$, can be determined using the equation from Hildebrand (1983),

$$M_{\text{gas}} = 100 \frac{F_v d^2}{\kappa_v B_v(T_{\text{dust}})}, \quad (3)$$

where $F_v$ is the continuum flux density at a frequency $v$, $d$ the distance (293 pc; Ortiz-León et al. 2018), $\kappa_v$ the dust opacity, $B_v(T_{dust})$ the Planck function evaluated at a dust temperature $T_{dust}$, and the factor 100 the assumed gas-to-dust mass ratio. Here, $T_{dust}$ is assumed to be 30 K, which is a typical dust temperature for protostellar envelopes at ~100 au scales (Whitney et al. 2003). For a measured $F_v = 0.396 \pm 0.079$ Jy in a 1.4'' diameter aperture (i.e., the same as was used in the spectral extraction) and a dust opacity of $\kappa_v = 1.84 \text{ cm}^{-2} \text{ g}^{-1}$ at a frequency of 340 GHz (Ossenkopf & Henning 1994), a gas mass of $M_{gas} = 0.11 \pm 0.02 \text{ M}_\odot$ is derived. The derived gas mass is consistent within a factor of 2 with other recent measurements (e.g., Tobin et al. 2018; Tychoniec et al. 2020). Assuming that this gas mass is predominantly in H₂ (i.e., taking a mean molecular weight $\mu_{H_2} = 2.8$ per hydrogen atom for gas composed of 71% hydrogen; Kauffmann et al. 2008), this results in $N_{H_2} = 4.7 \pm 0.9 \times 10^{55}$ molecules. This is in good agreement with $N_{H_2} = 6.6 \pm 0.9 \times 10^{55}$ molecules derived for IRAS 2A from C¹⁸O lines up to $J = 9-8$ with Herschel-HIFI (assuming a CO/H₂ ratio of $10^{-4}$ and rotational temperature of ~40 K; Yıldız et al. 2013). However, the derived radius of the SO₂ emitting area (~100 au) is smaller than the physical radius of the aperture (205 au for 1.4'' diameter aperture) used for computing $N_{H_2}$. The derived $N_{H_2}$ can be scaled to the derived emitting area of SO₂ assuming that the density scales as $n_H \propto R^{-p}$,

$$N_{H_2, T>100K} = N_{H_2} \left( \frac{R_{100K}}{R_{\text{ap}}} \right)^{3-p}, \quad (4)$$

where $R_{100K}$ is the derived size emitting area of SO₂ from the ALMA and MIRI-MRS maps (~100 au; see Sect. 3) and $R_{ap}$ is physical radius of the aperture (205 au for 1.4'' diameter aperture). Taking a density power-law index of 1.7 for IRAS 2A1 (Kristensen et al. 2012), this results in $N_{H_2, T>100K} = 1.8 \pm 0.4 \times 10^{55}$ molecules. This is almost equal to what is derived directly from the H₂ lines (1.7 ± 0.9 × 10⁵⁵ molecules; Appendix D), which suggests that the H₂ 0-0 lines themselves are likely sensitive to most of the $T > 100$ K gas despite the high rotational temperature ($T_{rot} = 356 \pm 41$ K).

The total number of SO₂ molecules derived from the emission of the vibrational ground state in the ALMA data can be directly compared to H₂. This results in an SO₂ abundance on the order of 1.0 ± 0.3 × 10⁻⁸ with respect to H₂, which is in agreement with other recent SO₂ abundance measurements in
low-mass Class 0 systems ($>6.6 \times 10^{-10}$; Artur de la Villarmois et al. 2023). The abundance of $\sim 10^{-8}$ implies that the SO₂ gas is not the dominant sulfur carrier in IRAS 2A since the cosmic [S/H] abundance is about $\sim 10^{-5}$ (Savage & Sembach 1996; Goicoechea et al. 2006). It is also means that SO₂ does not contain a significant amount of the total volatile sulfur budget in dense clouds (volatile [S/H] $\sim 10^{-7}$, i.e., the amount of sulfur that is not locked up in refractory formats; Woods et al. 2015; Kama et al. 2019). Furthermore, the derived SO₂ abundance is on the lower side compared with estimates of the SO₂ abundances in interstellar ices ($10^{-8}-10^{-7}$; Boogert et al. 1997, 2015; Zasowski et al. 2009) and cometary ices ($\sim 10^{-7}$; Altwegg et al. 2019; Rubin et al. 2019). In particular, SO₂ ice was also recently detected in absorption toward IRAS 2A itself with a very similar abundance with respect to H₂ as other low-mass protostars ($\sim 10^{-7}$; Rocha et al. 2024).

Although the gaseous SO₂ abundance is on the lower side compared to the ices, it suggests that the observed gaseous SO₂ could be sublimated from the ices in the central hot core since the sublimation temperature of SO₂ ($\sim 60$ K, E<sub>bin</sub> = 3010 K; Penteado et al. 2017) is similar to H₂O and many complex organics. It also naturally explains the compactness of the SO₂ emission in both the MIRI data and the ALMA data (see Figs. 3 and 6) with an extent very similar to that of H₂O, HDO, and complex organics (see Figs. B.1 and A.3). In fact, the emitting areas derived from the MIRI-MRS and ALMA integrated intensity maps ($\sim 100$ au in radius) agree well with that of a hot core based on the luminosity (i.e., R<sub>100K</sub> = 147 au, where R<sub>100K</sub> $\approx 15.4 \sqrt{L_{bol}/L_{\odot}}$; Bisschop et al. 2007; van't Hoff et al. 2022). The hot core origin is further supported by the line width in the ALMA data ($\sim 3.5$ km s<sup>-1</sup>), which is similar to that of lines from HDO and complex organics.

Another possibility for the SO₂ emission could be weak shocks on $\le 100$ au scales such as accretion shocks at the disk-envelope boundary (e.g., Sakai et al. 2014; Oya et al. 2019; Artur de la Villarmois et al. 2019, 2022). Recent shock models have suggested that SO₂ could be a good tracer of such accretion shocks (Miura et al. 2017; van Gelder et al. 2021). Abundances up to $10^{-8}-10^{-7}$ with respect to H₂ are easily reached in low-velocity accretion shocks, as long as a significant ultraviolet (UV) radiation field is present (i.e., stronger than the interstellar radiation field; van Gelder et al. 2021). Given the luminosity of IRAS 2A ($\sim 60-90 L_{\odot}$; Murillo et al. 2016; Karska et al. 2018), having a considerable UV field in the inner envelope out to 100 au scales is likely (see e.g., van Kempen et al. 2009; Yıldız et al. 2012, 2015). The emission morphology seen in the ALMA data (Fig. 6) shows that the emission is extended more in the direction of the outflow than along the disk. Furthermore, the line width of $\sim 3.5$ km s<sup>-1</sup> is consistent with those of the complex organics whereas accretion shocks are expected to show broader emission lines ($\sim 10$ km s<sup>-1</sup>; Oya et al. 2019; Artur de la Villarmois et al. 2019, 2022). A hot core origin is therefore a more likely explanation for the gaseous SO₂.

The SO₂ emission in IRAS 2A could be tracing similar components as the MIR SO₂ absorption that was detected toward multiple high-mass protostellar sources (e.g. Keane et al. 2001; Dungee et al. 2018; Nickerson et al. 2023). Typical temperatures of this SO₂ absorption are $\sim 100-300$ K (e.g., Dungee et al. 2018; Nickerson et al. 2023), which is very similar to the rotational temperature derived for IRAS 2A ($\sim 100$ K), though temperatures up to $\sim 700$ K have also been reported (Keane et al. 2001). The most common origin of the SO₂ absorption toward these high-mass sources is also suggested to be a hot core rather than a shock, similar to what SO₂ is tracing in IRAS 2A, with

the main difference that it is in absorption against the bright infrared continuum of the central high-mass protostar. The typical SO₂ abundances derived for these high-mass protostars are $\ge 10^{-7}$ with respect to H₂ (Keane et al. 2001; Dungee et al. 2018; Nickerson et al. 2023), which is an order of magnitude higher than what is derived for IRAS 2A.

It remains unknown whether the presence of ro-vibrational lines of SO₂ in IRAS 2A is unique or whether it is more common among low-mass protostars. One other low-mass protostar in the JOYS+ sample, NGC 1333 IRAS1A, shows emission of the SO₂ $v_3$ Q-branch (R- and P-branches are not detected) and similarly does not show any emission from the $v_1$ and $v_2$ bands, suggesting that infrared radiative pumping may also be responsible for the emission in the $v_3$ band of this source. IRAS 2A has a high luminosity ($\sim 60-90 L_{\odot}$; Murillo et al. 2016; Karska et al. 2018) and is suggested to currently be in a burst phase (e.g., Hsieh et al. 2019; van't Hoff et al. 2022), whereas IRAS 1A has a lower luminosity ($\sim 10 L_{\odot}$; Tobin et al. 2016). On the other hand, the low-mass Class 0 protostar IRAS 15398-3359 does not show gaseous SO₂ and has a lower luminosity of $\sim 1.5 L_{\odot}$ (Yang et al. 2018, 2022). Given the importance of infrared pumping in detecting the SO₂ lines, a high luminosity could be important. However, a larger sample of JWST/MIRI-MRS spectra of low-mass and high-mass embedded protostellar systems is needed to further investigate the importance of luminosity (and other source properties) on the presence of MIR SO₂ emission. Future MIRI-MRS observations from the JOYS+ program will cover a few other luminous low-mass protostellar systems (e.g., Serpens SMM1, $L_{bol} \sim 100 L_{\odot}$; Karska et al. 2018), as well as several high-mass protostellar systems (e.g., IRAS 18089-1732, $L_{bol} \sim 10^4 L_{\odot}$; Urquhart et al. 2018).

## 5. Conclusions

This paper presents one of the first medium-resolution MIR spectra and images taken with JWST/MIRI-MRS of a Class 0 protostellar source, NGC 1333 IRAS 2A, and presents the first detection of gas-phase SO₂ emission at MIR wavelengths. We analyzed the spectral lines of the $v_3$ asymmetric stretching mode of SO₂ using LTE slab models and find a rotational temperature of 92 ± 8 K. This is very similar to the rotational temperature of 104 ± 5 K derived from the pure rotational lines in the high-resolution ALMA data. Since the SO₂ emission in the MIRI-MRS data is optically thin, the column density could not be constrained accurately due to the degeneracy with the emitting area. However, the total number of molecules can be constrained and is predicted by LTE models of the MIRI-MRS data to be a factor of $\sim 2 \times 10^4$ higher than that derived from the ALMA data. Based on these results, our main conclusions are as follows:

- The $v_3$ asymmetric stretching mode of SO₂ detected around 7.35 µm with MIRI is not in LTE but rather radiatively pumped by a strong infrared radiation field scattered out to $\sim 100$ au distances. The discrepancy between the best-fit LTE models of the MIRI-MRS and ALMA data can be explained by a vibrational temperature ($\sim 200$ K) that is higher than the rotational temperature derived within the vibrational ground state (104 ± 5 K). The vibrational temperature is consistent with the brightness temperature of the continuum around 7.35 µm ($\sim 180$ K). The similarity in rotational temperatures suggests that MIRI-MRS and ALMA are in fact still tracing the same molecular gas.
- Assuming that ALMA is probing the total amount of SO₂, the abundance of gaseous SO₂ is estimated to be 1.0 ± 0.3 × 10⁻⁸ with respect to H₂, which is consistent with both a hot core and accretion shock origin. Based on the size of the emitting area (~100 au in radius) and the small line width of the SO₂ lines (~3.5 km s⁻¹) in the ALMA data, a hot core origin is suggested to be the most likely.

Our results show the importance of taking non-LTE effects into account when analyzing ro-vibrational lines at MIR wavelengths. The synergy between JWST probing the hot spots and the spatial and spectral resolution of ALMA proves to be crucial in determining the physical origin of molecules in embedded protostellar systems. Future astrochemical modeling with non-LTE effects included will be valuable for inferring the physics and the chemistry of SO₂ in the earliest phases of star formation and their connection to the sulfur depletion problem.

**Acknowledgements.** We would like to thank the anonymous referee for their constructive comments on the manuscript, and Valentin Christiaens, Matthias Samland, and Danny Gasman for valuable support with the MIRI-MRS data reduction. 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 program #1236. The following National and International Funding Agencies funded and supported the MIRI development: NASA; ESA; Belgian Science Policy Office (BELSPO); Centre Nationale d'Etudes Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum fur Luftund Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economía y Competitividad; The Netherlands Research School for Astronomy (NOVA); The Netherlands Organisation for Scientific Research (NWO); Science and Technology Facilities Council; Swiss Space Office; Swedish National Space Agency; and UK Space Agency. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2021.1.01578.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The PI acknowledges assistance from Aida Ahmadi from Allegro, the European ALMA Regional Center node in the Netherlands. M.v.G., E.v.D., L.F., H.L., K.S., and W.R. acknowledge support from ERC Advanced grant 101019751 MOLDISK, TOP-1 grant 614.001.751 from the Dutch Research Council (NWO), The Netherlands Research School for Astronomy (NOVA), the Danish National Research Foundation through the Center of Excellence "InterCat" (DNRF150), and DFG-grant 325594231, FOR 2634/2. The work of M.E.R. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. P.J.K. acknowledges support from the Science Foundation Ireland/Irish Research Council Pathway programme under Grant Number 21/PATH-S/9360. L.M. acknowledges the financial support of DAE and DST-SERB research grants (SRG/2021/002116 and MTR/2021/000864) from the Government of India.

## References

Altwegg, K., Balsiger, H., & Fuselier, S. A. 2019, ARA&A, 57, 113

Arce, H. G., Borkin, M. A., Goodman, A. A., Pineda, J. E., & Halle, M. W. 2010, ApJ, 715, 1170

Argyriou, I., Glasse, A., Law, D. R., et al. 2023, A&A, 675, A111

Artur de la Villarmois, E., Jørgensen, J. K., Kristensen, L. E., et al. 2019, A&A, 626, A71

Artur de la Villarmois, E., Guzmán, V. V., Jørgensen, J. K., et al. 2022, A&A, 667, A20

Artur de la Villarmois, E., Guzmán, V. V., Yang, Y. L., Zhang, Y., & Sakai, N. 2023, A&A, 678, A124

Belov, S. P., Tretyakov, M. Y., Kozin, I. N., et al. 1998, J. Mol. Spectrosc., 191, 17

Beuther, H., van Dishoeck, E. F., Tychoniec, L., et al. 2023, A&A, 673, A121

Bisschop, S. E., Jørgensen, J. K., van Dishoeck, E. F., & de Wachter, E. B. M. 2007, A&A, 465, 913

Blake, G. A., Sutton, E. C., Masson, C. R., & Phillips, T. G. 1987, ApJ, 315, 621

Boogert, A. C. A., Schutte, W. A., Helmich, F. P., Tielens, A. G. G. M., & Wooden, D. H. 1997, A&A, 317, 929

Boogert, A. C. A., Pontoppidan, K. M., Knez, C., et al. 2008, ApJ, 678, 985

Boogert, A. C. A., Gerakines, P. A., & Whittet, D. C. B. 2015, ARA&A, 53, 541

Boonman, A. M. S., van Dishoeck, E. F., Lahuis, F., et al. 2003, A&A, 399, 1047

Bosman, A. D., Bruderer, S., & van Dishoeck, E. F. 2017, A&A, 601, A36

Bottinelli, S., Ceccarelli, C., Williams, J. P., & Lefloch, B. 2007, A&A, 463, 601

Briggs, A. G. 1970, J. Chem. Education, 47, 391

Bruderer, S., Harsono, D., & van Dishoeck, E. F. 2015, A&A, 575, A94

Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2023, ASP Conf. Ser., 527, 583

Carr, J. S., & Najita, J. R. 2008, Science, 319, 1504

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

Christiaens, V., Gonzalez, C., Farkas, R., et al. 2023, J. Open Source Softw., 8, 4774

Codella, C., Maury, A. J., Gueth, F., et al. 2014, A&A, 563, L3

Codella, C., Bianchi, E., Podio, L., et al. 2021, A&A, 654, A52

Drozdovskaya, M. N., van Dishoeck, E. F., Jørgensen, J. K., et al. 2018, MNRAS, 476, 4949

Dungee, R., Boogert, A., DeWitt, C. N., et al. 2018, ApJ, 868, L10

Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, J. Mol. Spectrosc., 327, 95

Francis, L., van Gelder, M. L., van Dishoeck, E. F., et al. 2024, A&A, in press https://doi.org/10.1051/0004-6361/202348105

Garufi, A., Podio, L., Codella, C., et al. 2022, A&A, 658, A104

Giesen, C., Beuther, H., van Dishoeck, E. F., et al. 2023, A&A, 679, A108

Goicoechea, J. R., Pety, J., Gerin, M., et al. 2006, A&A, 456, 565

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

Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2022, J. Quant. Spec. Radiat. Transf., 277, 107949

Grant, S. L., van Dishoeck, E. F., Tabone, B., et al. 2023, ApJ, 947, L6

Greenfield, P., & Miller, T. 2016, Astron. Comput., 16, 41

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

Hildebrand, R. H. 1983, QJRAS, 24, 267

Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2019, ApJ, 884, 149

Jones, O. C., Álvarez-Márquez, J., Sloan, G. C., et al. 2023, MNRAS, 523, 2519

Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2005, ApJ, 632, 973

Jørgensen, J. K., Kuruwita, R. L., Harsono, D., et al. 2022, Nature, 606, 272

Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2004, A&A, 416, 603

Kama, M., Shorttle, O., Jermyn, A. S., et al. 2019, ApJ, 885, 114

Karska, A., Kaufman, M. J., Kristensen, L. E., et al. 2018, ApJS, 235, 30

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

Keane, J. V., Boonman, A. M. S., Tielens, A. G. G. M., & van Dishoeck, E. F. 2001, A&A, 376, L5

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

Kushwahaa, T., Drozdovskaya, M. N., Tychoniec, Ł., & Tabone, B. 2023, A&A, 672, A122

Labiano, A., Argyriou, I., Álvarez-Márquez, J., et al. 2021, A&A, 656, A57

Law, D. D., Morrison, J. E., Argyriou, I., et al. 2023, AJ, 166, 45

Lovas, F. J. 1985, J. Phys. Chem. Ref. Data, 14, 395

Maury, A. J., Belloche, A., André, P., et al. 2014, A&A, 563, L2

McClure, M. 2009, ApJ, 693, L81

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

McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127

Miura, H., Yamamoto, T., Nomura, H., et al. 2017, ApJ, 839, 47

Müller, H. S. P., & Brünken, S. 2005, J. Mol. Spectr., 232, 213

Müller, H. S. P., Thorwirth, S., Roth, D. A., & Winnewisser, G. 2001, A&A, 370, L49

Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, J. Mol. Struct., 742, 215

Murillo, N. M., van Dishoeck, E. F., Tobin, J. J., & Fedele, D. 2016, A&A, 592, A56

Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, A&A, 637, A39

Nickerson, S., Rangwala, N., Colgan, S. W. J., et al. 2023, ApJ, 945, 26

Öberg, K. I., Boogert, A. C. A., Pontoppidan, K. M., et al. 2011, ApJ, 740, 109

Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJ, 865, 73

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Oya, Y., López-Sepulcre, A., Sakai, N., et al. 2019, ApJ, 881, 112

Penteado, E. M., Walsh, C., & Cuppen, H. M. 2017, ApJ, 844, 71

Perotti, G., Christiaens, V., Henning, T., et al. 2023, Nature, 620, 516

Person, W. B., & Zerbi, G. 1982, Vibrational Intensities in Infrared and Raman Spectroscopy, Studies in Physical and Theoretical Chemistry (Amsterdam: Elsevier Scientific Publishing Company)

Pineau des Forêts, G., Roueff, E., Schilke, P., & Flower, D. R. 1993, MNRAS, 262, 915
Pontoppidan, K. M., Schöier, F. L., van Dishoeck, E. F., & Dartois, E. 2002, A&A, 393, 585

Rieke, G. H., Wright, G. S., Böker, T., et al. 2015, PASP, 127, 584

Rocha, W. R. M., van Dishoeck, E. F., Ressler, M. E., et al. 2024, A&A, in press, https://doi.org/10.1051/0004-6361/202348427

Rubin, M., Altwegg, K., Balsiger, H., et al. 2019, MNRAS, 489, 594

Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999, MNRAS, 306, 691

Sakai, N., Sakai, T., Hirota, T., et al. 2014, Nature, 507, 78

Salyk, C., Pontoppidan, K. M., Blake, G. A., Najita, J. R., & Carr, J. S. 2011, ApJ, 731, 130

Savage, B. D., & Sembach, K. R. 1996, ARA&A, 34, 279

Schutte, W. A., Boogert, A. C. A., Tielens, A. G. G. M., et al. 1999, A&A, 343, 966

Sonnentrucker, P., González-Alfonso, E., & Neufeld, D. A. 2007, ApJ, 671, L37

Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6

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

Taquet, V., López-Sepulcre, A., Ceccarelli, C., et al. 2015, ApJ, 804, 81

Taquet, V., Codella, C., De Simone, M., et al. 2020, A&A, 637, A63

Terwisscha van Scheltinga, J., Ligterink, N. F. W., Boogert, A. C. A., van Dishoeck, E. F., & Linnartz, H. 2018, A&A, 611, A35

Tobin, J. J., Dunham, M. M., Looney, L. W., et al. 2015, ApJ, 798, 61

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2016, ApJ, 818, 73

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2018, ApJ, 867, 43

Tsai, S.-M., Lee, E. K. H., Powell, D., et al. 2023, Nature, 617, 483

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

Tychoniec, Ł., Manara, C. F., Rosotti, G. P., et al. 2020, A&A, 640, A19

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

Urquhart, J. S., König, C., Giannetti, A., et al. 2018, MNRAS, 473, 1059

van der Tak, F. F. S., Lique, F., Faure, A., Black, J. H., & van Dishoeck, E. F. 2020, Atoms, 8, 15

van Dishoeck, E. F., Grant, S., Tabone, B., et al. 2023, Faraday Discuss., 245, 52

van Gelder, M. L., Tabone, B., van Dishoeck, E. F., & Godard, B. 2021, A&A, 653, A159

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

van't Hoff, M. L. R., Harsono, D., van Gelder, M. L., et al. 2022, ApJ, 924, 5

Wells, M., Pel, J. W., Glasse, A., et al. 2015, PASP, 127, 646

Whitney, B. A., Wood, K., Bjorkman, J. E., & Cohen, M. 2003, ApJ, 598, 1079

Wilson, T. L. 1999, Rep. Prog. Phys., 62, 143

Woods, P. M., Occhiogrosso, A., Viti, S., et al. 2015, MNRAS, 450, 1256

Wright, G. S., Wright, D., Goodson, G. B., et al. 2015, PASP, 127, 595

Wright, G. S., Rieke, G. H., Glasse, A., et al. 2023, PASP, 135, 048003

Yang, Y.-L., Green, J. D., Evans, Neal J., I., et al. 2018, ApJ, 860, 174

Yang, Y.-L., Green, J. D., Pontoppidan, K. M., et al. 2022, ApJ, 941, L13

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

Yıldız, U. A., Kristensen, L. E., van Dishoeck, E. F., et al. 2013, A&A, 556, A89

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

Zasowski, G., Kemper, F., Watson, D. M., et al. 2009, ApJ, 694, 459

<footnote n="1">¹ Leiden Observatory, Leiden University, PO Box 9513, 2300RA Leiden, The Netherlands
e-mail: vgelder@strw.leidenuniv.nl</footnote>

<footnote n="2">² Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA</footnote>

<footnote n="3">³ Max Planck Institut für Extraterrestrische Physik (MPE), Giessen-bachstrasse 1, 85748 Garching, Germany</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France</footnote>

<footnote n="5">⁵ Chalmers University of Technology, Department of Space, Earth and Environment, Onsala Space Observatory, 439 92 Onsala, Sweden</footnote>

<footnote n="6">⁶ European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany</footnote>

<footnote n="7">⁷ SETI Institute 189 Bernardo Avenue, 2nd Floor, Mountain View, CA 94043, USA</footnote>

<footnote n="8">⁸ Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="9">⁹ INAF-Osservatorio Astronomico di Capodimonte, Salita Moiariello 16, 80131 Napoli, Italy</footnote>

<footnote n="10">¹⁰ NASA Ames Research Center, Space Science and Astrobiology Division M.S. 245-6 Moffett Field, CA 94035, USA</footnote>

<footnote n="11">¹¹ Department of Experimental Physics, Maynooth University, Maynooth, Co Kildare, Ireland</footnote>

<footnote n="12">¹² UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="13">¹³ Bay Area Environmental Research Institute and NASA Ames Research Center, Moffett Field, CA 94035, USA</footnote>

<footnote n="14">¹⁴ Laboratory for Astrophysics, Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>

<footnote n="15">¹⁵ School of Earth and Planetary Sciences, National Institute of Science Education and Research, Jatni 752050, Odisha, India</footnote>

<footnote n="16">¹⁶ Homi Bhabha National Institute, Training School Complex, Anushaktinagar, Mumbai 400094, India</footnote>
## Appendix A: Additional ALMA figures and table

<table><caption>**Table A.1.** Transitions of SO₂ and ³⁴SO₂ covered in the 2021.1.01578.S ALMA program.</caption><thead><tr><th rowspan="2">Isotopologue</th><th colspan="2">Transition</th><th rowspan="2">Frequency<br>GHz</th><th rowspan="2">E<sub>up</sub><br>K</th><th rowspan="2">A<sub>ij</sub><br>s<sup>-1</sup></th><th rowspan="2">g<sub>up</sub></th><th rowspan="2">Detection</th></tr><tr><th>J<sub>K<sub>a</sub>,K<sub>c</sub></sub></th><th>J'<sub>K'<sub>a</sub>,K'<sub>c</sub></sub></th></tr></thead><tbody><tr><td rowspan="5">SO<sub>2</sub></td><td>8<sub>2,6</sub></td><td>7<sub>1,7</sub></td><td>334.673353</td><td>43.1</td><td>1.29 &times; 10<sup>-4</sup></td><td>17</td><td>Y</td></tr><tr><td>23<sub>3,21</sub></td><td>23<sub>2,22</sub></td><td>336.089228</td><td>276.0</td><td>2.67 &times; 10<sup>-4</sup></td><td>47</td><td>B</td></tr><tr><td>16<sub>7,9</sub></td><td>17<sub>6,12</sub></td><td>336.669581</td><td>245.1</td><td>5.84 &times; 10<sup>-5</sup></td><td>33</td><td>B</td></tr><tr><td>13<sub>2,12</sub></td><td>12<sub>1,11</sub></td><td>345.338538</td><td>93.0</td><td>2.38 &times; 10<sup>-4</sup></td><td>27</td><td>B</td></tr><tr><td>16<sub>4,12</sub></td><td>16<sub>3,13</sub></td><td>346.523878</td><td>164.5</td><td>3.43 &times; 10<sup>-4</sup></td><td>33</td><td>Y</td></tr><tr><td rowspan="5">³⁴SO<sub>2</sub></td><td>20<sub>8,12</sub></td><td>21<sub>7,15</sub></td><td>334.996573</td><td>334.0</td><td>6.19 &times; 10<sup>-5</sup></td><td>41</td><td>N</td></tr><tr><td>9<sub>4,6</sub></td><td>9<sub>3,7</sub></td><td>345.285620</td><td>79.2</td><td>2.88 &times; 10<sup>-4</sup></td><td>19</td><td>B</td></tr><tr><td>4<sub>4,0</sub></td><td>4<sub>3,1</sub></td><td>345.678787</td><td>47.1</td><td>1.31 &times; 10<sup>-4</sup></td><td>9</td><td>B</td></tr><tr><td>17<sub>4,14</sub></td><td>17<sub>3,15</sub></td><td>345.929349</td><td>178.5</td><td>3.64 &times; 10<sup>-4</sup></td><td>35</td><td>Y</td></tr><tr><td>28<sub>2,6</sub></td><td>28<sub>1,27</sub></td><td>347.483124</td><td>390.7</td><td>2.65 &times; 10<sup>-4</sup></td><td>57</td><td>N</td></tr></tbody></table>

**Notes.** The spectroscopic information of both SO₂ and ³⁴SO₂ are taken from CDMS (Müller et al. 2001, 2005; Endres et al. 2016). The entry of SO₂ was mostly based on the work of Müller et al. (2005) and the entry of ³⁴SO₂ is based on several spectroscopic works (e.g., Lovas 1985; Belov et al. 1998, for transitions at ALMA Band 7 frequencies). The last column indicates whether the specific transition was detected (Y), blended (B), or not detected (N).

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

**Fig. A.1.** Continuum-subtracted ALMA spectra (blue) for the SO₂ 8,2,6 – 7,1,7 (left, E<sub>up</sub> = 43 K), SO₂ 16,4,12 – 16,3,13 (middle, E<sub>up</sub> = 164 K), and ³⁴SO₂ 17,4,14 – 17,3,15 (right, E<sub>up</sub> = 178 K) transitions. The best-fit LTE slab models are overlaid in orange. All other transitions (see Table A.1) of SO₂ and ³⁴SO₂ are either blended with other molecular species or not detected.
<figure><img src="image_9.png" /><figcaption>**Fig. A.2.** χ² map for SO₂ derived from pure rotational lines in the ALMA data. The χ² values are inversely normalized by the minimum χ² and shown in color. The 1, 2, and 3σ confidence intervals are presented as the black contours. The best-fitting emitting radius, *R* (in au), is represented by the white contours and is calculated for each grid point by minimizing the χ² assuming a circular emitting area of πR². The χ² map indicates that the SO₂ emission has an excitation temperature of 104 ± 5 K and is marginally optically thick (τ ~ 0.1) with an emitting area of ~ 85 au in radius.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. A.3.** Integrated intensity maps of the HDO 3<sub>3.1</sub> – 4<sub>2.2</sub> (*E*<sub>up</sub> = 335 K), CH<sub>3</sub>OH 16<sub>1.15</sub> – 15<sub>2.14</sub> (middle, *E*<sub>up</sub> = 332 K), and CH<sub>3</sub>OCHO 27<sub>10.17</sub> – 26<sub>10.16</sub> (right, *E*<sub>up</sub> = 290 K) transitions in color. The images are integrated over [-2,2] km s<sup>-1</sup> with respect to the V<sub>lsr</sub> of 6.7 km s<sup>-1</sup>. The white vertical bar in the colorbar on top of each image indicates the 3σ threshold. The 0.875 mm continuum is overlaid as black contours. The main continuum peak is associated with IRAS2A1 and the secondary peak toward the south with IRAS2A2. The directions of the two outflows originating from IRAS2A1 and IRAS2A2 are indicated with the colored arrows in the middle panel (Tobin et al. 2015). The size of the beam is shown in the bottom right, and in the bottom left a scale bar is displayed.</figcaption></figure>
**Appendix B: Additional MIRI-MRS figures**

<figure><img src="image_11.png" /><figcaption>**Fig. B.1.** Integrated intensity map of the CO$_2$ $v_2 = 1-0$ Q-branch (top left), H$_2$O $v_2 = 1-0$ 4$_{2,3} - 5_{3,2}$ (top right), H$_2$ 0-0 S(1) (bottom left), and H$_2$ 0-0 S(7) (bottom right) observed with MIRI-MRS in color. The images are integrated over [14.9,15.0] µm for CO$_2$, [7.145-7.15] µm for H$_2$O, and [-0.01,0.01] µm with respect to the transitions of H$_2$. The extent of the continuum around the respective wavelengths is overlaid as white contours. A white scale bar is displayed in the bottom left of each panel, and the size of the PSF is presented as the filled white circle in the bottom right. The directions of the two outflows originating from IRAS2A1 and IRAS2A2 are indicated with the colored arrows in the top-left panel (Tobin et al. 2015). The H$_2$O emission is peaking on-source on the same scales as SO$_2$, whereas the CO$_2$ and H$_2$ 0-0 S(1) emissions are clearly peaking in the outflow toward the southwest. The H$_2$ 0-0 S(7) line peaks mostly on-source but also has a component in the southwestern outflow.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. B.2.** Spectrum of IRAS2A in blue centered around the *v*₃ band of SO₂ as observed by JWST/MIRI-MRS. The local continuum fit is shown in red and the best-fit LTE slab model in orange at a +5 mJy offset. In green, the SO₂-subtracted data are presented, revealing the 7.4 µm ice absorption feature that was hidden by the *P*-branch of the SO₂ *v*₃ band. The *Q*-branch is slightly overfitted, resulting in a subtraction residual. The ice absorption bands are further analyzed by Rocha et al. (2024).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. B.3.** JWST/MIRI-MRS continuum-subtracted spectrum (blue) centered on the *v*₂ band of SO₂ around 19 µm. No clear emission or absorption features of SO₂ are detected. Overlaid are LTE slab models obtained using the best-fit parameters derived from the *v*₃ band around 7.35 µm (red) and the pure rotational lines in the ALMA data (green). The MIRI slab model is scaled down by a factor of 10⁻³ for clarity. Both models are offset with respect to the data and are corrected for an extinction of A<sub>V</sub> = 55 mag (Rocha et al. 2024) using a modified version of the McClure (2009) extinction law (see Appendix C).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. B.4.** JWST/MIRI-MRS continuum-subtracted spectrum (blue) centered on the *v*₁ band of SO₂ around 8.5 µm. The two overlapping MIRI-MRS sub-bands (channel 2A and 2B) have been stitched for clarity. No clear emission or absorption features of SO₂ are detected. Overlaid are slab models obtained using the best-fit parameters derived from the *v*₃ band around 7.35 µm (red) and the pure rotational lines in the ALMA data (green). Both models are offset with respect to the data and are corrected for an extinction of A<sub>V</sub> = 55 mag (Rocha et al. 2024) using a modified version of the McClure (2009) extinction law (see Appendix C).</figcaption></figure>
## Appendix C: Extinction correction

A modified version of the extinction law of McClure (2009) was created for the extinction correction because the depths of the dominant absorption features (e.g., H₂O, CO₂, and silicates) did not match with those measured toward IRAS2A. The total extinction (in units of optical depth) toward IRAS2A ($\tau_{tot}(\lambda) = A_\lambda/1.086$) was therefore decomposed into two components,

$$ \tau_{tot}(\lambda) = \tau_{ice,silicate}(\lambda) + \tau_{ext}(\lambda), \qquad (C.1) $$

where $\tau_{ice,silicate}(\lambda)$ is the differential extinction caused by the ice and silicate absorption features and $\tau_{ext}(\lambda)$ is the absolute extinction.

The differential extinction $\tau_{ice,silicate}(\lambda)$ was obtained in a similar matter to what is used for typical ice analysis studies. The optical depth was computed with respect to a third-order polynomial fitted through obvious absorption-free wavelengths (i.e., 5.1, 5.3, 7.6, 22, and 24 µm; see the top panel of our Fig C.1; Rocha et al. 2024). The optical depth could then be calculated via

$$ \tau_{ice,silicate}(\lambda) = -\ln\left(\frac{F_{\lambda}}{F_{\text{cont}}}\right), \qquad (C.2) $$

where $F_\lambda$ is the measured flux with MIRI-MRS and $F_{cont}$ the continuum flux. The differential extinction in units of optical depth is presented in the bottom panel of Fig C.1.

The absolute extinction $\tau_{ext}$ was computed by fitting a power-law model to the extinction law of McClure (2009) for an extinction of $A_K > 1$. Only the wavelengths outside of the major absorption features were taken into account (i.e., 4.9-5.4, 7.2-7.3, and 28-29 µm). In this case, the power of the power law was also fitted. The extinction law was fitted for the case of $A_V = 1$ mag (assuming $A_K = A_V/7.75$ mag; McClure 2009) in units of optical depth ($\tau_{ext}(\lambda) = A_\lambda/1.086$). The fit to the extinction law is presented in Fig. C.2. Using this fit, the absolute extinction was derived to be

$$ \tau_{ext}(\lambda) = 0.085\lambda^{-0.25}A_V. \qquad (C.3) $$

Here, $A_V = 55$ was adopted based in the depth of the silicate absorption (Rocha et al. 2024), leading to $\tau_{ext} \sim 3$ around 5 µm and $\tau_{ext} \sim 2$ around 25 µm (see bottom panel of Fig C.1). The total extinction toward IRAS2A in units of optical depth is presented in the bottom panel of Fig. C.1.

<figure><img src="image_15.png" /><figcaption>**Fig. C.1.** Determining the extinction due to ice and silicate absorption. Top: Observed MIRI-MRS spectrum of IRAS2A (blue) with the global continuum estimate (red) overplotted. The continuum is based on a third-order polynomial fit to emission- and absorption-free wavelength ranges. The positions of the H₂ v = 0 rotational transitions are indicated with the vertical dashed gray lines. Bottom: Differential extinction in units of optical depth (blue) of the absorption features as a function of wavelength derived from the global continuum estimate (Eq. (C.2)). Overplotted is the absolute extinction (orange) computed using Eq. (C.3); see the main text for more details. The total extinction (green) is computed from the sum of the two (Eq. (C.1)).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Fig. C.2.** McClure (2009) extinction law for $A_V$ = 1 mag (assuming $A_K = A_V/7.75$ mag) in blue with the best-fit power-law model to the absorption-free regions (dots) overplotted in red. Both the extinction law and the best-fit model are plotted in units of optical depth for $A_V$ = 1 mag but are directly proportional to $A_V$.</figcaption></figure>

# Appendix D: H₂ analysis

## D.1. Fitting of the H₂ lines

The MIRI range covers the $\nu = 0$ pure rotational lines of H₂ from the S(0) line at 28.22 µm till the S(8) line at 5.05 µm. However, the S(0) line lies at the very edge of channel 4C ($\lambda > 24$ µm) where the sensitivity drops significantly and the flux calibration is very poor. Some weak line emission appears to be present (see Fig. D.1, bottom-middle panel) but since the peak intensity is equally strong as that of residual fringes present surrounding the line, it was excluded from the rotational diagram analysis. All other H₂ $\nu = 0$ transitions are detected except for the S(8) transition. The detected transitions were fitted with a simple Gaussian emission profile (see Fig. D.1). A flux calibration uncertainty of 5% was assumed. The H₂ line fluxes were corrected for total extinction based on a modified version of the extinction law of McClure (2009) introduced

in Appendix C The integrated line fluxes are reported in Table D.1.

## D.2. Fitting the rotational diagram

A rotational diagram was created following the formalism of Goldsmith & Langer (1999) and presented in Fig. D.2. The rotational diagram can be best fitted using a two-component model: a warm component with a rotational temperature of $T_{\text{warm}} = 356 \pm 41$ K and a column density of $N_{\text{warm}} = 5.6 \pm 2.9 \times 10^{23}$ cm$^{-2}$ and a hot component with $T_{\text{hot}} = 902 \pm 158$ K and a column density of $N_{\text{hot}} = 4.4 \pm 5.1 \times 10^{21}$ cm$^{-2}$. The total number of molecules can be constrained to $N_{\text{H}_2} = 1.7 \pm 0.9 \times 10^{55}$ molecules. An emitting area with a diameter of $1.4''$ is used (i.e., $R_{\text{source}} = 205$ au at a distance of 293 pc), equal to the size of the aperture at 7.35 µm. However, since the H₂ emission is assumed to be optically thin, the value of $N_{\text{tot}}$ does not depend on the assumed emitting area.

<table><caption>**Table D.1.** Integrated line fluxes of H₂ transitions.</caption><thead><tr><th>Transition</th><th>λ<br>μm</th><th>E<sub>up</sub><br>K</th><th>g<sub>u</sub></th><th>F<sub>line</sub><br>W m<sup>-2</sup> arcsec<sup>-2</sup></th><th>τ<sub>tot</sub><sup>(1)</sup></th><th>F<sub>line</sub> × exp (τ<sub>tot</sub>)<br>W m<sup>-2</sup> arcsec<sup>-2</sup></th><th>N<sub>u</sub><br>cm<sup>-2</sup></th></tr></thead><tbody><tr><td>0-0 S(7)</td><td>5.511</td><td>7196.7</td><td>57</td><td>4.5±1.4(-17)</td><td>2.68</td><td>4.3±1.3(-16)</td><td>3.2±1.6(18)</td></tr><tr><td>0-0 S(6)</td><td>6.109</td><td>5829.8</td><td>17</td><td>9.9±3.4(-18)</td><td>4.44</td><td>5.5±1.9(-16)</td><td>7.8±4.3(18)</td></tr><tr><td>0-0 S(5)</td><td>6.910</td><td>4586.1</td><td>45</td><td>4.5±0.4(-17)</td><td>4.25</td><td>2.1±0.2(-15)</td><td>6.5±1.9(19)</td></tr><tr><td>0-0 S(4)</td><td>8.025</td><td>3474.5</td><td>13</td><td>6.3±0.5(-17)</td><td>3.29</td><td>1.1±0.2(-15)</td><td>9.1±2.6(19)</td></tr><tr><td>0-0 S(3)</td><td>9.665</td><td>2503.7</td><td>33</td><td>3.1±1.3(-17)</td><td>6.39</td><td>1.2±0.5(-14)</td><td>3.2±2.0(21)</td></tr><tr><td>0-0 S(2)</td><td>12.279</td><td>1681.6</td><td>9</td><td>2.3±0.6(-17)</td><td>5.94</td><td>5.6±1.4(-15)</td><td>6.7±3.1(21)</td></tr><tr><td>0-0 S(1)</td><td>17.035</td><td>1015.1</td><td>21</td><td>1.8±0.4(-16)</td><td>4.08</td><td>7.0±1.4(-15)</td><td>6.8±2.7(22)</td></tr></tbody></table>

**Notes.** *a(b)* means $a \times 10^b$. <sup>(1)</sup>Total extinction in units of optical depth at the respective wavelength $\lambda$ computed using Eq. (C.1). The line optical depths of the H₂ transitions themselves are all small ($\tau_{\text{line}} \ll 1$).
<figure><img src="image_17.png" /><figcaption>**Fig. D.1.** Observed MIRI-MRS spectrum (blue) centered around $H_2 \nu = 0$ pure rotational lines from the S(7) line in the top-left panel to the S(0) line in the bottom-middle panel. The best-fit Gaussian model is presented in orange in each panel. All presented lines are considered detected except for the S(0) line since the flux is similar to that of several strong residual fringes located in the vicinity of the S(0) line.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. D.2.** Rotational diagram of $H_2$ derived from the MIRI-MRS data. The black data points indicate the measured fluxes that are corrected for extinction (see Table. D.1). The best-fit two-component model is presented as the solid green line, and the two individual components are shown as the dashed blue (warm) and orange (hot) lines. The best-fit parameters are displayed in the top right.</figcaption></figure>