© The Authors 2024

1 Introduction

Photodissociation regions (PDRs) are ideal targets for probing strong variations in the conditions of the interstellar medium (ISM). They appear at the borders of clouds that are illuminated by far-ultraviolet (FUV) photons originating from massive young stars that formed in the cloud, or from the interstellar radiation field (ISRF). The transition from the H II region to the PDR regime starts where the extreme ultraviolet (EUV, >13.6 eV) radiation is no longer sufficient to maintain the ionization and FUV radiation (<13.6 eV) drives the local physics instead. The PDR physics is driven by the balance between the FUV flux and the attenuation by the medium, resulting in sharp spatial transitions between ionized, atomic, and molecular hydrogen layers. When PDRs are located in nearby clouds, they offer the observational opportunity to spatially resolve the transitions between these layers in a single object (Hollenbach & Tielens 1997; Wolfire et al. 2022).

The Orion Bar is a frequently studied PDR (e.g., Tielens et al. 1993; Hogerheijde et al. 1995; van der Werf et al. 1996; Kassis et al. 2006; Pellegrini et al. 2009; Goicoechea et al. 2016) with a density of between several 10⁴ and 10⁶ cm⁻³, which is typical for high-density PDRs. It is subject to a high-intensity FUV radiation field with Go in the range of (2.2–7.1) × 10⁴ (Berné et al. 2022; Habart et al. 2024; Peeters et al. 2024), where G₀ represents a dimensionless parameter defining the FUV intensity with respect to the ISRF intensity defined by Habing (1968). The proximity (414 pc, Menten et al. 2007) and very high surface brightness of the Bar allow infrared observations at excellent signal-to-noise ratio (S/N) levels. Recent observations at high spatial resolution of the HCO⁺ J = 4−3 line (Goicoechea et al. 2016) and vibrationally excited H₂ lines (Habart et al. 2023) revealed the existence of many filamentary substructures at the surface of the Bar. Underlying these substructures, the Bar has a simple and approximately linear geometry, a perpendicular orientation with respect to the radiation field, and a nearly edge-on orientation with respect to the observer, meaning that it has served as an exemplary benchmark target for theoretical PDR models (e.g., Tielens et al. 1993; Jansen et al. 1995; Young Owl et al. 2000; Shaw et al. 2009; Joblin et al. 2018).

The PDRs4all Early Release Science program for JWST (ERS 1288; Berné et al. 2022) observed the Orion Bar with NIRCam and MIRI imaging, as well as integral field unit (IFU) spectroscopy with NIRSpec and MIRI medium-resolution spectroscopy (MRS). The imaging data have led to new insights in the detailed geometry of this PDR, and suggest a terrace-like structure with three consecutive dissociation fronts as described in Habart et al. (2024). The rich near-infrared (NIR) to mid-infrared (MIR) spectroscopic data produced by this program are discussed in a series of papers focusing on the NIRSpec IFU observations (Peeters et al. 2024), the aromatic infrared bands (AIBs) in both the NIRSpec and MIRI data (Chown et al. 2024), the spatial variations of dust properties (Elyajouri et al. 2024), and the spectra from a dense disk in the same field of view as the PDR (Berné et al. 2024). In the present work, we focus on the emission lines observed in the MIRI spectra, which cover the 4.9–28 µm range.

Emission lines in the MIR are key probes of the local conditions of the gas in PDRs and the surrounding environment. Constraints on the physical conditions of the gas can be obtained by direct observations of the MIR lines, and by comparing the intensity of the lines to theoretical models. These diagnostics can be better calibrated by studying these lines towards PDRs, where the ionized, neutral atomic, and molecular layers are spatially resolved. They can then be used for studies of distant star-forming galaxies, where PDRs contribute significantly to the emission (Wolfire et al. 1990; Farrah et al. 2007; Cormier et al. 2012), but where the ISM regimes cannot be spatially resolved. For example, the MIR fine-structure lines of different ionization stages of metals ([Fe II], [Ar III], [S IV], etc.) and the hydrogen recombination lines can be used to constrain the gas-phase abundances, excitation conditions (e.g., [Fe II] lines), and radiation field hardness (e.g., [Ne III] 15.5 µm / [Ne II] 12.8 µm) within the ionized gas (e.g., Verma et al. 2003; Bernard-Salas et al. 2009; Rubin et al. 2011; Koo et al. 2016). In the warm atomic gas, excitation of [Fe II] can provide diagnostics, while in the warm molecular gas close to the dissociation front in a PDR, the MIR contains rotational H₂ lines, which are key tracers. The ratios between pure-rotational H₂ lines constrain the kinetic gas temperature and density, as well as the local strength of the UV radiation field (Habart et al. 2005; Kaufman et al. 2006; Sheffer et al. 2011). An updated and more versatile “PDR toolbox” was developed in the context of the PDRs4All collaboration (Pound & Wolfire 2023), which enables users to derive constraints from a variety of diagnostics in a straightforward manner. It supports diagnostics in the MIRI wavelength range such as [Fe II] and [Ar III] lines, and the rotational H₂ lines. By providing the line list in the present work and the PDR toolbox to the community, the science products of the PDRs4All project will make future PDR studies based on MIR data more straightforward and consistent.

The MIRI IFU spectra that we focus on in the present article were extracted for five exemplary regions, consistent with the apertures originally defined for the analysis by Peeters et al. (2024). Each of these five apertures captures a specific regime in the PDR front: the ionized region, the atomic region, and three individual dissociation fronts. The scope of this work is limited to these five regions of the Orion Bar, with the main goal being to clarify which emission lines are expected in each regime. The data reduction and extraction of the spectra are described in Sect. 2. In Sect. 3, the identification, intensity measurements, and central wavelength measurements of the observed emission lines are described, and the results are presented in the main table of this work. Future work is expected to make full use of the spatial resolution of JWST, and the high S/N of these data so that analyses at single-spaxel resolution (0.2″) are possible. The provided line list will enable such efforts for both the PDRs4All Orion Bar data and future data of other objects dominated by PDR emission. In Sect. 4, we develop an analysis concentrating in particular on the H I recombination lines emitted in the ionized region, the pure-rotational H₂ lines of the three dissociation fronts, and a number of metal fine-structure lines. This initial analysis summarizes the nature of the gas conditions that produce the rich spectral contents of the MIRI MRS observations.

Fig. 1 Examples of flux differences in overlapping regions between MIRI bands. Blue and orange curves show individual segments, with the bands indicated at the top of each panel. These were extracted from the band cubes produced by the pipeline, in this example for the H II template aperture. The black line shows the stitched spectrum, which is offset with respect to the blue and orange lines due to the additive offsets applied to the 12 segments, which are cumulative with Ch 2 LONG used as a reference.

2 Data

2.1 Reduction

The MIRI MRS observations were processed from uncalibrated data through the official pipeline, using 1.12.5 of the jwst Python package. The pipeline and reference files from the Calibration Reference Data System (CRDS) have been continually updated since commissioning, and at the time we performed our latest data reduction we used the CRDS context jwst_1147.pmap. The wavelength calibration at the time of writing is accurate up to a few km s⁻¹ at short wavelengths, and about 30 km s⁻¹ at the longest wavelength (Argyriou et al. 2023). The spectral resolution ranges from R ~ 3500 to 1500 depending on the wavelength, or about 85–200 km s⁻¹ in terms of radial velocity (Labiano et al. 2021). Concerning the astrometric calibration, the pointing accuracy is about 0.45″ without target acquisition, and the typical accuracy of the assigned coordinate system is about 0.3″, where the main source of uncertainty is the guide star catalog (Patapis et al. 2024). For reference, the spaxel size varies between 0.2 and 0.3″. The spectrophotometric calibration was initially based on a single standard star (Gordon et al. 2022), but data from additional stars were recently introduced. The inclusion of these additional calibration data (since jwst_1094.pmap), led to a much improved matching of the continuum flux between the four channels of MIRI MRS, with only minor flux offsets in the overlap regions between the channels (see also Sect. 2.2 and Fig. 1). Together with support for time-dependent calibrations (since jwst version 1.11.0), the changes much improved the Channel 4 calibration in the case of the Orion Bar data, removing most local flux oscillations or broad artifacts so that the hot dust and telescope continuum resemble a smooth blackbody spectrum. As reported in the JWST user documentation¹, there may remain a 10% systematic uncertainty in the calibration, based on the deviations that are observed when using different calibration stars. There is also a spectral leak artifact, where second-order light at 6 µm arrives at a different part of the detector, causing a spurious signal at 12.2 µm in the spectra. We did not correct for this leak, as it is broad and does not interfere with our measurements of the emission line intensity.

We made additional corrections to the default pipeline steps to improve the quality of the final products. In the following section and in Fig. 1, we refer to the four channels of MIRI MRS as Ch 1, 2, 3, or 4, and to the band setting for each channel as SHORT (S), MEDIUM (M), or LONG (L). In the first stage of the pipeline (Detector1), which applies detector-level corrections to the exposures, we modified the threshold for the jump step to 3.0 standard deviations, instead of the default value of 4.0, to improve the cosmic ray detection. One of the corrections in the second stage of the spectroscopic pipeline (Spec2), is removing the fringes resulting from coherent reflections within the detectors (Argyriou et al. 2020). The default fringe correction (fringe) divides the images by a static extended source fringe flat and removes most of the fringe amplitude, while residual fringing remains primarily in Ch 2 LONG, which is known to have a poor fringe mitigation at this date. Therefore, we also applied the optional residual_fringe step of Spec2, which further reduces fringing by fitting sinusoidal functions at the detector level, and works with constrained spatial frequencies in order to avoid removing physical features.

The Spec3 stage then combines the multiple exposures (dither pattern and mosaic positions of both science and background), and builds a 3D cube as the final product. The dedicated background exposures of the program were reduced up to Spec2 using the same workflows described above, to obtain the 1D collapsed spectra produced through the extract_1d step at the end of this pipeline stage. The 1D background spectra were then used in the master_background step of Spec3, as recommended for bright extended sources (as opposed to a pixel-based background subtraction available in Spec2). We applied the outlier_detection step as it did not introduce new artifacts and resulted in an improved identification and removal of the remaining cosmic rays and bad pixels. Moreover, the tweakreg step was skipped and we took advantage of the MIRI images taken in parallel with the MRS observations to improve the world coordinate system (WCS) alignment. Because the parallel imaging was performed simultaneously, the telescope had an identical orientation during both the MIRI MRS observations and the MIRI imaging. The larger fields of view of the images allow for the use of the Gaia DR3 astrometric catalog (Gaia Collaboration 2023) for WCS correction purposes, by performing a source detection in the F770W images. Finally, we configured the cube_build step of Spec3 to produce the 12 individual data cubes corresponding to each of the four channels channels and three bands, to better correct for any offsets between the bands as described in the next section.

Fig. 2 Apertures from which the five template spectra were extracted, overlayed on two slices of the MIRI MRS data cube. The horizontal axis is aligned with the y-coordinate of the IFU and the 9×1 mosaic strip. The blue component is the total flux at 11.3 µm (104 MJy sr−1), which contains some continuum, but mainly traces the 11.3 µm AIB and peaks in the atomic region. The red component highlights the locations of the dissociation fronts using the intensity of the H2 S (1) line at 17.0 µm(10−3 erg s−1 cm−2 sr−1). The two-dimensional color legend shows where one component (red and blue) or both components (dark purple) are bright.

2.2 Extraction and additional corrections

We focus on five spectra extracted from the IFU data cubes, obtained by averaging the flux over apertures covering different spatial locations. The chosen apertures are defined by Peeters et al. (2024), and capture the five main zones representative of the physics observed in this section of the Orion Bar (Fig. 2). The resulting spectra are called the “template spectra” in this work as well as Peeters et al. (2024). We refer to the zones as “HII” for the HII region, “Atomic” for the atomic zone near the ionization front, and “DF1”, “DF2”, “DF3” for the “dissociation fronts”, the three H I/H₂ transitions that appear as bright filaments in maps of the H₂ emission. It should be noted that the S/N and resolution of the data are high enough to apply line diagnostics in a detailed spaxel-per-spaxel manner that spatially resolves the structures in the PDR. The scope of this work is focused on making a line list based on the template spectra, and showing a first analysis that demonstrates what future highly resolved analyses could include. For a detailed description of the contents of the five apertures and a discussion of the PDR structure based on the imaging and NIRSpec IFU data, we refer the reader to the respective papers of Habart et al. (2024) and Peeters et al. (2024).

Some remaining oscillations in the spectrum are an artifact of undersampling issues and the cube building method (Law et al. 2023). For the purposes of detecting individual AIB emission features, Chown et al. (2024) applied an additional empirical correction based on data of the calibration star 10 Lac from the CALSPEC archive (Bohlin et al. 2014). This further reduced the fringes, but also introduced artifacts that resemble lines or cause lines to exhibit double peaks. Because these artifacts interfere with the identification of real lines, we do not use the spectra of Chown et al. (2024) for this work.

Even with the absolute flux calibration, the many corrections in the pipeline, and the stacking of the spectra over the template apertures, there are still minor flux mismatches in the overlap regions of the 12 spectral segments produced by MIRI MRS. We measured each of the 11 jumps in the continuum, by computing the median for each pair of neighboring segments over the region where they overlap, and taking the difference. The typical values for the jumps are no more than 5%, or a few 10 MJy sr⁻¹. These jumps in the flux density do not affect the extraction of the line intensities, as we apply a local continuum measurement for each line individually (see Sect. 3.2). However, to improve the continuum for lines that are located in the overlap regions and to deliver higher quality template spectra, we perform a spectral stitching correction.

To match the continua over all 12 segments, we compute 12 absolute offsets by taking the cumulative sum of the 11 relative offsets calculated above. In this calculation, Ch 2 LONG is used as a reference segment, and we set its offset to zero by construction. For the other segments, the offsets are applied as an additive constant to the flux. While the individual jumps are rather small, the cumulative offsets do result in some corrections ≥100 MJy sr⁻¹, especially in Ch 1 and Ch 4. Three examples of the overlap regions and the stitching method are shown in Fig. 1. To merge the continuum-matched segments, a sliding weighted average is used for each overlap region for a gradual transition. The segment weights change linearly from 1 to 0 for the blue-ward segment, and from 0 to 1 for the redward segment, as the wavelengths are processed from the blue end to the red end of the overlap region. For the overlap regions between channels, the wavelength resolution of the longer channel is always lower. In this case, the grid of the short wavelength segment is used, and the long-wavelength segment is interpolated onto this grid before the sliding average is applied.

To determine if any multiplicative corrections might be needed because of potential calibration issues, we measured the intensities of lines that are present in the overlap regions, and inspected the ratio of the intensities obtained individually from each pair of segments. We selected the H II template spectrum for this task, since the emission lines in the overlap regions are brightest there. For the line at 5.711 µm (Fig. 1), the intensity ratio of Ch 1 M to Ch 1 S is 0.95. Similarly, we find a ratio of 1.02 for the line at 6.565 µm (Ch 1 L/Ch 1 M), and 0.93 at 8.760 µm (Ch 2 M/Ch 2 S). The Ch 3 L to Ch 4 S transition has a continuum that matches well without any correction (right panel of Fig. 1), but there are two lines present at 17.883 and 17.933 µm for which the ratios are 0.95 and 0.99. The last two intensity ratios are different, even though the lines are present in the same overlap region. These differences indicate that a mix of multiplicative and additive offsets could be needed, to match both the line intensity and the continuum level. Since not every overlap region has lines, and since a cumulative multiplicative correction across the segments would potentially inflate the fractional systematic uncertainties, we do not apply a multiplicative correction. Instead, we suggest that for a group of lines in one and the same segment, one should consider a fractional systematic uncertainty of a few percent which affects all those lines multiplicatively with the same factor. Taking the standard deviation of the segment-to-segment line intensity ratios measured above, we estimate that the uncertainty on this systematic factor is around 3%.

3 Lines

3.1 Identification

We have composed a list of the lines present in the five extracted spectra, omitting any features which we deem to be artifacts, or where local fluctuations in the spectral baseline make it hard to identify or measure a potential line. To distinguish features from artifacts, we compared the spectrum to observations of the 10 Lac calibration star, and any oscillations or peaks present in this observed spectrum (residual fringes and other artifact patterns) are considered suspicious and not a real feature of the Orion Bar. We also exclude the candidates of which the amplitude is similar to that of the local fluctuations and noise patterns, by measuring how this local noise affects the line measurement (see Sect. 3.2). We used a S /N > 3 as a rule of thumb to include a line in the list, but we made a few exceptions to complete a series of lines for a species, such as the H I lines listed in Sect. 3.1.1.

To identify the lines, our starting point was the line list originally created for the PDRs4all program (Berné et al. 2022; Peeters et al. 2024), available as a Science Enabling Product on the PDRs4All website². This line list was composed by combining results from the Cloudy ionization model (Ferland et al. 2017), and the Meudon PDR model (Le Petit et al. 2006). These models also predict the strength of the lines, and we used this information to filter out lines deemed too weak for detection, which is those with a predicted brightness lower than 10⁻⁷ erg s⁻¹ cm⁻² sr⁻¹. This line list was revised based on the spectra obtained as part of the observing program, once they were available. We match the simulated lines in the list with the central wavelengths we observe in the five template spectra. As an extra check for the accuracy of the identifications, we compared the measured intensities and their ratios (Sect. 3.2) to the intensities predicted by a Cloudy model (Sect. 4.1) or a PDR model (Sect. 4.3.3). The primary reason for performing this check was to verify the identity of lines for which the observations reveal a small wavelength shift compared to the observations (see also Sect. 3.3). For example, when the wavelength of a suspected H recombination line is slightly different from the theoretical value, the identification can be deemed more confident when the intensity ratios of the H lines match the theory. For the H₂ lines, the order of magnitude was checked against the PDR model included in the line list. A ratio mismatch was found in a few cases, for example due to the presence of overlapping lines. Overview plots of the spectrum and the line identifications are shown in Fig. 3. The wavelengths and transitions corresponding to the indicated lines were compiled in the main table of this work, which also lists the energy levels, Einstein coefficients, intensities with uncertainty estimates, and observed central wavelengths (Table A.1).

3.1.1 Hydrogen recombination lines

Roughly half of the identified lines are H I recombination lines that originate from the H II region between the star and the PDR (H II template), and the ionized layer that borders the PDR lying between the observer and the neutral gas (Atomic and DF templates). In the MIRI wavelength range, atomic hydrogen transitions with a lower principal quantum number n_l of 5 or higher can be observed, with the series of higher n_l having more but fainter lines. The brightest H I line is the 5-6 transition (shorthand for n = 5 ← n = 6 or Pfund α) at 7.46 µm. For the Humphreys series (n_l = 6), we observe all four bright lines, with the upper state n_u ranging from 10 (5.13 µm) to 7 (12.37 µm). All 15 available lines with n_l = 7 are observed from 7 to 23 (4.92 µm) to 7–8 (19.062 µm).

The remaining H I lines we observe have n_l ranging from 8 to 10, and are significantly fainter. Because of variations in the noise and artifacts throughout the spectra, some of the fainter lines are not detected or not reliably measured. We detect the n_l = 8 series from 8–29 (6.316 µm) to 8–10 (16.21 µm). The 8–12 (10.503 µm) transition was not measured because it is hidden in the wing of the extremely bright [S IV] 10.510 µm transition. The first transition of this series, 8–9 at 27.80 µm, is not detected because the current MIRI calibration is not accurate enough to distinguish it from the very high continuum flux at that wavelength. We measured the n_l = 9 series from 9–27 (8.308 µm) to 9–11 (22.34 µm), where the 9–10 line is not available as it is outside the MIRI MRS wavelength range. Finally, were report the n_l = 10 series from 10–15 (16.41 µm), up to 10–20 (12.16 µm). The 10–14 and 10–13 transitions are in the wavelength range but we do not list them, since they are too weak to be reliably measured, and because 10–13 (22.33 µm) overlaps with 9–11 at 22.34 µm while the spectral resolution is not sufficient to separate their contributions. Several transitions with higher n_u are still detectable, but their low intensity makes many of them too faint to measure at a S/N > 3, and we opt to omit the series beyond this point.

Because of the lower transition rates and more sparsely populated upper levels, some H I lines are faint enough to be affected by the data quality, typically when the intensity ≲10⁻⁵ erg s⁻¹ cm⁻² sr⁻¹. The intensity and uncertainty measurements are described in Sect. 3.2. Care has to be taken with 6–8, 7–11, and 8–16, which form a closely spaced trio near 7.5 µm. Similarly, 7–9, 8–23, 8–22, and 8–17 are in close proximity to other lines, and for such cases we employed a decomposition as explained in Sect. 3.2, with an example shown in Fig. 4.

3.1.2 Molecular hydrogen lines

The brightest H₂ lines are the pure-rotational lines (J → J − 2) in the vibrational ground state (v = 0). We observe all the pure-rotational transitions available in the MIRI wavelength range, which range from lower rotational quantum number J = 1 to J = 8, in other words all eight transitions denoted as 0–0 S (1) to 0–0 S(8). Only the S (0) line at 28.22 µm is missing. The calibrated wavelength range of MIRI MRS extends just far enough to include the wavelength of this line, but the continuum level is very high and the S/N and is too low to enable a detection. The intensity of these lines is several 10⁻³ erg s⁻¹ cm⁻² sr⁻¹, and these values are used in the H₂ excitation analysis presented in Sect. 4. We also detect pure rotational lines in the v = 1 vibrationally excited state, from 1–1 S (5) to 1–1 S (9). While transitions down to 1–1 S (1) lie within the available wavelength range, they are too faint to detect. The NIRSpec observations of this program show the continuation of the v = 0 and v = 1 ladders towards higher J, up to 0–0 S (19) and 1–1 S (17) (Peeters et al. 2024).

There are also a few v = 1−0 rovibrational transitions in the wavelength range that have the same upper level as the v = 1 − 1 pure rotational transitions mentioned above. For example, 1–1 S (6) (v = 1, J = 8 → v = 1, J = 6) has the same upper level as 1-0 O(10) (v = 1, J = 8 → v = 0, J = 10). We inspected the shape of the spectrum in detail, at the wavelengths of the 1− 0 O(10) to 1−0 O(13) transitions (5.0590, 5.6294, 6.3086, and 7.1267 µm, respectively). There is a tentative detection for each of these transitions upon inspection by eye, but the lines are not prominent enough for inclusion in Table A.1 due to the smaller Einstein coefficient. Another tentative detection, is the rovibrational line 4−3 O(8) at 5.0981 µm, while the other vibrationally excited series such as 3−2, 2−1, 3−3, or 2−2 are weak in the MIR, and rather expected at shorter wavelengths (NIRSpec data or Kaplan et al. 2017).

In work by Rosenthal et al. (2000), unresolved observations of the Orion Molecular Cloud outflow (OMC-1) with the Short Wavelength Spectrometer (SWS) of the Infrared Space Observatory (ISO), detected the bright pure-rotational emission driven by shocks, from 0–0 S (1) to S (25) and from 1–1 S (5) to S (17). Observations of the Orion Bar with ISO-SWS were obtained by F. Bertoldi (priv. comm.) and used by Habart et al. (2004) and Joblin et al. (2018), and these detected the pure-rotational lines from 0–0 S (0) to S (5). Previous deep spectroscopic observations in the NIR revealed a large set of rovibrational lines with v = 0 to v = 10 in the Orion Bar, using the Immersion Grating Infrared Spectrometer of McDonald Observatory Kaplan et al. (2017). Spatially resolved observations of the MIR pure-rotational lines were previously obtained through ground-based observations, where the 0–0 S (1), S (2), and S (4) lines were mapped at a 2″ resolution (Allers et al. 2005). The JWST data of our program provide the first maps of the Orion Bar for a complete set of H₂ lines, extending well past 0–0 S (5), at a resolution that resolves the individual dissociation fronts. An in-depth spatially dependent analysis of the H₂ excitation from the MIR to NIR will be presented by Sidhu et al. (in prep.).

Fig. 3 Overview of template spectra and lines, continued in Fig. 3. Before plotting, the segments were normalized individually, each divided by a constant factor, which is the average flux of that spectrum as integrated over the wavelength range of the plot. Then, multiplicative offsets using a factor of 1.5 were applied for clarity. Vertical lines and legend: Identified lines and species. Dashed black lines: Unidentified lines (UID). Sparse dotted black lines: End of a MIRI MRS band, where changes in the spectral resolution occur. We note a peculiar detail, namely broad absorption or emission or artifact features around 7.0 µm near the Aril line, which remain unidentified. The oscillations in the top panel are a good example of remaining fringes even after the fringe flat and residual fringe corrections in the pipeline.

Fig. 4 Examples of intensity determination for overlapping or nearby lines, in H II template spectrum. The x-axis has been shifted so that the zero point matches the wavelength in µm displayed in the top left corner of each panel. Blue curve: Fitted model consisting of a linear continuum and two or three Gaussian profiles. Dotted lines: central wavelengths of the Gaussian profiles. Upper left panel: [Ni II] 6.636 µm and H I 6.638 µm. Lower left panel: complex of three H I lines, with fitted wavelengths at 7.4955, 7.5030, and 7.5086 µm. Right panel: H I 5–6 and contributions by He I multiplets in its left wing.

3.1.3 Fine-structure lines

All the metal lines we detect are forbidden fine-structure lines. Based on ISO-SWS observations, Rosenthal et al. (2000) separated which lines originate from the foreground ionized and PDR region, and those from an outflow of OMC-1 excited by powerful shocks. Their Table 2 lists the lines for the Orion Bar PDR, and we detect the same set of lines with a few exceptions. The MIRI wavelength range does not extend to [Si II] 34.8 µm, but we do detect [S I] 25.2 µm in DF3, and the detailed spatial distribution of sulfur is discussed in a separate paper (Fuente et al. 2024). The newly detected and identified lines compared to these ISO data are [Ni II] 10.68 µm, [Ni III] at both 7.35 µm and 11.00 µm, and [Cl II] 14.37 µm. The Cl II line was also detected by Rubin et al. (2011) using the Spitzer Space Telescope. Most of these lines originate only from the ionized layer, while others come from both the ionized and atomic layers between the observer and the molecular part of the Orion Bar. This is discussed further in Sect. 4.2.

The brightest lines have intensities of several 10⁻² erg s⁻¹ cm⁻² sr⁻¹, are the forbidden lines of Ar ii, Ar iii, S iii, S iv, Ne ii, and Ne iii. They are measured without saturation issues and are present in all five apertures. Because of their brightness, they are frequently used as tracers for H ii regions, and as diagnostics for the ionic abundance ratios or the level of excitation for example, including in extragalactic objects (e.g., Rubin et al. 2007). The lines of Fe ii, Fe iii, Ni ii, Ni iii, Cl ii, and P iii have intensities down to 10⁻³ erg s⁻¹ cm⁻² sr⁻¹. Of these, the species with the weakest fine-structure lines are Ni ii and Nim, but they are still detected all five templates. The S/N for most of the weak lines is still > 10, meaning that the metal lines have a very high S/N in these data, and the uncertainty will be dominated by the systematics in the calibration and our spectral order matching factors (see Sect. 2). The exceptions are [Ni iii] 11.0023 µm and [S i] 25.25 µm for which the S /N is lower than 3 for some of the measurements. In Sect. 4.2, we use the brightest Ar, Ne, S, Fe, and Ni lines to trace the variation in radiation field hardness between the five template regions.

A caveat is that for [Ni iii], the wavelength as given by the Atomic Spectra Database (ASD, Kramida et al. 2022) of the National Institute of Standards and Technology (NIST) is 7.3492 µm, while the measured wavelength is 7.3524 µm, which differs by 0.0032 µm. A smaller difference exists for [Ni iii] 11.0023 µm, which is instead observed at 11.0035 µm. It is likely that this offset occurs because the atomic data for these transitions are not sufficiently accurate, since only a Ritz wavelength and no observed wavelength is listed on the ASD. The Ni ionization state is known to be similar to the Fe ionization state (Lutz et al. 1996) because of their very similar ionization potentials (see also Sect. 4.2 and Table 1). Considering that the Fe ii and Fe iii ions are both present, the Ni ii and Ni iii ions should also be present if Ni is abundant enough. When we assign [Ni iii] to the lines at 7.35 and 11.00 µm, the observed intensity for these lines is consistent with the expectations based on our theoretical models. In Sect. 4.2, we confirm the correspondence between the Fe iii Fe ii and Ni iii Ni ii ratios based on our measurements.

The measured wavelength difference for [Ni iii] 7.35 µm is 0.0035 µm, or 0.0032 µm if corrected for the velocity of the H ii region relative to the frame of the calibrated JWST data (see Sect. 3.3). This is much larger than the standard deviation of our H i line wavelength offset measurements of about 0.0003 µm. With the above information, our best empirical estimate for the [Ni iii] 7.35 µm wavelength is 7.3524 ± 0.0003 µm. For the other [Ni iii] line, with the same reasoning and velocity correction as above, the remaining offset is 0.0008 µm and our empirical wavelength estimate is 11 .0031 ± 0.0003 µm. For consistency, we refer to the corrected lines using the wavelength from the ASD in the rest of this work.

3.1.4 Helium recombination lines

While many He i lines are detected upon visual inspection of the H ii template spectrum, most of them are located in the blueward wing of bright H i recombination lines (e.g., in Fig. 4). In most of these cases, the brightness of the H i line is several orders of magnitude higher than the He i line, and we could not obtain reliable measurements. An example is the subtle contribution of¹ F^o ← ¹D at 6.769 µm to the wing of the H i 7–12 line. Some of these transitions are not listed in the ASD, but they are present in the output of the Cloudy model and in the “Atomic Line List” by van Hoof (2018). In the discussion below there are often many He I lines present near the same wavelength, and we mention those with the highest Einstein coefficient.

The brightest H I line (5–6, 7.4599 µm) in our spectrum is contaminated by lines at 7.4561 µm from the ³G ←³ Hº multiplet and the ¹G ←¹ Hº singlet, and a similar group of lines at 7.4538 µm from the ³Fº ← ³G multiplet and¹ Fº ← ¹G singlet. Performing a fit as demonstrated in the right panel of Fig. 4 results in a correction of about ~8% to the intensity of this H i line. There might be a similar contamination for the H I 6–8 line at 7.503 µm; in addition to the other H I 8–17 (7.495 µm) and 7–11 (7.508 µm) lines, the Cloudy model predicts possible contributions by He I multiplets at 7.499 µm and 7.505 µm at the 1% to 10% level. With the mixing of three H I lines and two He I multiplets, there is no straightforward way to extract the small He I contribution from this complex. For this reason we do not report the intensities of most He I lines, and we only use an empirical shape to correct the H I intensity of certain lines, as discussed above.

There are a few exceptions where the He i lines are sufficiently isolated to allow a measurement, as listed in Table A.1. The line at 6.7217 µm is part of a clearly defined double peak together with a weak H I line. It is identified as a He I line based on the Cloudy model, but the Atomic Line List by van Hoof (2018) lists many transitions around the same wavelength, with similar upper level energies, and it is unclear which would be the dominant contribution. Therefore, we note the transition as "multiple" in Table A.1, with an upper level energy close to all the listed upper levels. There is another clearly defined line at 7.4334 µm, where the ³D ← ³Fº multiplet is present, with no interference from other nearby lines, and at 10.88 µm the ³S ← ³Pº multiplet is measured.

3.1.5 Unidentified lines

Several lines for which the origin is unclear are marked as “UID” in Table A.1. We compared the observed wavelengths of these lines to the PDRs4All line list and the ASD for all elements with atomic number <30 and ionization up to stage iv. Some lines have no matches, while others match individual wavelengths in these line lists, many of which are highly excited lines of H₂ lines or atomic recombination lines. None of the candidates were convincing however, and we summarize the ruled-out species and transitions below, per template region type (ionized, atomic, dissociation front), and note which unidentified lines are the brightest in each region.

In the H II template spectrum, the brightest unidentified lines are at 5.240, 5.740, 5.596, and 7.878 µm (in order of decreasing brightness), with intensities of several 10⁻⁵ erg s⁻¹ cm⁻² sr⁻¹. The 5.240 and 7.878 µm lines only appear in this region and have no matches in any of our line catalogs. The 5.596 µm line is also exclusive to the H II template and is a somewhat close match to the recombination lines of Mg I at 5.5973 µm or C I at 5.5983 µm, but there is still a significant wavelength offset and no recombination lines of similar expected intensity are detected. The 5.740 µm line is visible in all templates and matches Li II 5.7406 µm, which is unlikely because of the same reason.

The weaker unidentified lines in the H II region (several 10⁻⁶ erg s⁻¹ cm⁻² sr⁻¹) are at 5.628, 5.658, 6.366, 6.498, 8.410, 10.386, and 14.333 µm. The line at 6.498 µm is located close to the 6.501 µm H I line. Close wavelength matches from the line lists include a He I recombination line at 6.4984 µm and a H₂ rotational line at 6.4999 µm (v = 3, J = 9 → 7), or a N I recombination line at 6.4985 µm. The identification as a vibrationally excited H₂ line is unlikely because the line is not observed for the DF templates. If the 6.498 µm line matches the N I 6.4985 µm recombination line, this could also explain another unidentified line at 10.386 µm, which could match N I 10.385 µ m. Since none of the other N I recombination lines are detected, the remaining option is a tentative detection of the He I 6.4984 µm line. Similarly the lines at 6.366 and 8.410 µm appear only in the ionized region, so coincidentally matching H₂ lines are unlikely candidates, because they do not appear in any of the molecular template spectra. The 6.366 µmline is close to [Ar III] 6.368 µm, but the latter has the same upper level as [Ar III] 8.99 µm, while the Einstein coefficient is smaller by a factor ~10⁴, and the observed flux ratio does not match this expected ratio. The 8.410 µm line matches very well in wavelength with C I 8.4095 µm, but again we do not see any C I recombination lines of similar excitation that are expected to be similarly strong in the MIR spectra. We note that C I and N I recombination lines have been observed in the PDRs4All NIR-Spec spectra (Peeters et al. 2024). For C I, the lines at 1.069 and 1.175 µmare measurable, but these have lower upper levels and an order of magnitude higher Einstein coefficients compared to those in the MIR, and the situation is analogous for N I. The line at 14.333 µm appears in all five template spectra, and we find no matches in the catalogs used.

In the Atomic template, the brightest lines are at 6.859, 6.498, 5.240, and 5.740 µm (about 10⁻⁵ erg s⁻¹ cm⁻² sr⁻¹), of which the last three are also present in the H II template as discussed above. The line at 6.859 µm is only present in the atomic region, and its wavelength matches the v = 8, J = 20 → 18 transition of H₂. As this line is not observed in the dissociation front templates, it is unlikely that this is the right identification. There is a N I recombination line nearby at 6.8560 µm, but as with several other lines mentioned above, the wavelength difference remains too large to serve as a convincing identification.

The DF1, DF2, and DF3 templates have a relatively strong feature at 18.651 µm, for which we find no clear identification. Its width is roughly twice that of the nearby S III line, so it is likely not a single atomic line. The dissociation fronts exhibit weak unidentified lines at 4.928, 4.968, 5.552, 5.634, 5.740, and 6.412 µm, of which 5.740 µm also appears in the H II and Atomic templates. For 5.634 µm, the 4−3 O(9) rovibrational transitions is a likely explanation (see also Sect. 3.1.2). Several others coincide in wavelength with a CO line of high excitation (e.g., v = 2, J = 20), or several highly excited H₂ lines in the line list. The lack of other highly rotationally excited H₂ or CO lines rules out those identifications. There is a near perfect wavelength match with an S I line at 6.4126 µm, but this is again an implausible recombination line.

At around 7.0 µm and 9.0 µm, several strong broad features resembling emission and absorption are present. Both of these appear near Ar lines and are particularly strong for the H II template. We compared spectra with a version of the data where no background subtraction was applied and confirmed that these features are not an artifact of the background subtraction. The data of the 10 Lac calibration star (see Sect. 2.2) were also inspected at this wavelength and show no sign of similar effects. At the time of writing it remains unclear if these particular unidentified features are another type of artifact, although it is peculiar that features of this type only appear at wavelengths near the Ar lines.

We note that certain molecules might be candidates in the DF templates. However, we only included H₂, HD, and CO in our search. The work by Zannese et al. (in prep.) investigates other molecular candidates including ${\text{CH}}_3^{+}$, which has already been detected in a disk with much denser gas that is present in our field of view (Berné et al. 2023). A more in-depth and physically motivated search for unidentified line candidates, is deferred to future studies that focus on specific spatial regions and specific groups of unidentified lines.

Fig. 5 Demonstration of uncertainty measurement described in Sect. 3.2. The line in this example is [S I] 25.25 µm, and the oscillations are an artifact of MIRI MRS Ch 4. Blue segments with dots: local linear continuum estimates, with the center indicated. Orange crosses: Blue dot value plus the average of the continuum subtracted spectrum as integrated over the continuum window. The deviation between the orange cross and the blue dot samples the deviation that occurs due to local noise patterns and the linear continuum approximation. The median absolute deviation is used as the uncertainty on the line intensity integral. Blue circle and error bar: Continuum and median absolute deviation of the above samples. Orange error bar: Line intensity and uncertainty, visualized approximately as an amplitude.

3.2 Intensity and uncertainty measurements

For each line, the intensity in each of the five templates is given in Table A.1. The extraction method used for most of the line intensities is a simple local integration over the spectrum, with an estimate of the continuum subtracted. To determine an appropriate integration window and continuum level, we first derive a value for the full width at half maximum (FWHM), using instrumental resolution curves for MIRI adapted from Labiano et al. (2021, Fig. 9). We evaluate the resolution curve at the central wavelength of the line λ, and denote the resulting FWHM as w(λ). After fitting a Gaussian profile to a few lines, we found that the fitted FWHM is very close to the instrumental FWHM, meaning that the lines are unresolved at the MIRI MRS velocity resolution of around 100 km s⁻¹.

The continuum is estimated as a local linear function, by taking the median value of the flux in the two nearby windows [λ − 4w(λ), λ − 2w(λ)] and [λ + 2w(λ), λ + 4w(λ)], and applying a linear interpolation between these two continuum levels. We found that this median-based approach worked best for the large number of lines in our data, as it is more robust against local outliers and does not suffer from overfitting, as opposed to fitting a polynomial which would require manual adjustment for individual lines. The estimated continuum levels were inspected by eye and proved reasonable for most lines. The line intensity (in erg s⁻¹ cm⁻² sr⁻¹) is then obtained by converting the flux to per-wavelength units (erg s/cm²/sr/µm), subtracting the local linear continuum, and integrating over the wavelength window [λ − 2w(λ), λ + 2w(λ)~]. For the measurement of the central wavelength, see Sect. 3.3.

Certain groups of lines are closely spaced, and some have problematic continuum determinations with the method above. For a total of around ten of these cases, we applied an individualized approach, fitting a functional model consisting of the sum of two or more Gaussians, on top of a first to third degree polynomial continuum. Three of these cases are shown in Fig. 4 to illustrate this method. The order of the polynomial for the continuum was larger than linear for the cases where significant local curvature due to AIB emission was present. The width of the fitting window and degree of the continuum were manually fine-tuned for each group of the overlapping lines. The fitted amplitude of each Gaussian component is then converted to an intensity by calculating the area. For these manually fit cases, the wavelengths were measured by fitting this model with variable central wavelengths for the Gaussian components, as opposed to the centroid method discussed in Sect. 3.3.

The amount of noise in the spectrum contributed by local variations in the spectral baseline depends on the noise features near each wavelength. The error bars on the data points of the spectrum do not suffice, as there can also be noise due to local artifacts (e.g., fringing, other oscillating patterns, or small spectral features). Therefore we use an empirical local uncertainty estimate for each line, performed as follows. First, we apply our linear continuum subtraction to around ten wavelength windows near the line, creating a set of samples of the local spectrum. The continuum subtracted samples are then integrated over the same interval width that was used for the line intensity, yielding the values ${\delta }_i={\displaystyle {\int }_i\left[F_{\lambda }(\lambda )-C_i(\lambda )\right]\text{d}\lambda }$, where i is one of the wavelength intervals, F_λ(λ) is the flux density per unit wavelength, and 𝒞i(λ) the local model for the continuum. With a perfect continuum subtraction, and a noiseless spectrum, each integral δ_i would be zero. Since there is noise, and since a linear continuum is not a perfect representation of the underlying continuum, every δ_i deviates from zero and is a quantity in the same units as the line intensity (erg s⁻¹ cm⁻² sr⁻¹). An example is shown in Fig. 5. We use the median of the absolute values |δ_i| as an estimate for the uncertainty on the line intensity integral (median absolute deviation). To visualize this uncertainty on the flux axis in Fig. 5, the line intensity and its uncertainty were converted to amplitude values, by dividing by the Gaussian normalization factor $\sqrt{2\pi }\sigma$, with $\sigma ={\left(2\sqrt{2\text{\hspace{0.17em}ln\hspace{0.17em}}2}\right)}^{-1}w(\lambda )$, the standard deviation derived from the MIRI MRS FWHM curve.

In addition to the local noise in the spectral baseline, we also note the calibration uncertainties between the segments which we estimated as roughly 3% (Sect. 2.2). We expect this systematic contribution to affect all lines within one and the same segment with the same factor. For applications where detailed uncertainty models are required, a covariance matrix can approximate this effect using the off-diagonal elements $C_{ij}=f_S^2{I}_i{I}_j$ for lines in the same segment. Here f_s is the fractional systematic uncertainty (~0.03), and I_i the intensity of line i. The diagonal elements are $C_{ii}={\sigma }_i^2+f_S^2{I}_i^2$, with σ_i the local uncertainty estimated as described above. For the parts of the analysis that use line ratios (Sect. 4), one can simply ignore the systematic calibration uncertainty for lines in the same segment, or add the fractional uncertainty to the numerator and denominator otherwise. The latter corresponds to a $\sqrt{2\times {0.03}^2}\approx 4\text{\%}$ relative systematic uncertainty on the ratio of two line intensities, which can be added in quadrature to the local uncertainty.

Fig. 6 Histogram of fractional wavelength deviation (Sect. 3.3). The H I and H2 lines originate from different gas layers, and the difference between the velocity averages is about 7.6 km s−1.

3.3 Observed wavelength offset

In addition to the intensity, we also measured the central wavelength of each line, so that the observed wavelength λ_obs can be compared to theoretical or experimental rest-frame wavelength λ_rest. To determine the center of each line, we calculate its centroid over a window spanning from λ_rest − w(λ) to λ_rest + w(λ), and the equation used is ${\lambda }_{\text{obs}}=\left[{{\displaystyle \int }}^{\text{}}\lambda F_{\lambda }(\lambda )\text{d}\lambda \right]/\left[{{\displaystyle \int }}^{\text{}}{F}_{\lambda }(\lambda )\text{d}\lambda \right]$. The measured wavelength shifts are given in a column of Table A.1, defined as Δλ_obs = λ_obs − λ_rest. The most precise atomic and molecular data are those for the Η i and H₂ lines, and a histogram of the relative wavelength shift Δλ_obs/λ_rest is shown in Fig. 6 for these lines. The average relative wavelength shift for the Η i lines is 11.3 km s⁻¹, and the standard deviation of this sample is 6.0km s⁻¹, while the H₂ lines are centered around l8.9 km s⁻¹ with a standard deviation of 5.4km s⁻¹. Dividing the standard deviations of the distributions by the square root of the number of lines yields the uncertainty on the measurement of the averages: 11.3 + 0.9 km s⁻¹ for H i and 18.9 + 1.5 km s⁻¹ for H₂. The difference between these two components is 7.6 + 1.7 km s⁻¹.

High-resolution spectroscopic observations of С and S recombination lines show that the emission lines originating from the Bar have an intrinsic width of ~2.5 km s⁻¹ and a red-shift corresponding to a radial velocity of v_LSR ≈ 10.6 km s⁻¹ (Goicoechea & Cuadrado 2021) with respect to the Local Standard of Rest (LSR). In the same study, the surrounding H ii region has a line width of ~20 km s⁻¹ and a blueshift of v_LSR ≈ −5 km s⁻¹, as measured using He recombination lines. The expected velocity difference between the H ii region and the H₂ component of the Bar is therefore about 15 km s⁻¹ ; the difference we observe between the two averages reported in the previous paragraph is 7.6 ± 1.7 km s⁻¹.

The spectroscopic pipeline³ adjusts the wavelength grid of the spectra, so that the final products are in the barycentric (heliocentric) reference frame. Hence the data still include the motion of the Sun, with respect to the LSR that is used to report v_LSR in the literature. Using the definition of the LSR from the “coordinates” module of Astropy, we find that a radial velocity of 0 km s⁻¹ in the heliocentric frame corresponds to v_LSR = −17.17 km s⁻¹ for the coordinates of the Orion Bar. Converting the velocities of our Η i and H₂ components, we find v_LSR(H i) = −5.9 km s⁻¹ and v_LSR(H₂) = 1.7 km s⁻¹. The velocity we observe for the ionized region is close to the literature result of −5 km s⁻¹ cited in the previous paragraph, while the results for the Bar differ by about 9 km s⁻¹. The velocity difference between the H ii and neutral layer is qualitatively similar to the previous observations, as the H ii region is also blueshifted with respect to the molecular component, but the size of the velocity difference we observe is about half the literature value. Therefore, the H₂ rotational lines likely probe a different part of the molecular region than the С and S recombination lines.

The wavelength offset with respect to previous observations reported above is smaller than the reported wavelength calibration uncertainty of 27km s⁻¹ at 28 µm, but similar to the 9km s⁻¹ accuracy that is expected near 5 µm (Argyriou et al. 2023). A shift of this size is smaller than one resolution element, considering that the spectral resolution of MIRI MRS varies between 3500 and 1500, from the shortest to the longest wavelengths respectively, corresponding to a radial velocity resolution of about 85–200 km s⁻¹ (Labiano et al. 2021). Our radial velocity measurements are more precise than the resolution, as the shape of the (instrumental) line profile is resolved by the detector pixels and the wavelength grid in the final products.

4 First analysis

4.1 Hydrogen recombination

We compare the observed Η i line ratios to recombination theory (Hummer & Storey 1987; Storey & Hummer 1995), using the recent Case В coefficients from Prozesky & Smits (2018). The assumption of Case В theory is that all Lyman transitions are optically thick, and absorbed in the same location where they were emitted by the gas (Baker & Menzel 1938). Practically, this is equivalent to ignoring transitions to the ground state, and case В approximates the emission well under typical nebular conditions (Osterbrock 1962). In Fig. 7, we show the observed intensity of lines with lower state n_l = 6 and n_l = 7 over the available range of upper states n_u, for the H ii template. To compare the observed intensity with the theoretical emissivity, we divide the respective values by that of the 5-6 line at 7.46 µm.

In other words the ratio I_l−u/I₅₋₆ is shown for the observations, while the theoretical data points are the emissivity coefficient ratios j_l−u/j₅₋₆. The lower panel of Fig. 7 shows the deviations from the theory, using the quantity (I_l−u/I₅₋₆) / (j_l−u/j₅₋₆), or in other words the observed line ratio divided by the theoretical line ratio.

Consistent with Peeters et al. (2024), we choose the emissivity coefficients from Prozesky & Smits (2018) with the electron temperature T_e = 10⁴ K and the electron density n_e = 10³ cm⁻³, and zero radiation field. These parameters work well to describe the observed emission ratios within the error bars. We also compared with different models from the same set, for example by changing n_e to 10² cm⁻³ or T_e to 5 × 10³ K, and the differences in the observed-to-theory ratio are minor compared to the error bars. There is a systematic offset however, where the data points for the H II template differ by a factor of 1.1 from Case B on average. For the other four template apertures, the observed-to-theory ratios have a similar scatter, but they are centered at 1.0 and do not show a significant systematic offset. The most likely explanation for the offset in the H II template is He I contamination, as these are strongest in the H II spectrum, and nearly all H I lines used in this analysis overlap with one or more He I lines according to our Cloudy-based line list. As discussed in Sect. 3.2 and Fig. 4, we were able to quantify this contamination for H I 5–6, where it is about 10%. For the other lines, the changes in the line profile are too subtle to separate the contamination. Since the 5–6 intensity was corrected, but the others were not, a similar contamination level of 10% would explain the offset by a factor of 1.1 we are seeing in Fig. 7.

There are no deviations that resemble the typical effects of dust attenuation, where a stronger attenuation is expected at shorter wavelength, and at 10 µm by the characteristic feature of silicate dust. Previous work has shown that the dust in the Orion Bar has a patchy attenuation profile, strongly varying in magnitude between different areas (Walmsley et al. 2000; Peeters et al. 2024). But the dust column affecting the ionized layer appears to have too little extinction in the MIR to be constrained by this set of lines.

Finally, we also compared the observations to a model developed for the planning and interpretation of the PDRs4All observations (see Berné et al. 2022), which is based on Cloudy, a code that models the gas as a function of the depth into the medium and computes the total emission (Ferland et al. 2017), instead of just providing emissivity coefficients. This is the same model from which we derived part of the theoretical line list, which we use in Sect. 3.1 to identify the lines. The model parameters were taken from Pellegrini et al. (2009) and Shaw et al. (2009), and can be briefly summarized as follows. The main assumptions for the radiation field are an illuminating star with a Kurucz type model with T_eff = 39 600 K for the illuminating star, for which the emission rate of ionizing photons is scaled to Q_LYC = 9.8 × 10⁴⁸ s⁻¹. The gas density assumes a constant pressure model, with an initial electronic density of n_0,e = 3160 cm⁻³.

The line intensity ratios predicted by the Cloudy model were also added to Fig. 7, and we find that this model overestimates the line ratios at the shortest wavelengths. Since the Cloudy model includes dust internal to the H II and the ionization front, but not the foreground dust, we test if a dust extinction correction would bring the model results more in line with the observations. The PDRs4All NIRSpec overview by Peeters et al. (2024) shows that magnitude of the foreground extinction varies between an A(V) of 0.9 and 1.9, when assuming the R(V)-parameterized average Milky Way curve by Gordon et al. (2023), evaluated at R(V) = 5.5. We use the same curve, and the extinction-corrected comparison is also shown in Fig. 7, for the maximal value of A(V) = 1.9. Even with this A(V), the deviation from the observations remains. The Cloudy model assumes a spherical shell geometry, and the way the observed emission is recorded in the model might not match the way we observe the ionized layer of the Orion Bar. The dust distribution and its presence in the HII region are also assumptions that could cause deviations. The Case B theory makes no geometrical or dust assumptions at all, and only provides emissivity coefficients which we use directly. A good match between Case B theory and our observations, could mean that the main contribution to the H I emission comes from an ionized gas layer subject to very little dust or geometry effects.

Fig. 7 Analysis of H I recombination line ratios. Left panel: transitions to lower levels 6 and 7 for H II template spectrum, comparing H I emission line ratios between observations and theory. Labels next to each data point indicate the transitions. The Cloudy model with foreground dust correction applied (Cloudy + FG dust) is shown assuming the largest value for the Orion Bar foreground (A(V) = 1.9), making it clear that the effect is small. Right panel: zoom-in to better show the deviations between the observations, the Case B model, and the Cloudy model. Bottom panels: observed line ratio divided by theoretical line ratio, for the three models as coded by the legend. Error bars: Uncertainty due to local noise (black) and with the 3% systematic calibration uncertainty added (gray).

4.2 Ionic line ratios

The emission line ratios of two ions of the same element can be used to probe the hardness of the radiation field. The spectra contain pairs of ionization stages for several elements. Using the brightest lines for each element, we show intensity ratios in Fig. 8 for the ions of Ar ([Ar III] 8.99 µm / [Ar II] 6.99 µm), Ne ([Ne III] 15.6 µm / [Ne II] 12.8 µm), S ([S IV] 10.5 µm /[S III] 18.7 µm), Fe ([Fe III] 22.9 µm / [Fe II] 26 µm) and Ni ([Ni III] 7.35 µm / [Ni II] 6.64 µm). To support the discussion below, we provide an overview of the observed species and their ionization potential (from the ASD) in Table 1. The table shows the photon energy threshold to ionize each species to the next stage by photoionization. For example, the energy required to produce S IV is the value listed for S III.

As shown in Fig. 8, there is a large difference of the line ratio profiles comparing Ar, Ne, and S versus Fe and Ni, as the ions originate from different layers along the line of sight. The Ar, Ne and S ions reside mainly in the ionized layer, since the photon energies required to produce Ar II (15.8 eV), Ne II (21.6 eV), and S III (23.3 eV) are greater than the H II threshold of 13.6 eV. The ratio of the Ar iii/Ar ii pair is significantly lower in the Atomic template, and this behavior is consistent for both the available Ar III lines. This is because the intensity of the [Ar ii] 6.99 µm line in the Atomic template is actually stronger than that in the H II template, as listed in Table A.1. Because of the ionization thresholds to produce both Ar ii (15.8 eV) and Ar iii (27.6 eV) are significantly lower than those of the Ne and S ion pairs, the Ar ions remain present in regions with softer radiation fields, closer to the atomic layer. The Ar ratios hence probe a wider range of physical conditions and are more complex to interpret.

The Fe and Ni ions have much lower ionization potentials, and the line ratios therefore have a very different behavior. Since the photon energies to produce Fe II (7.9 eV) and Ni II (7.6 eV) are below 13.6 eV, the denominator in the line ratios of Fig. 8 will have strong contributions by the atomic H layer, as these species remain ionized past the H ionization front. The measured intensities for Fe III and Ni III are at their highest in the H II template, while the intensities of Fe II and Ni II reach their highest value in the Atomic template. The energies to produce Fe III (16.3 eV) and Ni III (18.2 eV) are similar, and therefore the line ratios are expected to have a nearly identical behavior. This is demonstrated in Fig. 8, which further confirms that Ni III was identified correctly in Sect. 3.1.3.

Finally, we can use the ratios in Fig. 8 and the ionization potentials from Table 1 to discuss the hardness of the radiation field. All fine-structure lines (except for S I) have a maximal intensity in either the H II or Atomic template, and then get systematically weaker for DF1, DF2, and DF3 in that order, indicating that the column densities of all ions are lower for regions further away from the star. However, the relative abundances of the upper ionization stages (Ar III, Ne III, and S IV) compared to the stages below, are observed to be larger for the DF templates. Both the Atomic and DF template regions contain a part of the ionized layer bordering the PDR, but it appears that the radiation field affecting this ionized gas is harder at the locations of the DF templates. The ratios for Ne and S exhibit trends that are similar to each other, as expected from their similar ionization potentials in Table 1 and existing observations of H II regions (e.g., Giveon et al. 2002; Martín-Hernández et al. 2002). We also added a line ratio probing Ne iii/S IV to Fig. 8, as these species are sensitive to radiation of around 41.0 eV and 34.9 eV respectively, and hence this ratio probes a narrower energy range. This ratio has a slight downward trend, indicating that the radiation field hardness in the 35–41 eV range decreases mildly, while the overall hardness probed by the other line ratios over a wider energy range of roughly 20–40 eV increases strongly.

Fig. 8 Line intensity ratios for several ion pairs as a function of template position. The wavelengths are given in the legend in µm. Each curve is normalized to its value at the H II position. These ratios can be interpreted as indicators for radiation field hardness, and are sensitive to different photon energies (thresholds in Table 1), or different layers of the gas (Fe and Ni).

4.3 Molecular hydrogen

Molecular hydrogen plays an important role in the evolution of interstellar matter and star formation. Present in a variety of interstellar conditions (gas density and UV flux), H₂ is the most abundant molecule that efficiently cools the gas and may also heat it. Due to its high UV opacity, H₂ also shields itself and other molecules from photodissociation, and collisions with excited H₂ trigger endothermic chemistry (see e.g., Le Petit et al. 2006).

4.3.1 Excitation mechanisms and heating and cooling processes

Molecular hydrogen contributes to the gas thermal budget in the course of its excitation-deexcitation cycle that we summarize here. The H₂ excitation occurs mostly through radiation and collisions, with a minor contribution by the H₂ formation process with a rate that is ten times lower than for UV radiation, and so we do not consider the latter process here. Radiative excitation or pumping of H₂ takes place through strong absorption bands in the FUV (λ ≤ 100 nm). In most cases (about 9 of 10), following this UV excitation, the molecule deexcites to the ground electronic state where it cascades down through the rovibrational states. This radiative pumping thus provides both vibrational and rotational excitation of H₂ (Black & Dalgarno 1976; Black & van Dishoeck 1987).

On the other hand, rotational excitation of H₂ in its vibrational ground state (v = 0) occurs through inelastic collisions with gas species (mostly H, He and electrons). The molecule may then radiate away the energy taken up in the collision through the MIR lines, which contributes to gas cooling (see section 4.3.3). Collisional excitation of vibrational levels (v ≥ 1) is rare because it requires a high gas temperature (the typical energy between vibrational levels is 0.5 eV or ~6000 K). At high densities (≳10⁴ cm⁻³), collisional deexcitation of H₂ from radiatively pumped levels will lead to gas heating.

4.3.2 Excitation diagrams in the Orion Bar

The Orion Bar is known to have a large amount of warm H₂ gas, which can be traced by pure rotational H₂ lines (Parmar et al. 1991; Habart et al. 2004; Allers et al. 2005; Kaplan et al. 2021). As a first result for these detailed H₂ data in the MIR, we present estimations of the H₂ excitation temperature based on excitation diagrams of the pure-rotational lines measured in the five template spectra. The excitation diagram for the DF3 template is presented in Fig. 9, and we analyzed analogous excitation diagrams for the other four spectra. In this type of diagram, the column density of each upper level u is derived from the observed line intensities and known atomic parameters for the transitions. The quantity on the vertical axis is N_u/g_u = 4πI_ul/(g_uA_ul[E_u – E_l]), and plotted as a function of E_u/k. In these equations, I_ul is the measured integrated intensity of the line (erg s⁻¹ cm⁻² sr⁻¹), g_u is the statistical weight of the upper level, E_u and E_l are the energies of the upper and lower levels, A_ul is the Einstein coefficient of the transition, and k is the Boltzmann constant to express E_u as a temperature. Here we assumed the emission lines to be optically thin, given their low Einstein coefficients (Table A.1), so that an optical depth term is not used in the calculation of N_u/g_u (e.g., Goldsmith & Langer 1999).

To describe the H₂ excitation in Fig. 9, we fit two linear components to the data for which we assumed a 3% fractional uncertainty to account for the calibration error explained in Sect. 3.2. The resulting parameters are presented in Table 2, as derived from the fitted ordinate intercept and slope, based on the following linear prescription: ln(N_u/g_u) = b + aE_u = ln(N⁰) – E_u/T*. With i the index of the component (1 or 2), the fitted representative excitation temperature is $T_i^{\ast }=-a_i^{-1}$, and the parameter $N_i^0=e^{b_i}$ is a factor that scales with the column density. The total column density for each component is derived from the fit result as $N_i=\underset{u}{{{\displaystyle \sum }}^{\text{}}}{N}_{u,i}=e^{b_i}\underset{u}{{{\displaystyle \sum }}^{\text{}}}{g}_u\exp \left(E_u{a}_i\right)$, where the summation is $Z\left(T=-a_i^{-1}\right)$, the rotational partition function of H₂ evaluated at the fitted temperature. The equation used for Table 2 is then $N_i=e^{b_i}Z\left(-a_i^{-1}\right)$, where Z can be approximated as Z(T) = 0.02477T[1 – exp(1 – 6000/T)]⁻¹ (Herbst et al. 1996). While the results presented in Table 2 cover all five template spectra, the H II and Atomic regions have one or two missing lines, as can be seen from the masked out values in Table A.1. This results in larger uncertainties on the reported parameters for these two regions.

A more thorough H₂ analysis by Sidhu et al. (in prep.) will investigate the spatial variations in detail, by performing similar two-component fits on a spaxel-per-spaxel basis. This analysis will also include the H₂ lines observed with the NIRSpec data, to provide an in-depth characterization of the rovibrational excitation. For all five spectra, we assumed the same two-component shape for the excitation diagram, with an ortho-to-para ratio (OPR) of 3. The OPR measured by Kaplan et al. (2021) is 2.99 based on a large number of 1−0 rovibrational lines from 1.45 to 2.45 µm, for a similar warm molecular region of the Orion Bar PDR, somewhat further to the South-West with respect to our DF regions. The combined NIRSpec and MIRI spectroscopy allow for the mapping of local changes in the OPR, and preliminary results from the spatially dependent analysis by Sidhu et al. (in prep.) show that the OPR is very close to 3 in the areas of interest, and mostly constant over the field of view.

Although the diagram in Fig. 9 is curved, and suggests a continuous distribution of excitation temperature, we characterize the H₂ excitation with only two components. The $T_1^{\ast }$ component is fitted to the data of low excitation, which is rotational lines within v = 0 from S (1) to S (7), while the $T_2^{\ast }$ component is similarly fitted to rotational lines within the excited state v = 1, from S (5) to S (9). We also include the v = 0 S (8) line in the $T_2^{\ast }$ fit, as it appears to follow the linear trend, and therefore belongs to the high excitation component, indicating that radiative pumping is also important in the v = 0 level. In the following, we assume that $T_1^{\ast }$ is a good proxy to the gas temperature. While some curvature affects all transitions, causing mismatches to the linear model, the latter can still be used to quantify the average slope over the levels up to S (7). To model curved excitation diagrams in external galaxies, prescriptions based on an underlying power-law distribution for the temperature have been used (e.g., Togi & Smith 2016), which is appropriate when a collection of PDRs and shocks is considered. In our case, the temperature distribution is determined by the temperature gradient across the H₂ emission zone, and there is no direct argument that a power law distribution can approximate the resulting mix of temperatures. It is clear that these data probe shape of the excitation diagram very well, and could therefore be sensitive to the continuous range of gas temperatures. Developing a method to infer a suitable temperature distribution and reproduce the curvature of the diagram, would allow for more accurate gas mass estimations, but is outside the scope of this work. In the next section, we further our analysis by comparing the H₂ line intensity data to a PDR model.

The derived $T_1^{\ast }$ or gas temperature is in the range from 660 to 720 K, for all five regions. The gas temperature of the ionized and atomic layers is expected to be much higher, while the H₂ rotational temperature in the H II and Atomic apertures is similar to that of the DF apertures. Therefore, the H₂ emission may originate from a different part of the PDR, behind the H II and Atomic regions. Since the lines used for this component are located in different MIRI bands, we do not need to account for correlations induced by the systematic calibration uncertainty. The derived H₂ column density is between the orders 10²⁰ and 10²¹ cm⁻², with a significant increase for DF2 and DF3. Comparable gas temperatures are found for the three dissociation fronts, which are at different distances from the ionization front (0.02 pc for DF1, 0.03 pc for DF2, and 0.04 pc for DF3). This is compatible with our knowledge about the geometry of this region, where DF1, DF2, and DF3 are the edges of a terrace-like structure (Habart et al. 2024; Peeters et al. 2024). The observed temperature is somewhat higher compared to the estimations of 430 to 630 K given by Allers et al. (2005, Table 4), which used a different technique based on the S (1)/S (2) and S (2)/S (4) line ratios. We note that restricting our fit to the 0–0 S (1), S (2) and S (4) lines gives gas temperatures around 500 K, in better agreement with the Allers et al. (2005) results, in particular for their apertures C, D and F which are associated with dissociation fronts at a similar distance from the PDR edge, as are DF2 and DF3.

As the energy of the upper level increases, we expect a transition from collisional to radiatively pumped H₂ excitation. The $T_2^{\ast }$ component of our fits models the v = 1 lines from S (5) to S (9), which have much higher upper energy levels, making it unlikely that the excitation is dominated by thermal collisions for the physical conditions of the Orion Bar. The derived excitation temperature is around 2700 K, with a column density that is significantly lower, of the order 10¹⁷ to 10¹⁸ cm^–2. Again, DF2 and DF3 exhibit the highest column densities of excited H₂. The high rotational temperatures are compared to a model in the next section, where this behavior is discussed further. For this fit, there are some positive correlations between the line intensities due to the calibration uncertainty, since the 1–1 S (6) and S (7) lines are both measured in Ch 1 MEDIUM, while 0–0 S (8), 1–1 S (8), and 1–1 S (9) are in Ch 1 SHORT. From the covariance matrix suggested in Sect. 3.2, only the correlations between the three lines in Ch 1 SHORT are significant, with a correlation coefficient of around 0.6 at most. We therefore expect the correlations to have a minor effect on the fit results, given that the uncertainties for the affected lines are small, while their correlation coefficients are not very close to 1.

We note that for the results discussed above, we did not apply a dust attenuation correction to the line intensities from which the excitation diagrams were derived. Nevertheless, it can be seen in Fig. 9 that the S (3) is affected by the silicate dust absorption feature at 10 µm. As will be discussed in more detail in Sect. 4.3.4, we applied a simple attenuation correction to estimate the effect on the results. Briefly summarized, an attenuation correction results in equivalent rotational temperatures, while for the derived column densities there could be a significant increase of around a factor ~3 (right panel of Table 2). However, the latter factor depends heavily on the relative strength of the silicate feature in the assumed attenuation law.

Fig. 9 H2 excitation diagram for DF3, with linear fits and comparison to Meudon model. Black crosses: column density of the upper level populations derived from the data, versus the energy of the upper level with respect to the ground state. Error bars are shown where they exceed 5% (includes 3% systematic uncertainty). Red circles: column densities derived from our Meudon PDR model. Blue line: fit including lines up to 0–0 S (7), excluding the 0–0 S (3) line (gray) as it is subject to extinction by the 10 µm silicate feature. Orange line: fit including 0–0 S (8) and 1–1 S (5) to 1–1 S (9). Temperatures derived from the slope of each line are shown in the legend.

Fig. 10 Top panel: structure of the HI/H2 transition predicted by Meudon model as function of depth. Densities of H and H2, with vertical line that indicates the transition defined as n(H) = n(H2). The local emissivity of each gas layer is shown for four representative H2 lines. Bottom panel: temperature and cumulative emissivity versus model depth. Black line and left y-axis: gas temperature for each layer in the model. Colorful lines and right y-axis: cumulative of the emissivity contributions by each layer (weighted by length scale). By inspecting the A(V) ranges over which the cumulative emissivity curves rise the most, we can see the typical range of gas temperatures from which the lines originate. For the 1–1 S (7) line for example, the observed excitation temperature is ~2700 K, while around 75% of the emission originates from gas cooler than 1000 K, due to UV pumping.

4.3.3 Comparison to PDR model

To interpret the observed H₂ emission, we compare our line flux data to the results of an existing Orion Bar model, computed with version 1.5.4 of the Meudon PDR Code (Le Petit et al. 2006), using parameters from the literature (Joblin et al. 2018). This model is isobaric, with a thermal pressure P_th = 2.8 × 10⁸K cm⁻³, an ISRF with scaling factor G₀ = 3.1 × 10⁴, and a geometrical factor of ~2 that corresponds to an inclination angle of 60°. To model processes driven by dust such as H₂ formation and photoelectric heating, a dust size distribution of the MRN type (Mathis et al. 1977) is used with a power law exponent of -3.5 and a grain size range of 3 to 300 nm, and optical properties corresponding to a mixture of graphite (30%) and silicate (70%) grains. For the attenuation of the radiation field, we used an extinction curve according to the analytical model by Fitzpatrick & Massa (1990), assuming R(V) = 5.62 and N(H)/E(B – V) = 1.05 × 10⁻²² cm⁻². While a newer version of the Meudon code was used compared to Joblin et al. (2018), most of the above parameters were adopted from the literature without further adjustments, to compare the existing knowledge to the new data. Updating the parameters by means of fits and model grids will result in a better match to the observations, but this task is outside the scope of this work, and is part of the work by Meshaka et al. (in prep.).

Using the depth-dependent output of the Meudon model, we examine the gas temperature and the emissivity of several lines in each gas layer of the model, as shown in Fig. 10. For each H₂ line, the most significant emissivity contributions originate from a continuous range of depths and temperatures, as illustrated by the cumulative emissivity curves. The highly excited H₂ lines of v = 0 (S (5) and above) peak in the first hot zones of the H/H₂ transition, and in the model, the 0–0 S (7) emission originates from the gas layer where the temperature curve has a shallow gradient as a function of depth, decreasing from around 1100 to 800 K. A similar plateau in the emissivity curve is observed for S (5), and S (6) (not shown). While the temperature decreases, the 0–0 S (7) emissivity stays mostly constant (yellow curve in top panel of Fig. 10), suggesting an important contribution of nonthermal, radiative excitation. On the other hand, the gas temperature in the model roughly matches the slope in the excitation diagram, so a mix of thermal collisions and radiative pumping is likely for S (5) to S (7). This occurs at a depth where n(H) is still larger than n(H₂), and there is little emissivity in these lines past the H I/H₂ transition defined as n(H) = n(H₂) (vertical line in Fig. 10).

To follow the behavior of the less excited lines (0–0 S (4) and below) we show the 0–0 S (1) emissivity profile: we see that it peaks around the H I/H₂ transition and corresponds to a gas temperature around 450 K. The curve for S(2) peaks at a similar location and temperature, while S(3) and S(4) (not shown) have a similar profile with a peak at the H I/H₂ transition, where the gas temperature is around 600 K. This is in good agreement with the $T_1^{\ast }$ temperatures we found in Sect. 4.3.2 and supports a thermal, collisional origin for these lines. Using the cumulative emissivity curve (red curve in bottom panel of Fig. 10), we determine the depths where the 10%, 50%, and 90% quantiles of the emission are reached. We find that there are significant contributions to S (1) from gas ranging from around 320 to 960 K, while the 50% quantile corresponds to around 500 K. This explains why the curvature of the diagram seems to smoothly continue from the high-excitation lines (driven by a mix of collisional and radiative excitation) to the low-excitation lines, even when the latter are mostly thermally driven.

Finally the emissivity profile of the highly excited 1–1 S (7) line has two maxima at A(V) ≃ 0.6 and 1.2. The first peak at A(V) ≃ 0.6 is probably due to the competing effects of the rising H₂ abundance (Fig. 10 upper panel) and the decreasing radiative pumping. The second peak points at a collisional contribution at the H I/H₂ front, which then decreases along with the temperature, in a steeper fashion than 0–0 S (1) because the collisional excitation of 1–1 S (7) requires higher temperatures. Analogous to the paragraph above, the 10% and 90% emission quantiles reveal a large temperature range from which 1–1 S (7) originates, with significant contributions by gas roughly between 500 and 2600 K. The above double-peaked emissivity profile and temperature range apply to all four v = 1 lines studied in this work. While there is a minor contribution from the 2600 K gas, which is close to the measured rotational temperature of around 2700 K, the bulk of the emission (~70%) originates from gas cooler than 1000 K. Therefore, collisional excitation cannot be the main population process for the upper levels in the v = 1 vibrational state, and pumping by UV photons is expected to be the dominant excitation mechanism.

With the temperature structure of the model described above, we now compare the observed (DF3) and theoretical excitation diagrams (black crosses and red circles in Fig. 9); we note that the model intensities were not rescaled to the observations. For the cold component of the gas, the S (2) and S (4) lines match well, while there is a minor mismatch with S (1) that indicates a slightly different curvature. There is a large but expected mismatch for S (3), due to the 10 µm silicate extinction (see Sect. 4.3.4). The S (5) to S (8) lines on the other hand are strongly overestimated: the level populations based on the line emission from the model are up to a factor of 3 higher than the values derived from the observations. This deviation, in the transition region between the two linear components of the excitation diagram, is observed for all five template spectra. Continuing to the v = 1 pure-rotational transitions, the model does match the data. This means that the PDR model correctly matches the most strongly UV-excited layer, but needs parameter adjustments to model the transition into the layers where thermal H₂ excitation becomes dominant. Despite the fact that we use literature values for the parameters instead of values optimized through fitting, this fixed model still qualitatively reproduces the transition in the slope near 0-0 S(8).

Future adjustments to the model parameters will need to slightly increase the amount of cool gas at ~500 K, and more strongly reduce the amount of hot gas at ~900 K, steepening the S (1) to S (8) part of the curve while keeping the behavior of the v = 1 levels unchanged. Changes to the model, for the Orion Bar specifically, are under development by Meshaka et al. (in prep.). They replace the isotropic ISRF with scale factor G₀ of the model presented here, by a synthetic stellar spectrum representative of Theta 1 Orionis C. This paper will explore alternate recipes for specific processes in the PDR model, such as the excitation of H₂ that is newly formed on dust grains. In combination with these changes, they will optimize the thermal pressure and distance with respect to the star for DF1, DF2, and DF3 individually. Through these efforts, our data have the potential to reveal how well the state-of-the-art PDR codes model the excitation, heating and cooling mechanisms related to H₂. Alongside measurements of the AIB emission in the same spectra, these data could have important implications for understanding the thermal balance in PDRs, because the small dust grains that carry the AIB also contribute significantly to the UV extinction and gas heating via the photoelectric effect (Bakes & Tielens 1994).

4.3.4 Dust attenuation of S(3)

We applied no dust attenuation correction for the above analysis, as we expected the extinction in the mid-IR to be small and similar for the different lines, thereby having only a minor effect on the derived results for H₂. Despite this expectation, the 0–0 S (3) line had to be excluded from the above analysis, as the attenuation due to the 10 µm silicate dust feature appears to be rather strong. We measured the decrement of the S (3) data with respect to an interpolated value using a local linear fit, and find values that correspond to 0.8 mag for HII, 1.1 mag for Atomic, 1.2 mag for DF1, 0.7 mag for DF2, and 0.9 mag for DF3 (or factors of 0.49, 0.35, 0.35, 0.51, and 0.44 respectively). We denote these decrements as a difference between extinction values: A(S(3)) – A(other).

As an estimate for the typical strength of the 10 µm silicate feature, we assume that the attenuation curve is given by the R(V)-dependent average Milky Way extinction curve of Gordon et al. (2023) evaluated for R(V) = 5.5, in a screen geometry. We note that we deliberately distinguish between the words “attenuation” and “extinction” in this discussion. As the scattering albedo at the wavelength of the S(3) line is small, the difference between extinction and attenuation will depend on the geometry of the dust relative to the emitting gas. For the assumed curve, A(λ)/A(V) is 0.095 at the wavelength of the S (3) line, while this ratio for the other MIR H₂ lines is about 0.045. The observed S (3) decrement can hence be written as A(S (3)) – A(other) ≈ 0.095A(V) – 0.045A(V) = 0.05A(V). Using this to extrapolate the observed decrement of S (3) to the V-band extinction, yields A(V) values of roughly 15 (HII), 23 (Atomic), 23 (DF1), 14 (DF2), and 18 (DF3), for the extinction that affects the emission originating from the H₂.

A recent version of the “H₂ tool” available in the PDR toolbox (Pound & Wolfire 2023), supports simultaneous fitting of A(V) and an H₂ excitation diagram. Under the same assumptions for the attenuation model, the simultaneous fit yields equivalent results for A(V) as the more ad-hoc estimates presented above. To clarify the effect of attenuation on the H₂ analysis of Sect. 4.3.2, we applied the attenuation correction to the line intensities using the above assumptions, set up analogous excitation diagrams, and applied the same fitting method. The results presented in the right panel of Table 2, show that the effect on the derived temperature is minimal, while the derived column densities log N_i(H₂) increase by ~0.3 to 0.5, or a factor of ~2 to 3.

The results of this attenuation correction depend strongly on two assumptions. Firstly, the correction for the H₂ column density presented in Table 2 is completely dependent on the assumed A(S (3))/A(other) ratio. For larger values of this ratio, the same S (3) decrement can be reproduced with a smaller value of A(other), resulting in a smaller correction for the H₂ column. Secondly, the extrapolated A(V) values are highly dependent on both the silicate feature strength, and the typical V-band to MIR extinction ratio. More precisely when the ratio of A(S (3)) – A(other)]/A(V) is higher, then the derived A(V) will be smaller. Therefore, deviations from other published values are to be expected, and the A(V) values presented in this section should only be interpreted as a scaling factor for the S(3) decrement. Therefore, corrected H₂ column densities presented in Table 2 only serve as an example of the typical impact of the attenuation correction.

To set up a more realistic attenuation curve, several aspects of the observed system need to be considered. According to the diagram of the PDR structure by Peeters et al. (2024, their Fig. 14), the foreground extinction by the ionized layer contributes A(V) ~ 0.9–1.9 mag and the atomic layer contributes A(V) ~ 10–12 mag, in front of DF1. In the analysis of the previous paragraph, we used a simple extrapolation to estimate the V-band extinction, based on the measured attenuation of S(3) relative to the other H₂ lines and an extinction curve, which assumes a screen geometry for the dust. However, part of the dust is spatially mixed with the emitting H₂ gas, which affects estimates of the “internal” extinction in the PDR that corresponds to the integrated dust column of the atomic layer and the layer from which the H₂ emission originates. The observed attenuation might thus differ significantly from the total extinction by the dust column along the line of sight, even in the absence of scattered light. To estimate the dust column, and connect the attenuation at different wavelengths more accurately, a suitable attenuation curve model needs to be set up, which is done by making assumptions about the spatial distribution of the emitting gas and the dust. Two formalisms for the geometry were explored by Peeters et al. (2024), to determine the internal PDR extinction based on H₂ lines with the same upper level (1−0 S (1) at 2.1 µm and 1−0 O(5) at 3.2 µm). In their “foreground only” case, all the dust is assumed to be in front of the H₂ emission (i.e. in the atomic region only), and the estimate for DF1 is A(V) ~ 9. In their “intermingled” case, the dust is perfectly mixed with the H₂ emission and the estimated dust column increases to A(V) ~ 37. This illustrates the large differences for the derived dust column, which result from geometrical assumptions. Meshaka et al. (in prep.) will present a more thorough analysis of the dust attenuation, making use of both the near-IR and mid-IR H₂ transitions that have the same upper level.

5 Conclusions

The MIRI MRS data of PDRs4All allow spatially resolved studies of the MIR lines and the diagnostics derived from those lines, for the conditions of the ISM in the Orion Bar PDR. With the updated JWST pipeline and reference files, the calibration has been greatly improved since the completion of the PDRs4All observing program. From the improved data, we extracted five template spectra that represent the H II region, the atomic region near the ionization front, and the three dissociation fronts in the footprint of the IFU mosaic. We provide an overview of the lines present in these regions, with measurements of the intensity, uncertainty, and central wavelength. Thanks to the combination of spectral resolution and depth, we reveal a rich inventory of H I recombination lines, H₂ pure-rotational lines, and metal fine-structure lines. This inventory will enable more in-depth studies that will make full use of the spatial resolution of these IFU data, such as maps of the H₂ excitation by Sidhu et al. (in prep.), or the sulfur abundance by Fuente et al. (2024). In addition, new constraints based on these lines will be applied to further improve theoretical PDR models, such as those based on the Meudon code (e.g., Meshaka et al., in preparation).

An initial analysis is presented, summarizing the behavior of the gas from which the lines originate. The H I recombination lines originating from the ionized layer bordering the PDR are found to have ratios that closely match the emissivity coefficients from Case B recombination theory. At MIR wavelengths, the expected dust extinction has minimal effects on these lines. However, in the H II region, there is a contamination of up to 10% due to He I recombination lines. For the strongest lines, the spectral resolution and S/N of these data are sufficient to separate the contaminating component.

Further diagnostics of the ionized layer are provided by the fine-structure lines of several ionization stages of Ne, P, S, Cl, Ar, Fe, and Ni. We provide a brief overview of the ionization conditions by showing how the ion ratios vary for Ar III/Ar II, Ne III/Ne II, S IV/S III, Fe III/Fe II, and Ni III/Ni II. The metals in the H II layer appear more strongly ionized for the regions located in the foreground of the dissociation fronts, meaning that the radiation field is harder at these locations. We note that the Fe III/Fe II and Ni III/Ni II ratios behave very similarly to each other, but very differently from the other line ratios, as the Fe II and Ni II lines contain significant contributions from the neutral PDR gas.

We used the rotational H₂ lines to set up excitation diagrams and find similar excitation conditions for all three dissociation fronts. The smooth curvature in the trend of the excitation diagrams shows the departure from the linear (single temperature) approximation, and probes a mix of two effects: first, the temperature distribution along the line of sight, and second, the gradual transition from thermal collisions to radiative pumping as the dominant excitation mechanism. We compared the diagrams to a theoretical PDR model for which the parameters were set to values from the literature. This model qualitatively reproduces the data, and exhibits a curved excitation diagram where the transition in the slope of the diagram occurs at the same excitation energy as that derived from observations The main quantitative discrepancy is that the 0–0 S (5) to S (8) intensities are strongly overestimated, calling for parameter adjustments that reduce the amount of warm H₂ at ~900 K. On the other hand, the thermally driven 0–0 S (1) to S (4) lines of low excitation and the 1–1 S (5) to S (9) lines of vibrationally excited H₂ driven by UV pumping are well approximated by the model without parameter adjustments.

Acknowledgements

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 #1288, and the specific observations analyzed can be accessed via https://doi.org/10.17909/pg4c-1737. D.V.D.P. acknowledges support for program #1288 provided by NASA through a grant from the Space Telescope Science Institute. C.B. is grateful for an appointment at NASA Ames Research Center through the San José State University Research Foundation (80NSSC22M0107). JRG thanks the Spanish MCINN for funding support under grant PID2019-106110GB-100. T.O. is supported by JSPS Bilateral Program, Grant Number 120219939. MGW was supported in part by JWST Theory grant JWST-AR-01557.001-A. A.P. would like to acknowledge financial support from Department of Science and Technology - SERB via Core Research Grant (DST-CRG) grant (SERB-CRG/2021/000907), Institutes of Eminence (IoE) incentive grant, BHU (incentive/2021-22/32439), Banaras Hindu University, Varanasi and thanks the Inter-University Centre for Astronomy and Astrophysics, Pune for associateship. H.Z. acknowledges support from the Swedish Research Council (contract No 2020-03437) J.H. is sponsored (in part) by the Chinese Academy of Sciences (CAS), through a grant to the CAS South America Center for Astronomy (CASSACA) in Santiago, Chile. M.B. acknowledges DST INSPIRE Faculty fellowship and thanks The Inter-University Centre for Astronomy and Astrophysics for visiting associateship. A.F. is grateful to the European Research Council (ERC) for funding under the Advanced Grant project SUL4LIFE, grant agreement No101096293 and to the Spanish MICINN for funding support from PID2019-106235GB-I00. Work by Y.O. and M.R. is carried out within the Collaborative Research Centre 956, subproject C1, funded by the Deutsche Forschungsgemeinschaft (DFG) - project ID 184018867.

Appendix A Additional table

All identifications and measurements for the lines are listed in Table A.1.

References

All Tables

All Figures

Fig. 1 Examples of flux differences in overlapping regions between MIRI bands. Blue and orange curves show individual segments, with the bands indicated at the top of each panel. These were extracted from the band cubes produced by the pipeline, in this example for the H II template aperture. The black line shows the stitched spectrum, which is offset with respect to the blue and orange lines due to the additive offsets applied to the 12 segments, which are cumulative with Ch 2 LONG used as a reference.

In the text

Fig. 2 Apertures from which the five template spectra were extracted, overlayed on two slices of the MIRI MRS data cube. The horizontal axis is aligned with the y-coordinate of the IFU and the 9×1 mosaic strip. The blue component is the total flux at 11.3 µm (104 MJy sr−1), which contains some continuum, but mainly traces the 11.3 µm AIB and peaks in the atomic region. The red component highlights the locations of the dissociation fronts using the intensity of the H2 S (1) line at 17.0 µm(10−3 erg s−1 cm−2 sr−1). The two-dimensional color legend shows where one component (red and blue) or both components (dark purple) are bright.

In the text

Fig. 3 Overview of template spectra and lines, continued in Fig. 3. Before plotting, the segments were normalized individually, each divided by a constant factor, which is the average flux of that spectrum as integrated over the wavelength range of the plot. Then, multiplicative offsets using a factor of 1.5 were applied for clarity. Vertical lines and legend: Identified lines and species. Dashed black lines: Unidentified lines (UID). Sparse dotted black lines: End of a MIRI MRS band, where changes in the spectral resolution occur. We note a peculiar detail, namely broad absorption or emission or artifact features around 7.0 µm near the Aril line, which remain unidentified. The oscillations in the top panel are a good example of remaining fringes even after the fringe flat and residual fringe corrections in the pipeline.

In the text

Fig. 4 Examples of intensity determination for overlapping or nearby lines, in H II template spectrum. The x-axis has been shifted so that the zero point matches the wavelength in µm displayed in the top left corner of each panel. Blue curve: Fitted model consisting of a linear continuum and two or three Gaussian profiles. Dotted lines: central wavelengths of the Gaussian profiles. Upper left panel: [Ni II] 6.636 µm and H I 6.638 µm. Lower left panel: complex of three H I lines, with fitted wavelengths at 7.4955, 7.5030, and 7.5086 µm. Right panel: H I 5–6 and contributions by He I multiplets in its left wing.

In the text

Fig. 5 Demonstration of uncertainty measurement described in Sect. 3.2. The line in this example is [S I] 25.25 µm, and the oscillations are an artifact of MIRI MRS Ch 4. Blue segments with dots: local linear continuum estimates, with the center indicated. Orange crosses: Blue dot value plus the average of the continuum subtracted spectrum as integrated over the continuum window. The deviation between the orange cross and the blue dot samples the deviation that occurs due to local noise patterns and the linear continuum approximation. The median absolute deviation is used as the uncertainty on the line intensity integral. Blue circle and error bar: Continuum and median absolute deviation of the above samples. Orange error bar: Line intensity and uncertainty, visualized approximately as an amplitude.

In the text

	Fig. 6 Histogram of fractional wavelength deviation (Sect. 3.3). The H I and H2 lines originate from different gas layers, and the difference between the velocity averages is about 7.6 km s−1.
In the text

Fig. 7 Analysis of H I recombination line ratios. Left panel: transitions to lower levels 6 and 7 for H II template spectrum, comparing H I emission line ratios between observations and theory. Labels next to each data point indicate the transitions. The Cloudy model with foreground dust correction applied (Cloudy + FG dust) is shown assuming the largest value for the Orion Bar foreground (A(V) = 1.9), making it clear that the effect is small. Right panel: zoom-in to better show the deviations between the observations, the Case B model, and the Cloudy model. Bottom panels: observed line ratio divided by theoretical line ratio, for the three models as coded by the legend. Error bars: Uncertainty due to local noise (black) and with the 3% systematic calibration uncertainty added (gray).

In the text

	Fig. 8 Line intensity ratios for several ion pairs as a function of template position. The wavelengths are given in the legend in µm. Each curve is normalized to its value at the H II position. These ratios can be interpreted as indicators for radiation field hardness, and are sensitive to different photon energies (thresholds in Table 1), or different layers of the gas (Fe and Ni).
In the text

Fig. 9 H2 excitation diagram for DF3, with linear fits and comparison to Meudon model. Black crosses: column density of the upper level populations derived from the data, versus the energy of the upper level with respect to the ground state. Error bars are shown where they exceed 5% (includes 3% systematic uncertainty). Red circles: column densities derived from our Meudon PDR model. Blue line: fit including lines up to 0–0 S (7), excluding the 0–0 S (3) line (gray) as it is subject to extinction by the 10 µm silicate feature. Orange line: fit including 0–0 S (8) and 1–1 S (5) to 1–1 S (9). Temperatures derived from the slope of each line are shown in the legend.

In the text

Fig. 10 Top panel: structure of the HI/H2 transition predicted by Meudon model as function of depth. Densities of H and H2, with vertical line that indicates the transition defined as n(H) = n(H2). The local emissivity of each gas layer is shown for four representative H2 lines. Bottom panel: temperature and cumulative emissivity versus model depth. Black line and left y-axis: gas temperature for each layer in the model. Colorful lines and right y-axis: cumulative of the emissivity contributions by each layer (weighted by length scale). By inspecting the A(V) ranges over which the cumulative emissivity curves rise the most, we can see the typical range of gas temperatures from which the lines originate. For the 1–1 S (7) line for example, the observed excitation temperature is ~2700 K, while around 75% of the emission originates from gas cooler than 1000 K, due to UV pumping.

In the text