© The Authors 2024

1 Introduction

The processing of the interstellar medium (ISM) is an important driver of the evolution of galactic ecosystems. Massive stars profoundly shape the ISM, on small and large scales, by their radiation, stellar winds and at the end of their lives as super-novae. The concerted action of many massive stars leads to the formation of various phases of the ISM (e.g., Field et al. 1969; McKee & Ostriker 1977; Wolfire et al. 2003; Haffner et al. 2009): cold diffuse clouds – the cold neutral medium (CNM) –, a warmer intercloud phase – the warm neutral medium (WNM) and warm ionized medium (WIM) –, compact HII regions, and photodissociation regions (PDRs). While the CNM, WNM, and WIM fill large volumes, compact HII regions and PDRs are found either deeply embedded in or on the surfaces of dense molecular clouds in the vicinity of massive stars. Only in the latter case can they be readily observed at optical wavelengths. The study of HII regions and PDRs reveals insights into the fate of molecular clouds, the birth sites of stars, which may be compressed or destroyed. This allows estimates of the efficiency with which new stars are formed (for a review see Wolfire et al. 2022).

The extreme ultraviolet (EUV, E > 13.6 eV) radiation of massive stars creates an HII region around the star, where hydrogen (ionization potential IP = 13.6 eV) is ionized. The most massive stars can ionize helium (IP = 24.6 eV) and doubly ionize carbon (IP = 24.4 eV). The HII region borders on a molecular cloud and the far-ultraviolet radiation (FUV, 6 eV < E < 13.6 eV) of the star creates a PDR on the surface of the molecular cloud, where hydrogen is neutral, but carbon is still ionized (IP = 11.3 eV). The radiation being attenuated in the cloud, hydrogen transitions from atomic to molecular, and carbon transitions from singly ionized to neutral to being locked up in CO. The electron density in the cloud controls the chemistry, for instance, the formation of many molecules, and the coupling of matter to magnetic fields. The ionization fraction, x_e ≈ n_e/n_H, decreases from the PDR surface (x_e ~ 10⁻⁴) to the shielded molecular layers (x_e ~ 10⁻⁸) (e.g., Caselli et al. 1998; Maret & Bergin 2007; Goicoechea et al. 2009).

In the past, many studies have used the far-infrared ²P_3/2 -²P_1/2 fine-structure line of ionized carbon at 158 µm to study the surface layers of PDRs and the dynamics of the ISM (in Orion Stacey et al. 1993; Herrmann et al. 1997; Goicoechea et al. 2015b; Pabst et al. 2019, 2020, 2021, 2022). Recent developments of broadband receivers have made it possible to observe multiple radio recombination lines (RRLs) simultaneously and at angular and spectral resolutions comparable to that of the [C II] 158 µm line. RRLs are emitted when a free electron recombines with an ion and cascades down through Rydberg states to the ground state. While hydrogen and helium RRLs stem from the HII region, carbon RRLs are bright at the C⁺/C/CO transition in the PDR (Natta et al. 1994; Wyrowski et al. 1997b,a, 2000). While the [C II] line is often moderately optically thick toward PDRs, the carbon RRLs remain optically thin and can be used to infer electron densities and temperatures, which were previously largely unknown. In bright regions such as the Orion Bar additional RRLs of sulfur (IP = 10.4 eV) can be detected (Goicoechea & Cuadrado 2021).

The most prominent PDR for our eyes is the PDR created by the massive Trapezium stars in the Orion Nebula (M42). The Orion Nebula is the nearest site of massive star formation (d ≃ 414 pc Menten et al. 2007; Kounkel et al. 2018; Großschedl et al. 2018) and therefore studied readily at any wavelength in great detail. Besides the central nebula, also called the Huygens Region (e.g., O’Dell et al. 2020), the Orion Nebula complex encompasses the two H II regions M43 and NGC 1973, 1975, and 1977. The Orion Nebula complex hosts four star-forming dense molecular cores, of which the Orion Molecular Core 1 (OMC1) is strongly irradiated by the Trapezium stars in front, of which the O7V star θ¹ Ori C is the most massive. OMC1 contains several smaller regions of interest, among which the prototypical edge-on PDR of the Orion Bar and the Beckmann-Neugebauer/Kleinman-Low (BN/KL) object with massive star formation and strong molecular outflows (e.g., Gómez et al. 2005; Bally et al. 2011, 2017; Morris et al. 2016). Enshrouding the Trapezium stars with associated HII region in front of OMC1, several layers of ionized (emitting in tracers of ionized gas such as [O III], [O II], [O I], [S III], and [N II] optical forbidden lines) and neutral gas (emitting in [C II] and the optical [O I] forbidden lines, and absorbing in the H I hyperfine-structure line) have been detected throughout the years (e.g., van der Werf & Goss 1989; Wen & O’Dell 1993; O’Dell et al. 1993; Abel et al. 2004, 2006, 2019). On a scale of 4 pc, the Trapezium stars and the O9.5IV star θ² Ori A have created and are illuminating an expanding bubble filled with hot gas, the Extended Orion Nebula (EON) with the foreground neutral Veil Shell (Güdel et al. 2008; Pabst et al. 2019).

The H II region of M43, excited by the B0.5V star NU Ori, is surrounded by an expanding half-shell of neutral gas (Pabst et al. 2020), which abuts the OMC2. The H II region of NGC 1973, 1975, and 1977 also drives the expansion of a neutral shell (Pabst et al. 2020). NGC 1977 borders on the molecular core OMC3, where the stellar radiation of the B1V star 42 Ori creates an edge-on PDR. However, the physical conditions in the various parts of the cloud complex are not well-constrained. Also the dynamic relation between the ionized, neutral and molecular phases of the H II region with adjacent PDR is not understood in many cases, the different gas phases often being observed at slightly different velocities. In both cases, RRLs are the link from the ionized to the neutral gas (hydrogen and helium RRLs), and from the neutral to the molecular phase (carbon RRLs).

We here present a broadband study of mm-wave RRLs toward the Orion Nebula complex, including M43 and NGC 1973, 1975, and 1977, in a series of observations of 10 representative positions using the 40 m Yebes antenna. Earlier studies of carbon RRLs focussed on the brightest part of the Orion Nebula (Balick et al. 1974; Wyrowski et al. 1997a; Tsivilev 2014; Salas et al. 2019; Cuadrado et al. 2019), while we include fainter regions that are representative for the extended outskirts of the Orion Nebula complex. Our fainter observed positions are situated in dynamic regions of interest, previously studied in other tracers, such as the EON southeast of the Orion Bar, the border between OMC1 and M43, M43 itself, and OMC3, that are bright enough to be detectable in RRLs with state-of-the-art telescopes (Salas et al. 2021 include also OMC2 and OMC3). In particular, these regions have been studied in detail in [C II] line emission in Pabst et al. (2020, 2021, 2022) and our selection is based on the results of these studies. We discuss the dynamics of the PDR gas compared to the H II gas and the turbulent motions. Furthermore, we aim to determine the physical conditions in the respective cloud environment from a comparison of the carbon RRLs with the [C II] lines using nonlocal thermal equilibrium (non-LTE) models.

This paper is organized as follows. In Sect. 2, we summarize the observations used in this paper. Section 3 compares the RRL, [C II] and CO line properties. Section 4 contains a discussion of the implications of the results obtained in Sect. 3 and we derive physical conditions (electron temperatures and densities) from the observed lines. We conclude with a summary of our results in Sect. 5.

2 Observations

2.1 Yebes RRL observations

We used the 40 m radiotelescope at Yebes Observatory (Guadalajara, Spain) to observe ten representative positions throughout the Orion Nebula complex (summarized in Table 1) in the Q band (31.3–50.6 GHz; project 22A012, PIs: C. H. M. Pabst, J. R. Goicoechea). We employ the new Nanocosmos high-electron-mobility transistor (HEMT) Q band receiver and fast Fourier transform spectrometers, covering 18 GHz of instantaneous bandwidth per polarization at a spectral resolution of 38 kHz (i.e., a typical velocity resolution of 0.3 km s⁻¹) (Tercero et al. 2021). The half-power beam width (HPBW) in the Q band ranges from 36″ to 54″. We observed the 10 positions in position switching mode with a reference position at (Δα, Δδ) = (−800″, 1200″) relative to the central position of the nebula, θ¹ Ori C, (RA, Dec) = (5h35^m16.47^s, −5°23^m22.90^s).

We reduced the data using the GILDAS software¹. The antenna temperature $T_{\text{A}}^{*}$ was converted to the main-beam temperature by $T_{\text{mb}}=T_{\text{A}}^{*}{F}_{\text{eff}}/{\eta }_{\text{mb}}$ with the main-beam efficiency η_mb between 0.49 and 0.67 (higher at lower frequencies²) and the forward efficiency F_eff = 0.967. The root-mean-square (rms) noise in the main-beam temperature achieved outside the line window is usually between 4 and 12 mK per velocity channel (0.37 km s⁻¹) in the C59α line and between 11 and 23 mK per velocity channel (0.23 km s⁻¹) in the C51α line, after total observing times of 8–18 h (depending on position; on-source times in position-switching mode and resulting rms for the C59α and C51α RRLs are given in Table 1), varying with telescope conditions.

In order to detect faint features in the spectra, we stacked RRLs of same quantum leap An and similar principal quantum number n weighted by the rms noise, as is traditionally done for RRL spectra (Balser et al. 2006; Emig et al. 2019). Individual spectra, i.e. before stacking, of positions toward the molecular cores in OMC1, BNKL-PDR and TRAP-PDR (cf. Sect. 3 for a summary of the observed positions), have additional features from molecular emission, which we discard in the present study. We stacked the n = 51–59 α transitions, the n = 64–74 β transitions, the n = 76–84 γ transitions, excluding n = 74 and n = 75 due to bad baselines, and n = 80 (blended with 56α line) and n = 83 (blended with 91δ line), the n = 81–92 δ transitions, excluding n = 82 (blended with 88є line) and n = 91 due to bad baselines, and the n = 87–99 є transitions, excluding n = 88 (blended with 82δ line). We subtracted standing wave features before stacking in some spectral chunks. The resulting velocity resolution after stacking is 0.36km s⁻¹ with a rms noise per channel of 1–2 mK.

2.2 SOFIA [C II] observations

We complement the RRL observations with existing observations of the [C II] 158 µm line toward the Orion Nebula complex (Pabst et al. 2019). These velocity-resolved observations (covering the velocity range from −90 to 95 km s⁻¹ at a velocity resolution of 0.3 km s⁻¹) were obtained as a square-degree sized, fully sampled map by the upgraded German Receiver At Terahertz Frequencies (upGREAT, Risacher et al. 2016) onboard the Stratospheric Observatory for Infrared Astronomy (SOFIA, Young et al. 2012) with a native angular resolution of 14.1″. Details on the observing strategy and data reduction are given in Higgins et al. (2021). The resulting velocity range also contains the three hyperfine components of the [¹³C II] line. We convolve the map to a beam of 45″ to compare with the Yebes observations.

2.3 IRAM CO(2–1) observations

We also make use of ¹²CO J =2–1 (230.5 GHz) and ¹³CO J = 2–1 (220.4 GHz) line maps taken with the IRAM 30 m radiotelescope (Pico Veleta, Spain) at a native angular resolution of 10.7″. The central region (1º × 0.8º) around OMC1 was originally mapped in 2008 with the HERA receiver array (Berné et al. 2014). These fully sampled CO maps were enlarged using the EMIR receiver and FFTS backends, being part of the Large Program “Dynamic and Radiative Feedback of Massive Stars”. Goicoechea et al. (2020) provide details on the observing strategy and on how the old HERA and new EMIR CO maps were merged. Here we present the final maps, that include NGC 1973, 1975, and 1977, the last observations having been taken in October 2020. The resulting maps have a velocity resolution of 0.4 km s⁻¹. We convolve the maps to a beam of 45″ to compare with the Yebes observations.

3 Results

Table 1 gives the definitions of the positions observed in RRLs with the Yebes telescope with coordinates and a classification of the respective exciting star. Figures 1 and 2 show the locations of those positions on the SOFIA [C II] and IRAM ¹²CO(2–1) and¹³ CO(2–1) maps of the Orion Nebula complex at an angular resolution of 45″. We observed four positions toward the bright OMC1, BAR-INSIDE, BNKL-PDR, TRAP-PDR, and EAST-PDR. The position BAR-INSIDE covers the CO-bright part of the Orion Bar and the background molecular cloud southeast of the Orion Bar within our beam (average of 45″). BNKL-PDR is aimed at the BN/KL object and its surface PDR. TRAP-PDR is pointed at the Trapezium cluster and covers the bright PDR/clump west of the Trapezium stars that is the surface of the Orion South cloud. EAST-PDR covers the bright PDR east of OMC1, toward the (optically) Dark Lane. In addition, we observed six positions toward fainter regions outside OMC1. Pointing ORI-P1 lies on the northeastern border of OMC1 at the far side of the Dark Lane, ORI-P2, ORI-P3, and ORI-P4 are situated in M43, ORI-P5 lies southeast of the Orion Bar in the cavity of the EON, and ORI-P6 lies toward OMC3. Appendix A contains zoom-ins toward OMC1, M43, and OMC3, that show the structure averaged within each beam on the SOFIA [C II] intensity map and the Spitɀer 8 µm intensity map in their original resolution.

Figures 3 and 4 show the stacked spectra of the observed carbon and helium RRLs and hydrogen RRLs, respectively, of the α, β, γ, δ, and e transitions (with mean $\overline{n}=55$ for the α transitions, $\overline{n}=69$ for the β transitions, $\overline{n}\approx 80$ for the γ transitions, $\overline{n}\approx 87$ for the δ transitions, and $\overline{n}\approx 93$ for the є transitions). The velocity separation due to different atomic weight of the Rydberg atoms between helium RRLs and carbon RRLs is 27.4 km s⁻¹, and between hydrogen RRLs and carbon RRLs the velocity separation is 149.4 km s⁻¹. Hence, the helium RRLs appear as a broad component with the carbon RRLs on its shoulder, while the carbon and helium RRLs are well-separated from the hydrogen RRLs.

In the brightest positions (all four positions in OMC1) we detect the α, β, γ, δ, and є transitions of the hydrogen, helium and carbon RRLs. In ORI-P1 and ORI-P5 we detect the hydrogen, helium and carbon RRLs in higher transitions than the α transition. ORI-P2 has significant emission of hydrogen and carbon RRLs, but no helium RRL emission. ORI-P3 has weak carbon RRLs, but strong hydrogen RRL emission and we detect the α helium RRL. ORI-P4 has weak hydrogen and carbon RRL emission, but no significant helium RRL. In ORI-P6 we detect only the strong carbon RRL and no hydrogen or helium RRL emission. Tables 2, 3 and 4 give the observed integrated signal-to-noise ratios from Gaussian line fits for the stacked carbon, hydrogen, and helium RRLs, respectively. Extrapolating from the tabulated signal-to-noise ratios, we would be able to detect CRRLs after stacking up to ∆n = 7, HRRLs up to ∆n = 37 and HeRRLs up to ∆n = 10 in the brightest position TRAP-PDR. Tables 5, 6, and 7 give the results of the Gaussian fits for the carbon, hydrogen and helium α RRLs. Appendices B, C, and D contain tables giving the results of Gaussian fits to the detected β, γ, δ, and є RRLs of carbon, hydrogen, and helium.

Figure 5 shows the stacked Cnα RRL, [¹²C II] and ¹³CO(2– 1) lines toward the observed positions in the Orion A molecular cloud. We use the¹³ CO(2–1) line for the analysis and discussion throughout the paper as it traces the depth of the CO-emitting layer in the PDR, while the ¹²CO(2–1) line is optically thick and the C¹⁸O layer traces only the more FUV-shielded and densest parts of the molecular cloud. Figure 6 shows the [¹²C II] and [¹³C II] lines comparatively, while Fig. 7 shows the ¹²CO(2– 1) and¹³ CO(2–1) spectra in each observed position. In the following sections we describe the velocity structure of the Cnα RRL, [C II] and ¹³CO emission lines and the intensity ratios of the carbon RRLs and the [¹²C II] line. We also compare the velocity profiles of the carbon RRLs with those of the hydrogen RRLs. Fig. 8 shows the stacked Cnα RRLs and the stacked Hnα RRLs together, the velocities centered on their respective rest frequency.

Fig. 1 Observed positions (except positions in OMC1) as white circles with a FWHM of 45″ on [12C II], 12CO(2–1), and 13CO(2–1) maps at 45″. The black rectangle outlines the region of OMC1, shown as a zoom-in in Fig. 2. The red circle outlines M43. The yellow and orange stars mark the positions of θ1 Ori C and θ2 Ori A, respectively, the pink star marks NU Ori, the purple star marks 42 Orionis.

Fig. 2 Zoom into OMC1. Observed positions in OMC1 as turquois circles with a FWHM of 45″ on [12C II], [13C II], 12CO(2–1), and 13CO(2–1) maps at 45″. The yellow and orange stars mark the positions of θ1 Ori C and θ2 Ori A, respectively.

Fig. 3 Stacked Cnα, Cnβ, Cnγ, Cnδ, and Cnє radio-recombination lines with Henα, Henβ, Henγ, Henδ, and Henє radio-recombination lines toward the ten targeted positions. In the position BAR-INSIDE the component 7.5 km s−1 blue of the main component may be a radio-recombination line from sulfur, while detection of the sulfur RRL in other positions is hampered by the carbon RRL component corresponding to the Veil (cf. Sect. 4.1). Spectra in in the upper top are offset from the C/Henє spectrum in steps of 0.025 K, sand pectra in the lower two rows are offset in steps of 0.01 K.

3.1 Velocity structure

In order to obtain line parameters (peak temperature, peak velocity, line width) we perform Gaussian fits to the observed stacked α RRLs, the [C II] lines and the ¹³CO lines. Results of the Gaussian fits are given in Tables 8–7. Where we fit multiple components, we append a number to the name of the position. The position BAR-INSIDE is special in this regard, as the main components of the α CRRL, [¹²C II] and¹³ CO lines seems to consist of two close Gaussian components, rather than one. However, the [¹³ C II] line is too weak to extract two components and we opt to use a one-component fit for the velocity structure and line intensity analysis, where we need the [¹³C II] line. The peak velocity of the [¹³C II] line very nearly matches the peak velocity of the one-component fit to the [¹²C II] line in BAR-INSIDE, and we assume that the emission is weighted toward the higher-column density region of the CO-bright Orion Bar within the beam. Furthermore, Kabanovic et al. (in prep.) observed a velocity shift between the peaks of the [¹²C II] and [¹³C II] line in OMC3, and we follow them in performing a fit taking into account self-absorption in the [C II] line in ORI-P6, using the recipe of Guevara et al. (2020). We give the peak velocity of the lines relative to the respective rest frequency in the local-standard-of-rest (LSR) frame. While the statistical errors (given in the tables) are small, we estimate that the systematic errors of the measured intensities with the Yebes 40 m telescope are about 10%.

Studying the obtained line parameters, the FIR [C II] lines tend to be somewhat broader than the mm-wave Cnα lines, apart from ORI-P3 and ORI-P4 (associated with M43). The [¹²C II] line in PDRs tends to be marginally optically thick and is therefore expected to be opacity broadened. However, as observed by Ossenkopf et al. (2013), opacity broadening in itself is not sufficient to explain the difference in line width between the [¹²C II] and the [¹³C II] lines. Comparing the line widths of the [¹³C II] lines with the Cnα lines, reveals different behaviors: Some [¹³C II] lines are much narrower than the Cnα lines (e.g., BNKL-PDR, TRAP-PDR, ORI-P1-2), while others are similar (e.g., BAR-INSIDE, EAST-PDR-2, and ORI-P6). In BNKL-PDR, ORI-P1, ORI-P3, and ORI-P5 line parameters of the [¹³C II] line are very different when using the F=2–1 component compared to when adding all three hyperflne components. Due to the faintness of the two far hyperflne satellites (F =1–0 and F=1–1), we opt to fit only the F=2–1 satellite, apart from BNKL-PDR, where the broad line wings heavily contaminate the F=2–1 satellite and the F=1–0 and F=1–1 satellites are strong enough to be detected individually. We observe velocity shifts between the [¹²C II] and the [¹³C II] line in ORI-P2, ORI-P3 and ORI-P5.

Also the ¹²CO line is optically thick and thus opacity broadened, while the ¹³CO line is usually optically thin. Adopting T_p(¹³CO)/T_P(¹²CO) ≃ 1 − exp [−τ (¹³CO)], the optical depth of the ¹³CO line is estimated to be less than 0.7 in our sample. Hence, we use the ¹³CO lines to extract velocity information. The ¹³CO lines are somewhat narrower than the Cnα lines in most positions, a notable exception being the second component of ORI-P1 (ORI-P1-2). In ORI-P3 and ORI-P4 the main¹³ CO lines seem to trace the molecular background (cold, with narrow lines), while the Cnα lines in these positions are much broader. Many ¹³CO lines have wings or a second component. Most notably, the ¹³CO lines in BNKL-PDR have very broad wings from high-velocity gas in BN/KL outflows, that are also present to a lesser degree in the [C II] line, but not pronounced in the Cnα line. In TRAP-PDR, ORI-P1-2, and ORI-P5 the ¹²CO and¹³ CO lines are observed at somewhat different velocities, which might be an opacity effect.

In BAR-INSIDE, the Cnα line has a second faint component at υ_LSR = 3.0 ± 0.2 km s⁻¹ besides the main component at 10.6 km s⁻¹, but the [C II] line has a second fainter component at −2.3 ± 0.2 km s⁻¹. The Cnα line in TRAP-PDR has a faint component at 0.9 ± 0.4 km s⁻¹ besides the main component at 9.2 km s⁻¹, while the [C II] line has a weak component at 0.5 ± 0.1 km s⁻¹. The lines in EAST-PDR and ORI-P1 have two components, separated by 6 km s⁻¹ and 3 km s⁻¹, respectively. The lines in ORI-P2, ORI-P4, ORI-P5 and ORI-P6 are approximately fitted with a single emission component, with a peak velocity increasing from 8.7 km s⁻¹ to 11.4 km s⁻¹ .In ORI-P3 neither the peak velocity (8.2 km s⁻¹) nor the line width of the Cnα line match the profile of the [¹²C II] line, that has two peaks (at 5.7 km s⁻¹ and 10.3 km s⁻¹) and a broad component between the two peaks. The [¹³C II] line is almost coincident with the red component of the [¹²C II] line. The ¹³CO line in ORI-P3 exhibits non-Gaussian wings.

In summary, the peak velocities of the [C II], ¹³CO, and Cnα lines differ slightly in all positions, except in BNKL-PDR where they are all similar. Figure 9 shows the peak velocities obtained from Gaussian fits to the lines. We observe five different behaviors (or classes) in the layering of the peak velocities of the main components, from lower to higher LSR velocity (excluding ORI-P3):

1. Velocity layering [C II]-Cnα-CO (EAST-PDR-1, ORI-P6);

2. Velocity layering CO-Cnα/[C II] (EAST-PDR-2, ORI-P1-2, ORI-P2);

3. Velocity layering CO-[C II]-Cnα (BAR-INSIDE, TRAP-PDR, ORI-P1-1);

4. Velocity layering [C II]/Cnα-CO (ORI-P4);

5. Velocity layering [C II]/CO-Cnα (ORI-P5).

Comparing the Cnα and Hnα RRLs toward the observed positions, we note that toward OMC1 the neutral gas and the ionized gas move at very different velocities, as has been known for a long time and repeatedly studied in increasing detail (e.g., Balick et al. 1974; O’Dell et al. 2009, 2020; Goicoechea et al. 2015b; Abel et al. 2019). This behavior is not linked to the classification above. The ionized gas is blue-shifted by 14.8 km s⁻¹ toward the Orion Bar and by 12.0 km s⁻¹ toward BNKL-PDR and TRAP-PDR. Toward EAST-PDR the velocity shift between the main component of the neutral gas and the ionized gas is 13.6 km s⁻¹. Also in ORI-P5, southeast of the Orion Bar, the velocity shift is 14.8 km s⁻¹. Opposed to that, the ionized gas is blueshifted by only 1.8 km s⁻¹ toward ORI-P3, the center of M43. Here, however, the CRRL is blueshifted from the [C II] line that traces the background PDR and we do not observe a CRRL component that corresponds to the background velocity. In ORI-P1 the ionized gas is blueshifted by 3.4 km s⁻¹ with respect to the brighter neutral component and by 6.6 km s⁻¹ with respect to the fainter component. Toward ORI-P6, OMC3, we do not detect ionized gas.

We do not observe systematic velocity differences between the Cnα, Cnβ and Cnγ lines. Also in the Hnα, Hnβ and Hnγ lines there is no systematic ordering of the peak velocities. The velocity shifts are usually less than 1 km s⁻¹ and do not increase with quantum leap Δn, but scatter about zero velocity shift.

Fig. 4 Stacked Hnα, Hnβ, Hnγ, Hnδ, and Hnє radio-recombination lines toward the ten targeted positions. Spectra in in the top row are offset from the Hnє spectrum in steps of 0.15 K, and spectra in the lower two rows are offset in steps of 0.01 K.

Fig. 5 Stacked Cnα RRL, [12C II], and 13CO(2–1) lines toward the 10 targeted positions. The 13CO(2–1) line is scaled by the factor indicated in light gray in each panel.

Fig. 6 [12C II] and [13C II] F=2–1 lines toward the 10 targeted positions.

Fig. 7 12CO(2–1) and 13CO(2–1) lines toward the 10 targeted positions.

Fig. 8 Stacked Cnα (blended with Henα) with respect to the Cnα rest frequency and stacked Hnα RRLs with respect to the Hnα rest frequency toward the ten targeted positions.

Fig. 9 Peak velocities in the observed positions of the [C II], Cnα, and 13CO lines obtained from Gaussian fits. The background shading indicates the velocity-ordering class, defined in the main text, the position and component belongs to: red is class 1, green is class 2, blue is class 3, yellow is class 4, and purple is class 5.

Fig. 10 Intensity ratios of Cnα over Cnβ and Cnγ, respectively (left axis) and [12C II] over Cnα (right axis).

3.2 Intensity ratios

Figure 10 shows the intensity ratios of the carbon RRLs and the intensity ratio of the [C II] line over the Cnα line (intensity units are K km s⁻¹). Of all positions, BNKL-PDR and TRAP-PDR have the lowest [C II]/Cnα intensity ratio with 1.1 × 10³ and 1.3 × 10³, respectively. ORI-P1 also has a low [C II]/Cnα intensity ratio with 1.6 × 10³ The majority of the regions has [C II]/Cnα ratios between 2 × 10³ and 4 × 10³. Only in ORI-P4 and ORI-P5 the ratios are elevated at about 5.5 × 10³ and 9.1 × 10³, respectively, but the errors are large. The blue component of EAST-PDR has an elevated [C II]/Cnα ratio compared to the other PDRs in OMC1 and much elevated CRRL intensity ratios, but again the errors are large.

The mean of the Cnα/Cnβ intensity ratios of all positions is 3.6 ± 1.2, without the outlier EAST-PDR-1 with a ratio of 9 ± 4. ORI-P1-2 and ORI-P2 have elevated ratios of 5.5 ± 1.8 and 5.4 ± 1.1, respectively, while ORI-P5 has a lower ratio of 1.5 ± 1.0.

The Cnα/Cnγ intensity ratios in OMC1 lie close to the mean ratio of 8.0 ± 0.7. ORI-P1-1 and ORI-P2 are outliers with Cnα/Cnγ ratios of 12 ± 5 and 10 ± 4, respectively, but due to their faintness the errors are also largest. In ORI-P1-2 the Cnγ intensity appears to be equal to the Cnβ intensity (Cnα/Cnγ intensity ratio of 6 ± 3), but the errors in the fits allow for the ‘normal’ intensity ordering. The double-peaked structure of ORI-P1 is less pronounced in the Cnβ line, but very pronounced in the Cnγ line. ORI-P6 has a somewhat lower ratio of 6.4 ± 1.8.

Figure 11 shows the integrated intensity ratios of the hydrogen RRLs, where detected. We observe little variation in the HRRL line ratios between the positions. Mean values with standard deviations are 3.24 ± 0.06 for the Hnα/Hnβ intensity ratio, 7.1 ± 0.5 for the Hnα/Hnγ intensity ratio, 12.6 ± 1.3 for the Hnα/Hnδ intensity ratio, and 16.9 ± 1.1 for the Hnα/Hnє intensity ratio. Where the line are fainter and errors are larger (ORI-P2, ORI-P5), we observe the most deviation from the mean values.

Fig. 11 Intensity ratios of Hnα over Hnβ, Hnγ, and Hnδ, respectively.

4 Discussion

4.1 Velocities of the [C II], Cnα and CO-emitting gas

In order to interpret the velocity structure of our observations, we first run representative PDR models. Figure 12 shows the stratification of emission layers in a PDR as computed by the Meudon PDR code (Le Petit et al. 2006) for an incident radiation field of G₀ = 100 (appropriate for regions far from OMC1) and three different densities (n = 10³, 10⁴, 10⁵ cm⁻³), including ¹³C fractionation. We include the emission lines and abundances relevant for this discussion. Typically, the [C II] line is bright at the surface of a PDR, where the gas is warmest, but extends to the cooler gas at A_V ~ 2. The C54α line, which we observe with the Yebes 40 m telescope, peaks somewhat deeper into the PDR, as the optical depth is weighted (∝ T^−2.5) toward cooler C⁺-rich gas (A_V ≃ 1–2). Beyond A_V ≃ 2 carbon is locked up in CO and low-J CO rotational lines are excited. In lower-density models (e.g., n = 10³ cm⁻³), the regions of bright emission of the [C II] and C54α line are slightly broader in A_V space than at higher density (e.g. n = 10⁵ cm⁻³), while CO emission sets on closer to the surface, as the C⁺/C/CO transition is pulled toward the lower A_v for higher densities.

Considering the typical stratification of a PDR, as expected from the models, the velocity layering expected for a face-on PDR, with a photoevaporative flow from the surface, is [C II]-Cnα-CO (blue to red, class 1 in our classification). Photoevaporative flows from PDRs are expected to have velocities of about 1 km s⁻¹ (Bertoldi 1989; Störzer & Hollenbach 1999; Bron et al. 2018; Maillard et al. 2021). However, the layering expected for face-on PDRs (such as TRAP-PDR and ORI-P5) we encounter only in two edge-on PDRs, ORI-P4 and ORI-P6. For perfectly edge-on PDRs, we do not expect a radial velocity shift. ORI-P1 likely has one face-on and one edge-on component (cf. Sect. 4.2), while ORI-P2 on the edge of M43 is likely observed edge-on and ORI-P3 in the center of M43 face-on. With the geometries of the cloud in the positions we observe it is unclear what causes the observed velocity differences between the lines. In most positions, we can only speculate about the origins. The velocity ordering of class 1 (EAST-PDR-1, ORI-P6) can be explained by assuming that the illuminating star is located somewhat in front of the molecular cloud, while on the other hand the class-2 behavior (EAST-PDR-2, ORI-P1-2, ORI-P2) would mean that the cloud is situated closer to the observer than the illuminating star.

Furthermore, it is likely that the cloud parts move with respect to one another due to dynamics induced by stellar feedback from the most massive stars at earlier evolutionary stages, among which fossil outflows (Kavak et al. 2022b), but also gravitational infall (Hacar et al. 2017). We may thus observe the remnants of these earlier dynamics as velocity differences between the layers of the PDR. Also the dynamics currently induced by the winds and radiation of the massive stars in the Orion Nebula complex (Pabst et al. 2019, 2020) and protostellar jets of less massive objects (e.g., Méndez-Delgado et al. 2021; Kavak et al. 2022a) can possibly lead to shear in the cloud parts due to the irregular shapes of the clouds and their relative orientation with respect to the stars. It is also possible that velocity gradients within our beam get weighted differently in the different emission lines. This might be the case in particular in BAR-INSIDE, where we average emission from the bright Orion Bar with emission from the background molecular cloud southeast of the Orion Bar, which is known to be slightly blue-shifted with respect to the Orion Bar (Tauber et al. 1994; Cuadrado et al. 2017).

The observed velocity differences between the [¹²C II] and the Cnα lines are less than 1 km s⁻¹, except in ORI-P3, with a median of the absolute values of all positions of 0.25 km s⁻¹. In ORI-P3 the observed Cnα line does not have line parameters that match either the foreground or the background neutral component, but it could also be a blend of multiple components that we are unable to identify due to the low S/N ratio. The observed velocity differences between the¹³ CO and the Cnα lines scatter between 0 and 2 km s⁻¹ with a median of 0.75 km s⁻¹. The spectrum of ORI-P3 does have a weak ¹³CO component with a peak velocity matching the Cnα line, but with very different line widths; it might also be a wing of the main component. The main ¹³CO component corresponds to the main [¹²C II] component of the background PDR.

The lines in BAR-INSIDE show slight deviations from Gaussianity, which can be leveraged by fitting the line with two components. The Orion Bar is the prototypical edge-on PDR close to the Trapezium stars. Our beam averages emission from the narrow Orion Bar and emission from the background molecular cloud southeast of the Orion Bar. The line is dominated by the bright edge-on PDR. Also in other positions the velocity shifts may be caused by averaging multiple components with varying strengths in the beam (without leading to much deviation from Gaussianity, however). Given our relatively large beam and the systematic uncertainties in performing Gaussian fits to the observed lines, we refrain from drawing strong conclusions from the small (but significant by eye inspection) observed peak velocity differences between the different tracers.

The broad wings in BNKL-PDR are the result of molecular outflows and shocks produced by proto-stellar activity in the molecular core, likely the dynamical decay of a multiple system (Gómez et al. 2005; Bally et al. 2011), visible in both molecular lines (Rosenthal et al. 2000; González-Alfonso et al. 2002; Goicoechea et al. 2015a; Bally et al. 2017) and the [C II] line (Goicoechea et al. 2015b; Morris et al. 2016). Relatively little C⁺ is produced in shocks (Hollenbach & McKee 1989), hence the [C II] line and the CRRLs are relatively unaffected and mostly trace the PDR gas (main peak) and the FUV-illuminated surface of the shocked gas (wings). In contrast, the [O I] 63 µm line and CO lines are heavily affected by the shocks and outflows (Stacey et al. 1993; Berné et al. 2014), consisting of a PDR component with moderate temperature and a hot component likely produced from a J-type shock (Rosenthal et al. 2000; González-Alfonso et al. 2002; Goicoechea et al. 2015a).

Position TRAP-PDR lies in the region that has been extensively studied in optical wavelengths (e.g., O’Dell et al. 1993, 2020; Abel et al. 2019). The velocity components now known as Orion’s Veil were first detected in H I 21 cm absorption (van der Werf & Goss 1989). In [C II] and CRRL emission, we observe only the PDR behind the Trapezium stars (υ_LSR ≃ 9.1 km s⁻¹) and the Veil component B (υ_LSR ≃ 0.5 km s⁻¹). The CRRL Veil component, however, seems to lie on a broader line (wing), while the [C II] Veil component is approximately Gaussian.

In EAST-PDR and ORI-P1 we observe two pronounced velocity components. While the weaker (blue) component of EAST-PDR probably corresponds to the Dark Lane in the foreground of the molecular cloud with its main PDR, the position ORI-P1 lies close or on the border between OMC1 and M43. Possibly this part of the cloud is a dense wall, through which the hot plasma created by the Trapezium stars cannot escape. The presence of two velocity components may give rise to the assumption that this part of the cloud is created by a cloud-cloud collision between the clouds containing OMC1 and the cloud containing M43. However, the molecular cores OMC1 and OMC2 close to M43 both belong to the continuous structure of the Integral-Shaped Filament (ISF), the dense molecular spine of the Orion Nebula complex.

In each position in OMC1, we observe a faint blue-shifted component in the [C II] line at small positive or negative velocities, that most likely corresponds to the Veil in positions other than BAR-INSIDE. In the Cnα line toward BAR-INSIDE we observe a weak second component, blue-shifted from the main component by 7.5 ± 0.5 km s⁻¹, that does not have a corresponding [C II] component. While the velocity observed here does not match exactly the velocity of a sulfur RRL (the velocity shift between C and S should be −8.6 km s⁻¹), an earlier study of the Orion Bar edge has reported the detection of sulfur RRLs toward this position (Goicoechea & Cuadrado 2021). The intensity ratio between the main component and the blue-shifted component is about 10 and the line widths are similar. Hence, we tentatively conclude that this component is indeed a sulfur RRL with a similar intensity ratio compared to the carbon RRL even somewhat deeper inside the cloud. We do not detect such a faint component without [C II] counterpart in the other positions: In TRAP-PDR and BNKL-PDR the sulfur RRL may be hidden in the component corresponding to the [C II] Veil component at υ_LSR = 0.5 km s⁻¹ and υ_LSR = 1.7 km s⁻¹, respectively, while in EAST-PDR the sulfur RRL may be hidden in the blue strong component (EAST-PDR-1) of the carbon RRL. In positions outside OMC1, the S/N ratio is not sufficient to detect the potential sulfur lines with ~ 10 times lower intensity in the Cnα spectra.

Fig. 12 PDR model outputs from the Meudon PDR code for an incident radiation field of G0 = 100 and three different densities (n = 103, 104, 105 cm−3). The upper panel shows the normalized intensity, assuming LTE for the CRRL emission, the lower panel shows the abundances of the relevant species.

4.2 Physical conditions in the neutral gas

The comparison of the carbon RRLs with the [C II] line emission allows for a first-order estimate of the physical conditions in the neutral, not yet fully molecular PDR, in particular the electron density n_e and the electron temperature T_e. Since the CRRL intensities and the [C II] intensities scale differently with density and temperature, the two intensities can be used to derive the physical conditions in the gas, assuming that the emissions stem from the same regions in the PDR. This latter assumption is not exactly valid as the CRRLs come from closer to the molecular gas, offsetting the derived values from the true physical conditions. The necessary theoretical framework to compute CRRL intensities given physical conditions was developed by Natta et al. (1994); Smirnov et al. (1995); Tsivilev (2014); Salgado et al. (2017a,b); Salas et al. (2021). The non-LTE excitation model computes the Cn integrated intensity for Δn = 1, 2, 3, 4, 5, and the [¹²C II] and [¹³C II] integrated intensities, taking into account optical depth effects, on a grid of T_e = 20–1000 K, n_e = 0.01– 200 cm⁻³, and C⁺ column densities log₁₀ N(C⁺ [cm⁻²]) = 17.5– 19.0. For n we take the mean principal quantum number of the stacked lines. The model assumes that hydrogen nuclei are equally distributed in atomic hydrogen and molecular hydrogen (i.e., the hydrogen nuclei fractions are ƒ(HI) = ƒ(H₂) = 0.5 with ƒ(HI) = n_HI/n_H and $f\left({\text{H}}_2\right)=2n_{{\text{H}}_2}/n_{\text{H}}$) and that the hydrogen nuclei density is given by n_H = n_e/𝒜_C with 𝒜_C = 1.4 × 10⁻⁴ the carbon gas-phase abundance in Orion (Sofia et al. 2004), i.e. all electrons stem from carbon ionization. We also assume that the gas kinetic temperature equals the electron temperature. We interpolate the model grid using a spline interpolation.

The peak temperature of the [¹²C II] line for an extended source is given by $$T_{\text{P}}=\left(J\left(T_{\text{ex}}\right)-J\left(T_{\text{bg}}\right)\right)\left[1-{\text{e}}^{-{\tau }_{[\text{CII}]}}\right]$$(1) $$=\left(\frac{91.2 \text{K}}{e^{91.2 \text{K}/T_{\text{ex}}}-1}-\frac{91.2 \text{K}}{{\text{e}}^{91.2 \text{K}/T_{\text{bg}}}-1}\right)\left[1-{\text{e}}^{-{\tau }_{[C\text{II}]}}\right],$$(2)

where 91.2 K is the energy level separation of the [¹²C II] fine-structure transition, T_ex is the excitation temperature of the [¹²C II] transition, T_bg is the temperature of a background permeating the [C II] emitting region, and τ_{[C II]} is the optical depth of the transition (Goldsmith et al. 2012). The background temperature T_bg generally arises from the cosmic microwave background (CMB) and the dust FIR field, external or internal to the [C II] emitting region. In the case of the [C II] line, the CMB is unimportant and the dust FIR field with temperature T_d is attenuated by the dust optical depth τ_d, such that $J\left(T_{\text{bg}}\right)=J\left(T_{\text{d}}\right)\left(1-{\text{e}}^{-{\tau }_{\text{d}}}\right)$. In our case, the warm dust is intermingled with the [C II] emitting region in the PDR. In OMC1, the continuum brightness temperature at 158 µm T_c = J(T_bg) has been observed to range between 1 and 7 K, thus also being negligible (Goicoechea et al. 2015b). The [¹³C II] peak temperature can be computed from the same equation, where the optical depth is scaled with the ¹²C/¹³C iso-topic ratio. In Orion, we assume a ¹²C/¹³C the isotopic ratio of 67 (Boreiko & Betz 1996).

The [¹²C II] excitation temperature is related to the kinetic gas temperature T_kin by $${\text{e}}^{91.2 \text{K}/T_{\text{ex}}}=\left(1+\frac{n_{\text{cr}}}{n}\right){\text{e}}^{91.2 \text{K}/T_{\text{kin}}},$$(3)

where n_cr is the critical density of the transition and n the gas density (cf. Eq. (13) in Goldsmith et al. 2012, with G = 0). If there are multiple collision partners, that collisionally excite the transition (electrons, hydrogen atoms, and hydrogen molecules), we define the gas density as the hydrogen nuclei density $n=n_{\text{H}}+2n_{{\text{H}}_2}$, and $$\frac{n_{\text{cr}}}{n}=\beta \left({\tau }_{[\text{CII}]}\right)\frac{A_{\text{ul}}}{C_{\text{ul}}},$$(4)

where $\beta \left({\tau }_{[\text{CII}]}\right)=\left(1-e^{-{\tau }_{[\text{CII}]}}\right)/{\tau }_{[\text{CII}]}$ is the 158 µm photon escape probability, A_ul is the Einstein coefficient for spontaneous radiative de-excitation, and $C_{\text{ul}}={\gamma }_{\text{ul},\text{e}}{n}_{\text{e}}+{\gamma }_{\text{ul},\text{H}}{n}_{\text{H}}+{\gamma }_{\text{ul},{\text{H}}_2}{n}_{{\text{H}}_2}$ is the combined collision de-excitation rates with coefficients γ_i· (see Goldsmith et al. 2012, for parametric expressions).

The peak temperature of a high-frequency carbon RRL (v > 2 GHz), under the assumption that the line is optically thin and there is no background, is given by (cf. Salas et al. 2019) $$T_{\text{P}}={\tau }_{\text{CRRL}}{b}_{\text{u}}{T}_{\text{e}}\text{, }$$(5)

where τ_CRRL is the optical depth of the CRRL in LTE and b_u = N_u/N_u,LTE is the departure coefficient of the upper level of the CRRL transition, defined as the ratio of the level population N_u over the equilibrium level population N_u,LTE (Shaver 1975). Departure coefficients were computed using the equations of Salgado et al. (2017a). The peak temperature of the CRRL scales approximately with $n_{\text{e}}^2{T}_{\text{e}}^{-1.5}$ (Salgado et al. 2017b). The line-integrated optical depth is given by $${\tau }_{\text{CRRL }}\Delta v=1.069\times {10}^7\Delta nM{T}_{\text{e}}^{-2.5}{\text{e}}^{{\chi }_n}{n}_{\text{e}}N\left({\text{C}}^{+}\right)\text{Hz},$$(6)

with the oscillator strength of the transition M, ${\chi }_n=157800n^{-2}{T}_{\text{e}}^{-1}$, and the C⁺ column density N(C⁺) (Menzel 1968; Salgado et al. 2017b; Salas et al. 2019).

We use the [¹²C II], [¹³C II] F=2–1 and the Cnα intensities to constrain the gas density and temperature in the PDR in a Bayesian analysis. The errors are estimated using the 0.68 confidence interval, i.e. at 1σ. In the model we use the line width of the [¹³C II] F=2–1 line to compute the [C II] intensities. The model computes the opacity broadening of the optically thick [¹²C II] line.

We make use of the dust spectral energy distributions (SEDs) by Pabst et al. (2019), to compute the dust properties in our positions. They combine the far-IR and submillimeter Herschel/PACS and SPIRE bands to obtain the dust temperature T_d and dust optical depth τ_{d,160 µm} at 160 µm in a modified blackbody fit at 36″ resolution. We extract dust properties with a beam of 45″. To calculate dust properties in BNKL-PDR and TRAP-PDR, where the SPIRE bands are saturated, we use a PACS-only SED (also saturated in a few pixels around BNKL-PDR). Since the FUV-heated dust FIR radiation in Orion is dominated by emission in the PACS bands, especially in warmer regions, values do not vary greatly between the PACS+SPIRE and the PACS-only SEDs. Table 11 summarizes the dust properties toward the observed positions. We note that the derived dust parameters render the J(T_bg) term in Eq. 2 negligible: The equivalent continuum brightness temperature T_c ranges between 10 and 20 K. To estimate the total gas column density in the PDR, we convert the dust optical depth to a hydrogen column density, N_H = N(HI) + 2N(H₂) using (Weingartner & Draine 2001; Pabst et al. 2020), $$N_{\text{H}}\simeq 6\times {10}^{24} {\text{cm}}^{-2}{\tau }_{\text{d},160\mu \text{m}}.$$(7)

Table 12 summarizes the results of the CRRL Bayesian analysis. Figure 13 illustrates how the observed intensities constrain T_e and n_e. The sharp breaks in some of the plots are artifacts from the interpolation of the sparsely sampled original grid, on which the departure coefficients were computed. In general, the observations constrain the physical conditions well. We note that the temperature within a PDR really decreases from the surface, which affects the [C II] and Cnα intensities differently. Hence, the inferred n_e and T_e are an ill-defined average over the [C II] and Cnα emitting layers along the line of sight and within our beam. With the physical conditions in our sample, T_e is predominantly set by the [C II] line intensity, while the Cnα intensity largely defines n_e. Hence, in general, we expect the derived temperature to reflect more the temperature at the PDR surface, where the [C II] emission is strongest, while the density will reflect more the deeper denser layers, where the CRRL emission is strongest. On the other hand, CRRLs will be enhanced compared to the [C II] intensity in clumpy PDRs (increased n_e). How these different conditions affect our type of line intensity analysis is difficult to judge as it strongly depends on the adopted PDR model with its detailed incorporated physical and chemical processes (e.g., constant density versus constant pressure, or homogeneous versus clumpy). More investigation is needed in this regard, which is beyond the scope of this paper.

The results in TRAP-PDR are not to be trusted, since the [¹³C II] intensity is much weaker than the [¹²C II] intensity scaled with the ¹²C/¹³C isotopic ratio (assumed to be 67; earlier studies used a ¹²C/¹³C isotopic ratio of 43 at θ¹ Ori C, Stacey et al. 1993, and references therein). The results in BAR-INSIDE, where the bright PDR lines seem to consist of two close Gaussian components, are most likely weighted toward the high-column density region of the Orion Bar rather than the background molecular cloud southeast of the Orion Bar.

Positions in OMC1 (except for TRAP-PDR) have high electron densities (ranging from 20 to 40 cm⁻³) and moderate temperatures (ranging from 210 to 240 K). Although ORI-P1 is close to OMC1, the electron densities and temperatures of the two components are much less than in OMC1, 2 and 8 cm⁻³ and 100 and 150 K. Possibly we observe ORI-P1-1 edge-on (though moving away from Trapezium), and ORI-P1-2 face-on (same systemic velocity as Trapezium), judging from the C⁺ column density. This is also reflected in the relative strength of the ¹²CO and ¹³CO lines: The blue component appears to be cooler and with higher column density, while the red component appears warmer with lower column density. This two-layer model, with a warm layer behind a cooler layer, for the Dark Lane, in which ORI-P1 lies, was first proposed by Herrmann et al. (1997), derived from their [C II] and [O I] 63 µm and 145 µm fine-structure FIR line observations. The gas density in the warmer layer derived in this work, assuming that all electrons result from carbon ionization and a carbon gas-phase abundance of 𝒜_C ≃ 1.4 × 10⁻⁴ (Sofia et al. 2004), is on the order of 10⁵ cm⁻³, in agreement with the assumption of Herrmann et al. (1997) for the molecular component, while the gas density in the cooler layer is somewhat lower with about 3 × 10⁴ cm⁻³.

In M43 the electron densities are much lower, 2–3 cm⁻³, and the temperatures are 140 K. The non-detection of a CRRL component from the background molecular cloud in ORI-P3 suggests a low density and/or high temperature and/or a low column density, but from the [C II] line properties we derive an optical depth of 1.2 and an excitation temperature of 130 K, values consistent with a moderately dense, moderately illuminated cloud surface with a significant C⁺ column density. On the west border of M43, where the H II region presses against the molecular ridge, Herrmann et al. (1997) derive a density of 3 × 10⁵ cm⁻³ for the molecular gas. Judging from the [O I] 63 µm and 145 µm intensities the density toward the center of M43 must be somewhat lower, but still within the typical range of PDRs (Herrmann et al. 1997). This is in agreement with the gas densities calculated from the electron densities derived in this work, n_H ≃ 2 × 10⁴ cm⁻³.

ORI-P5, associated with the background PDR in the EON, has an electron density of 2 cm⁻³ and a temperature of 180 K. The electron density in ORI-P6 is the lowest of all the sources, 1 cm⁻³, and the region is cool with 130 K. ORI-P6 lies in an edge-on PDR, giving rise to the high C⁺ optical depth and column density.

Our findings indicate that the PDRs in OMC1 have gas densities above 10⁵ cm⁻³ (n_H ≃ n_e/𝒜_C ≃ 1–3 × 10⁵ cm⁻³), while PDRs outside OMC1 have gas densities lower than 10⁵ cm⁻³. The PDRs in OMC1 are known to have high densities, and our estimates are in good agreement with earlier estimates (Tielens & Hollenbach 1985; Wolfire et al. 1990; Stacey et al. 1993; Herrmann et al. 1997; Goicoechea et al. 2015b; Cuadrado et al. 2019). Our estimates for the BNKL-PDR yield a somewhat higher [C II] excitation temperature than the average position toward BN/KL of Morris et al. (2016) with similar optical depth. Our temperature estimates in OMC1 are consistent with earlier findings from [C II] observations (Goicoechea et al. 2015b), but lower than results from CRRL analysis at the edge of the Bar (Cuadrado et al. 2019). The density estimates for M43 are somewhat higher than in Pabst et al. (2022), while the estimate for OMC3 is consistent. The derived gas temperatures for M43 and OMC3 are somewhat higher than estimates in Pabst et al. (2022). Our position ORI-P6 in OMC3 is equivalent to the position EFF5 of Salas et al. (2021), who derive lower temperatures (T_ex[C II] ≃ 63 ± 6 K and T ≃ 55 ± 2 K), but a similar electron density from C110α, C109α, and C102α lines.

In highly irradiated regions the ratio of the C⁺ column density divided by 𝒜_C (i.e., the hydrogen nucleus column density in the gas containing ionized carbon), inferred from the previous analysis, and the total H nucleus column density, derived from the dust optical depth, is low (N(C⁺)/𝒜_C/N_H below 0.25). In M43 and OMC3 the ratio is significantly elevated in comparison (between 0.3 and 0.6). This is in agreement with the higher density inferred for the highly irradiated regions in OMC1 and the EON, where the [C II] emission is more concentrated on the PDR surface (cf. Fig 12).

There are some issues with the adopted approach above. The position of TRAP-PDR gives special riddles. Here, the [¹³C II] line is weaker than expected assuming the standard isotopic ratio of carbon and relative strengths of the hyperfine components. Fractionation, reducing the ¹³C⁺ abundance by locking up carbon-13 in ¹³CO, should be a minor effect, according to our PDR models. We observe a curious velocity shift between the ¹²CO and ¹³CO line in TRAP-PDR – the ¹³CO line is blue-shifted , which we do not observe at all in other positions. Both CO lines are blue-shifted from the [C II] line in this position. We also observe velocity shifts between the [¹²C II] and [¹³C II] lines in some positions. While ORI-P3 has several components with possible self-absorption that might explain the apparent shift, the reason for the shift in ORI-P2 is not clear. It is curious that we do not observe a CRRL component from the background PDR in ORI-P3. It is unlikely that the density is too low or the temperature too high, since we observe both [C II] and CO as fairly strong lines. The [¹³C II] line in ORI-P5 is noisy, but still we observe a extreme shift, while the [¹²C II] profile seems symmetric.

Ideally one would like to use the [¹³C II] line intensity together with the Cnα, β, and γ intensities to constrain the physical properties from the line ratios of these optically thin lines, as Cuadrado et al. (2019) do for the edge of the Orion Bar. However, the predicted intensity ratios for the CRRLs of different Δn are very close to each other and the comparison with the observed intensity ratio therefore very sensitive to the observed intensity ratio. The Cnβ/Cnα ratio only varies from 0.22 to 0.28 in the temperature-density range of (T_e, n_e) = (20– 1000 K, 0.01–100 cm⁻³). Owing to the low S/N ratio, the errors in our observations are actually larger (due to baseline issues, multiple components within the line of sight and the helium RRL being blended with the CRRL) and often observed values lie outside this range. For the positions in OMC1, the errors in the ratios are less than 10%, but the inferred temperatures and densities lie outside the grid. In ORI-P6, where the CRRLs are strong, as well, and there is no ionized gas, the ratios also exhibit high values that indicate higher temperatures and/or densities than the grid range, which we deem unlikely considering the quiescent nature of this region and previous estimates (Kabanovic et al., in prep.).

Fig. 13 Observed CRRL, [12C II], and [13C II] intensities constrain electron temperature Te and electron density ne (orange, C55α; green, [12C II]; blue, [13C II]). The shaded area gives the 1σ error of each observation. Each panel uses the respective model grid with the column density inferred from the Bayesian analysis.

Fig. 14 Thermal pressure in the neutral gas derived in this work. Dashed lines indicates the thermal pressure in the adjacent H II region taken from the literature.

4.3 The pressure equilibrium

Figure 14 shows the thermal pressures computed from the physical conditions in the neutral PDR gas derived in this work, P_th/k_B = n_e/𝒜_CT_e, compared to the thermal pressure in the ionized gas, in the observed positions. The ionized gas in the Orion Nebula was analyzed by O’Dell & Harris (2010), the ionized gas in M43 was analyzed by Simón-Díaz et al. (2011), and physical conditions of the ionized gas in NGC 1977 were estimated by Pabst et al. (2020).

The total pressure in a PDR is the sum of several pressure terms: thermal, turbulent, magnetic, radiation and line pressure (cf. discussion in Sect. 4.1 in Pabst et al. 2020). We should usually expect the PDR to be governed by equipartition of the thermal, turbulent and magnetic pressures, while radiation pressure and line pressure are generally unimportant in Orion. However, as we shall see below (Sect. 4.6), at the spatial scales we probe the turbulent pressure is usually somewhat larger than the thermal pressure in our PDRs. The magnetic pressure in OMC1 is on the order of the thermal pressure (Chuss et al. 2019).

With the physical conditions derived in this work, the two gas phases, neutral and ionized, are in approximate pressure equilibrium in OMC1, albeit values in OMC1 are somewhat lower than earlier findings, that indicate high thermal pressures at the edge of the Orion Bar (p_th/k_B ~ 10⁸K cm⁻³, Cuadrado et al. 2019; Goicoechea et al. 2019, from mm-wave CRRLs and molecular tracers). CRRLs may be enhanced in dense clumps, while the [C II] line arises in the inter-clump medium. Other high-density tracers, such as high-J CO lines, also imply much higher pressures (Joblin et al. 2018). Pellegrini et al. (2009) find a maximum of p_th/k_B ≃ 8 × 10⁷ K cm⁻³ (close to the transition from H to H₂) from modeling various emission lines from the H⁺/H/H₂ transition in the Orion Bar. The thermal pressures in the PDR gas of the Orion Bar, M43 and OMC3 are consistent with those cited in Pabst et al. (2020), derived from analysis of the [C II] line alone. The thermal pressure of M43 is short of the thermal pressure of the ionized gas by a factor of 3, possibly reflecting the equipartition of the thermal, turbulent and magnetic pressure in the PDR. The thermal pressure in the EON, measured in ORI-P5, is higher than the thermal pressure in the ionized gas (n ~ 10² cm⁻³, O’Dell & Harris 2010), suggesting that the ionized gas is insignificant in driving the expansion of the Veil Shell.

Overall, RRLs are useful tracers of the physical conditions and pressures in the respective gas components. Higher angular resolution observations will allow us to resolve finer structures and potentially clumps in the PDRs to fully characterize the properties of the different gas phases.

4.4 Physical conditions in the ionized gas

If ionized gas can freely expand against a background of neutral (PDR) gas, the ionization front moves through the neutral gas at a velocity of ${\upsilon }_{\text{if}}~\frac{c_{\text{s},\text{I}}^2}{2c_{\text{s},\text{II}}}$, where c_S,I is the sound speed of the neutral gas and C_s,II is the sound speed of the ionized gas. Depending on the temperature of the neutral gas, υ_if ~ 0.3–0.5 km s⁻¹. The ionized gas will be accelerated away from the background along the pressure gradient reaching values of υ_i ~ 2c_s,II for free expansion (Spitzer 1968; Kahn 1969).

In OMC1 the ionized gas is enshrouded in a multilayer system of ionized and neutral components (e.g., Abel et al. 2019; O’Dell et al. 2020) in front of the background PDR. The ionized gas traced by our hydrogen RRLs streams away from the background PDR vary rapidly (υ_i ≃ 12–15 km s⁻¹) due to the strong pressure gradient from the interface to the dilute foreground interstellar medium. The expansion velocity is on the order of the sound speed in the ionized gas, υ_i ~ C_s,II. We note that the peak velocity of the HRRL emitting gas does not match any of the components Abel et al. (2019) identify, but O’Dell et al. (2020) comment that the H41α emission of Goicoechea et al. (2015b) should arise from the [O III] emitting layer, which is shifted from the PDR velocity by υ_evap,[OIII] = 9 ± 3 km s⁻¹ (their observations). The main component of our observed HRRLs is a few km s⁻¹ bluer than the [O III] component, and it is unclear why. Especially the position TRAP-PDR samples part of the Orion-South cloud, where several velocity systems have been identified (O’Dell et al. 2021a,b). This is however not clearly reflected in our HRRL spectrum in TRAP-PDR, where the line is still very nearly Gaussian. All HRRL seem to be broadened by other effects than thermal motion of the atoms alone (cf. Sect. 4.6 for a discussion on extra line broadening).

We note that we do not observe [C II] line and CRRL components from the ionized gas toward OMC1 (nor in ORI-P1). In principle, ionized gas emits in the [C II] line and CRRLs, but the radiation from the close Trapezium cluster is energetic enough to doubly ionize carbon and the emission from singly ionized carbon will be weak. The EON, to which our position ORI-P5 corresponds, contains largely hot plasma and little ionized gas, and we do not observe the weak [C II] line and CRRL components from the ionized gas either. In M43 and NGC 1977, in contrast, the illuminating stars are less massive and therefore their radiation is less energetic, and a significant fraction of carbon is expected to be singly ionized. We observe those [C II] lines from the ionized gas as broad components (Δυ_FWHM ~ 5– 10 km s⁻¹) in spectra toward the center of M43 and NGC 1977 (Pabst et al. 2017, 2020).

In M43 the ionized gas is contained within an expanding shell. The HRRL in ORI-P3, the center of M43, is blueshifted from the CRRL by 1.5 km s⁻¹. However, the CRRL in ORI-P3 is blue-shifted by 2 km s⁻¹ from the background PDR traced by [C II] and CO. Possibly this CRRL also stems from the ionized gas, where carbon is singly ionized. The ionized gas is streaming away from the background PDR with 3.5 km s⁻¹. The HRRL is significantly broader than the CRRL due to the different atomic weight (20.5 km s⁻¹ compared to 7.4 km s⁻¹), but the CRRL is broader than in the other positions. In the rim of the M43 shell, in positions ORI-P2 and ORI-P4, the ionized gas is blue-shifted from the [C II] and CO-traced PDR gas by 3 and 5 km s⁻¹ , respectively, and the CRRLs peak velocities and line widths are close to the peak velocities and line widths of the [C II] line.

Plotting the peak intensities of the individual hydrogen transitions in the observed positions reveals that the ionized gas is in local thermal equilibrium (LTE). This is characterized by the intensities being proportional to frequency. The nearly constant ratios between the different orders of RRLs suggest that the gas has similar electron temperatures in all positions. The ionized gas in the Huygens Region, however, is denser (n_e ≃ 5 × 10³ cm⁻³ close to θ¹ Ori C, O’Dell & Harris 2010) than the ionized gas in M43 (n_e ≃ 500 cm⁻³, Simón-Díaz et al. 2011).

Figure 15 shows the integrated intensity ratios of the helium RRLs over the hydrogen RRLs. Using I_Henα/I_Hnα = n(He⁺)/n(H⁺) (e.g., Roshi et al. 2017), where we assume the H and He RRLs coincide spatially, we derive a helium ioniza-tion fraction. The Henα/Hnα ratio is high in OMC1 (0.070 ± 0.002, 0.0824 ± 0.0009, 0.0826 ± 0.0006, 0.084 ± 0.001 for BAR-INSIDE, BNKL-PDR, TRAP-PDR, and EAST-PDR, respectively) and significantly lower in the center of M43 (0.014 ± 0.003). ORI-P1 and ORI-P5 take intermediate values (0.051 ± 0.003 and 0.049 ± 0.006). In ORI-P2 we do not observe the He lines, meaning a He⁺/H⁺ ratio close to zero.

In OMC1 the ionizing radiation is energetic enough to ionize helium (IP = 24.6 eV) and to doubly ionize carbon (IP = 24.4 eV). In M43 there still is some ionized helium, but only about a quarter of the fraction in OMC1. ORI-P1 and ORI-P5 are still subject to radiation from the Trapezium stars and have a fair amount of ionized helium. The He RRL in ORI-P2, associated with M43, is too weak to be detected, but the noise level is consistent with the helium fraction in ORI-P3.

Fig. 15 Intensity ratio of Henα over Hnα RRLs.

Fig. 16 Stacked X IInα RRLs toward BNKL-PDR (upper panel) and TRAP-PDR (lower panel) in the velocity frame of the C+ RRL, together with the carbon and helium RRLs (scaled to the intensity of the X IInα RRL) in the velocity frame of the helium RRL for reference.

4.5 Detection of heavy-ion RRLs toward BNKL-PDR and TRAP-PDR

In the H II region of OMC1, the Huygens Region, carbon is doubly ionized and we do not observe the ionized gas in [C II] or CRRLs. Toward the positions BNKL-PDR and TRAP-PDR, the helium RRLs are strongest and, given the sensitivity of our observations, we can potentially detect α RRLs from C⁺ and/or O⁺, that arise in regions where helium is ionized. The first detection of RRLs of C⁺ and/or O⁺ (dubbed X II) toward OMC1 has been reported recently by Liu et al. (2023), by matching the theoretical frequencies with highly sensitive multiband radio observations. Given the reported line width (Δυ_FWHM > 10 km s⁻¹) and the velocity offset between the same C⁺ and O⁺ RRL transitions of ~3.5 km s⁻¹ , it is not possible to distinguish between the two ions. Using the frequencies Liu et al. (2023) give for RRLs of C⁺, we stack the X II nα lines for n = 81, 82, 83, 85, 86, 88, 89, 90 toward BNKL-PDR and n = 81, 82, 83, 85, 86, 88, 89 toward TRAP-PDR, where we apply a standing-wave baseline correction and resample to a velocity resolution of 2.0 km s⁻¹. We exclude spectra that contain molecular lines and other, earlier identified, RRLs, or where a standing-wave baseline correction is not possible. We obtain a rms noise of 1.8 mK toward BNKL-PDR and 2.4 mK toward TRAP-PDR. Figure 16 shows the resulting stacked X IInα line together with the carbon and helium RRLs (scaled to the intensity of the X IInα line) in the velocity frame of the helium RRL for reference.

From a Gaussian fit of the X IInα RRL, we obtain a peak temperature of T_mb ≃ 9.4 ± 0.8 mK for BNKL-PDR and T_mb ≃ 6.7 ± 0.8 mK for TRAP-PDR. The peak velocity and the line width depend critically on the baseline correction, especially for TRAP-PDR, where the line is somewhat weaker and baselines are worse. We tentatively obtain, assuming X=C, a peak velocity of υ_LSR ≃ −5 ± 1 km s⁻¹ and a line width of Δυ_FWHM ≃ 21 ± 2 km s⁻¹ for BNKL-PDR. For TRAP-PDR we tentatively report υ_LSR ≃ −1 ± 2 km s⁻¹ and Δυ_FWHM ≃ 27 ± 4km s⁻¹. If X=O, the peak velocity is redder by 3.5 km s⁻¹ .In the integrated intensity we thus detect the X IInα RRL in BNKL-PDR at a 17σ level, and in TRAP-PDR at 12σ significance.

The X IInα RRLs closely resemble those of helium. Assuming X=O would shift the peak velocities further from the H and He RRL peak velocities. The ionization potential of C⁺ is very close to that of helium, while O⁺ has a much higher ionization potential at 35 eV, hence it is likely that C⁺ and O⁺ RRLs arise in different ionized layers, while C⁺ and He RRLs may stem from the same layer. Among the many layers Abel et al. (2019); O’Dell et al. (2020) identify in the Huygens Region, they identify a high-ionization component in optical [O III] lines at υ_LSR ≃ 0 km s⁻¹, attributed to the main ionization front (MIF). This may indicate that O⁺ RRLs arise in the same layer in the MIF, that cannot be distinguished from a layer of lower-ionization gas, where the C⁺ and He RRLs arise. Our beam toward TRAP-PDR includes parts of their NE (northeast) region and the SW (southwest) dark cloud, while our BNKL-PDR region lies somewhat to the north of the region of their observations. O’Dell et al. (2020) also identify a high-ionization component in optical [O III] lines at υ_LSR ≃ 9 km s⁻¹, that they attribute to shocked gas moving into the MIF. We do not see that component in our observations; possibly this component is too localized to contribute significantly in our beam.

While the peak velocities and line widths we obtain are in agreement with line properties reported by Liu et al. (2023), our line intensities seem to be somewhat smaller (a factor of ~3), while our velocity resolution is a factor of 1.5 better, but our beam is 1.5 times larger. From the intensity ratio with the Hnα intensities, assuming that the intensity ratio equals the density ratio if they stem from the same region, we obtain n(X⁺⁺)/n(H⁺) ≃ 3.8 ± 0.4 × 10⁻³ in BNKL-PDR and n(X⁺⁺)/n(H⁺) ≃ 2.6 ± 0.4 × 10⁻³ in TRAP-PDR, which is 8 and 5 times larger, respectively, than the expected abundance ratio in Orion of C⁺ and O⁺ together, ~5 × 10⁻⁴ (Sofia et al. 2004). As Liu et al. (2023) observe, dielectronic recombination may be important in enhancing the intensity of heavy-ion RRLs. Detailed photoionization modeling of the excitation conditions would shed light on the region where RRLs from C⁺ and O⁺ arise, but this is beyond the scope of this paper.

4.6 Line broadening

In this section we discuss the observed line widths in both the [C II] and the hydrogen, helium and carbon RRLs. While RRLs can be broadened by a number of physical processes, this does not apply to the [C II] line. Also, some of these processes are expected to only apply to higher quantum numbers (lower frequencies) than those of our observations. Collisional and radiation broadening lead to the line having a Voigt profile in general, a convolution of a Gaussian profile, produced by turbulence, with a Lorentzian profile (Salgado et al. 2017b). However, our observations cover only quantum numbers n < 100, hence we expect the profiles to be dominated by Doppler broadening. In fact, the lines can be fitted very well by Gaussian profiles.

A discussion of the extra line broadening (in addition to thermal broadening) in optical lines can be found in O’Dell et al. (2023). Our observed radio and far-infrared lines are usually broader than the thermal line width, while preserving a Gaussian profile. The extra line broadening observed in optically thin lines can either stem from turbulence or Alfvén waves, according to O’Dell et al. (2023). We follow along the lines of that discussion. While our rather large beam width can be problematic in comparing results to those from higher spatial resolution optical spectra, we note that our high velocity resolution does make it technically possible to distinguish separate velocity components averaged in our beam (as long as the separation is large enough compared to the line width), as we clearly observe in BAR-INSIDE (averaging the Orion Bar and the molecular cloud southeast of the Orion Bar), where the [C II], Cnα, and ¹³CO lines deviate from a single-component Gaussian profile. Especially in OMC1, where velocity shifts occur on rather small spatial scales in certain regions (O’Dell et al. 2021a), this is critical to check. However, we emphasize that the lines in our sample are surprisingly Gaussian, in contradiction with this given. It seems that our beam avoids these regions, e.g. toward TRAP-PDR we do not cover much of the Orion-South cloud, where spatial and spectral complexity is a known property (also O’Dell et al. 2021b). We expect more quiescent, more extended regions to be less affected by small-scale velocity shifts. However, we can estimate the effect on the line widths of large-scale (turbulent) motions averaged within our beam with an average beam size of 45″ (0.09 pc at the distance of Orion) for the PDR gas.

To estimate this effect, we compare the [¹³C II] and ¹³CO lines at 16″ resolution with those at 45″, that we use in the discussion, and the C51α line at 36″ with the C59α line at 54″. In the [¹³C II] line this line-broadening effect seems to be on the order of at most 20%, with the exception of BAR-INSIDE, where the two fits give a difference of 36%. Most accurately the effect of the beam size can be estimated for the ¹³CO line, as the S/N ratio is usually high. Here we see that the line width is smaller for the 16″ beam by at most 20%, but in the brightest sources it is smaller by less than 10%. For the Cnα lines the effect is on the order of 10 to 20%. For 10% narrower lines, the turbulent pressures we compute in Sect. 4.6.4 thus would be lower by 19%, while for 20% narrower lines the turbulent pressures would be lower by 36%. Where the turbulent Mach number (cf. Sect. 4.6.4) is close to or slightly higher than 1, turbulence would then be subsonic rather than sonic or slightly supersonic. Running ahead of the discussion, however, our main conclusions would remain unchanged by this modification, but we have to keep this cautionary note in mind when drawing conclusions.

In the following, we need the equation for the total line width. The thermal line width (FWHM) for an ideal gas is given by $$\Delta {\upsilon }_{\text{th }}^2=\frac{8\ln 2k_{\text{B}}T}{m},$$(8)

with T being the gas temperature and m the mass of the atom or molecule. The total line width, in the presence of turbulence, is then given by $$\Delta {\upsilon }^2=\Delta {\upsilon }_{\text{th }}^2+\Delta {\upsilon }_{\text{turb }}^2\text{. }$$(9)

The Alfvén line width, that can contribute to the extra line broadening, is given by (Henney et al. 2005) $$\Delta {\upsilon }_{\text{A}}^2=\frac{B^2}{4\pi \rho },$$(10)

where B is the magnetic field strength and ρ is the ion mass density of the gas. In addition, the Alfvénic Mach number 𝓜_A is defined as the ratio of the line width (standard deviation) over the Alfvén velocity. The spectral resolution of the data contributes to the observed line width, but in our case it is negligible: With a smallest observed line width of 1.3 km s⁻¹ in the ¹³CO line, the velocity resolution of these observations (0.4 km s⁻¹; 0.3 km s⁻¹ for the [¹³C II] lines, and 0.36 km s⁻¹ for the CRRLs) contributes at most 5% to the measured line width; in the vast majority of the ¹³CO, [¹³C II], and Cnα observations, this instrumental broadening is at most 1%.

4.6.1 Line broadening in the neutral gas

For carbon in the PDR at a temperature of T = 300 K, the thermal line width is Δυ_FWHM ≈ 1.1 km s⁻¹. The extra line-broadening component in the CRRLs toward OMC1 is therefore between 2.5 km s⁻¹ and 4.8 km s⁻¹. If this was due to Alfvén waves, the magnetic field can have a maximum value of 23 µG, assuming a density of 3 × 10⁵ cm⁻³. This is lower than the magnetic field observed in the Veil (B_los ≃ −50–75 µG, Troland et al. 2016) and much weaker than the magnetic field observed in the Orion Bar (B ≃ 300 µG, Chuss et al. 2019).

If we attributed the CRRL in ORI-P3 to the neutral gas, as we do in all other positions, the line width in ORI-P3 would indicate extreme turbulence in the CRRL emitting gas, which is however not reflected in the [C II] lines. The line width is consistent, though, with the thermal line width of gas at the temperature of ionized gas with a small contribution from turbulence. The CRRL emission in ORI-P3 would then stem from an ionized layer closer to the molecular background compared to the H and He RRLs. It is extremely puzzling, though, that we do not observe a component that stems from the PDR. In ORI-P2 and ORI-P4 the extra line broadening components are 3.1 and 5.1 km s⁻¹, respectively, which results in a maximum magnetic field strength estimate of 5 µG if Alfvénic.

The CRRL in ORI-P6 is likely also turbulence dominated. Assuming T = 100 K, the thermal line width is Δυ_FWHM ≈ 0.6 km s⁻¹ and the turbulent line width would be 1.8 km s⁻¹. The extra line broadening component is smaller than in OMC1, OMC3 being a more quiescent region. If the extra line broadening was dominated by Alfvén waves, the magnetic field in ORI-P6 could have a maximum value of 1.5 µG, assuming a density of 10⁴ cm⁻³.

Assuming these recent estimates of the magnetic field strength in OMC1 to be accurate and applicable to our positions, we conclude that the extra line broadening is highly sub-Alfvénic in the neutral gas (𝓜_A << 1).

4.6.2 Line broadening in the ionized gas

For hydrogen at a temperature of T = 10⁴ K, the thermal line width is Δυ_FWHM ≈ 21.4 km s⁻¹. For helium in thermal equilibrium with the hydrogen gas, the thermal line width is Δυ_FWHM ≈ 10.7 km s⁻¹.

If we assume that hydrogen and helium are in thermal equilibrium and that the turbulence is equal in the hydrogen and helium gas, we can solve the two equations for the total line widths of the hydrogen and helium lines, $\Delta {\upsilon }^2=\Delta {\upsilon }_{\text{th }}^2+\Delta {\upsilon }_{\text{turb }}^2$ (Eq. (9)), for the gas temperature and the turbulent velocity. This way, we derive a temperature of the ionized gas and a turbulent velocity of 1.1 × 10⁴ K and 14.0 km s⁻¹ toward the Orion Bar, 0.99 × 10⁴ K and 12.5 km s⁻¹ toward BNKL-PDR, 1.0 × 10⁴ K and 14.2 km s⁻¹ toward TRAP-PDR, 0.77 × 10⁴K and 21.0 km s⁻¹ toward EAST-PDR. In ORI-P3, however, the observed line widths are in good agreement with sole thermal broadening in a gas with T ≈ 0.9 × 10⁴ K.

These assumptions may not be exactly valid. We observe the centroid velocities of the helium RRLs to be more blueshifted with respect to the hydrogen RRLs in the positions toward OMC1 (including ORI-P1) and ORI-P5 in the EON. The biggest blueshift is observed in ORI-P5 with −3.4 km s⁻¹. Also in ORI-P3 we observe a significant blueshift of −2.1 km s⁻¹. It can be expected from dynamic modeling of the Huygens region that helium RRLs arise from gas closer to the stars (since its ion-ization potential is higher than that of hydrogen, “the High Ionization Zone”) and hence be more blueshifted (Henney 2003; O’Dell et al. 2023). On the other hand, according to recent observations with the JWST, the He Strömgren sphere is observed to be slightly smaller than the H Strömgren sphere (5 × 10⁻³ pc, Peeters et al. 2024). As the gas density will decrease as the gas expands into the Huygens Region, and the recombination lines will be strongly weighted with density, this demonstrates that the ionized gas is rapidly accelerated once it leaves the ionization front to about the sound speed over 5 × 10⁻³ pc – as indicated by the velocity of the HRRL – and that further acceleration only adds about 3.4km s⁻¹ to this, as indicated by the HeRRL. Further dynamical modeling taking the full geometry of M42 HII region and the stellar bubble into account may be very insightful in the dynamical feedback of this region.

In the ionized gas of OMC1, the magnetic field would need to be rather strong to produce the extra line broadening, about 400 µG. In the light of this, the assumption that the extra line broadening is only due to Alfén waves becomes unreasonable, as the frozen-in magnetic field should rather increase in strength from the ionized gas to the denser PDR. Rather, turbulence in the ionized gas is slightly super-Alfvénic (𝓜_A ≳ 1): the Alfvén speed in the ionized gas is at most about 15 km s⁻¹, while the extra line broadening we observe in OMC1 is between 12.5 and 21.0 km s⁻¹.

4.6.3 Line broandening in the molecular gas

The ¹³CO line in BNKL-PDR is characterized by very broad wings, that are due to the strong outflows from this region. The [C II] line toward this position also shows deviation from Gaus-sianity, with broad wings, albeit less than the ¹³CO lines. In contrast, the RRLs seem to be even less affected by this line-broadening effect, although we note that this is hard to judge as the carbon RRL is blended with the helium RRL and a blue weaker component at υ_LSR ≃ 1 km s⁻¹ , corresponding with a [C II] component. Likely the ¹³CO and [C II] lines consist of broad wings from shocked gas and a narrow component from the PDR, that emits in the carbon RRL.

The temperature of the molecular gas can be estimated from the peak temperature of the optically thick ¹²CO line by $$T_{\text{P}}\left({}^{12}\text{CO}\right)=J\left(T_{\text{ex}}\right)-J\left(T_{\text{bg}}\right),$$(11)

with the equivalent brightness temperature of a black body at T $J(T)=\left(hv/k_{\text{B}}\right)/\left(e^{hv/k_{\text{B}}T}-1\right)$, the excitation temperature of the ¹²CO line Tex and the cosmic background temperature T_bg = 2.7 K. The excitation temperature of the CO line is a good lower limit on the kinetic temperature in the CO gas. For the ¹²CO(J = 2–1) line, hν/k_B = 11.07 K. If hν/k_BT << 1, T_P ≃ T_ex ≤ T_kin. In this manner, we obtain estimates for the gas temperature in the CO gas of 100–130 K in OMC1, and 40 to 60 K in positions outside OMC1. We use these temperatures to determine the thermal line widths of the ¹³CO lines.

In molecular gas with a temperature of T = 100 K, the thermal line width of ¹³CO is Δυ_FWHM ≈ 0.2 km s⁻¹. Hence, the line widths we observe here are heavily dominated by the extra line broadening. Generally the line widths are somewhat smaller than those of the corresponding [C II] lines, between 1 and 2km s⁻¹ (the exception being BNKL-PDR, of course). Assuming an electron abundance of x_e ~ 10⁻⁵ and a density of 3 × 10⁵ cm⁻³, the magnetic field in the CO gas can be at most 3 µG to produce a corresponding Alfvén velocity. If the density increases from the surface layer of the PDR to the molecular gas, the magnetic field can be higher. Berné et al. (2014) conclude that the magnetic field is important in stabilizing the molecular cloud against gravitational collapse, but cite higher values of the magnetic field strength (50 to 250 µG Abel et al. 2004; Brogan et al. 2005). Hence, in the molecular gas, turbulence is highly sub-Alfvénic (𝓜_A << 1).

The ¹³CO lines in OMC1 are not consistently broader than those in regions more distance from the Trapezium stars, which suggests that the cloud is intrinsically turbulent and not much energy is transmitted deep into the cloud. The only exceptions are the observed lines in ORI-P6 that are noticeably narrower than in the other positions. OMC3 is a much more quiescent region than OMC1, M43 and the EON, with still significant star formation (Megeath et al. 2012, 2016, see Fig. 11 in Pabst et al. 2021). This suggests that at least some turbulence is caused by the presence of massive stars, but our sample is too small to draw definite conclusions. Berné et al. (2014) conclude that only 0.01% of the radiative energy injected by stars (E_★ ≃ 5 × 10⁵¹ erg) is transferred into the molecular cloud as turbulent kinetic energy, E_k ≃ 3.6 × 10⁴⁷ erg. A large amount, 50%, of the injected wind energy (E_w ≃ 5 × 10⁴⁸ erg) is transferred into the large-scale expansion of the Veil Shell, E_kin ≃ 2.5 × 10⁴⁸ erg (Pabst et al. 2020).

4.6.4 Connecting the dots

Comparing the turbulent line widths in Table 13, we conclude that the turbulent line widths of each gas phase are about equal in BAR-INSIDE, EAST-PDR-2, ORI-P1-1, and ORI-P6. In the other positions/components, the line widths in the different phases vary a great deal with respect to one another. In BNKL-PDR this can be due to the fact the molecular outflows and associated shocks affects mostly the molecular lines and lines close to the molecular phase. In TRAP-PDR, the [¹³C II] line has somewhat lower turbulent line width than the Cnα and ¹³CO line. In ORI-P1-2, the turbulent line widths increase from the neutral to the molecular phase, while in ORI-P1-1, ORI-P5 and ORI-P6, the turbulent line widths decrease going from the neutral to the molecular phase. ORI-P2 and ORI-P4 show unordered behavior, with the turbulent line width in the Cnα line largest, which may be due to the faintness of the CRRLs. In ORI-P3(-2), the Cnα line may stem not from the neutral gas, but from the ionized gas, as discussed earlier.

We can also calculate the turbulent Mach numbers 𝓜_turb of the different layers. The one-dimensional turbulent Mach number 𝓜_turb,z is the ratio of the velocity dispersion σ_turb over the sound speed. The total Mach number is given by 𝓜_turb = $\sqrt{3}{ℳ}_{\text{turb },\text{z}}$ if the turbulence is isotropic (Orkisz et al. 2017). For interstellar gas the sound speed is given by $c_{\text{s}}\simeq \sqrt{\gamma p_{\text{th}}/\rho }=\sqrt{\gamma k_{\text{B}}T/m}$, where p_th is the thermal pressure, ρ is the mass density of the gas, γ is the adiabatic index, which is 5/3 for a monatomic ideal gas and 7/5 for a diatomic ideal gas, T is the temperature, and m is the atomic or molecular weight of the gas.

In ionized gas with a temperature of T = 10⁴ K, the sound speed is about c_s ≃ 14 km s⁻¹. In neutral atomic hydrogen at a temperature of T = 300K, the sound speed is approximately c_s ≃ 2.0 km s⁻¹, while at T = 100K it is c_s ≃ 1.2 km s⁻¹. If hydrogen is molecular with T = 100 K, c_s ≃ 0.8 km s⁻¹. If half of the hydrogen nuclei are in atomic gas and half are locked up in molecular hydrogen, the sound speed can be computed to be effectively c_s ≃ 1.0 km s⁻¹. In molecular hydrogen gas with T = 30 K, the sound speed is c_s ≃ 0.4 km s⁻¹.

We assume that the hydrogen nuclei fractions in atomic and molecular form, respectively, are each 0.5 in the neutral gas (i.e., ƒ(HI) = ƒ(H₂) = 0.5). In BAR-INSIDE, the turbulence in the neutral gas is slightly subsonic. In BNKL-PDR, the turbulence is subsonic in [¹³C II], but supersonic in Cnα, while in TRAP-PDR and EAST-PDR it is slightly supersonic. In the positions outside OMC1 turbulence in the neutral gas is supersonic. In ORI-P6 the turbulent Mach number is lowest. In the molecular gas Mach numbers are significantly higher than in the neutral gas, all positions are supersonic (in BNKL-PDR the high Mach number is caused by the molecular outflows). Plotting the turbulent line widths and Mach numbers, respectively, of the three PDR tracers ([¹³C II], Cnα, ¹³CO) against each other does not reveal an overall correlation of the turbulent properties between the different tracers/PDR layers, nor do Pearson correlation coefficients. In general, we observe that the Mach numbers increase from the neutral to the molecular gas, albeit by very different factors. In the ionized gas in OMC1, turbulence is slightly subsonic with Mach numbers between 0.7 and 1.0.

Table 14 summarizes the turbulent pressures and the ratio of the turbulent pressure over the thermal pressure in the observed PDRs, derived from the Bayesian analysis of our line intensities. For the ¹³CO lines we assume the same density as for the [C II]/CRRL layer and the excitation temperature of the ¹²CO line as a lower limit for the gas temperature. The turbulent pressure, $p_{\text{turb }}=\rho {\sigma }_{\text{turb }}^2$, in the [C II]/CRRL layer is usually on the order of the thermal pressure, albeit somewhat larger within a factor of 3. Only in ORI-P2 and ORI-P5 the turbulent pressure in [¹³C II] is 6–7 times larger than the thermal pressure, and in ORI-P4 the turbulent pressure in Cnα is 13 times larger than the thermal pressure. In the molecular gas, the turbulent pressure dominates over the thermal pressure with factors between 1.9 (ORI-P6) and 28 (ORI-P1-2) and is usually within a factor of 2 of the turbulent pressure in the CRRL layer (except in ORI-P4, where it is somewhat lower).

Turning to the ionized gas, assuming a density of 5 × 10³ cm⁻³ for the ionized gas (O’Dell & Harris 2010), the turbulent pressure in the ionized gas is on the order of 5 × 10⁷ K cm⁻³, thus on the order of the thermal pressure in the ionized gas and the turbulent pressure in the PDRs in OMC1. In M43, the turbulent pressure in the ionized gas is negligible, the Hnα and Henα lines are in first order consistent with thermal broadening alone (cf. Sect. 4.6.2). In NGC 1977, Pabst et al. (2020) report a gas density of the ionized gas of n ≃ 40 cm⁻³ and a width of the [C II] line of Δυ_FWHM = 11.6 km s⁻¹. With these givens, assuming a gas temperature of T ~ 10⁴ K, the turbulent pressure can be computed to be p_turb ≃ 2 × 10⁵ K cm⁻³, which is a factor of 5 lower than the thermal pressure of the ionized gas and the turbulent pressure in the PDR of OMC3.

Estimates from [C II] observations in Pabst et al. (2020) suggest that the thermal, turbulent and magnetic pressures in the Veil Shell and the PDR of M43 are in approximate equiparti-tion, while the thermal pressure is an order of magnitude larger than the turbulent and magnetic pressure in the Orion Bar and the PDR of NGC 1977. The radiation pressure in the Veil Shell and M43 is about equal to the other pressure terms, while it is an order of magnitude lower in the Orion Bar and the PDR of NGC 1977. The Veil Shell and M43 are very dynamic regions that expand rather rapidly (with 13 km s⁻¹ and 6 km s⁻¹, respectively), while NGC 1977 expands more slowly (with 1.5 km s⁻¹). Possibly the still rapid expansion of the Veil Shell and M43 injects turbulence into the interstellar medium, while the shell of NGC 1977 interacts more quiescently in a Spitzer-type expansion and has aged sufficiently to have let turbulence be dissipated away. The time scale for turbulent dissipation may be estimated by τ ~ l/c_s (Kolmogorov 1941; Elmegreen & Scalo 2004); if the characteristic length scale is l ~ 1 pc and the sound speed c_s ~ 1 km s⁻¹, the dissipation timescale is τ ~ 10⁶yr, which is on the order of the age of NGC 1977, 0.4 × 10⁶ yr. The Orion Bar is heavily irradiated by the nearby Trapezium stars and the high turbulence may be a result of flows caused by the strong radiation field and pressure gradients in this region.

This same behavior is reflected in our samples, where the ratios of the turbulent pressure over the thermal pressure in positions in OMC1 is smaller than in other positions. We note that we derive rather low temperatures of the PDRs in OMC1 compared to earlier studies: despite exhibiting high turbulent pressure, the thermal pressure in the PDRs in OMC1 may be still dominating the total pressure in the PDRs. The turbulent pressure we derive in M43 (ORI-P2 and ORI-P4) is an order of magnitude higher than in Pabst et al. (2020) and larger than the thermal pressure. Those two positions lie on the edges of M43, though, while Pabst et al. (2020) look at the central background PDR. Possibly, turbulence is increased where the gas can freely stream away. We derive a ten times higher turbulent pressure in OMC3 than Pabst et al. (2020) derive for the PDR associated with the large shell of NGC 1977, rendering the thermal pressure and the turbulent pressure in OMC3 in approximate equipartition.

In conclusion, at the spatial scales our observations probe we observe that turbulence in the Orion A molecular cloud is highly supersonic in the molecular gas, but sub-Alfvénic, varying between slightly sub- and supersonic and sub- and super-Alfvénic in the neutral gas, but slightly subsonic and highly super-Alfvénic in the ionized gas. The turbulent pressure dominates the pressure balance in the neutral and molecular gas, but plays a minor role compared to the thermal pressure in the ionized gas.

The data presented here may help settle the flow of turbulent energy in regions of massive star formation. Since the gas phases border each other, it is likely that turbulence is transmitted from one gas phase into the other. In our case, the extra line broadening is larger in the ionized gas by an order of magnitude, hence it is unlikely that it stems from the turbulence of the molecular gas. We cannot exclude the possibility of Alfvén waves being responsible for the extra line broadening in the ionized gas, but turbulent motions from gas flows are likely to contribute in this setting. How magnetic fields and turbulence are passed on between the phases of the ISM is of yet unknown. Another obvious contribution to the line width is the dynamic structure of the observed regions, blister H II regions (OMC1) and thermally expanding H II regions (M43). However, lines are nearly Gaussian, not giving away much hints on the dynamic origin of the line width.

An alternative way to study the flow of turbulent energy in star-forming regions, is to statistically analyze the full data cubes using structure functions. However, the scale at which turbulence is important is comparable to the spatial resolution of our [C II] and CO data cubes (for the ionized gas Arthur et al. (2016) gives a turbulent scale of 22″; O’Dell (2001) gives an average turbulent scale of 18″, depending on emission line), while the larger-scale morphology of the region becomes important at our resolutions and the radial velocity distribution in each spectrum is complicated. For this approach higher-resolution data cubes of lines from neutral atomic PDR gas are needed, which may be available with the next generation of far-infrared and radio telescopes.

It is generally acknowledged that molecular clouds are strongly turbulent, stabilizing the cloud against collapse, usually interpreted as the result of a turbulent cascade from cloud scales to filament scales (e.g., Hennebelle & Falgarone 2012; Federrath 2016; Padoan et al. 2016). Gas flows caused by local energy injection in the interstellar medium are believed to be turbulent, while also gravitational collapse can induce turbulence (Heigl et al. 2020). On galactic scales, gravity may be the dominant source of turbulence in molecular clouds (Krumholz & Burkhart 2016). If the line broadening is caused by turbulence, there may be two reasons for turbulence in the molecular gas: either that the ¹³CO lines are broadened by turbulence that is transmitted from the ionized gas through the neutral gas to the molecular gas, or else that the molecular clouds are intrinsically turbulent which is transmitted to the neutral gas and exacerbated by photoevapora-tion on the surface. In the case of the Orion B molecular cloud, the compression of the molecular gas caused by H II regions leads to the formation of filaments and thus increased star formation (Gaudel et al. 2023), while turbulence in the molecular gas is highly supersonic (Orkisz et al. 2017).

5 Conclusion

We have analyzed multiple hydrogen, helium and carbon mm-wave RRLs observed with the Yebes 40 m telescope toward 10 representative positions in the Orion Nebula complex, including M43 and OMC3. Comparing with the [C II] 158 µm line as a PDR surface tracer and the ¹³CO (2–1) line as a tracer of the molecular gas deeper in the PDR, we observe slight velocity shifts between the main components of these two tracers and the carbon RRLs, that cannot always be attributed to geometry, but may indicate intrinsic shear. The hydrogen RRL emitting ionized gas in OMC1 is streaming away from the PDR surface in the background at a velocity of 12–15 km s⁻¹, while in M43 the streaming velocity is 3.5 km s⁻¹.

Combining the carbon RRL intensities with the intensities of the [¹²C II] main component and its [¹³C II] F=2–1 satellite, we infer electron temperatures and electron densities. This approach has its shortcomings, as the three lines originate from slightly different PDR layers, which is reflected in different observed line parameters (peak velocity and line width). We obtain robust estimates, that allow us to compare the thermal pressure in the observed PDRs with other pressure terms (turbulent, magnetic) and the pressure in the ionized gas. Besides in the position southeast of the Orion Bar, where the expansion of the Veil is driven by the contained hot plasma, the PDR and the ionized gas are in approximate pressure equilibrium.

In OMC1, the RRLs show that helium is ionized, while in M43 only a quarter of the amount of helium is ionized. This reflects the difference in radiation field, of a O7V star in OMC1 and a B0.5V star in M43. In OMC3, we do not observe ionized gas. Toward two positions in OMC1, BNKL-PDR and TRAP-PDR, we detect very faint lines that can be attributed to RRLs of C⁺ and/or O⁺ (X IInα RRLs), stemming from the ionized layers in OMC1.

The observed lines, whether H, He, or C RRLs, are significantly broader than thermal broadening would allow. The only exception is M43, where we attribute the observed carbon RRL to the ionized gas instead of the background PDR due to its peak velocity and line width. Here, the line widths are consistent with thermal broadening alone. In the other positions we explore the possibility that the extra line broadening is due to Alfvén waves. We conclude that the extra line broadening is super-Alfvénic and subsonic in the ionized gas, Alfvénic and varying between slightly sub- and supersonic in the neutral gas at the H/H₂ transition, but sub-Alfvénic and supersonic in the molecular gas. Hence, we conclude that the extra line broadening is mainly due to interstellar turbulence that is non-Alfvénic, while the precise origin of such turbulence remains unknown.

By observing RRLs of different elements and in different ionization states, future telescopes such as the SKA will be able to resolve much finer structures, which will result in a refinement of the estimates of the physical conditions in exquisite detail, in Orion as well in more distant massive star forming regions.

Acknowledgements

We thank the anonymous referee for valuable comments that helped improve the manuscript. We thank C.R. O’Dell for valuable discussions on the properties of OMC1. This work is based in part on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA was jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Insti-tut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. CHMP and JRG thank the Spanish MCINN for funding support under grant PID2019-106110GB-100. CHMP thanks the Dutch NWO for funding support through a Rubicon grant no. 13257.

Appendix A Observed positions on SOFIA/upGREAT [C II] and Spitzer/IRAC4 8 µm maps

Figures A.1 to A.3 show the average beam (FWHM = 45″) toward the observed positions on maps of the SOFIA/upGREAT [C II] intensity and the Spitzer/IRAC4 8 µm intensity converted to [C II] intensity using the relations provided by Pabst et al. (2022), zoomed in to discern the structure that is averaged in each beam, in their original resolution (SOFIA: FWHM = 16″, Spitzer: FWHM = 1.9″). As the IRAC4 8 µm image reveals, each beam contains finer structure, that is also lost in the upGREAT [C II] image. Future telescopes (such as the SKA and far-infrared heterodyne receivers on large telescopes) can provide the spatial resolution necessary to resolve the structure in velocity-resolved observations, which will lead to much refined estimates of the physical conditions in the observed regions.

Fig. A.1 Zoom into OMC1. Left: SOFIA/upGREAT [C II] intensity at 16″. Right: Spitzer/IRAC4 8 µm intensity converted to [C II] intensity at 1.9″. Observed positions in the region of OMC1 (BAR-INSIDE, BNKL-PDR, TRAP-PDR, EAST-PDR, ORI-P1, and ORI-P5) are indicated as purple circles, that represent the average 45″ beam size of our Yebes 40 m telescope observations.

Fig. A.2 Zoom into M43. Left: SOFIA/upGREAT [C II] intensity at 16″. Right: Spitzer/IRAC4 8 µm intensity converted to [C II] intensity at 1.9″. Observed positions in the region of M43 (ORI-P2, ORI-P3, and ORI-P4) are indicated as purple circles, that represent the average 45″ beam size of our Yebes 40 m telescope observations.

Fig. A.3 Zoom into OMC3. Left: SOFIA/upGREAT [C II] intensity at 16″. Right: Spitzer/IRAC4 8 µm intensity converted to [C II] intensity at 1.9″. Observed positions in the region of OMC3 (ORI-P6) are indicated as purple circles, that represent the average 45″ beam size of our Yebes 40 m telescope observations.

Appendix B Line fits to Cnβ, γ, δ and ϵ RRLs

Tables B.1, B.2, B.3, and B.4 give the fit parameters of the Cnβ, γ, δ and ϵ RRLs toward the ten targeted positions.

Appendix C Line fits to Hnβ, γ, δ and ϵ RRLs

Tables C.1, C.2, C.3, and C.4 give the fit parameters of the Hnβ γ, δ and ϵ RRLs toward the ten targeted positions.

Appendix D Line fits to Henβ, γ, δ and ϵ RRLs

Tables D.1, D.2, D.3, and D.4 give the fit parameters of the Henβ, γ, δ and ϵ RRLs toward the ten targeted positions.

References

All Tables

All Figures

	Fig. 1 Observed positions (except positions in OMC1) as white circles with a FWHM of 45″ on [12C II], 12CO(2–1), and 13CO(2–1) maps at 45″. The black rectangle outlines the region of OMC1, shown as a zoom-in in Fig. 2. The red circle outlines M43. The yellow and orange stars mark the positions of θ1 Ori C and θ2 Ori A, respectively, the pink star marks NU Ori, the purple star marks 42 Orionis.
In the text

	Fig. 2 Zoom into OMC1. Observed positions in OMC1 as turquois circles with a FWHM of 45″ on [12C II], [13C II], 12CO(2–1), and 13CO(2–1) maps at 45″. The yellow and orange stars mark the positions of θ1 Ori C and θ2 Ori A, respectively.
In the text

Fig. 3 Stacked Cnα, Cnβ, Cnγ, Cnδ, and Cnє radio-recombination lines with Henα, Henβ, Henγ, Henδ, and Henє radio-recombination lines toward the ten targeted positions. In the position BAR-INSIDE the component 7.5 km s−1 blue of the main component may be a radio-recombination line from sulfur, while detection of the sulfur RRL in other positions is hampered by the carbon RRL component corresponding to the Veil (cf. Sect. 4.1). Spectra in in the upper top are offset from the C/Henє spectrum in steps of 0.025 K, sand pectra in the lower two rows are offset in steps of 0.01 K.

In the text

	Fig. 4 Stacked Hnα, Hnβ, Hnγ, Hnδ, and Hnє radio-recombination lines toward the ten targeted positions. Spectra in in the top row are offset from the Hnє spectrum in steps of 0.15 K, and spectra in the lower two rows are offset in steps of 0.01 K.
In the text

	Fig. 5 Stacked Cnα RRL, [12C II], and 13CO(2–1) lines toward the 10 targeted positions. The 13CO(2–1) line is scaled by the factor indicated in light gray in each panel.
In the text

	Fig. 6 [12C II] and [13C II] F=2–1 lines toward the 10 targeted positions.
In the text

	Fig. 7 12CO(2–1) and 13CO(2–1) lines toward the 10 targeted positions.
In the text

	Fig. 8 Stacked Cnα (blended with Henα) with respect to the Cnα rest frequency and stacked Hnα RRLs with respect to the Hnα rest frequency toward the ten targeted positions.
In the text

	Fig. 9 Peak velocities in the observed positions of the [C II], Cnα, and 13CO lines obtained from Gaussian fits. The background shading indicates the velocity-ordering class, defined in the main text, the position and component belongs to: red is class 1, green is class 2, blue is class 3, yellow is class 4, and purple is class 5.
In the text

	Fig. 10 Intensity ratios of Cnα over Cnβ and Cnγ, respectively (left axis) and [12C II] over Cnα (right axis).
In the text

	Fig. 11 Intensity ratios of Hnα over Hnβ, Hnγ, and Hnδ, respectively.
In the text

	Fig. 12 PDR model outputs from the Meudon PDR code for an incident radiation field of G0 = 100 and three different densities (n = 103, 104, 105 cm−3). The upper panel shows the normalized intensity, assuming LTE for the CRRL emission, the lower panel shows the abundances of the relevant species.
In the text

	Fig. 13 Observed CRRL, [12C II], and [13C II] intensities constrain electron temperature Te and electron density ne (orange, C55α; green, [12C II]; blue, [13C II]). The shaded area gives the 1σ error of each observation. Each panel uses the respective model grid with the column density inferred from the Bayesian analysis.
In the text

	Fig. 14 Thermal pressure in the neutral gas derived in this work. Dashed lines indicates the thermal pressure in the adjacent H II region taken from the literature.
In the text

	Fig. 15 Intensity ratio of Henα over Hnα RRLs.
In the text

	Fig. 16 Stacked X IInα RRLs toward BNKL-PDR (upper panel) and TRAP-PDR (lower panel) in the velocity frame of the C+ RRL, together with the carbon and helium RRLs (scaled to the intensity of the X IInα RRL) in the velocity frame of the helium RRL for reference.
In the text

	Fig. A.1 Zoom into OMC1. Left: SOFIA/upGREAT [C II] intensity at 16″. Right: Spitzer/IRAC4 8 µm intensity converted to [C II] intensity at 1.9″. Observed positions in the region of OMC1 (BAR-INSIDE, BNKL-PDR, TRAP-PDR, EAST-PDR, ORI-P1, and ORI-P5) are indicated as purple circles, that represent the average 45″ beam size of our Yebes 40 m telescope observations.
In the text

	Fig. A.2 Zoom into M43. Left: SOFIA/upGREAT [C II] intensity at 16″. Right: Spitzer/IRAC4 8 µm intensity converted to [C II] intensity at 1.9″. Observed positions in the region of M43 (ORI-P2, ORI-P3, and ORI-P4) are indicated as purple circles, that represent the average 45″ beam size of our Yebes 40 m telescope observations.
In the text

	Fig. A.3 Zoom into OMC3. Left: SOFIA/upGREAT [C II] intensity at 16″. Right: Spitzer/IRAC4 8 µm intensity converted to [C II] intensity at 1.9″. Observed positions in the region of OMC3 (ORI-P6) are indicated as purple circles, that represent the average 45″ beam size of our Yebes 40 m telescope observations.
In the text