© The Authors 2026

1 Introduction

Explosive events triggered by mergers of protostars generate transitory infrared (IR) line luminosity bursts and inject substantial mechanical energy (∼10⁴⁷⁻⁴⁸ erg) into the natal molecular cloud, serving as an early source of mechanical feedback in clustered star-forming regions (e.g., Bally & Zinnecker 2005). The resulting bow-shock tips, or “bullets,” propagate through the cloud, generating shocks that cool via molecular and atomic line emission, depending on their velocity, magnetic field, gas density, and ambient ultraviolet (UV) radiation. Stellar mergers must occur relatively frequently in clustered star-forming regions (perhaps once per century in our Galaxy; e.g., Zapata et al. 2020, 2023) and they are therefore likely to be common in luminous galaxies. The Orion BN/KL outflow (e.g., Kwan & Scoville 1976; Beckwith et al. 1978; Snell et al. 1984) is the nearest example of an “explosive outflow” (e.g., Gómez et al. 2005; Zapata et al. 2009; Goddi et al. 2011; Bally et al. 2011, 2020). Its distinctive kinematics and diverse gas cooling lines (atomic, H₂, CO, H₂O, and OH; e.g., Werner et al. 1984; Rosenthal et al. 2000; González-Alfonso et al. 2002; Zapata et al. 2009; Peng et al. 2012; Goicoechea et al. 2015a) indicate the presence of shocks, in which the mechanical energy of the explosion is converted into gas heating and compression as the outflow impacts the swept-up gas and the ambient molecular cloud. However, the nature of these shock(s), ‘jump” (J-type), “continuous” (C-type) or mixed, as well as the role of UV radiation, whether shock-generated or external, remain uncertain. On the one hand, UV-irradiated shocks are best traced by velocity-resolved observations of the [O I] 63 and 145 µm and [C II] 158 µm cooling lines, but only a few sources have been mapped in the wings of all three lines. In particular, the scarcity of [O I] 145 µm data limits constraints on the physical conditions of the [O I]-emitting gas across velocities. On the other hand, only recently have shock models (e.g., Hollenbach & McKee 1979, 1989; Draine 1980; Kaufman & Neufeld 1996) included realistic descriptions of UV-irradiated shocks (Lesaffre et al. 2013; Melnick & Kaufman 2015; Godard et al. 2019, 2024; Lehmann et al. 2020, 2022; Kristensen et al. 2023).

Here, we present the highest so far angular resolution [O I] 63 and 145 µm maps of the BN/KL outflow (6″ and 13″, respectively) carried out at sub-km s⁻¹ resolution and complemented by sensitive line maps of [C II] 158 µm, a tracer of UV radiation and stellar feedback (e.g., Tielens & Hollenbach 1985; Stacey et al. 1993; Goicoechea et al. 2015b; Pabst et al. 2019), and OH, a tracer of the outflow kinematics (e.g., Melnick et al. 1987, 1990; Goicoechea et al. 2006a, 2015a). We obtained these observations with the Stratospheric Observatory for Infrared Astronomy (SOFIA; Young et al. 2012) using GREAT¹, providing a legacy dataset for studying the properties of shocked gas.

The paper is organized as follows. In Sect. 2, we provide details on the outflow. In Sect. 3, we describe the airborne observations. In Sect. 4, we present the main observational findings, while in Sect. 5, we detail the physical conditions as a function of velocity in the outflow. Finally, Sect. 6 discusses the nature of the [O I]-emitting shocks based on comparisons with state-of-the-art shock models and places our results in the context of FIR line emission from other massive star-forming regions.

2 The Orion BN/KL explosive outflow

Embedded in the heart of the Orion Molecular Cloud-1 core (OMC-1), just behind the Orion Nebula (M42), the Becklin–Neugebauer/Kleinmann–Low (BN/KL) region is the closest (∼414 pc; Menten et al. 2007) high-mass star-forming region (Genzel & Stutzki 1989; O’Dell 2001; Bally 2008). In addition to the two well-known runaway stars BN and source I, which are the most massive objects in this region (e.g., Lonsdale et al. 1982; Scoville et al. 1983; Menten & Reid 1995; Bally et al. 2020), the field hosts the first identified molecular outflow, exhibiting high-velocity CO emission, with line wings extending to over ±100 km s⁻¹, together with a wide-angle, shock-excited H₂ outflow (e.g., Kwan & Scoville 1976; Beckwith et al. 1978; Snell et al. 1984). Observations of this outflow have long provided key benchmarks for the development of shock models (Draine & Roberge 1982; Neufeld & Dalgarno 1989; Hollenbach & McKee 1989). Today we know that this is an explosive outflow produced by a dynamical decay event (e.g., a stellar merger), which triggered the acceleration of BN and source I, about 500 yr ago (e.g., Gómez et al. 2005; Zapata et al. 2009; Bally et al. 2011; Goddi et al. 2011; Bally et al. 2020). The center of the explosion is located between the current positions of these stars.

ALMA CO observations reveal an approximately spherically symmetric “Hubble–Lemaître flow,” that is, radial velocities scaling with the projected distance from the center. The flow is composed of over a hundred narrow bow-shock wakes, or “fingers,” within ∼1^′ (∼0.1 pc) of the center (Bally et al. 2017). Due to the isotropic distribution of the flow, the red- and blueshifted fingers appear to overlap when projected onto the plane of the sky, a characteristic of explosive outflows. The bullets show shock-excited H₂ v = 1–0 S (1) (2.12 µm) and forbidden [Fe II] 1.64 µm emission at their tips, together with H₂ and lower velocity CO emission in their wakes (e.g., Zapata et al. 2009; Kristensen et al. 2007, 2008; Nissen et al. 2012; Bally et al. 2015; Youngblood et al. 2018). The H₂ emission shows a two-lobe northwest–southeast (NW–SE) orientation (e.g., Allen & Burton 1993; McCaughrean & Mac Low 1997), likely because the other part of the outflow – or its IR emission – is partially obscured by the dense dust ridge along OMC-1, which extends toward the northeast (NE). The outflow contains ∼8 M_⊙ of molecular gas (Snell et al. 1984). About half of this mass is at expansion velocities below 20 km s⁻¹. The other half belongs to the high-velocity outflow, which exhibits some of the brightest H₂ emissions in the sky. Peak 1, located ∼30″ northwest of BN, is the strongest H₂-emitting region of the outflow (see Fig. 1).

The total H₂ line luminosity across the outflow is 120 ± 60 L_⊙ (Rosenthal et al. 2000). Observations of H₂ lines toward Peak 1 reveal excitation temperatures rising from ∼600 K in the v = 0 rotational transitions, to ∼2500 K in the vibrationally excited transitions (Rosenthal et al. 2000; Youngblood et al. 2018). Observations suggest a more excited component at ∼5000 K, likely caused by H₂ formation pumping after shock dissociation and reformation (Geballe et al. 2017). Peak 1 also exhibits excited CO and H₂O emission associated with warm, ∼500 K, and hot, ∼2500 K molecular gas (González-Alfonso et al. 2002; Goicoechea et al. 2015a). In addition, low angular and spectral resolution FIR observations revealed bright [O I] fine-structure line emission (e.g., Werner et al. 1984; Herrmann et al. 1997; Lerate et al. 2006; Goicoechea et al. 2015a) while Herschel/HIFI enabled the first velocity-resolved [C II]158 µm maps of the region (Goicoechea et al. 2015b; Morris et al. 2016).

The BN/KL outflow also shows intense OH emission, with $L_{\text{OH}}/L_{{\text{H}}_2\text{O}}\simeq 0.4$, suggesting ongoing H₂O photodissociation (Goicoechea et al. 2015a). The low-lying OH rotational lines fall at FIR wavelengths, where the dust continuum is very intense. These lines require very high densities to be collisionally excited (n_cr ≳ 10⁸⁻⁹ cm⁻³) and are typically optically thick. As gas densities are typically lower, their excitation temperatures (T_ex) often fall below the continuum temperature, producing “P-Cygni” profiles in molecular outflows (first detected toward BN/KL; Betz & Boreiko 1989; Melnick et al. 1990; Goicoechea et al. 2006a).

Fig. 1 BN/KL outflow and the Trapezium cluster observed with JWST/NIRCam H2 F212N (McCaughrean & Pearson 2023). The dashed square marks the field of view mapped with SOFIA/GREAT. Cyan contours show the [O I] 63 µm redshifted line-wing emission over vLSR = 25–30 km s−1, from 15 to 90 K km s−1 in steps of 15 K km s−1.

3 SOFIA airborne observations and data reduction

We mapped multiple FIR lines with GREAT during eight SOFIA flights between January 2014 and February 2021. Apart from the flights taking place in in spring 2021 (observing cycle 8), which we operated out of Cologne-Bonn (CGN, Germany) during the pandemic, we carried out all flights from Palmdale, California. We conducted all observations during GREAT Consortium time (project IDs: 83_0004 (cycle 1), 83_0428 (cycle 4), and 83_0630 (cycles 8–9). The flights accumulated 5 hours of observing time.

In 2014, we carried out the first observations of the OH ²Π_3/2 J=5/2-3/2 line (2514.3 GHz ≃ 119.2 µm) with the single-pixel detector of GREAT (Heyminck et al. 2012). In 2021, we extended the map using the corresponding channel of the multi-color 4GREAT receiver (Duran et al. 2021). We observed the [O I] ³P₁–³P₂ line (4744.8 GHz ≃ 63 µm) with the 7-pixel high-frequency array (HFA) of upGREAT (Risacher et al. 2018). We simultaneously observed the [C II] ²P_3/2–²P_1/2 line (1900.5 GHz ≃ 158 µm) and the [O I] ³P₀–³P₁ line (2060.1 GHz ≃ 145 µm) with the two polarization-split upGREAT low-frequency sub-arrays (LFA): the 7-pixel V-polarization array tuned to the [O I] line and the H-polarization array to the [C II] line. This configuration processed 21 signals in parallel. We employed Fast Fourier Transform spectrometers, updated from Klein et al. (2012), which provided 32k channels across the ∼0.5–4 GHz intermediate-frequency bands. Table A.1 summarizes the instrument configurations.

The instrument and observatory performance (flight altitude, precipitable water vapor) determined the actual observing mode and source coverage (map sizes), while previous Herschel spectroscopic maps of Goicoechea et al. (2015a) guided them. We reference all observations to the nominal source position RA = 05^h35^m14.3^s, Dec = −05°22′33.7″ (J2000) in the hot core region. We performed the upGREAT observations of both [O I] lines and [C II] in fast total power on-the-fly slewing mode (0.4 sec integration time per dump), using 3″ sampling in both RA and Dec We observed the source twice, scanning once in Dec and once in RA, with the fine sampling driven by the 6.3″ beam of the HFA. The map size (144″×126″) refers to the on-sky coverage of the central array pixels, while the hexagonal array’s outer pixels provide a larger coverage. We selected a “clean” on-sky reference position at ∆RA, ∆Dec=(−1700″, +900″). To achieve the most uniform sampling, we tilted the array axis by 19.1° relative to the scanning direction (Risacher et al. 2016a). We observed the OH ²Π_1/2 (1834.7 GHz ≃ 163.4 µm) line during the LFA commissioning in 2016, using chopped (1 Hz) on-the-fly mode with 6″ spatial sampling. The reference OH line frequencies (Table A.1) correspond to the line-strength-weighted frequency of each hyperfine-structure triplet, which are spectrally unresolved. We observed the OH 163.4 µm line of the Λ-doublet in the upper sideband, thereby avoiding to blend with the 163.1 µm line of the doublet in our dual sideband receiver.

We applied a large chop throw of 240″ toward negative RA offsets. Observing the OH ²Π_3/2 119.2 µm line proved most challenging. At this frequency, the shortage of local oscillator power needed to pump the HEB mixers forced us to couple the signal via a Martin–Puplett interferometer (Heyminck et al. 2012); this, in turn, limited the receiving bandwidth of the sky signal to ∼1.6 GHz. Thus, the expected velocity coverage of the outflow barely fits into the effective bandpass (∼190 km s⁻¹). This required pointed (raster) observations, in double-beam chopped mode, for optimal system stability (wobbler throw 300″, at −20° counter-clockwise against positive RA). In 2014, an exploratory 5×5 raster on a 10″ grid was observed (with the single-pixel M-channel of GREAT), showing the emission to be compact; in 2021 with 4GREAT we added a 5×5 raster on a 6″ grid.

The GREAT consortium derived the half-power main beam sizes and efficiencies from planet observations. Telescope operators established the pointing on nearby optical reference stars with an accuracy of 1–2″. At the beginning of each flight series, we aligned the instrument’s optical axis with these imagers by observing planets. Although the co-alignment between the central pixels is better than 1–2″ (and the data header accounts for it), we generally tracked the highest frequency channel in a given configuration – typically the central HFA pixel in most flights. We amplitude-calibrated the raw data with the KOSMA kalibrate software (Guan et al. 2012). We corrected the data for atmospheric extinction and calibrated them in T_mb². When SOFIA flew between 12.2–13.4 km in altitude, the atmospheric transmission remained smooth near the target velocities.

We processed the calibrated data with GILDAS³. We inspected the data to remove noisy spectra and other artifacts (e.g., spikes from onboard interference). We smoothed the spectra to a channel resolution of 0.15 km s⁻¹ (0.3 km s⁻¹ for the OH 119 µm line) and removed first-order baselines. The spectra were then gridded into a data cube through a convolution with a Gaussian kernel of a full width at half maximum (FWHM) of ∼1/3 of the telescope beamwidth at the different frequencies. The 1σ rms noise in the maps are ∼3.5 K ([O I] 63 µm), ∼2.5 K ([O I] 145 µm), and ∼2.0 K ([C II] 158 µm), and ∼4.3 K (OH 163 µm) per velocity channel (as provided by the NOISE routine). To determine the noise over a given velocity bin, we used ${\sigma }_{\text{bin}}={\sigma }_{\text{ch}}\delta {\upsilon }_{\text{ch}}\sqrt{N_{\text{ch}}}$, where σ_ch is the rms per channel, δv_ch the channel width, and N_ch the number of channels in that bin. For specific applications (e.g., line intensity ratio maps), we convolved the data cubes with Gaussian kernels to a common angular resolution, corresponding to FWHMs of 13″ or 15″.

4 Results

Figure 1 shows the core of OMC-1 observed with JWST/NIRCam F212N (McCaughrean & Pearson 2023), dominated by H₂ v = 1–0 S (1) line emission at 2.12 µm. The dashed square marks the field of view mapped with SOFIA/GREAT, encompassing the wide-angle BN/KL outflow, including the bright Peak 1 and Peak 2 regions. The cyan contours in Fig. 1 show the [O I] 63 µm redshifted line-wing emission at ∼20 km s⁻¹ relative to the cloud systemic velocity of ≃8–9 km s⁻¹ (hereafter v_LSR,0; Bally et al. 1987; Berné et al. 2014). This [O I] emission peaks close to BN (Sect. 4.3) and follows the H₂ emission from the outflow. The field also includes the massive stars of the Trapezium cluster, in the foreground, located at ∼0.2 pc in front of the cloud. Their far-UV (FUV; 6 < E < 13.6 eV) radiation, dominated by θ¹ Ori C star, illuminates all surfaces of OMC-1.

Fig. 2 Total line intensity maps (in erg s−1 cm−2 sr−1 = “cgs”) integrated over the complete line profile. The beam size is indicated in the bottom-right corner of each panel except for H2, where it is too small to display (0.1″). The cyan triangle shows a position near the Trapezium, where the [O I] 63µm line intensity peaks (Table D.1). The H2 image refers to the JWST/NIRCam F212N image (McCaughrean & Pearson 2023). The CO J = 10–9 map was obtained with Herschel/HIFI at 20″ resolution (Goicoechea et al. 2019).

Fig. 3 Total line intensity ratio maps (derived from line intensities in erg s−1 cm−2 sr−1). (a) [O I] 63/145 at a common angular resolution of 13″, (b) [O I] 63/[C II] 158 at 15″, and (c) [O I] 145/ [C II]158 at 15″. The dashed box shows the (100″×80″) area used to extract both the line luminosities in the outflow region (Table 2) and the line intensities as a function of velocity (see Figs. 7 and 8). Contours show the FIR 79 µm-continuum obtained with Herschel/PACS, from 2.5 to 10 (103 Jy), centered at the position of the hot core (Goicoechea et al. 2015a).

4.1 Integrated line intensity maps

Figure 2 shows [O I] 63 145 µm, [C II] 158 µm, and OH 163 µm total line intensity maps (integrated over their complete line profile) in units of erg s⁻¹ cm⁻² sr⁻¹, which are used to compare the emitted line luminosity and line cooling budget. In these maps we mark the positions of key stars (BN and source I) and environments (Peak 1, Peak 2, and the Trapezium region), which we study throughout the manuscript. The [O I] 63 µm emission is the brightest FIR atomic fine-structure line in the region and roughly follows the H₂ emission along the wide-angle outflow. Still, the line intensity peaks southwest of the Trapezium stars, hereafter the “Trapezium position,” at (+20″, −50″) from the map center (see details in Sect. 4.2). This region is part of the large-scale, face-on photodissociation region (PDR), with n_H ≳ 10⁵ cm⁻³ (e.g., Herrmann et al. 1997; Goicoechea et al. 2015b, 2019), at the interfaces between OMC-1 and the foreground H II region, which is photoionized by θ¹ Ori C. This PDR component is located closer to the observer than the BN/KL outflow and produces bright narrow line emission – often called the “spike” component – at nearly all positions in the field.

The [O I] 145 µm line is strikingly bright and as bright as the [C II] 158 µm line, displaying a similar spatial distribution as the [O I] 63 µm emission. As observed before (e.g., Goicoechea et al. 2015b; Morris et al. 2016), the [C II] 158 µm integrated intensity map does not follow the H₂ emission. Instead, [C II] 158 µm shows strong emission close to the Trapezium, along with lower line intensity emission from the FUV-irradiated surface of OMC-1. The OH ²Π_1/2 J=3/2-1/2 line emission at 163 µm is confined to the outflow, peaking toward the inner hot core region, with only very faint emission from the quiescent cloud and PDR. The spatial distribution of the OH 163 µm emission resembles that of the high-J CO emission (e.g., Peng et al. 2012; Goicoechea et al. 2015a, 2019).

Figure 3 shows the atomic fine-structure line intensity ratio maps (integrated over the entire line profile). Typically, these intensity ratios provide diagnostics of the gas density, temperature, and FUV radiation field strength in the warm gas associated with PDRs and shocks. Although these maps do not reveal the velocity structure of the gas, they display a different spatial distribution compared to the absolute line intensity maps (Fig. 2). The [O I] 63/[C II] 158 intensity ratio map shows a clear maximum (≃20) toward Peak 1 (Fig. 3b). The [O I] 145/[C II] 158 map exhibits a similar spatial distribution, with a high ratio (≳1) toward Peak 1 (Fig. 3c). Given the much higher critical densities for collisional excitation of the [O I] 63 and 145 µm lines compared to that of the [C II] 158 µm line (see Table A.1), these peaks trace the highest density gas (assuming similar FUV illumination). In general, the [O I] 63/[C II] 158 µm line intensity ratios toward Peak 1 are higher than those observed toward protostellar outflows (e.g., Liseau et al. 2006) and PDRs (e.g., Bernard-Salas et al. 2012). Figure 3 (contours) also shows the FIR continuum map at 79 µm (Goicoechea et al. 2015b), tracing the extended warm dust emission from the hot core region. The continuum emission is fainter toward the Peak 1 and Peak 2. The [O I] 63/145 ratio shows its lowest values (≃10) in a region roughly orthogonal to the wide-angle outflow (Fig. 3a), following the quiescent cloud, the “extended ridge”. Radiative transfer models (Sect. 5.2) indicate that the low intensity ratios are consistent with the [O I] 63 µm line being optically thick, and with the gas having lower density than in the post-shock gas.

Fig. 4 Velocity-resolved spectra at representative positions, all from maps convolved to 15″ (except CO 10–9, with a beam of 20″). Left: Complete spectra, with offsets in arcseconds given in parenthesis; continuum levels of [C II] 158 µm and [O I] 145 µm shifted for clarity. The OH 163 µm emission lines (yellow) are scaled by a factor of three or five. Right: zoom on faint line wings and foreground line features.

4.2 Velocity-resolved far-IR line profiles

4.2.1 Line wings

The GREAT maps enable us to resolve the velocity structure of the [O I] 63 and 145 µm emission and to compare it with other FIR lines: [C II] 158 µm, OH, and high-J CO. Figure 4 presents such a comparison for four representative positions: Peak 1, the hot core region, Peak 2, and the Trapezium position. The first three positions lie on the outflow region and show line wing emission. As is apparent from these spectra, the OH and CO lines exhibit similarly broad line profiles, with ∆v_FWHM ≳ 30 km s⁻¹. These lines reveal prominent wings, with full widths at zero intensity (FWZI) exceeding 150 km s⁻¹ (see Fig. 5 and Kwan & Scoville 1976; Zapata et al. 2009; Peng et al. 2012, for additional velocity-resolved CO studies). These line profiles differ markedly from those of the FIR atomic fine-structure lines (Fig. 4), which consist of a bright narrow component, ∆v_FWHM ≃ 4 km s⁻¹, together with moderately broad wings. This suggests that the shocks producing the atomic fine-structure emission differ in nature or location from those driving the highest velocity CO and OH emission.

The FWZI of the [O I] 63 and 145 µm lines toward Peak 1 is ∼85 km s⁻¹ and ∼62 km s⁻¹, respectively. The [O I] FWZIs exceed that of the [C II] 158 µm line (≲ 50 km s⁻¹). This implies that the highest velocity component either lacks FUV illumination or that C⁺ is rapidly converted into other species in the post-shock gas (see also Morris et al. 2016). The [O I] 145 µm line is optically thin and is not affected by foreground absorptions (see Sect. 4.2.3). Toward Peak 1, the [O I] 145 µm line profile can be fitted with two Gaussians of similar intensity but different line FWHM: a narrow one (∆v_FWHM ≃ 3.6 ± 0.1 km s⁻¹ centered at v_LSR ≃ +8.2 ± 0.1 km s⁻¹), consistent with gas in the quiscent face-on PDR, and a broad one (∆v_FWHM ≃ 19.6 ± 0.3 km s⁻¹ at v_LSR ≃ +7.4 ± 0.1 km s⁻¹), implying supersonic velocity dispersions associated with shocked gas.

In Appendix B, we show the velocity channel maps of the observed FIR lines. As in Fig. 1, the redshifted [O I] 63 and 145 µm emission follows the H₂ and CO emission along the two lobes of the wide-angle outflow, whereas the red-shifted [C II] 158 µm emission is weak toward Peak 1 and Peak 2 (Fig. B.1). The [O I] 63 and 145 µm emission from the blueshifted wing peaks around BN and Peak 1. At v_LSR,0, the atomic lines show extended emission in the SW–NE direction, orthogonal to the wide-angle outflow, mostly tracing the illuminated rim of OMC-1’s ridge.

Assuming that most of the narrow-line emission arises from the face-on PDR at the rims of OMC-1, as suggested by its spatial distribution, we estimate that ∼65%, ∼50%, and ∼15% of the observed [O I] 63, [O I] 145, and [C II] 158 µm emissions, respectively, originate from shocked gas. Nonetheless, some of the narrow line emission may originate from low-velocity outflows or from shocks oblique to the line of sight, so the above percentages should be considered lower limits.

4.2.2 [O I] 63 µm narrow-line “spike” emission

For bright, optically thick [O I] 63 µm emission, the continuum-subtracted peak main-beam temperature, T_P = J(T_ex) − J(T_dust), where J(T) is the Planck-corrected radiation temperature, yields the excitation temperature, T_ex, of the line, arising here from the extended PDR in front of the outflow: $$T_{\text{ex}}([\text{OI}]\text{\hspace{0.17em}63\hspace{0.17em}μm})=\frac{227.7}{\ln \left(1+\frac{227.7}{T_{\text{P}}([\text{OI}]\text{\hspace{0.17em}63\hspace{0.17em}μm})+J(T_{\text{dust}})}\right)}\text{\hspace{0.17em}K},$$(1)

where ∆E/k_B = (E_u − E_l)/k_B = 227.7 K, and J(T_dust) = T_c is the continuum brightness temperature at 63 µm. In most observed positions, T_P([OI] 63 µm) ≃ 100–120 K, and J(T_dust) ≃ 5–100 K (from the 63 µm continuum levels of Werner et al. 1984), where the highest value refers to the hot core region. These T_P([OI] 63 µm) values translate to T_ex ≳ 200 K, which is a lower limit to the gas temperature (e.g., Fig. C.2).

Far from the outflow, the atomic fine-structure lines are narrow and lack high-velocity wings. This is exemplified by the spectra toward the Trapezium position (Fig. 4), where these lines show Gaussian profiles and line FWHM of about 5 km s⁻¹ (Table D.1). Still, this particular position – resembling a cavity structure – shows remarkably strong [O I] 63 µm emission, with a peak brightness temperature of T_peak = 330 K (at the native 6″ resolution). This implies T_k ≥ T_ex ≃ 440 K (with T_c ≃ 5 K, and assuming a beam filling of one), consistent with warm-to-hot gas in a dense PDR facing the Trapezium. This excited PDR is located near the region of maximum CH⁺ J = 1–0 emission in OMC-1 (position #3 of Goicoechea et al. 2019), a signature of strong FUV irradiation and elevated temperatures. JWST images reveal this and similar cavities illuminated by FUV radiation from the Trapezium; e.g., Fig. 4 of McCaughrean & Pearson (2023) and Fig. 4 of Habart et al. (2024). These “crenellations” (O’Dell et al. 2015) at OMC-1’s illuminated surface were likely carved by protostellar outflows (see also Kavak et al. 2022).

4.2.3 Foreground [O I] 63 µm line absorption

The [O I] 63 µm line is a ground-state transition. The strong dust continuum at 63 µm from the central BN/KL region (Werner et al. 1976; Cernicharo et al. 2006; Goicoechea et al. 2015b), together with the presence of low excitation gas along the line of sight, produces absorption features in the blueshifted [O I] 63 µm line wing, where T_ex < T_c. In this region, the [O I] 63 µm line profile shows a narrow, a few km s⁻¹ FWHM, absorption component at v_LSR ≃ +6 km s⁻¹, which is typically associated with gas in the hot core itself (e.g., Blake et al. 1987; Tercero et al. 2010). In addition, it shows blueshifted absorptions, at v_LSR ≃ −1 km s⁻¹ (toward Peak 1) and at v_LSR ≃ +2 km s⁻¹ (toward the hot core and Peak 2). These blueshifted absorptions coincide with [C II] 158 µm emission peaks (left panels of Fig. 4) known to arise from foreground material in Orion’s Veil (Goicoechea et al. 2015b; Pabst et al. 2019), including the “Northern Dark Lane” and the “Dark Bay”, which are apparent in optical images (O’Dell 2001). Because of the lower density in these components (n_H ≈ 10³ cm⁻³; Pabst et al. 2020), the fainter 158 µm continuum, and the much lower critical density of the [C II] 158 µm transition, leading to T_ex > T_c, these foreground gas layers appear as [C II] 158 µm emission components. Away from the hot core region, the [O I] 63 µm line profiles show no absorption components.

4.2.4 OH P-Cygni profiles from the high-velocity outflow

Figure 4 (left panel, yellow histograms) display the OH ²Π_1/2 J=3/2-1/2 line profile (one of the Λ-doublets of the excited OH 163 µm transition). These observations have much higher angular and spectral resolution than previous observations with the Kuiper Airborne observatory (KAO), the Infrared Space Telescope (ISO), and Herschel (Melnick et al. 1987, 1990; Goicoechea et al. 2006a, 2015a). The OH 163 µm line profile toward Peak 1 shows a very broad profile, with a FWZI of ≳ 150 km s⁻¹ and ∆v_FWHM ≃ 40 km s⁻¹, thus comparable to the CO lines (e.g., Kwan & Scoville 1976; Snell et al. 1984; Zapata et al. 2009; Peng et al. 2012; Bally et al. 2017). As opposed to the atomic fine-structure lines, the OH 163 µm emission does not show a narrow line component. Indeed, the line is very faint toward the Trapezium position (Fig. 4), implying that the broad-line emission originates in shocked gas and regions that are subject to strong radiative pumping, but not from the quiescent cloud gas. This rotationally excited OH line arises from the first excited state (at E_u/k = 270 K) of the 2Π_1/2 ladder. Given its high excitation requirements (n_cr ≃ 10⁹ cm⁻³), the upper level can only be populated by collisions in very dense, warm gas, and/or by absorption of FIR dust photons in the ²Π_1/2–²Π_3/2 cross-ladder transition at ∼53 µm (e.g., Goicoechea & Cernicharo 2002; Goicoechea et al. 2006a). Owing to the different dust continuum levels, the latter process contributes to the hot core region, while the former likely dominates in Peak 1 and Peak 2, with possible lateral (i.e., not along the line of sight) FIR illumination from the hot core. The broad OH 163 µm line profiles trace the outflow at all velocities, and the line intensity tends to peak toward the hot core (see also Fig. B.1), where the effect of FIR pumping is strongest.

Figure 5 compares the [O I] 63 µm, CO J = 10–9, and OH ²Π_3/2 J=5/2-3/2 (one of the Λ-doublets of the OH 119 µm ground-state transition) line profiles toward the hot core position (roughly the FIR dust continuum peak). The [O I] 63 µm spectrum shows the three absorption components at v_LSR ≃ +6, +2, and −1 km s⁻¹. The first component is consistent with gas in the hot core, while the other two originate from foreground gas associated with Orion’s Veil. Betz & Boreiko (1989) provided the first sub-km s⁻¹ resolution observations of the OH 119 µm rotational line with the KAO. These observations revealed a self-absorption profile, characterized by absorption on the blueshifted side and emission on the redshifted side. The SOFIA telescope delivers a factor of three better angular resolution (11.3″ at 119 µm), which provides less angular dilution of the FIR continuum source. The resulting OH 119 µm line shows a more prominent P-Cygni profile (the line shows absorption and emission of similar strength), with a FWZI of ≃160 km s⁻¹, much like CO lines. Compared to Betz & Boreiko (1989), SOFIA observations show that the entire blueshifted wing is absorbed from v_LSR ≃ −80 to +8.5 km s⁻¹, while the redshifted emission extends up to v_LSR ≃ +80 km s⁻¹, with a velocity coverage similar to that of the CO J = 10–9 line. Clearly, OH is present in the high-velocity outflowing gas, with terminal (maximum) velocities of about 80 km s⁻¹ – the fastest gas producing absorption projected along the line of sight. A closer examination of the OH 119 µm P-Cygni profile reveals that the main absorption dip occurs at v_LSR ≃ +6 km s⁻¹, which corresponds to the hot core absorption (for a zoomed-in spectrum, see Fig. B.3). We leave a detailed analysis of the OH excitation to a forthcoming paper.

Fig. 5 P-Cygni profile of the OH 119 µm line toward the hot core region (the FIR continuum peak), compared to other line profiles. See Fig. B.3 for a zoom on the low-intensity features.

4.3 A possible outflow or wind bow shock around BN

Figure 6 (left) shows a map of the [O I] 63 µm emission in the v_LSR range from +20 to +25 km s⁻¹. In addition to the red-shifted extended emission from the explosive outflow, we detect an emission peak near BN, located slightly northwest (NW) of this young massive star (∼8–13 M_⊙; e.g., Rodríguez et al. 2005). This emission does not coincide with the more extended, so-called low-velocity “18 km s⁻¹ outflow,” which expands in the SW–NE direction (Genzel et al. 1981; Plambeck et al. 1982) and is likely driven by Source I (Beuther & Nissen 2008; Wright et al. 2022). The new [O I] 63 µm component arises from a broad emission feature, with a line FWHM of 15 ± 1 km s⁻¹, centered at v_LSR ≈ 20 km s⁻¹ (blue spectrum in Fig. 6). This emission is confined to the vicinity of BN, indicating either a compact outflow or a knot or shell of shock-excited gas surrounding the star. It envelops the IR nebulosity associated with BN (see Shuping et al. 2004; Bally et al. 2020) but has not been reported in other line tracers, likely owing to the nature of this shocked gas and the complex substructure of the region in terms of emission, velocity, and extinction.

The centroid of the [O I] 63 µm line feature is compatible with the LSR velocity measured from emission lines associated with BN (v_LSR ≃ 23 km s⁻¹; Gómez et al. 2008; Goddi et al. 2011; Plambeck et al. 2013). The apparent extent of this [O I] 63 µm feature is < 10″ and remains undetected at 145 µm. The blue- and redshifted [O I] 63 µm emission around the LSR velocity of BN shows a similar spatial distribution, suggesting that, if it originates from an outflow, the outflow is spatially unresolved. Given the uncertain properties of BN’s circumstellar disk, particularly its geometry (Jiang et al. 2005; Beuther et al. 2010; Bally et al. 2020), it remains unclear whether this structure traces a disk-driven outflow or another form of shocked gas.

Interestingly, the [O I] 63 µm intensity peak is shifted to the NW of BN by ≲3″, and is roughly aligned with the proper motion of this runaway star. Thus, this [O I] 63 µm emission may correspond to an unresolved wind bow shock produced by the interaction of BN’s supersonic motion and wind with the ambient cloud gas (Baranov et al. 1971; Wilkin 1996). Adopting estimates of the wind mass-loss rate (Ṁ_w), wind velocity (v_w), and BN’s velocity (v_∗) relative to the surrounding cloud of density n₀ (Scoville et al. 1983; Bally et al. 2020), the bow-shock radius (r_bs) would be several hundred au, or ≲ 1″ (Tan 2004). However, the wind mass-loss rate and velocity remain uncertain. Since $r_{\text{bs}}\propto {(}^{{\stackrel{\cdot }{M}}_w}{v}_{*}^{-1}$, larger wind parameters would yield a correspondingly larger r_bs (for fixed n₀ and v_∗), potentially consistent with the apparent size of the [O I] 63 µm feature.

5 Analysis

In the following, we restrict our analysis to a 100″ × 80″ (0.2 pc × 0.16 pc) region around the outflow, indicated with a dashed square in Fig. 3. Specifically, we excluded the Trapezium region. In addition, we used the spectra of Peak 1 as a template to determine the physical conditions of the post-shock gas.

5.1 Velocity-resolved line intensity ratios

Figure 7 shows the FIR line intensity ratios as a function of velocity toward the Peak 1 position (filled squares) obtained from maps convolved to a common angular resolution of 15″. The lowest [O I] 63/145 line intensity ratio, 11.4 ± 0.1, occurs at v_LSR,0, corresponding to the LSR velocity of the quiescent gas in OMC-1 (dashed magenta line in Fig. 7), where the narrow line component is dominated by emission from the extended PDR in front of the outflow. This intensity ratio is comparable to that observed in the Orion Bar (e.g., Bernard-Salas et al. 2012). The [O I] 63/145 intensity ratio then increases in the line wings, reaching ∼30 at redshifted velocities. The [O I] 63/145 intensity ratio shows an asymmetry between the blue- and redshifted wings. Foreground [O I] 63 µm absorption at v_LSR between roughly −5 and +7 km s⁻¹ contributes to this asymmetry, but the blueshifted wing shows smaller ratios at higher, more negative velocities. These lower ratios suggest a correspondingly lower density in the blueshifted post-shock gas (Sect. 5.2).

Figure 7 also shows the velocity dependence of the [O I] 145/[C II] 158, [O I] 63/[C II] 158, and OH 163/[C II] 158 line intensity ratios. All these ratios reach a minimum at the LSR velocity of the face-on PDR, v_LSR,0, and gradually increase in the line wings, indicating a lower contribution of C⁺ in the highest velocity gas. The [O I] 145 µm line is particularly bright, with [O I] 145/[C II] 158 intensity ratios of ≃2–3 in the line wings. The OH 163/[C II] 158 line intensity ratio increases from ≃0.1 at the line center to 4 ± 1 in the wings, thus showing a strong dependence on radial velocity. The redshifted line wing shows very high [O I] 63/[C II] 158 intensity ratios, up to ≃50–100, much higher than the expected ratios in PDRs with n_H ≈ 10⁵ cm⁻³ (e.g., Kaufman et al. 1999). For reference, in the prototypical dense PDR, the Orion Bar, at ∼2^′ south of the Trapezium, these line intensity ratios are much lower, [O I] 63/145 ≃ 9, [O I] 63/[C II] 158 ≃ 8, [O I] 145/[C II] 158 ≃ 0.9, and OH 163/[C II] 158 ≃ 0.03 (e.g., Bernard-Salas et al. 2012; Joblin et al. 2018).

The empty squares in Fig. 7 show the FIR line intensity ratios as a function of LSR velocity, averaged over the 100″ × 80″ region, thus representing the average conditions in the outflow. These ratios are similar to those toward Peak 1, with the main differences being a slightly higher [O I] 63/145 intensity ratio in the blueshifted wing and slightly lower [O I] 145/[C II] 158, [O I] 63/[C II] 158, and OH 163/[C II] 158 ratios. This reflects the higher relative contribution of the [C II] 158 µm emission⁴ in the SE part of the outflow (see Figs. 3b and 3c). Figure B.2 shows these results in the form of intensity-ratio velocity channel maps. These maps illustrate the same trends described above: the [O I] 63/[C II] 158 and [O I] 145/[C II] 158 ratios in the redshifted wing peak toward Peaks 1 and 2, following the H₂ outflow. These regions coincide with enhanced [O I] 63/145 intensity ratios.

Table 1 summarizes the intensities of the atomic fine-structure lines toward Peak 1, measured in three velocity intervals of the line profile: the redshifted wing, the blueshifted wing, and the narrow component. To increase the S/N of the line-intensity measurements and to cover the relevant outflow velocity range, we derived the wing intensities using velocity bins of 25 km s⁻¹ in the wings. We adopted the [+21, +46] and [−30, −5] km s⁻¹ LSR velocity intervals for the redshifted and blueshifted wings, respectively. That is, radial velocities of 13 to 38 km s⁻¹ with respect to v_LSR,0. This choice excludes the velocity ranges of the [O I] 63 and [C II] 158 µm lines that are affected by foreground absorption and emission in the blueshifted wing (see Fig. 4). We also subtracted the contribution of the [¹³C II] F = 2–1 hyperfine-structure line in the redshifted wing using multi-Gaussian fitting. In Sect. 6.2, we use these line-wing intensities to constrain the nature of the shocks that give rise to the FIR atomic line emission.

Fig. 6 Possible compact outflow or wind bow shock around BN star. Left: spatial distribution of the redshifted [O I] 63 µm emission in the vLSR range from +20 to +25 km s−1. This map (in RA and Dec) shows extended emission following the wide-angle H2 outflow and a possible spatially unresolved outflow around BN, revealed as a distinct broad spectral component. Right: [O I] 63 µm spectra toward Source I (red), BN (blue), and Peak 1 (black). Only the blue spectrum shows the broad spectral feature around the LSR velocity of BN (dashed green line).

Fig. 7 Line intensity ratios as a function of velocity in a 15″ beam toward Peak 1 (filled squares) and averaged over the 100″ × 80″ region of the BN/KL outflow (empty squares), binned in intervals of 5 km s−1 from maps convolved to a common resolution of 15″. The dashed magenta line marks the LSR velocity of the quiescent gas in OMC-1.

5.2 Post-shock gas physical conditions and [O I]63 145 µm excitation across different outflow velocity components

The observed [O I] 63/[O I] 145 µm intensity ratios (I₆₃/I₁₄₅) in the line wings are high, reaching up to ∼30 (e.g., Table 1). This is close to the high-temperature limit for optically thin emission under local thermodynamic equilibrium (LTE⁵) conditions (T_ex = T_k). This would imply N(O) < a few 10¹⁸ cm⁻² and high densities, meaning n_H ≫ n_cr, with n_cr,63µm ≃ 5 × 10⁵ cm⁻³ and n_cr,145µm ≃ 6 × 10⁶ cm⁻³ for optically thin emission. However, LTE is rarely fulfilled, and even for a three-level system such as the ³P_J levels of oxygen, subtle non-LTE (NLTE), line opacity, and FIR pumping effects can dominate the [O I] emission, thus, the resulting I₆₃/I₁₄₅ ratio (Monteiro & Flower 1987; Elitzur & Asensio Ramos 2006; Goicoechea et al. 2009; Goldsmith 2019).

To estimate the beam-averaged physical conditions (n_H, T_k) and the atomic oxygen column densities, N(O), in the post-shock gas across the different velocity components, we used a nonlocal NLTE radiative transfer model (Goicoechea et al. 2006b, 2009) with up-to-date O–H₂ and O–H inelastic collision rate coefficients (Lique et al. 2018). This nonlocal approach accounts for [O I] 63 µm line opacity and line trapping, as well as [O I] 145 µm suprathermal excitation (T_ex > T_k) and population inversions, which “local” single escape probability methods fail to correctly capture. To match the observed [O I] 145 µm line FWHM (broad component), we adopted a turbulent line width of 20 km s⁻¹. Following Goicoechea et al. (2009), we also examined the impact of FIR radiative pumping on the [O I] fine-structure level populations. Appendix C provides details and explores the excitation conditions over a wide parameter space.

Figure 8 presents the results of a model grid displayed in the velocity-dependent parameter space I₆₃/I₁₄₅ (v) versus I₁₄₅(v). The light grey to black curves show radiative transfer models for four atomic oxygen column densities, from 2×10¹⁶ to 2×10¹⁹ cm⁻², increasing from left to right, at two gas temperatures, 200 and 500 K. Tick marks denote log n_H (in cm⁻³) and are separated by one decade along curves of constant N(O). In general, as the opacity of the [O I] 63 µm line increases, the intensity ratio I₆₃/I₁₄₅ decreases roughly as ∼ 1/τ₆₃, while the [O I] 145 µm line is rarely optically thick. In dense gas, intensity ratios I₆₃/I₁₄₅ ≲ 15 typically imply optically thick [O I] 63 µm emission. For the adopted parameters, the gas density also determines the [O I] excitation regime. For n_H ≳ n_crit, the I₆₃/I₁₄₅ ratio increases with n_H owing to a significant increase in the [O I] 63 µm excitation temperature and a decrease in τ₆₃ toward the optically thin regime. In this regime, the ratio also increases with decreasing gas temperature. The colored squares in Fig. 8 represent the observed values in LSR velocity bins of 5 km s⁻¹ toward Peak 1 (filled squares) and averaged over the BN/KL outflow (open squares). Red squares correspond the redshifted line wing, in the range [+10, +25] km s⁻¹, while blue squares correspond to the blueshifted line wing, [−20, 0] km s⁻¹. Green squares in Fig. 8 correspond to the quiescent cloud gas at OMC-1 velocities, v_LSR = [+5, +10] km s⁻¹. The green squares are consistent with large column densities of atomic oxygen, N(O) ≃ 10¹⁹ cm⁻², optically thick [O I] 63 µm emission, n_H of several 10⁵ cm⁻³, and temperatures T_k ≃ 200–400 K – thus compatible with the observed T_P([O I]63 µm) values (Sect. 4.2) –, typical of a dense PDR with a thermal pressure, P_th/k, of 10⁷ to 10⁸ K cm⁻³. This agrees with previous observations of the large-scale PDRs at the surfaces of OMC-1 (Herrmann et al. 1997; Goicoechea et al. 2015b, 2019). An increase in velocity within the line wings (red and blue squares) is accompanied by a rise in the [O I] 63/145 intensity ratio, consistent with higher gas densities and less opaque [O I] 63 µm emission. The lower I₁₄₅(v) intensities as velocities increase simply reflect smaller N(O) column densities. Toward Peak 1, the [O I] 63 and 145 µm line-wing emission is consistent with n_H of several 10⁵ to 10⁶ cm⁻³, T_k ≲ 500 K, and N(O) of a few 10¹⁷ cm⁻². Figure 8 also indicates that the post-shock gas density is slightly lower in the blueshifted gas. When averaged over the entire outflow, the I₆₃/I₁₄₅(v) ratio is highest at the largest velocities, where the H₂ and O column densities are lowest. The higher ratios imply n_H up to several 10⁶ cm⁻³, consistent with gas compression in the post-shock gas.

Fig. 8 [O I] 63/145 µm intensity ratio versus [O I] 145 µm intensity. The light grey to black curves show a grid of radiative transfer models for four values of N(O): 2×1016, 2×1017, 2×1018, and 2×1019 cm−2, from left to right, at 200 and 500 K. Tick marks denote log nH (in cm−3) and are separated by one decade along curves of constant N(O). The colored squares show the observed values in LSR velocity bins of 5 km s−1 toward Peak 1 (filled squares) and averaged over the BN/KL outflow (open squares). Red squares correspond to the redshifted line wing, while blue squares correspond to the blueshifted line wing. Green squares correspond to OMC-1 velocities (see Tables D.2 and D.3).

6 Discussion

6.1 Peak 1 as a laboratory for shock modeling

Early models of Peak 1 have required two types of shocks to explain the observations (e.g., Chernoff et al. 1982; Draine & Roberge 1982; Neufeld & Dalgarno 1989; Haas et al. 1991; Kaufman & Neufeld 1996; Le Bourlot et al. 2002; Esplugues et al. 2014). The accepted picture was that a high-velocity molecular flow (the “plateau”) produces nondissociative C-type shocks (with shock velocities v_S ≳ 25 km s⁻¹) as the flow interacts with the ambient molecular cloud, generating high temperatures (above 1000 K). These shocks lead to bright H₂ and high-J CO emission along with large H₂O/CO > 1 abundance ratios. Newer C-type shock models included the effects of modest external UV illumination (Melnick & Kaufman 2015) to explain the unusual detection of narrow-line O₂ toward Peak 1 (Chen et al. 2014). In addition, fast, v_S = 70–80 km s⁻¹, dissociative, J-type shocks were proposed to occur where a high-velocity wind impacts the expanding swept-up material (Hollenbach 1985), producing bright atomic fine-structure line emission. However, BN/KL’s explosive outflow is driven by bullets moving through the molecular cloud (e.g., Bally et al. 2015, 2020), rather than by stellar jets or winds, as in protostellar outflows. In Sect. 6.2 we explore the nature of the shocks driving the atomic fine-structure emission.

Fig. 9 FIR atomic line intensities as a function of shock velocity and preshock gas density (from Lehmann et al. 2022). The right panels show externally irradiated shocks. The horizontal blue and red lines indicate the observed line-wing intensities toward Peak 1 (Table 1).

6.2 Nature of the [O I]–emitting shocks

Figure 9 compares the intensities of the [O I] 63 µm, [O I] 145 µm, and [C II] 158 µm lines observed in the red and blue wings toward Peak 1 (Table 1) with those predicted by the Paris–Durham shock code⁶ (Godard et al. 2019) along the direction perpendicular (face-on) to the shock front. In this comparison, we assume an emission beam filling factor of f = 1. The models are from Lehmann et al. (2022), which account for magnetized molecular shocks with velocities in the range v_S = 5–80 km s⁻¹, preshock gas densities in the range n_H,0 = 10²–10⁶ cm⁻³, enabling us to include an external FUV field $(G_0^{\text{ext}})$, in units of the Mathis field (Mathis et al. 1983).

The strength of the preshock magnetic field prior to the explosion is uncertain. Dust polarimetric observations of the plane-of-sky magnetic field in the quiescent cloud indicate values of a few hundred µG (Houde et al. 2009; Guerra et al. 2021), consistent with line-of-sight field strengths derived from CN Zeeman measurements (Crutcher et al. 1999; Crutcher 1999). Thus, all models in Fig. 9 adopt a reference magnetic parameter b = 1, where the preshock transverse magnetic field strength is given by B = b (n_H,0/cm⁻³)^1/2 µG. For this range of parameters, shocks can be of type C, C^∗, or CJ at v_S < 25 km ⁻¹, or J-type at v_S > 25 km⁻¹ (Lehmann et al. 2020, 2022). The models assume either $G_0^{\text{ext}}=0$ (left panels) or $G_0^{\text{ext}}={10}^3$ (right panels). Accurately determining the flux of external FUV photons reaching each bow-shock structure–whether locally self-generated, generated in nearby shocks, emitted by massive protostars within the region (BN, source I, and source n; Shuping et al. 2004), or originating from the Trapezium stars–is challenging, as it depends on both the porosity of the medium and the properties of the shocks and the protostars. The representative value of $G_0^{\text{ext}}={10}^3$ lies within the very rough estimates from FIR continuum observations toward Peak 1 (Goicoechea et al. 2015a).

In the absence of external irradiation $(G_0^{\text{ext}}=0)$, the predicted [O I] line intensities converge toward a plateau at high shock velocities, the level of which is primarily set by the preshock gas density. This is consistent with the dissociative shock models of Hollenbach & McKee (1989). Shocks with sufficiently high velocities generate a local FUV radiation field (e.g., Lehmann et al. 2020), whose intensity scales as $$G_0^{\text{shock}}={10}^3\left(\frac{n_{\text{H},0}}{{10}^4\text{\hspace{0.17em}cm}{\text{}}^{-3}}\right){\left(\frac{{\upsilon }_S}{80\text{\hspace{0.17em}km\hspace{0.17em}s}{\text{}}^{-1}}\right)}^3,$$(2)

for 30 < v_S ≲ 100 km s⁻¹ (Lehmann et al. 2020; Godard et al. 2024), and thus strongly depend on v_S . At low shock velocities, models with $G_0^{\text{ext}}=0$ underpredict the observed [C II] 158 µm intensity by several orders of magnitude. Marginally reproducing the [C II]158 µm intensity would require very fast shocks, v_S ≳ 80 km s⁻¹. However, this is in contradiction with the ∆v_FWHM and FWZI values of the observed fine-structure lines (Sect. 4.2), assuming that the flow is roughly directed along the line of sight. Therefore, a significant contribution from locally self-generated UV radiation is unlikely.

External FUV radiation, which leads to carbon photoionization and CO photodissociation, significantly enhances the intensities of atomic fine-structure lines (right panels of Fig. 9). Still, owing to the enhanced importance of C⁺ recombination and chemical destruction at high gas densities, the [C II] 158 µm intensity does not increase monotonically with density.

Under external irradiation, the observations no longer constrain a unique set of shock parameters, but are instead consistent with a family of solutions involving J-type shocks with v_S ≃ 30–40 km s⁻¹ and n_H,0 ≳ a few 10⁴ cm⁻³. Models with $G_0^{\text{ext}}={10}^3$ reproduce the [C II] 158 µm intensity of the redshifted wing, tracing gas moving into OMC-1, whereas the blueshifted wing emission, tracing gas moving toward the Trapezium, appears to require higher values of $G_0^{\text{ext}}$. In the following, we focus on the redshifted wing emission.

In these irradiated shocks, the predicted I₆₃/I₁₅₈ line intensity ratio primarily scales with the preshock gas density and shows only a weak dependence on the shock velocity for v_S > 30 km s⁻¹ (see Fig. 10). The above family of solutions predicts I₆₃/I₁₅₈ ≳ 60, consistent with the redshifted wing emission toward Peak 1, which shows intensity ratios up to ∼50–100 (Fig. 7). Given the inferred post-shock densities from the [O I] excitation analysis in Sect. 5.2, n_H of several 10⁵ to 10⁶ cm⁻³, the resulting compression factors, F = n_H/n_H,0 ≳ 30, also support magnetized J-type shocks, where F ∼ v_S /b (Godard et al. 2024). The compressed magnetic field would have a strength of several mG, similar to the values observed toward BN/KL (Pattle et al. 2017; Guerra et al. 2021). All in all, the [O I] and [C II] line-wing emission is best reproduced by externally irradiated, dissociative J-type shocks with velocities of 30–40 km s⁻¹, preshock densities of a few 10⁴ cm⁻³, and compression factors consistent with post-shock densities derived from the [O I] excitation analysis.

A more detailed interpretation would require the geometry, time dependence, and spatial variations of b and $G_0^{\text{ext}}$ across the outflow to be fully accounted for. Still, the results shown in Fig. 9 provide a general interpretative framework: dense bullets, with proper motions up to several hundred km s⁻¹, accelerated by the explosion event (i.e., Doi et al. 2002; Bally et al. 2015, 2017), interact with the ambient molecular cloud, generating bow-shock structures seen in [Fe II] knots and H₂ (Fig. 1). The highest bullet velocities exceed the observed FWZIs of the CO and OH lines, and imply molecular gas dissociation and the production of significant FUV radiation in many bow-shock tips (Eq. (2)). Within this framework, surrounding fast shocks generate sufficient external FUV radiation to illuminate the inner layers of the bow-shock structures. The working surfaces of these bow shocks drive dissociative J-type shocks with velocities of a few tens of km s⁻¹, producing the observed [O I] emission.

Consistently with this scenario, the observed H I 21 cm line emission from atomic hydrogen (produced after H₂ (photo)dissociation) extends only up to v_LSR ≤ 31 km s⁻¹ (see Fig. 20 of van der Werf et al. 2013) and broadly resembles the [O I] 63 µm emission, both in its spatial distribution and in its lack of the highest velocity emission seen in CO or OH. The broader CO and OH line profiles indicate that the shocks responsible for their emission differ from those driving the atomic fine-structure emission. For instance, Melnick & Kaufman (2015) concluded that a C-shock with v_S = 23 km s⁻¹ propagating into a more shielded molecular clump (preshock density of 8 × 10⁴ cm⁻³ and $G_0^{\text{ext}}=1$) explains the detection of O₂ (Chen et al. 2014).

We conclude that, in contrast to what was found using previous models (e.g., Hollenbach & McKee 1989; Haas et al. 1991), fast shocks with v_S ≃ 70–80 km, s⁻¹ are not required to explain the observed atomic fine-structure line intensities and profiles. Instead, the observed lines are consistent with lower velocity, dissociative J-type shocks illuminated by external FUV radiation, with $G_0^{\text{ext}}\simeq {10}^2-{10}^3$, depending on the shock parameters (e.g., Kristensen et al. 2023). This radiation arises from surrounding fast shocks and possibly from nearby massive (proto)stars, enhancing the [C II] 158 µm emission, particularly in blueshifted gas moving toward the cluster.

Fig. 10 [O I] 63/[C II] 158 µm line intensity ratio as a function of shock velocity and preshock gas density in an externally irradiated shock with $G_0^{\text{ext}}={10}^3$ (from Lehmann et al. 2022). The horizontal blue and red lines indicate the observed line-wing intensities toward Peak 1.

6.3 Mass-loss rate of the explosive outflow in Orion BN/KL

The energy loss from the outflow-driven shock provides a measure of the kinetic (mechanical) luminosity of the outflow, with $L_{\text{cool}}=(1-f_m)L_{\text{kin}}=(1-f_m)\frac{1}{2}{\stackrel{\cdot }{M}}_0{v}_S^2$. In this expression⁷, L_cool is the total shock cooling rate and (1 − f_m) denotes the fraction of the shock mechanical energy translated into internal excitation. In general, gas cooling proceeds through IR, H I, and optical line emission, L_cool ≈ L_IR + L_HI + L_optical, where $L_{\text{IR}}=L_{\text{atoms}}+L_{{\text{H}}_2}+L_{\text{molec}}$ is the contribution from atomic fine-structure lines, H₂ lines, and other molecular lines (CO, H₂O, OH, etc.). L_HI and L_optical are only relevant in dissociative J-type shocks (Lehmann et al. 2020). L_HI represents the hydrogen-line cooling emission (e.g., Lyα, Lyβ, and two-photon continuum) and L_optical refers to “optical” atomic forbidden-line emission, which arises from electronic transitions in the visible and NIR. From the above relation⁷, $${\stackrel{\cdot }{M}}_0=1.2\cdot {10}^{-2}\frac{{(1-f_m)}^{-1}}{f_{\text{IR}}}\left(\frac{L_{\text{IR}}}{{\text{L}}_{\odot }}\right){\left(\frac{{\upsilon }_S}{\text{km\hspace{0.17em}s}{\text{}}^{-1}}\right)}^{-2}{\text{\hspace{0.17em}M}}_{\odot }\text{\hspace{0.17em}yr}{\text{}}^{-1},$$(3)

where f_IR is the fraction of the total cooling due to IR emission lines. In general, f_IR ≈ 0.5–0.8 in C-type shocks (e.g., Kaufman et al. 1999; Karska et al. 2025) and f_IR ≈ 0.5 in J-type shocks (e.g., Lehmann et al. 2020). Equation (3) has clear limitations when different types of shocks coexist, or when part of the emitted line luminosity arises from FUV-heated gas. In the following analysis, we assume that v_S ≃ 30–40 km s⁻¹ represents the velocity range of the possible shock types in BN/KL’s outflow (if more than one), whether dissociative or nondissociative. Table 2 summarizes the IR line luminosities, integrated over the entire profile, across the outflow (dashed square in Fig. 3). Inserting the outflow IR line luminosity, L_IR ≃ 320 L_⊙, Eq. (3) yields a mass-loss rate of Ṁ₀ ≃ (9.1 ± 2.6) × 10⁻³ M_⊙ yr⁻¹, using f_IR = 0.5 (Fig. 3 of Lehmann et al. 2020) and (1 − f_m) = 0.75 (Kaufman & Neufeld 1996). Even if the derived Ṁ₀ is an approximate value, the estimated mass-loss rate of BN/KL’s outflow is higher than that of steady protostellar outflows. Given a dynamical timescale of t_dyn ≃ 500 yr (e.g., Zapata et al. 2009), the derived outflow mass is M = t_dyn Ṁ₀ ≃ (3.3–5.9) M_⊙. This estimate agrees with, yet may be more accurate than, mass estimates obtained from column densities computed under the LTE assumption (Snell et al. 1984; Peng et al. 2012). Still, future theoretical work should revisit the validity of Eq. (3) in the complex scenario described in Sect. 6.2, considering its dependence on multiple shocks and $G_0^{\text{ext}}$.

6.4 Mass-loss rate from the putative BN outflow

Here, we assume that the spatially unresolved [O I] 63 µm emission detected around BN – both spatially and in LSR velocity (see Sect. 4.3) – arises from shocked gas produced by a compact outflow launched by this star. When [O I] 63 µm is the dominant gas coolant in a J-type “wind” shock–consistent with detection of this outflow only in [O I] 63 µm–Hollenbach (1985) showed that the line luminosity scales with the outflow mass-loss rate as Ṁ (M_⊙ yr⁻¹) ≃ 10⁻⁴ L₆₃ (L_⊙), provided that the flux of material into the shock satisfies n_H,0 v_S ≲ 10¹² cm⁻² s⁻¹. In this regime, the particle flux is related to the [O I] 63 µm line intensity by n_H,0 v_S ≃ 10¹³ I₆₃, with n_H,0, v_S , and I₆₃ in cgs units.

The broad [O I] 63 µm line emission around BN is characterized by I₆₃ ≃ 5×10⁻² erg s⁻¹ cm⁻² sr⁻¹ and L₆₃ ≃ 0.6 L_⊙ in ≲ 4000 au. Thus, its contribution to the total L₆₃ of the entire map is negligible. Using the relations above, we obtain, n_H,0 v_S ≃ 5×10¹¹ cm⁻² s⁻¹ and Ṁ ≃ 6 × 10⁻⁵ M_⊙ yr⁻¹, which is typical of accreting high-mass protostars (e.g., Beuther et al. 2002; Leurini et al. 2015). Assuming further that the shock velocity is comparable to the [O I] 63 µm line FWHM of the spectral component around BN, ∆v_FWHM ≃ 15 km s⁻¹ (Sect. 4.3), we derive a preshock gas density of several 10⁵ cm⁻³. These results indicate that the kinematics, density, and mass-loss rate associated with the compact [O I] 63 µm emission around BN are consistent with a disk-accretion–driven outflow originating from BN. Further observations of additional tracers at high angular resolution will be needed to confirm the outflow versus wind bow shock scenario, or to constrain alternative interpretations.

6.5 Comparison with other high-mass star-forming regions

To investigate whether explosive outflows produce distinct spectroscopic features, here we compare the observed FIR [O I], [C II], and CO line luminosities of the BN/KL outflow with those measured in the central regions of other, more distant high-mass star-forming regions where similar observations have been carried out (e.g., Leurini et al. 2015). In the BN/KL outflow region, the FIR atomic line luminosity is dominated by [O I] 63 µm⁸, which is ∼15 times brighter than the [C II] 158 µm emission. The total [O I] 63 µm and 145 µm line luminosity is remarkably high, 86.5 L_⊙, comparable to the total H₂ and CO luminosity (Rosenthal et al. 2000; Goicoechea et al. 2015a). This yields L_[OI]/L_CO ≃ 0.9 and $L_{[\text{OI}]}/L_{{\text{H}}_2}\simeq 0.7$. That is, about 25% of the total IR gas cooling is through [O I] emission lines (Table 2). Furthermore, the contribution to L_[OI] from the broad spectral component, neglecting the narrow spike [O I] emission, is 55 L_⊙.

Table 3 compares the FIR line luminosities toward the central regions of several high-mass star-forming regions, on scales of ∼0.2 pc comparable to our BN/KL map. Interestingly, interferometric observations of the central regions of DR21 and G5.89−0.39 reveal CO outflows exhibiting fingers and Hubble-Lemaître flow kinematics, consistent with the presence of explosive outflows (e.g., Zapata et al. 2013, 2020; Guzmán Ccolque et al. 2024). These other regions exhibit fractional luminosities comparable to those of BN/KL outflow, both displaying high L_OI/L_bol and $L_{\text{CO}}^{\text{FIR}}/L_{\text{bol}}$ ratios, ≳10⁻⁴. In contrast, W3-IRS5, a proto-Trapezium system (e.g., Rodón et al. 2008), and NGC7538-IRS1, a massive star-forming core hosting an ultracompact H II region, a late O–to–early B star (Sandell & Sievers 2004), exhibit lower L_OI/L_bol and $L_{\text{CO}}^{\text{FIR}}/L_{\text{bol}}$ luminosities. Pending more robust statistics, these high fractional luminosities may represent a distinctive signature of explosive outflows.

Still, the BN/KL outflow stands out in FIR [O I] and CO line surface luminosities, per unit area, and the I₆₃/I₁₅₈ intensity ratios exceed those in most regions, even when integrating over the full line profile, suggesting a denser gas environment.

7 Summary and conclusions

In this paper, we present velocity-resolved, sub-km s^−1- resolution maps of the [O I] 63 and 145 µm fine-structure emission lines along the wide-angle explosive outflow in Orion BN/KL. We complement this data with new [C II] 158 µm and FIR OH line maps, all obtained with the GREAT heterodyne receivers on board SOFIA. Our main conclusions are as follows:

• We separated the quiescent and the outflow (broad [O I] line-wing) components, with the latter following the shock-excited H₂ emission and the former primarily tracing the foreground, face-on PDR at the UV-illuminated rims of OMC-1. The [O I] 63 µm line has a FWZI of ∼85 km s⁻¹. The FWZI of the [C II] 158 µm line is smaller, ≲ 50 km s⁻¹, implying that either the highest velocity gas lacks FUV illumination or (more likely) C⁺ is quickly converted into other species;

• The broad [O I] 145 µm component (∆v_FWHM ≃ 20 km s⁻¹) is narrower than the corresponding CO and OH emission (∆v_FWHM ≳ 30 km s⁻¹), which also exhibit larger FWZIs (≳ 150 km s⁻¹). This suggests that the shocks producing the [O I] emission differ either in their nature or in exact location from those driving the highest velocity CO and OH emission;

• The [O I] 63/145 and [O I] 63/[C II] 158 intensity ratios attain remarkably high values in the line wings (20–30 and 40–60, respectively), exceeding those observed in PDRs and typical protostellar outflows. They are consistent with dense (n_H ≳ 10⁵–10⁶ cm⁻³) and warm (T ≲ 500 K) post-shock gas;

• The fine-structure line-wing emission is consistent with originating in externally irradiated, magnetized, dissociative J-type shocks with v_S ≃ 30–40 km s⁻¹ and preshock densities of a few 10⁴ cm⁻³. External FUV radiation arises from surrounding fast shocks and possibly from massive (proto)stars in the region if the medium is sufficiently porous;

• The [O I] 63 and 145 µm luminosity along the outflow region is significantly high, amounting to 86.5 L_⊙, with 55 L_⊙ contributed by the broad component alone. This is comparable to the H₂ and CO line luminosities, yielding L_[OI]/L_CO ≃ 0.9 and $L_{[\text{OI}]}/L_{{\text{H}}_2}\simeq 0.7$. These values imply an outflow mass-loss rate of (9.1 ± 2.6) × 10⁻³ M_⊙ yr⁻¹ and a mass of M ≃ (3.3–5.9) M_⊙;

• We report the detection of broad [O I] 63 µm emission, slightly offset from BN but aligned with the direction of the star’s proper motion. This L₆₃ ≃ 0.6 L_⊙ emission could arise from an unresolved outflow (corresponding to a mass-loss rate of several 10⁻⁵ M_⊙ yr⁻¹) or a wind bow shock produced by the supersonic motion of this wind-blowing runaway B star.

The FIR [C II] and [O I] line intensities are among the brightest tracers of interstellar shocks and PDRs, probing the thermal budget and the presence of FUV radiation. In addition, sub-km s⁻¹-resolution maps reveal the gas dynamics, allowing us to quantify the contributions of radiative and mechanical feedback. Unfortunately, given the cessation of SOFIA operations, very-high-spectral-resolution mapping has been carried out for only a few regions of the sky. Therefore, a FIR heterodyne mission (airborne or space) would be the only telescope to provide this essential missing information.

Acknowledgements

We thank our referee for a constructive report that allowed us to improve the presentation of our work. This work is based on observations made with SOFIA (NASA/DLR). SOFIA was jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under the DLR contract 50 OK 0901 and 50 OK 1301 to the University of Stuttgart. GREAT was a development by the MPI für Radioastronomie and the KOSMA/Universität zu Köln, in cooperation with the DLR Institut für Optische Sensorsysteme. GREAT development was financed by the participating institutes, by DLR under Grants 50 OK 1102, 1103, and 1104, and within the Collaborative Research Centre 956, funded by the Deutsche Forschungsgemeinschaft (DFG). JRG, MGSM, and MZ thank the Spanish MICIU for funding support under grant PID2023-146667NB-I00. MZ acknowledges the JdC Postdoctoral Fellowship JDC2024-054658-I, funded by MICIU/AEI/10.13039/501100011033 and by the ESF+. BG and AG acknowledge the support from the Programme National “Physique et Chimie du Milieu Interstellaire” (PCMI) of CNRS/INSU with INC/INP funded by CEA and CNES.

References

Appendix A Observed lines, instrument configurations, and observing modes

Table A.1 summarizes the instrument configurations, observing parameters, and observed FIR line properties.

The observing modes are indicated: total power (TP) or chopped (ch) on-the-fly (OTF: data is acquired while slewing the telescope), or raster mode (the telescope stays pointed while integrating). The size of a given map is quoted for the central array pixel only; the step between dumps (along the scanning direction) and between adjacent OTF slews is given or the grid used for pointed raster observations.

Appendix B Complementary observational figures

This appendix presents additional figures that aim to better understand the velocity structure of BN/KL outflow.

Appendix B.1 Velocity channel maps

Figure B.1 shows velocity-channel maps (in K km s⁻¹ units) of the FIR lines observed with SOFIA and the CO J = 10–9 line observed with Herschel/HIFI (Goicoechea et al. 2019).

Fig. B.1 Velocity channel maps at their native angular resolution and in K km s−1 intensity units, from vLSR = −10 to +25 km s−1 in bins of 5 km s−1 [v, v + 5], with v indicated in the bottom-right corner of each panel. The map at v = +5 km s−1 shows the emission at the systemic velocity of the quiescent gas in OMC-1, vLSR,0 ≃ 8–9 km s−1, dominated by extended emission from a face-on PDR. The main reference positions discussed in the text are illustrated with symbols (see Fig. 2).

Figure B.2 shows channel maps of line intensity ratios for the different FIR atomic lines (with intensities in cgs units).

Fig. B.2 Velocity-resolved line intensity ratio channel maps (same velocity bins as in Fig. B.1) convolved at a common angular resolution of 15″, with the integrated line intensities in erg s−1 cm−2 sr−1. The map at v = +5 km s−1 shows the emission at the systemic velocity of the quiescent gas.

Appendix B.2 Line profiles toward the FIR continuum peak

Figure B.3 (right panel) shows a zoom on the low-intensity features of several FIR lines toward the FIR continuum peak position (the hot core region), with the main [O I] 63 µm absorption and [C II] 158 µm emission features labelled.

Fig. B.3 OH 119 µm P-Cygni profile toward the hot core (FIR continuum peak), compared to other line profiles. The OH 119 µm absorption dip agrees with [O I] 63 µm absorption at +6 km s−1 (hot core velocities). The [C II] 158 µm emission and [O I] 63 µm absorption at +2 km s−1 corresponds to a gas structure in the foreground.

Appendix C Nonlocal NLTE [O I] radiative transfer models

In our [O I] excitation and radiative transfer models (Goicoechea et al. 2006b, 2009), we adopted the fine-structure O–H₂ and O–H collisional rate coefficients calculated by Lique et al. (2018) in the 10–1000 K temperature range. As observed, the FIR [O I] fine-structure line emission in the outflow arises from molecular gas with some atomic H (van der Werf et al. 2013). Hence, we define the gas density as n_H = n(H) + 2n(H₂), and adopted n(H₂)/n(H) ≃ 4 in the models. We verified that this choice does not affect our main conclusions and results. We ran spherical models with uniform gas densities (n_H = 10⁴–10⁸ cm⁻³), temperatures (T_k = 200 and 500 K, based on the observed T_mb([O I] 63 µm) line peaks), and a fixed velocity dispersion, σ_turb = 8.5 km s⁻¹ (i.e., ∆v_turb,FWHM ≃ 20 km s⁻¹), matching the observed [O I] 145 µm broad component toward Peak 1.

Figure C.1 shows the results of the grid, displaying I₁₄₅ as a function of n_H (left panels) for different atomic oxygen column densities: 2 × 10¹⁷, 2 × 10¹⁸, and 2 × 10¹⁹ cm⁻², and gas temperatures. The right panels shows the resulting I₆₃/I₁₄₅ line intensity ratio as a function of n_H. Continuos curves refer to models without external FIR illumination. The red and cyan colored areas show the I₆₃/I₁₄₅ ratios measured toward Peak 1 in the red- and blueshifted wings, respectively (extracted from Table 1). These ratios are consistent with the presence of dense (n_H ≃ several 10⁵ to 10⁶ cm⁻³) and warm (T ≲ 500 K) post-shock gas, with N(O) ≃ a few 10¹⁸ cm⁻².

In these plots, the opacity of the [O I] 63 µm line increases as N(O) increases, from optically thin to thick emission. That is, from high to low densities–from right to left in the figures–as collisional excitation becomes less important and the low energy levels become more populated. This is reflected in very sub-thermal [O I] 63 µm excitation (T_{ex,63 µm} ≪ T_k) when n_H ≪ n_{cr,63 µm} (Fig. C.2, left panels). Figure C.2 shows the average T_{ex,63 µm} and T_{ex,145 µm} in the models, but we note that due to radiative transfer effects—such as line trapping—the local T_ex values vary with position, even under uniform physical conditions. The excitation of the [O I] 145 µm line is more complex and is governed by the spontaneous radiative decay of the upper level (³P₀), which is roughly five times lower than that of the lower level (³P₁). In addition, the collisional de-excitation rate coefficient for O–H₂ collisions of the ³P₀–³P₁ (145 µm) transition is significantly smaller than that of the ³P₁–³P₂ (63 µm) transition (e.g., Monteiro & Flower 1987; Lique et al. 2018). These properties lead to [O I] 145 µm population inversions (T_{ex,145 µm} < 0 K) at low densities (n_H < n_{cr,145 µm}) and suprathermal emission (T_{ex,145 µm} > T_k) at high densities (previously discussed by Liseau et al. 2006; Elitzur & Asensio Ramos 2006; Goicoechea et al. 2009; Goldsmith 2019). These features are evident in the right panels of Fig. C.2.

To investigate the role of FIR pumping in the [O I] excitation, the dotted curves in Figs. C.2 and C.1 represent the same models but adding external FIR continuum illumination. We model the continuum background as the sum of the cosmic microwave background and a modified blackbody, η B(T_d)(1 − exp(−τ_d,λ)), at a dust color temperature of T_d = 100 K. In this expression, η is the filling factor of the background FIR continuum emission, which we set to 0.5. The dust continuum emission is assumed to be optically thick at 100 µm, with τ_d,λ = 7 [100/λ(µm)]², and represents the typical FIR dust continuum emission close to the hot core region (see Goicoechea et al. 2015a). The main role of this FIR continuum is to produce [O I] 63 µm line absorption (when T_ex < T_cont) at low gas densities, n_H below a few 10⁴ cm⁻³, as the [O I] 63 µm line becomes optically thick, for N(O) ≳ 10¹⁸ cm⁻².

In general, FIR pumping increases the population of the ³P₁ level (in this case increasing the excitation temperature of the [O I] 63 µm line; see dotted curves in Fig. C.2, left), thereby reducing the population inversion of the ³P₀–³P₁ transition at [O I] 145 µm. At low to moderate gas densities, the [O I] 63 µm line intensity decreases, whereas the intensity of the [O I] 145 µm remains nearly unchanged. Consequently, the I₆₃/I₁₄₅ intensity ratio is slightly lower (dotted curves in Fig. C.1, right) than in the nonilluminated case. Hence, the inferred post-shock gas density may be slightly higher.

At high densities, however, the main effect of FIR pumping is to slightly reduce T_{ex,145 µm}, although it remains suprathermal (Fig. C.2, right), leading to a slightly higher I₆₃/I₁₄₅ intensity ratio (Fig. C.1, right).

Fig. C.1 Nonlocal and NLTE [O I] excitation and radiative transfer models: Predicted [O I] 145 µm line intensity (left) and [O I] 63/145 µm intensity ratio in cgs units (right) as a function of total gas density. The black, grey, and blue curves represent different column densities of atomic oxygen. Dotted curves represent models with external FIR illumination (see text). In the left panels, the red and cyan dashed lines show the [O I]145 µm intensities measured toward Peak 1 in the red- and blueshifted wings, respectively. In the right panels, the red and cyan dashed colored areas show the intensity ratios (with errors) measured toward Peak 1 in the red- and blueshifted wings, respectively (see Table 1).

Fig. C.2 Average excitation temperature of the [O I] 63 µm and [O I] 145 µm lines as a function of gas density for two gas temperatures, 200 and 500 K. The black, grey, and blue curves represent different column densities of atomic oxygen, indicated in the right panel. Dotted curves represent models with external FIR illumination (see text). The increase in the [O I] 63 µm excitation temperature at low densities reflects the role of FIR pumping by a modified black body at Td = 100 K. This effect reduces the population inversion of the transition at 145µm at low densities.

Appendix D Line intensity tables

In this appendix we tabulate the FIR line intensities toward the Trapezium position and Peak 1, as well as the average line intensities over the entire outflow region.

All Tables

All Figures

	Fig. 1 BN/KL outflow and the Trapezium cluster observed with JWST/NIRCam H2 F212N (McCaughrean & Pearson 2023). The dashed square marks the field of view mapped with SOFIA/GREAT. Cyan contours show the [O I] 63 µm redshifted line-wing emission over vLSR = 25–30 km s−1, from 15 to 90 K km s−1 in steps of 15 K km s−1.
In the text

Fig. 2 Total line intensity maps (in erg s−1 cm−2 sr−1 = “cgs”) integrated over the complete line profile. The beam size is indicated in the bottom-right corner of each panel except for H2, where it is too small to display (0.1″). The cyan triangle shows a position near the Trapezium, where the [O I] 63µm line intensity peaks (Table D.1). The H2 image refers to the JWST/NIRCam F212N image (McCaughrean & Pearson 2023). The CO J = 10–9 map was obtained with Herschel/HIFI at 20″ resolution (Goicoechea et al. 2019).

In the text

Fig. 3 Total line intensity ratio maps (derived from line intensities in erg s−1 cm−2 sr−1). (a) [O I] 63/145 at a common angular resolution of 13″, (b) [O I] 63/[C II] 158 at 15″, and (c) [O I] 145/ [C II]158 at 15″. The dashed box shows the (100″×80″) area used to extract both the line luminosities in the outflow region (Table 2) and the line intensities as a function of velocity (see Figs. 7 and 8). Contours show the FIR 79 µm-continuum obtained with Herschel/PACS, from 2.5 to 10 (103 Jy), centered at the position of the hot core (Goicoechea et al. 2015a).

In the text

Fig. 4 Velocity-resolved spectra at representative positions, all from maps convolved to 15″ (except CO 10–9, with a beam of 20″). Left: Complete spectra, with offsets in arcseconds given in parenthesis; continuum levels of [C II] 158 µm and [O I] 145 µm shifted for clarity. The OH 163 µm emission lines (yellow) are scaled by a factor of three or five. Right: zoom on faint line wings and foreground line features.

In the text

	Fig. 5 P-Cygni profile of the OH 119 µm line toward the hot core region (the FIR continuum peak), compared to other line profiles. See Fig. B.3 for a zoom on the low-intensity features.
In the text

Fig. 6 Possible compact outflow or wind bow shock around BN star. Left: spatial distribution of the redshifted [O I] 63 µm emission in the vLSR range from +20 to +25 km s−1. This map (in RA and Dec) shows extended emission following the wide-angle H2 outflow and a possible spatially unresolved outflow around BN, revealed as a distinct broad spectral component. Right: [O I] 63 µm spectra toward Source I (red), BN (blue), and Peak 1 (black). Only the blue spectrum shows the broad spectral feature around the LSR velocity of BN (dashed green line).

In the text

	Fig. 7 Line intensity ratios as a function of velocity in a 15″ beam toward Peak 1 (filled squares) and averaged over the 100″ × 80″ region of the BN/KL outflow (empty squares), binned in intervals of 5 km s−1 from maps convolved to a common resolution of 15″. The dashed magenta line marks the LSR velocity of the quiescent gas in OMC-1.
In the text

Fig. 8 [O I] 63/145 µm intensity ratio versus [O I] 145 µm intensity. The light grey to black curves show a grid of radiative transfer models for four values of N(O): 2×1016, 2×1017, 2×1018, and 2×1019 cm−2, from left to right, at 200 and 500 K. Tick marks denote log nH (in cm−3) and are separated by one decade along curves of constant N(O). The colored squares show the observed values in LSR velocity bins of 5 km s−1 toward Peak 1 (filled squares) and averaged over the BN/KL outflow (open squares). Red squares correspond to the redshifted line wing, while blue squares correspond to the blueshifted line wing. Green squares correspond to OMC-1 velocities (see Tables D.2 and D.3).

In the text

	Fig. 9 FIR atomic line intensities as a function of shock velocity and preshock gas density (from Lehmann et al. 2022). The right panels show externally irradiated shocks. The horizontal blue and red lines indicate the observed line-wing intensities toward Peak 1 (Table 1).
In the text

	Fig. 10 [O I] 63/[C II] 158 µm line intensity ratio as a function of shock velocity and preshock gas density in an externally irradiated shock with $G_0^{\text{ext}}={10}^3$ (from Lehmann et al. 2022). The horizontal blue and red lines indicate the observed line-wing intensities toward Peak 1.
In the text

Fig. B.1 Velocity channel maps at their native angular resolution and in K km s−1 intensity units, from vLSR = −10 to +25 km s−1 in bins of 5 km s−1 [v, v + 5], with v indicated in the bottom-right corner of each panel. The map at v = +5 km s−1 shows the emission at the systemic velocity of the quiescent gas in OMC-1, vLSR,0 ≃ 8–9 km s−1, dominated by extended emission from a face-on PDR. The main reference positions discussed in the text are illustrated with symbols (see Fig. 2).

In the text

	Fig. B.2 Velocity-resolved line intensity ratio channel maps (same velocity bins as in Fig. B.1) convolved at a common angular resolution of 15″, with the integrated line intensities in erg s−1 cm−2 sr−1. The map at v = +5 km s−1 shows the emission at the systemic velocity of the quiescent gas.
In the text

	Fig. B.3 OH 119 µm P-Cygni profile toward the hot core (FIR continuum peak), compared to other line profiles. The OH 119 µm absorption dip agrees with [O I] 63 µm absorption at +6 km s−1 (hot core velocities). The [C II] 158 µm emission and [O I] 63 µm absorption at +2 km s−1 corresponds to a gas structure in the foreground.
In the text

Fig. C.1 Nonlocal and NLTE [O I] excitation and radiative transfer models: Predicted [O I] 145 µm line intensity (left) and [O I] 63/145 µm intensity ratio in cgs units (right) as a function of total gas density. The black, grey, and blue curves represent different column densities of atomic oxygen. Dotted curves represent models with external FIR illumination (see text). In the left panels, the red and cyan dashed lines show the [O I]145 µm intensities measured toward Peak 1 in the red- and blueshifted wings, respectively. In the right panels, the red and cyan dashed colored areas show the intensity ratios (with errors) measured toward Peak 1 in the red- and blueshifted wings, respectively (see Table 1).

In the text

Fig. C.2 Average excitation temperature of the [O I] 63 µm and [O I] 145 µm lines as a function of gas density for two gas temperatures, 200 and 500 K. The black, grey, and blue curves represent different column densities of atomic oxygen, indicated in the right panel. Dotted curves represent models with external FIR illumination (see text). The increase in the [O I] 63 µm excitation temperature at low densities reflects the role of FIR pumping by a modified black body at Td = 100 K. This effect reduces the population inversion of the transition at 145µm at low densities.

In the text