© ESO 2020

1 Introduction

Stellar feedback, that is injection of energy and momentum from stars, is one of the most important input parameters in cosmological models that simulate and explain the evolution of our universe. Even small variations in this crucial parameter can lead to drastic changes in the results. Too much stellar feedback of early stars disrupts the ambient gas in too early a stage to form more stars and, eventually, planetary systems that allow for the formation of life. Too little feedback fails to prevent the interstellar gas from gravitational collapse to dense, cold clumps. Stellar feedback is often quantified as the star-formation rate (SFR), the rate at which (mostly low-mass) stars are formed. Studies on the interaction of massive stars with their environment generally focus on the effects of supernova explosions injecting mechanical energy into their environment and the radiative interaction leading to the ionization and thermal expansion of the gas. The explosion of a massive star ejects some 10 M_⊙ at velocities of some 10 000 km s⁻¹, injecting about 10⁵¹ erg into the surrounding medium. The hot plasma created by the reverse shock drives the expansion of the supernova remnant. The concerted effects of many supernovae in an OB association create superbubbles that expand perpendicular to the plane and may break open, releasing the hot plasma and any entrained colder cloud material into the lower halo (McKee & Ostriker 1977; McCray & Kafatos 1987; Mac Low & McCray 1988; Norman & Ikeuchi 1989; Hopkins et al. 2012).

However, there is evidence that also stellar winds have a profound impact on the interstellar medium (ISM). Observations and models have long attested to the importance of stellar winds from massive stars as sources of mechanical energy that can have profound influence on the direct environment (Castor et al. 1975; Weaver et al. 1977; Wareing et al. 2018; Pabst et al. 2019). These stellar winds drive a strong shock into its surroundings, which sweeps up ambient gas into dense shells. At the same time, the reverse shock stops the stellar wind, creating a hot tenuous plasma. Eventually, the swept up shells break open, venting their hot plasma into the surrounding medium.

Feedback by ionization of the gas surrounding a massive star will also lead to the disruption of molecular clouds (Williams & McKee 1997; Freyer et al. 2006; Dale et al. 2013). Initially, this will be through the more or less spherical expansion driven by the high pressure of the H II region (Spitzer 1978). Once this expanding bubble breaks open into the surrounding low density material, a champagne flow will be set up (Bedijn & Tenorio-Tagle 1981), rapidly removing material from the molecular cloud.

Bubble structures are ubiquitous in the ISM (e.g., Churchwell et al. 2006). Most of these are caused by expanding H II regions (Walch et al. 2013; Ochsendorf et al. 2014). However, due to the clumpy structure of the ISM, and molecular clouds in particular, the expanding shock fronts can be highly irregular. Depending on the morphology of the shock front, that is the formation ofa shell and more or less massive clumps, massive star formation can be either triggered or hindered (e.g., Walch et al. 2012).

The relative importance of these three different feedback processes (SNe, stellar winds, thermal expanding H II regions) is controversial. Given the peculiar motion of massive stars, SN explosions may occur at relatively large distances from their natal clouds and therefore have little effect on the cloud as the energy is expended in rejuvenating plasma from previous SN explosions, sweeping up tenuous intercloud material and transporting it to the superbubble walls (McKee & Ostriker 1977; Mac Low & McCray 1988; Norman & Ikeuchi 1989; Ochsendorf et al. 2015). While the mechanical energy output in terms of stellar winds is much less than from SN explosions, they will act directly on the natal cloud.

The Orion molecular cloud is the nearest region of massive star formation and has been studied in a wide range of wavelengths. The dense molecular cloud is arranged in an integral-shaped filament (ISF), which has fragmented into four cores, OMC1, 2, 3, and 4, that are all active sites of star formation. At the front side of the most massive core, OMC1, the massive O7V star θ¹ Ori C has ionized the well-known H II region M 42 (NGC 1976), that is known as the Orion Nebula. The Orion Nebula, appearing as the middle “star” in the sword of Orion at a distance of d ≃ 414 pc from us (Menten et al. 2007), is one of the most picturesque structures of our universe. The strong stellar wind from this star has created a bubble filled with hot plasma that is rapidly expanding into the lower-density gas located at the front of the core (Güdel et al. 2008; Pabst et al. 2019). In addition, two slightly less massive stars, B0.5V-type NU Ori and B1V-type 42 Orionis, have created their own ionized gas bubbles, M 43 and NGC 1977. These stars are not expected to have strong stellar winds and likely the thermal pressure of the ionized gas dominates the expansion of these two regions.

The current picture is as follows¹: the Trapezium cluster including θ¹ Ori C, the most massive stars in the Orion Nebula complex, is situated at the surface of its natal molecular cloud, the OMC1 region of the Orion A molecular cloud. Presumably, the Trapezium stars live in a valley of the molecular cloud, having swept away the covering cloud layers (see O’Dell et al. 2009 for a detailed discussion of the structure of the inner Orion Nebula). In the immediate environment of the Trapezium stars, the gas is warm and ionized, constituting a dense H II region, the Huygens Region. It is confined towards the back by the dense molecular cloud core and this is exemplified by the prominent ionization front to the south, the Orion Bar, a narrow and very dense structure Tielens et al. (1993); Goicoechea et al. (2016). To the west of the Trapezium stars, embedded in the molecular gas of OMC1, one finds two dense cores with violent emission features, strong outflows and shocked molecular gas, that are the sites of active intermediate- and high-mass star formation, the Becklin-Neugebauer/Kleinmann-Low (BN/KL) region and Orion South (S).

Further outward, the gas is still irradiated by the central Trapezium stars, granting us the beautiful images of the Orion Nebula. The outward regions were dubbed the Extended Orion Nebula (EON, Güdel et al. 2008). It is, in fact, a closed circular structure of about 5 pc in diameter surrounding the bright Huygens Region, the latter being offset from its center. In the line of sight near the Trapezium stars in the Huygens Region, Orion’s Veil (O’Dell & Yusef-Zadeh 2000) is observed as a layer of (atomic and molecular) foreground gas at 1− 3 pc distance from the blister of ionized gas that forms the bright Orion Nebula. H I and optical/UV absorption lines reveal multiple velocity components, some with velocities that link them to the expanding bubble part of the Veil² seen across the EON (Van der Werf et al. 2013; O’Dell 2018; Abel et al. 2019).

In the vicinity of the Orion Nebula, we find De Mairan’s Nebula (NGC 1982 or M 43) just to the north of the Huygens Region, and the Running-Man Nebula (NGC 1973, 1975, and 1977) further up north. M 43 hosts the central B0.5V star NU Ori and is shielded from ionizing radiation from the Trapezium cluster by the Dark Bay (O’Dell & Harris 2010; Simón-Díaz et al. 2011). NGC 1973, 1975, and 1977 possess as brightest star B1V star 42 Orionis in NGC 1977, which is the main ionizing source of that region (Peterson & Megeath 2008). The structure surrounding NGC 1973, 1975, and 1977 is dominated by the dynamics induced by 42 Orionis, hence, we will denote it by NGC 1977 only in the following.

The hot plasma filling stellar wind bubbles and the ionized gas of H II regions are both largely transparent for far-UV radiation. These non-ionizing photons are instead absorbed in the surrounding neutral gas, creating a warm layer of gas, the photodissociation region (PDR), which cools through atomic fine-structure lines (Hollenbach & Tielens 1999). The [C II] 158 μm fine-structure line of ionized carbon is the dominant far-infrared (FIR) cooling line of warm, intermediate density gas (T ~50−300 K, n ~ 10³−10⁴ cm⁻³). It can carry up to 2% of the total FIR emission of the ISM, most of the FIR intensity arising from re-radiation of UV photons by interstellar dust grains. Velocity-resolved line observations provide a unique tool for the study of gas dynamics and kinematics. Velocity-resolved [C II] and [¹³ C II] observations towards the Huygens Region/OMC1, an area of about 7.5′ × 11.5′, were obtained and analyzed by Goicoechea et al. (2015). Here, we use a large-scale, 60′ × 80′, study of velocity-resolved [C II] emission from the Orion Nebula M 42, M 43 and NGC 1973, 1975 and 1977³.

This paper is organized as follows. In Sect. 2, we review the observations that we used in the present study. In Sect. 3, we discuss the morphology of the shells associated with M 42, M 43, and NGC 1977 and derive gas masses and expansion velocities. Section 4 contains a discussion of the results presented in Sect. 3. We compare theobserved shell kinematics with analytical models. We conclude with a summary of our results in Sect. 5.

2 Observations

2.1 [C II] observations

Velocity-resolved [C II] line observations towards the Orion Nebula complex, covering M 42, M 43 and NGC 1973, 1975, and 1977, were obtained during 13 flights in November 2016 and February 2017 using the 14-pixel high-spectral-resolution heterodyne array of the German Receiver for Astronomy at Terahertz Frequencies (upGREAT⁴, Risacher et al. 2016) onboard the Stratospheric Observatory for Infrared Astronomy (SOFIA). We produced a fully-sampled map of a 1.2 square-degree-sized area at a angular resolution of 16′′. The full map region was observed in the array on-the-fly (OTF) mode in 78 square tiles, each 435.6′′ wide. Each tile consists of 84 scan lines separated by 5.2′′, covered once in both x and y direction. Each scan line consists of 84 dumps of 0.3 s, resulting in a root-mean-square noise of T_mb ≃ 1.14 K per pixel at a spectral resolution of 0.3 km s⁻¹. The original data at a native spectral resolution of 0.04 km s⁻¹ were rebinned to 0.3 km s⁻¹ channels to increase the signal-to-noise ratio (S/N). Ninety per cent of the total 2.2 million spectra required no post-processing, while most problematic spectra could be recovered using a spline baselining approach. A catalogue of splines is generated from data containing no astronomical signal. These splines could then be scaled to the astronomical data and more effectively remove the baselines than a polynomial fit. For a detailed description of the observing strategy and data reductions steps see Higgins et al. (in prep.).

2.2 CO 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) presented the on-the-fly (OTF) mapping and data reduction strategies. In order to cover the same areas mapped by us in the [C II] line, we started to enlarge these CO maps using the new EMIR receiver and FFTS backends. These fully-sampled maps are part of the Large Program “Dynamic and Radiative Feedback of Massive Stars” (see the observing strategy and calibration in Goicoechea et al. 2020).

Goicoechea et al. (2020) provide details on how the old HERA and new EMIR CO maps were merged. Line intensities were converted from antenna temperature () scale to main-beam temperature (T_mb) using appropriate beam and forward efficiencies. Finally, and in order to properly compare with the velocity-resolved [C II] maps, we smoothed the CO(2–1) data to an angular resolution of 16′′. The typical root-mean-square noise level in the CO(2–1) map is 0.16 K in 0.4 km s⁻¹ velocity resolution channels.

2.3 Dust maps

In order to estimate the mass of the [C II]-traced gas in the expanding shells, we make use of FIR photometry obtained by Herschel. Lombardi et al. (2014) present a dust spectral energy distribution (SED) fit from Herschel/Photoconductor Array Camera and Spectrometer (PACS, Poglitsch et al. 2010) and Spectral and Photometric Imaging Receiver (SPIRE, Griffin et al. 2010) photometricimages across a vast region of the Orion molecular cloud. However, they comment that the short-wavelength PACS 70 μm might be optically thick towards the dense molecular cores, hence they exclude it. Since we expect the expanding shell, we are mostly interested in, to mainly consist of optically thin, warm dust, we use the shorter-wavelength bands of PACS at 70, 100, and 160 μm, tracing the warm dust; the longer-wavelength SPIRE bands at 250, 350, and 500 μm, also included in the fit, are dominated by emission from cold (background) dust. We let the dust temperature T_d and the dust optical depth τ₁₆₀ be free parameters, and fit SEDs using a modified blackbody for fixed grain emissivity index β: (1)

We convolve and re-grid the PACS and SPIRE maps to the spatial resolution of the SPIRE 500 μm image, that is 36′′, at a pixel size of 14′′. SED fits to individual pixels are shown in Fig. A.1. We note that β = 0 results in the least residual in the PACS bands, but according to Hollenbach et al. (1991), at least β = 1 should be used. Often (e.g., in Goicoechea et al. (2015) for OMC1) β = 2 is adopted. The resulting dust temperature and dust optical depth vary considerably with β. Decreasing β from 2 to 1, decreases the dust optical depth by half, as does decreasing β from 1 to 0. The longer-wavelength SPIRE bands can only be fitted with β = 1−2. Since we use the conversion factors from τ to gas mass of Li & Draine (2001), we revert to β = 2 as suggested by their models. Lombardi et al. (2014) employ the β map obtained by Planck at 5′ resolution, which amounts to β ≃ 1.6. We note that τ₁₆₀ is biased towards the warm dust, which is beneficial for our purposes. We could have constrained the SED fits to the three PACS bands, using β = 2. This leaves the dust temperature and the dust optical depth in the shells mostly unchanged. The only exception is the southern part of the limb-brightened Veil shell, where the dust optical depth turns out to be 30% lower. In addition, the dust optical depth towards the molecular background is reduced by half. From the dust optical depth τ₁₆₀, we can compute the gas column density: (2)

where we have used a gas-to-dust mass ratio of 100 and assumed a theoretical absorption coefficient κ₁₆₀ ≃ 10.5 cm² g⁻¹ appropriate for R_V = 5.5 (Weingartner & Draine 2001).

To estimate the contribution of emission from very small grains (VSGs) in the 70 μm band, we compare the latter with the Spitzer/Multiband Imaging Photometer (MIPS) 24 μm image. VSGs in PDRs are stochastically heated and obtain temperatures that are higher than that of larger grains. In the shells of the Orion Veil, M 43, and NGC 1977, their contribution is small. In the bubble interiors, however, the emission from hot dust tends to dominate: dust in these H II region becomes very warm due to absorption of ionizing photons and resonantly trapped Lyα photons and radiates predominantly at 24 μm (Salgado et al. 2016). For dust in the PDR, the observed flux at wavelength shorter than 24 μm is due to emission by fluctuating grains and PAH molecules and we defer this analysis to a future study.

Fig. 1 Overview of the Orion Nebula (M 42), M 43, and NGC 1973, 1975, and 1977 in different wavelenghts. Left: DSS2 Hα emission (ESO Archive). Center: [C II] line-integrated emission (Pabst et al. 2019). Right: line-integrated 12CO(2–1) emission (Berné et al. 2014; Goicoechea et al. 2020). Hα emission stems from the ionized gas (T ~ 104 K), the [C II] lineis emitted by mostly neutral gas (T ~ 100 K), whereas CO traces the molecular gas (T ~ 30 K). M 42 and M 43 have a substantial amount of molecular gas in the background, while NGC 1977 seems to be devoid of it (the coverage of the CO map is not sufficient, however, see discussion of the expansion characteristics below). All three regions comprise ionized gas within their neutral limb-brightened shells.

2.4 Hα observations

We make use of three different Hα observations: the Very Large Telescope (VLT)/Multi Unit Spectroscopic Explorer (MUSE) image taken of the Huygens Region (Weilbacher et al. 2015), the image taken by the Wide Field Imager (WFI) on the European Southern Observatory (ESO) telescope at La Silla of the surrounding EON (Da Rio et al. 2009), and a ESO/Digitized Sky Survey 2 (DSS-2) image (red band) covering the entire area observed in [C II]. The DSS-2 image is saturated in the inner EON, which is basically the coverage of the WFI image. We use the MUSE image to calibrate the WFI image in units of MJy sr⁻¹ with a log fit of the correlation. The units of the MUSE observations are given as 10⁻¹² erg s⁻¹ cm⁻³, which we convert to MJy sr⁻¹ with the central wavelength λ ≃ 656.3 nm and a pixel size of 0.2′′. We in turn use the thus referenced WFI image to calibrate the DSS-2 image with a fit of the form y = a(1 − exp(−bx)), where x, y are the two intensities scaled by the fit parameters. To calculate the surface brightness from the spectral brightness we use a , that is I = I_λΔλ = I_νΔν.

We only use the thus calibrated DSS-2 image in NGC 1977 for quantitative analysis. The Hα surface brightness in the EON is largely due to scattered light from the bright Huygens Region (O’Dell & Harris 2010). Also Hα emission in M 43 has to be corrected for a contribution of scattered light from the Huygens Region (Simón-Díaz et al. 2011).

We correct for extinction towards NGC 1977, using R_V = 5.5, suitable for Orion. The reddening is E(B−V) = 0.08 (Knyazeva & Kharitonov 1998). For M 43, the reddening is E(B−V) = 0.49 (Megier et al. 2005). Extinction towards M 42 was studied by O’Dell & Yusef-Zadeh (2000), for example, and more recently by Weilbacher et al. (2015).

3 Analysis

3.1 Global morphology

From the comparison of the three gas tracers in Fig. 1, we see that [C II] emission stems from different structures than either Hα, tracing the ionized gas, or ¹²CO(2–1), tracing the cold molecular gas. In fact, from Fig. 1 in Pabst et al. (2019) we conclude that it stems from the same regions as the warm-dust emission and PAH emission. We can quantify this by correlation plots of the respective tracers (Pabst et al., in prep.). We clearly see bubble-like structures in [C II], that are limb-brightened towards the edges. The southern red circle in Fig. 2 indicates the outlines of the bubble of M 42. We note that the Trapezium stars, causing the bubble structure, are offset from the bubble center; in fact, they seem to be located at the northern edge of the bubble. We estimate a bubble radius of r ≃ 2.7 pc. However, when taking into account the offset of the Trapezium stars, the gas at the southern edge of the bubble is some 4 pc from the stars. The geometry of the bubble indicates that there is a significant density gradient from the north, where the Trapezium stars are located, to the south, the ambient gas there being much more dilute.

In M 43, the star NU Ori is located at the center of the surrounding bubble. 42 Orionis in NGC 1977 is slightly offset from the respective bubble center (cf. Fig. 2). Both bubbles are filled with ionized gas as observed in Hα emission (cf. Fig. 1) as well as radio emission (Subrahmanyan et al. 2001). In the main part of this section, we derive the mass, physical conditions and the expansion velocity of the expanding shells associated with M 42, M 43, and NGC 1977.

3.2 The expanding Veil shell – M 42

3.2.1 Geometry, mass and physical conditions

The prominent shell structure of the Orion Nebula is surprisingly symmetric, although it is offset from the Trapezium cluster. From the pvdiagrams (cf. Sect. 3.2.2 and Appendix C), we estimate its geometric center at (Δα, Δδ) =(−520′′,−550′′) and its radius with r =1360′′ = 0.38°, which corresponds to r ≃ 2.7 pc at the distance of the Orion Nebula; the geometric center of the bubble has a projected distance of 1.5 pc from the Trapezium cluster. Hence, the gas in the southern edge of the bubble is about 4 pc away from the Trapezium stars, whereas to the north the bubble outline is more ellipsoid and the gas in the shell is at only 0.5 pc distance.

The limb-brightened shell of the expanding Veil bubble is mainly seen in the velocity range (with respect to the Local Standard of Rest (LSR)) v_LSR = 5−8 km s⁻¹. The emission from the bubble itself can be found down to v_LSR ≃−7 km s⁻¹ (cf. Sect. 3.2.2). The bright main component, that originates from the surface of the background molecular cloud, lies at v_LSR ≃ 8 km s⁻¹; in the bright Huygens Region significant emission extends up to v_LSR ~ 15 km s⁻¹. Here, we also detect the [¹³C II] F = 2−1 line in individual pixels, corresponding to gas moving at v_LSR ~ 8 km s⁻¹ (see previous detections of [¹³C II] lines in Boreiko & Betz 1996; Ossenkopf et al. 2013; Goicoechea et al. 2015).

Figure 3 shows the average [C II] spectrum towards the Veil shell without OMC1 and the ISF. The [¹³C II] F = 2 − 1 line, one of the three [¹³C II] fine-structure lines, that is shifted by 11.2 km s⁻¹ with respect to the [¹²C II] line, is marginally detected. We estimate the detection significance over the integrated line from the fit errors at 5σ. From the [¹³C II] F = 2 − 1 line, we can compute the [C II] optical depth τ_[C II] and the excitation temperature T_ex:

where is the isotopic ratio for Orion (Langer & Penzias 1990) and 0.625 is the relative strength of the [¹³C II] F = 2 − 1 line (Ossenkopf et al. 2013); is the continuum brightness temperature of the dust background (Goicoechea et al. 2015). In general, in the limb-brightened shells. The measured peak temperatures yield τ_[C II] ≃ 0.1 and T_ex ≃ 144 K for the averaged spectrum. We obtain an average C⁺ column density of . These results have to be interpreted with caution, since we average over a large region with very different conditions (G₀, n) and the signal is likely dominated by the bright eastern Rim of the Veil shell (see Appendix B for a discussion thereof). Also, baseline removal for the [¹³C II] line is problematic. From the dust optical depth, we estimate a much larger column density towards the Veil shell, .

If we assume the [C II] emission from the limb-brightened shell to be (marginally) optical thick, that is τ_[C II] ≳ 1, we can alsoestimate the excitation temperature from the peak temperature of the [¹²C II] spectra by Eq. (4). We obtain excitation temperatures of T_ex ≃ 100 K in the Eastern Rim (and similar gas temperatures), and T_ex ≃ 50 K in the far edge of the shell. Here, the density presumably is somewhat lower, resulting in a higher gas temperature. The gas (or kinetic) temperature is given by: (5)

where n_cr ≃ 3 × 10³ cm⁻³ is the critical density for C⁺-H collisions (Goldsmith et al. 2012; Pabst et al. 2017). With the density estimates below, the gas temperature in the southern Veil shell is T_gas ≃ 70−125 K.

From the dust SEDs, Pabst et al. (2019) estimate for the mass of the gas in the Veil shell after accounting for projection effects. This is in good agreement with their mass estimate from the [¹³C II] line in the brightest parts of the shell. The bright eastern part of the shell, however, might not be representative of the entire limb-brightened shell of the Veil shell. From a line cut through this Eastern Rim (Fig. B.4), we estimate a density of n ≃ 9 × 10³ cm⁻³ from the spatial separation of the peaks in [C II] and CO emission (see Appendix B), assuming an A_V of 2 for the C⁺ /C/CO transition and N_H∕A_V ≃ 2 × 10²¹ cm⁻² (Bohlin et al. 1978). Although continuously connected to the limb-brightened Veil shell in space and velocity, this Eastern Rim might be a carved-out structure in the background molecular cloud, as has been argued from the presence of foreground scattered light from the Trapezium stars (O’Dell & Goss 2009). We note that the [C II] emission from the expanding shell connected to the Eastern Rim is confined to within the Rim, supporting the view that the Eastern Rim is a static confinement of the expanding bubble.

Morphologically, the Eastern Rim fits well in with the rest of the Veil shell, following its curvature well and forming one, seemingly coherent structure. In contrast, there is no connection to the ISF that spans up the dense core of the Orion molecular cloud in the submillimeter maps (cf. Fig. 1). It is then hard to conceive that this part of the nebula is a fortuitous coincidence of the expanding Veil shell encountering a background or foreground molecular cloud structure. On the other hand, the channel maps do not reveal the shift of the shell with increasing velocity in the Eastern Rim that are the signature of an expanding shell in the other parts of the Veil shell (cf. Fig. 7). Excluding the area associated directly with the dense PDR around the Trapezium(indicated by the circle in Extended Data Fig. 8 in Pabst et al. 2019), analysis of the dust SED yields M ≃ 2600 M_⊙ for the mass of the Veil shell. If we also exclude the Eastern Rim, this drops to M ≃ 1500 M_⊙.

From the visual extinction towards the Trapezium cluster, A_V ≃ 1.8 (O’Dell & Goss 2009), corresponding to a hydrogen column density of N_H ≃ 3.6 × 10²¹ cm⁻², the mass of the half-shell can be estimated as M ≃ 2πr²N_Hμm_H, where μ = 1.4 is the mean molecular weight. With r ≃ 2.5 pc, this gives M ≃ 1600 M_⊙. From H₂ and carbon line observations towards the Trapezium cluster, Abel et al. (2016) estimate a density of n ≃2.5 × 10³ cm⁻³ and a thickness of the large-scale Veil component B of d ≃ 0.4 pc, resulting in a visual extinction of A_V ≃ 1.5. In a new study of optical emission and absorption lines plus [C II] and H I line observations, Abel et al. (2019) determine n ≃1.6 × 10³ cm⁻³ for the density in the Veil towards the Trapezium cluster, their Component III(B). However, the central part of the Veil might not be representative of the entire Veil shell, since the [C II] emission from the Veil in the vicinity of the Trapezium stars is less distinct than in the farther parts of the Veil shell, suggesting somewhat lower column density in the central part. The visual extinction, as seen by MUSE (Weilbacher et al. 2015), varies across the central Orion Nebula and decreases to only A_V ≃ 0.7 to the west of the Trapezium cluster. Using this value reduces the mass estimate of the Veil shell to M ~ 600 M_⊙.

X-ray observations of the hot plasma within the EON indicate varying extinction due to the foreground Veil shell. Güdel et al. (2008) derive an absorbing hydrogen column density of N_H ≃ 4 × 10²⁰ cm⁻² towards the northern EON and N_H ≤ 10²⁰ cm⁻² towards the southern X-ray emitting region. Assuming again T_ex ≃ 50 K, these column densities correspond to [C II] peak temperatures of T_P ≃ 2 K and T_P ≃ 0.5 K, respectively,which would be below the noise level of our observations. Yet, we do detect the [C II] Veil shell towards these regions (cf. Figs. C.1 and C.2). From the lack of observable X-rays towards the eastern EON and the more prominent [C II] Veil shell in this region, we suspect that the absorbing column of neutral gas is higher in the eastern Veil shell.

With the MUSE estimate for the column density, we can estimate the density and check for consistency with our [C II] observations. With a shell thickness of d ≃ 0.3 pc and a visual extinction ofA_V ≃ 0.7, we estimate a density in the Veil shell of n ≃ 1.5 × 10³ cm⁻³. We will use this density estimate in the following, but we note that we consider this a lower limit as the density estimate from a typical dust optical depth, τ₁₆₀ ≃ 2 × 10⁻³, in the southern Veil shell and a typical line of sight, l ≃ r∕2, gives n ≃ 4 × 10³ cm⁻³. However, from the lack of detectable wide-spread CO emission from the Veil shell, also in the southern limb-brightened shell, we can set an upper limit for the visual extinction of A_V ~ 2, equivalent to a gas column of N_H ~ 4 × 10²¹ cm⁻² (Goicoechea et al. 2020).

With a shell thickness of d ≃ 0.3 pc and a density of n ≃ 1.5 × 10³ cm⁻³, the shell has a radial C⁺ column density of . Adopting this latter value and an excitation temperature of T_ex = 50 K, the calculated [C II] optical depth, given by: (6)

in the shell seen face-on is τ_[C II] ≃ 0.4, with a line width of . The expected peak temperature then, (7)

is T_P ≃ 6 K, which is in reasonably good agreement with the peak temperatures of the narrow expanding components in spectra in Fig. 6.

We can also calculate the mass of the [C II]-emitting gas from the luminosity of the shell, L_[C II] ≃ 170 L_⊙ in the limb-brightened Veil shell, assuming that the line is effectively optically thin. By (8)

where is the carbon gas-phase abundance (Sofia et al. 2004), A ≃ 2.3 × 10⁻⁶ s⁻¹ the Einstein coefficient for spontaneous emission, ΔE∕k_B ≃ 91.2 K the level separation of the two levels with statistical weights g_l = 2 and g_u = 4, we obtain M ≃ 680 M_⊙, assuming an excitation temperature of T_ex = 50 K. As we expect the line to be (marginally) optically thick in the limb-brightened shell, τ_[C II] ≃ 1−2, this value constitutes a lower limit of the shell mass. Correcting with this [C II] optical depth would yield a mass estimate of at least M ≃ 1100 M_⊙. For further analysis we will use the mass estimate from the dust optical depth, where we have subtracted the mass of the Eastern Rim, M ≃ 1500 M_⊙. This mass estimate is robust when considering only the PACS bands in the SED, reducing contamination by cold background gas. We recognize that the column density of the shell is highly variable and is lowest in the foreground shell. This is as expected given the strong pressure gradient from the molecular cloud surface towards us (cf. Sect. 4.1).

We judge that the large-scale arc-structured [C II] line emission mainly stems from the outer dense neutral shell indeed, rather than from the contained ionized gas. Evidence for this is the rather small line width where there is a significant signal, which is comparable to the line width of the main component from the neutral surface layer of the molecular cloud (Δv_FWHM ≃ 3−4 km s⁻¹, cf. Fig. 6). If the [C II] line originated from the hot (T ~ 10⁴ K) ionized gas, we would expect it to be broader, Δv_FWHM ~ 10 km s⁻¹, as observed in the H II region bordering the Horsehead Nebula (Pabst et al. 2017). Where the signal from the expanding shell is weak, we indeed observe such broad lines, hence part of the shell might still be ionized (see Sect. 4.4 for a discussion of the ionization structure of wind-blown bubbles). The limb-brightened shell is also visible in Hα, indicating that the surface is ionized. This is also illustrated by the line cut in Fig. 4, where the onset of the shell is marked by a steep increase in all tracers except CO. The relative lack of CO indicates low extinction of FUV photons within the southern shell and gives an upper limit on the density in this part, n ≲3 × 10³cm⁻³. From the correlation with other (surface) tracers ([C II], 8 μm, 160 μm) along the line cut we conclude that the shell has a corrugated surface structure, as there are multiple emission peaks in these tracers behind the onset of the shell.

The broader lines (Δv_FWHM > 5 km s⁻¹) we observe where the shell arc is barely visible, agree with an estimation of the expected [C II] emission from ionized gas in the shell. Here, we have (Pabst et al. 2017): (9)

which is, with an electron density of n_e ~ 50−100 cm⁻³ (O’Dell & Harris 2010) and a line of sight of l ~ 0.3 pc, what we find in the spectra from the pv diagram in Fig. 6. We do not expect much [C II] emission from the inner region of the EON, since the dominant ionization stage of carbon in the vicinity of a O7 star is C²⁺.

Fig. 2 [C II] line-integrated intensity from the Orion Nebula (M 42), M 43, and NGC 1977. The red rectangles 1 and 2 indicate the positions of the pv diagrams shown in Figs. 5 (horizontal) and 8 (vertical). Rectangle 3 indicates the position of the pv diagram shown in Fig. 17. The red circles delineate the approximate extent of the expanding bubbles of gas that are the Veil shell, M 43, and NGC 1977, and within which the respective masses and luminosities are computed. The light-red dashed rectangle circumscribes the area of the eastern bright arm shown in Figs. B.1 and B.2. The black line indicates the position of the line cut in Fig. 4. The stars mark the position of θ1 Ori C (yellow), θ2 Ori A (orange), NU Ori (pink), and 42 Orionis (purple). Unless noted otherwise, all coordinate offsets are given with respect to the position of θ1 Ori C, (α, δ) = (5h35′ 16.46′′, − 5°23′22.85′′).

Fig. 3 Average [C II] spectrum towards the Veil shell without OMC1 and ISF. The inlaid panel shows the residual of the spectrum, in the systemic velocity of the [13C II] F = 2 − 1 line, after subtracting the [12C II] fit. The [12C II] line can be fitted by two Gaussians with TP ≃ 11.2 ± 0.3 K, vP ≃ 8.8 ± 0.1 km s−1, ΔvFWHM ≃ 5.0 ± 0.1 km s−1 and TP≃ 2.3 ± 0.1 K, vP ≃ 4.0 ± 0.4 km s−1, ΔvFWHM ≃ 9.3 ± 0.5 km s−1; the [13C II] F = 2 − 1 component is fitted by a Gaussian with TP ≃ 0.11 ± 0.02 K, vP ≃ 8.2 ± 0.3 km s−1, ΔvFWHM ≃ 2.1 ± 0.6 km s−1. The main line could also be fitted with three Gaussian components, but we compare the respective combined peak temperatures of the [12C II] and [13C II] line for a first-order estimate of the excitation conditions.

Fig. 4 Line cut through the southwest of the Veil shell, indicated in Fig. 2. The onset of the shell is marked by a steep increase in all tracers, except CO, at Δδ ≃ −1400′′. With a line angle of 30° with respect to the vertical axis, this is at a distance of 3.25 pc from the Trapezium stars.

3.2.2 Expansion velocity

Persistent arc structures observed in position-velocity (pv) diagrams are the signature of bubbles that form in the ISM. By the aid of such diagrams one can estimate the expansion velocity of the associated dense shell. The pv diagrams for the Veil shell reveal the presence of a half shell. Figure 5 shows one such pv diagram as an example. Other cuts are shown in Appendix C. Pabst et al. (2019) estimate an expansion velocity of v_exp ≃ 13 km s⁻¹ for this half shell. This value is consistent with Gaussian fits to single (averaged) spectra taken from the pv diagram where this is possible (Fig. 6).

It is important to note that the bubble shell traced by [C II] emission expands in one direction only, which is towards us. Towards theback, the gas is confined by the background molecular cloud hindering rapid expansion. The prominent 8 km s⁻¹ emission originates from this dense background gas.

An alternative way to estimate the expansion velocity of the bubble is from channel maps of this region. With increasing v_LSR, the limb-brightened shell filament of M 42 is observed to move outward, away from the Trapezium stars, as can be seen from Fig. 7. The projected geometry of the filaments in different channels provides an estimate of v_exp. However, this only works well in the bright eastern arm of the shell, the so-called Rim, which is comparatively narrow (cf. Appendix B, Figs. B.1 and B.2). In other parts of the limb-brightened shell, the south and west, the filaments in different velocities do not line up consistently. To the north of the Trapezium stars, we do not see outward expansion at all. From these alternative methods (fitting Gaussians to single spectra, outlines of the limb-brightened shell) we find values agreeing with those computed above from the pv diagram, v_exp ≃ 10−15 km s⁻¹. Also the majority of pv diagrams shown in Figs. C.1 and C.2 are consistent with an overall expansion velocity of v_exp ≃ 13 km s⁻¹.

Figure 8 compares the velocity structures of [C II], ¹²CO(2–1), ¹³CO(2–1), and H I along a vertical cut at Δα =167′′ − 242′′. This figure illustrates the complex dynamics in the region of the Orion Nebula and contrasts the various kinematic structures in the different tracers. It gives an impression of the variety of emission structures and the potential synergies between the four data sets. The cut shows arc structures close to the central region (Huygens Region) that coincide with structures detected in H I and discussed in Van der Werf et al. (2013). Also, the shell surrounding M 43 is distinctly visible (cf. Fig. 12). We note that to the north the bubble-like arc of the Veil shell is disrupted by apparently violent gas dynamics in that region, even far to the west of the Huygens Region. Closer examination of pv diagrams in this region reveal multiple arcs indicative of bubble structures (see Appendix C for all pv diagrams through the Veil shell). South of OMC1, the observed difference in velocity (~ 1 km s⁻¹) between the [C II] background and the molecular cloud as traced by CO is typical for the advection flow through the PDR (Tielens 2010, Chap. 12). From our analysis, we conclude that the systematic uncertainty of the mass estimate are of the order of 50%, systematic uncertainties of the estimates of the extent of the shell and of the expansion velocity are about 30%.

Fig. 5 [C II] pv diagram along horizontal position cut 2 indicated in Fig. 2 (Δδ = −907′′ − (−831′′)). The lower panel shows the same cut with the arc structure for an expansion velocity of 13 km s−1 on a background velocity of 8 km s−1 (red dashed lines).

Fig. 6 [C II] spectra in the Veil shell taken from pv diagram in Fig. 5 with Gaussian fits. Spectra at Δα = −362′′ to Δα = −1268′′ are averaged over 151′′ × 75.5′′ in order toimprove the signal-to-noise ratio of the faint shell component to be fitted, others are averaged over 75.5′′ × 75.5′′. Coordinate offsets indicate the lower left-hand corner of the rectangle over which the spectra are averaged; Δδ = −907′′. The fit parameters are given in Table 1.

Fig. 7 Three-color image of velocity channels towards M 42 and M 43. Blue is the velocity channel vLSR =4−5 km s−1, green vLSR = 6−7 km s−1, and red vLSR = 9−10 km s−1. With increasing vLSR, the limb-brightened Veil shell filament is displaced outward, away from the bubble center. The gas of M 43 that is expanding towards us can be observed in the blue channel; the limb-brightened shell of M 43 has higher vLSR.

Fig. 8 PV diagrams along vertical position cut 1 indicated in Fig. 2 (Δα =167′′ − 242′′) in [C II], 12CO(2–1), 13CO(2–1) (Berné et al. 2014; Goicoechea et al. 2020), and H I (Van der Werf et al. 2013).

3.3 M 43

3.3.1 Geometry, mass and physical conditions

The nebula M 43 hosts the central star NU Ori, which is located approximately at the geometrical center of the limb-brightened shell. The shell radius is r ≃150′′ ≃ 0.3 pc with a thickness of Δr ≃ 0.05 pc. Figure 9 shows a three-color image of M 43. Hα stems from the center of M 43, constrained by a thin shell of [C II] and CO in the east and north. To the south lies the massive bulk of OMC1. The H II region of M 43 appears as a region of higher dust temperature, the shell is distinctly visible in [C II], IRAC 8 μm emission, FIR emission from warm dust (cf. Fig. 1 in Pabst et al. 2019), and CO. The limb-brightened shell reveals the structure of a PDR with 8 μm PAH emission, [C II] emission, and 160 μm warm dust emission originating at the inner, illuminated side while the CO emission originate from deeper within the shell (Fig. 10). M 43 is located just eastwards of the ridge of the molecular ISF, whose illuminated, [C II] emitting outliers form the background against which the half-shell expands.

From our dust SEDs, we derive a mass of M ≃ 210 M_⊙. However, this mass is likely contaminated by the mass of the molecular background, whose surface is visible in [C II] (cf. Fig. 8). Using only the PACS bands in the SED, reduces the mass estimate to M ≃ 110 M_⊙. From the [C II] luminosity, L_[C II] ≃ 24 L_⊙ and assuming T_ex ≃ 90 K (see below), the mass is estimated to be M ≃ 55 M_⊙.

We note that the limb-brightened shell of M 43 consists of two segments. The eastern shell lies closer to the central star and therefore is probably denser. We can derive the density of the shell by the distance between the [C II] and the CO peak in a line cut through that region (Fig. 10a), Δd ≃20′′ ≃ 0.04 pc between the [C II] peak at Δα ≃310′′ and the CO peak at Δα ≃330′′. Assuming a typical PDR structure with ΔA_V ≃ 2 between the PDR front, traced by [C II], and the CO peak, we arrive at a density of n ≃3 × 10⁴ cm⁻³. In the northern shell (Fig. 10b) we derive a density of n ≃8 × 10³ cm⁻³ with Δd ≃80′′ ≃ 0.16 pc between the [C II] peak at Δδ ≃560′′ and the CO peak at Δδ ≃640′′. With the shell extent above, that is , and the latter density estimate, we compute a mass of the expanding hemisphere of M ≃ 7 M_⊙, which is significantly lower than the mass derived from the dust optical depth and the [C II] luminosity. Given the potential contamination of these mass estimates by background molecular cloud material, we elected to go with the latter estimate for the gas mass.

When averaging over the eastern limb-brightened shell of M 43, we detect the [¹³C II] F = 2 − 1 line at a 4σ level, as shown in Fig. 11. From this⁵ we obtain an average and an excitation temperature of . This translates into a C⁺ column density of along the line of sight. Assuming that most of the material is located in the limb-brightened shell, we can estimate a density in the shell with the assumed length of the line of sight of l ~ r∕2. With r ≃ 0.3 pc, this gives n ≃ 3 × 10⁴ cm⁻³, which is about the same as the previous estimate from the [C II]-CO peak separation. With this density we compute a gas temperature of T_gas ≃ 90 K.

If we assume the [C II] emission from the limb-brightened shell to be (marginally) optical thick, we can estimate the excitation temperature from the [C II] peak temperature. This gives an average of T_ex ≃ 110 K, and hence a slightly higher gas temperaturethan from the [¹³C II] estimate.

Fig. 9 Three-color image of M 43. Red is Hα emission from the ionized gas, green is line-integrated CO(2–1) emission from molecular gas, and blue is the [C II] line-integrated intensity, tracing the neutral gas. The line cuts in Fig. 10 are taken along lines a and b.

Fig. 10 Line cuts along lines a and b, indicated in Fig. 9. The arrows indicate the peaks in [C II] emission (green) and CO emission (blue) that correspond to the shell.

Fig. 11 Average [C II] spectrum towards the eastern shell of M 43. The inlaid panel shows the residual of the spectrum, in the systemic velocity of the [13C II] F = 2 − 1 line, after subtracting the [12C II] fit. The [12C II] line can be fitted by a Gaussian with TP ≃ 46.2 ± 0.3 K, vP ≃ 9.6 ± 0.1 km s−1, ΔvFWHM ≃ 4.6 ± 0.1 km s−1; the [13C II] F = 2 − 1 component is fitted by a Gaussian with TP ≃ 1.1 ± 0.2 K, vP ≃ 8.6 ± 0.3 km s−1, ΔvFWHM ≃ 4.8 ± 0.8 km s−1. We note that a better fit of the [12C II] line can be obtained by fitting two components with TP ≃ 30.9 ± 1.3 K, vP ≃ 9.6 ± 0.1 km s−1, ΔvFWHM ≃ 3.3 ± 0.1 km s−1 and TP≃ 17.5 ± 1.3 K, vP ≃ 9.8 ± 0.1 km s−1, ΔvFWHM ≃ 6.8 ± 0.2 km s−1, but both fits leave significant residuals.

Fig. 12 [C II] pv diagram along vertical position cut 1 indicated in Fig. 2 (Δα =167′′ − 242′′), zoomed in to M 43 (cf. Fig. 8). The right panel shows the same cut with the arc structure for an expansion velocity of 6 km s−1 on a background velocity of 8 km s−1 (red dashed lines).

3.3.2 Expansion velocity

The [C II] pv diagram running through M 43 clearly exhibits a half shell structure (Fig. 12). From the spectra taken towards this region (Fig. 12), we measure an expansion velocity of v_exp ≃ 6 km s⁻¹. The bubble only expands towards us, away from the background molecular cloud. We do not see an expanding CO counterpart (cf. Fig. 8). The pv diagram is consistent with the radius of the shell of r ≃150′′ ≃ 0.3 pc and the expansion velocity of v_exp ≃ 6 km s⁻¹, both values are much less than those derived for the Trapezium wind-blown bubble. Also, in Fig. 12, we observe [C II] emission from within the shell arc, presumably stemming from the expanding ionized gas.

The eastern shell of M 43, that lies closer to the shell center, moves outward at a velocity of about v_exp ≃ 2.5 km s⁻¹ (estimated from channel maps, see above and Appendix B). This corroborates the thesis that this is a denser region compared to other structures in M 43. Again, the systematic uncertainty of the mass estimate is of the order of 50%, systematic uncertainties of the estimates of the extent of the shell and of the expansion velocity are about 30%.

3.4 NGC1977

3.4.1 Geometry, mass and physical conditions

As for M 42 and M 43, the shell structure of NGC 1977 is almost spherical, but again it is offset from the central star(s). We estimate its geometric center at (Δα, Δδ) =(195′′,155′′) (from 42 Orionis, (Δα, Δδ) = (−100.4′′,1984.7′′) from θ¹ Ori C) and its outer radius with r = 800′′ =0.22°, which corresponds to r ≃ 1.6 pc, assuming a distance of 414 pc (Menten et al. 2007). Hence, the projected distance of the bubble center from 42 Orionis is 0.5 pc. The thickness of the limb-brightened shell is Δr = 7′≃ 0.8 pc, somewhat thicker than the shell of M 42, but more dilute as well. Figure 14 shows a three-color image of NGC 1977. Visible Hα stems from the center of NGC 1977, constituting an H II region, surrounded by a shell of [C II] and constrained by a bulk CO cloud (OMC3) in the southwest.

The limb-brightened shell is seen in [C II] emission in the velocity range v_LSR = 9−15 km s⁻¹. The brightest [C II] component stems from the molecular core OMC3, which hosts a PDR at its surface towards 42 Orionis. The [C II] emission associated with this PDR is analyzed in detail by Kabanovic et al. (in prep.). With the FUV luminosity given by Kim et al. (2016), we compute an incident FUV intensity of G₀ ≃ 100 at the shell surface.

From our τ₁₆₀ map, we estimate the mass of the shell, where we exclude the region of OMC3. We obtain M ≃ 700 M_⊙. With the [C II] luminosity L_[C II] ≃ 140 L_⊙, we derive a mass of the [C II]-emitting gas of M ≃ 540 M_⊙ for T_ex = 50 K, somewhat lower than the mass derived from the dust opacity, but still in good agreement⁶.

Figure 15 shows the average spectrum towards the [C II]-bright shell of NGC 1977, excluding OMC3. Here, the detection of the [¹³C II] F = 2 − 1 line is marginal with σ ≃ 3. As the [C II] optical depth and excitation temperature are very sensitive to the exact [¹³C II] peak temperature, it is hard to derive reliable values with the uncertainties at hand (fit uncertainties and baseline ripples). With the results of a Gaussian fit as given in the caption of Fig. 15, we obtain τ_[C II] ≃ 0.7 ± 0.3 and . The resulting C⁺ column density then is . If we assume a column length of r∕2 ≃ 0.5 pc in the limb-brightened shell, we estimate a density of n ≃ 2.5 × 10³ cm⁻³. This results in an estimate of the gas temperature of T_gas ≃ 180 K, which is somewhat higher than expected for the rather moderate radiation field in NGC 1977.

The mean dust optical depth from NGC 1977 (without OMC3) is τ₁₆₀ ≃ 2.4 × 10⁻³. From this, with r∕2 ≃ 0.5 pc, we estimate a gas density of n ≃ 1 × 10⁴ cm⁻³ in the shell. With the average T_ex ≃ 70 K (from the dust optical depth and the peak temperature of the [C II]), this gives T_gas ≃ 90 K. From the dust optical depth and the [C II] peak temperature, we calculate τ_[C II] ~ 3, so [C II] emission is optically thick. This is in disagreement with the result from the [¹³C II] F = 2 − 1 line. However, the result obtained here seems to be more credible, as it yields a gas temperature that is more in line with a previous [C II] study of the Horsehead Nebula at similar impacting radiation field (Pabst et al. 2017). PDR models predict a surface temperature of ~ 150 K (Kaufman et al. 2006; Pound & Wolfire 2008), but the [C II]-emitting layer is expected to be somewhat cooler.

In the region where the spectrum in Fig. 18 is taken, we have an average τ₁₆₀ ≃ 8.6 × 10⁻⁴, corresponding to a hydrogen column density of N_H ≃ 5 × 10²¹ cm⁻². From Hα emission, we can estimate the electron density in the ionized gas that is contained in the shell. We obtain n_e ≃ 40 cm⁻³. Hence, with a column length of l ~ 2 pc the ionized gas is only a fraction of the total column density. Most of the dust emission within the column stems from the denser expanding shell. With a shell thickness of d ≃ 0.8 pc, appropriate for the limb-brightened shell, we estimate a gas density of n ≃ 2 × 10³ cm⁻³. From the column density⁷ (assuming that all carbon is ionized) and the spectrum we calculate an excitation temperature of T_ex ≃ 37 K and a [C II] optical depth of τ_[C II] ≃ 1. This corresponds to a gas temperature of T_gas ≃ 60 K, which is somewhat lower than in the limb-brightened shell as calculated above. However, a slight overestimation of the dust optical depth of the expanding material would lead to a significantly increased estimate of the gas temperature.

Towards the center of the H II region, we observe a very broad component in the [C II] spectrum (Fig. 16). This component likely stems from the ionized gas. The expected [C II] intensity from ionized gas at a temperature of T ~ 10⁴ K with an electron density of n_e ≃ 40 cm⁻³ and a line of sight of l ≃ 2 pc matches the observed intensity very well. We note that the line is broader than what one would expect from thermal broadening alone. The additional broadening might be due to enhanced turbulence and the expansion movement of the gas.

Fig. 13 [C II] spectra along position cut 1 (cf. Fig. 12) in M 43 with Gaussian fits. Each spectrum is averaged over 75.5′′ × 25.2′′. The expanding shell has a distinguished emission peak. Additionally to the main component and the shell component, we fit a third, broad, component, that possibly stems from the ionized gas within the shell. Coordinate offsets indicate the lower left-hand corner of the rectangle over which the spectra are averaged; Δα =242′′. The fit parameters are given in Table 2.

Fig. 14 Three-color image of NGC 1977. Red is Hα emission from the ionized gas, green is line-integrated CO(2–1) emission from molecular gas (coverage only in OMC3), and blue is the [C II] line-integrated intensity, tracing the neutral gas.

Fig. 15 Average [C II] spectrum towards [C II]-bright shell of NGC 1977 (without OMC3). The inlaid panel shows the (smoothed) residual of the spectrum, in the systemic velocity of the [13C II] F = 2 − 1 line, after subtracting the [12C II] fit. The [12C II] component can be approximated by a Gaussian with TP ≃ 16.5 ± 0.1 K, vP ≃ 12.0 ± 0.1 km s−1, ΔvFWHM ≃ 3.2 ± 0.1 km s−1; the [13C II] F = 2 − 1 component is fitted by a Gaussian with TP ≃ 0.21 ± 0.03 K, vP ≃ 11.8 ± 0.2 km s−1, ΔvFWHM ≃ 2.1 ± 0.5 km s−1.

Fig. 16 [C II] spectrum towards NGC 1977, averaged over circular area with radius 125′′ centered at 42 Orionis ((Δα, Δδ) = (−100.4′′,1984.7′′)), with Gaussian fits. The fit parameters are given in Table 3.

Fig. 17 [C II] pv diagram of NGC 1977 along position cut 3, indicated in Fig. 2 (Δδ =2085′′ − 2259′′). The lower panel shows the same cut with the arc structure for an expansion velocity of ± 1.5 km s−1 (red dashed lines).

Fig. 18 [C II] spectrum along pv diagram in NGC 1977 (Fig. 17) with Gaussian fits, averaged over 201′′ × 201′′. Coordinate offsets indicate the lower left-hand corner of the square over which the spectrum is averaged. The two narrow components stem from the expanding shell, while the broad component originates in the ionized gas that is contained within the shell. Since the bubble expands in both directions, the peak-velocity difference is twice the expansion velocity, vexp ≃ 1.3 km s−1. The fit parameters are given in Table 4.

3.4.2 Expansion velocity

Figure 17 shows a [C II] pv diagram through the center of the bubble associated with NGC 1977. We recognize evidence of expansion but the arc structure is much fainter than in M 42 and M 43 and disrupted. As opposed to M 42 and M 43, there is no background molecular cloud constraining the expanding gas: we see the bubble expanding in two directions.

Since the [C II] emission from NGC 1977 is fainter than that of M 42, we have opted to determine the expansion velocity from spectra towards this region (Fig. 18), resulting in v_exp ≃ 1.3 km s⁻¹. This is consistent with the pv diagram shown in Fig. 17, from which we obtain the bubble radius r ≃ 500′′ ≃ 1.0 pc with the expansion velocity of v_exp ≃ 1.5 km s⁻¹. As for the Veil shell and M 43, the systematic uncertainty of the mass estimate are of the order of 50%, systematic uncertainties of the estimates of the extent of the shell and of the expansion velocity are about 30%.

4 Discussion

OMC1, including the Trapezium cluster, is the nearest site of active intermediate- and high-mass star formation (e.g., Hillenbrand 1997; Megeath et al. 2012, 2016; Da Rio et al. 2014; Großschedl et al. 2019 for discussions of the stellar content). Most of the stars formed and forming are low-mass stars, but also some high-mass stars, that have profound impact on the characteristics and evolution of their environment. Outside the Huygens Region, a significant massive (triple) star is the O9.5IV star θ² Ori A, located just to the southeast of the Orion Bar. The radiation of this star dominates the ionization structure of the gas towards the south-east of the HuygensRegion (O’Dell et al. 2017). It also possesses strong winds. The most dominant star in the Orion Nebula is the O7V star θ¹ Ori C, the most massive Trapezium star. It is itself a binary, with possibly a third companion (Lehmann et al. 2010). While its exact peculiar velocity is controversial, the θ¹ Ori C is plowing away from the molecular cloud, the site of its birth, and will have travelled some 25 pc before it explodes as a supernova (O’Dell et al. 2009; Kraus et al. 2009; Pabst et al. 2019).

While the small, 0.5 pc sized Huygens Region has been extensively studied, studies on the much fainter EON are less numerous. O’Dell & Harris (2010) determine the temperature of the ionized gas within the large EON H II region to be T ~8.3 × 10³ K, while electron densities decrease from n_e ~ 3000 cm⁻³ at the Trapezium stars to about n_e ~ 30 cm⁻³ 20′ away, roughly as d⁻², with d the projected distance. Güdel et al. (2008) report that the EON exhibits significant X-ray emission, emanating from hot (T ~ 2 × 10⁶ K), dilute (n_e ~ 0.1−0.5 cm⁻³) gas. While, in principle, ionization followed by thermal expansion of the H II region can create a bubble of ~ 2 pc size, the hot gas is the tell-tale signature that the stellar wind of massive stars, rather than stellar radiation, is driving the expansion and forming the EON cavity. As Güdel et al. (2008) note, the emission characteristics in conjunction with its structure and young age render it unlikely that the bubble is a supernova remnant. While the observed morphology is in qualitative agreement with simple models for stellar winds from massive stars (Weaver et al. 1977), the observed plasma temperature is lower than expected from refined models and suggests that mass loading of the hot plasma has been important (Arthur 2012).

Recently, the inner shocked wind bubble surrounding θ¹ Ori C has been identified in optical line observations (Abel et al. 2019). This inner shock will heat the gas in the EON that drives the expansion of a larger, outer, shell. These same observations, however, suggest that the thus heated gas is only free to escape through the south-west of the Huygens Region.

The limb-brightened edge of the Veil shell, the dense shell associated with M 42, is readily observed edge-on. It is a closed, surprisingly spherically symmetric structure enveloping the inner Huygens Region and the EON (Pabst et al. 2019), and confining the hot X-ray emitting and ionized gas observed by Güdel et al. (2008) and O’Dell & Harris (2010) in the foreground. Likewise, the H II regions of M 43 and NGC 1977 are surrounded by rather dense shells. The limb-brightened shell of M 43 exhibits a PDR-like layered structure (cf. Fig. 10). At the southern edge of the shell surrounding NGC 1977 one encounters OMC3, which also possesses a PDR-like structure, irradiated by 42 Orionis (Kabanovic et al., in prep.).

For our analysis we adopt a distance of 414 ± 7 pc (Menten et al. 2007) towards the Orion Nebula complex, although more recent results suggest somewhat lower values, 388 ± 5 pc (Kounkel et al. 2017), which is in agreement with Gaia DR2 results (Großschedl et al. 2018). In view of other uncertainties in the analysis, we consider this a minor source of uncertainty.

4.1 The pressure balance

Table 5 summarizes the physical conditions in the H II regions and the limb-brightened PDR shells in M 42, M 43, and NGC 1977. Table 6 summarizes the various pressure terms in the total pressure p_tot = p_thermal + p_magnetic + p_turb + p_lines + p_rad in the PDRs of M 42 (Orion Bar and Veil shell), M 43, and NGC 1977.

The thermal pressure, p_thermal = nkT, in the Orion Bar is controversial. Constant-pressure and H II region models of atomic and molecular lines of Pellegrini et al. (2009) and observations of molecular hydrogen by Allers et al. (2005) (n ≃ 10⁵ cm⁻³, T_gas ≃ 500 K) indicate a pressure of p∕k_B ≃ 5 × 10⁷ K cm⁻³. Observations of [C II] and [O I] emission in the Orion Bar suggest similar values (n ≳ 10⁵ cm⁻³, T_gas ≳ 300 K) (Bernard-Salas et al. 2012; Goicoechea et al. 2015). High-J CO observations indicate a higher gas pressure within the Bar, 3 × 10⁸ K cm⁻³ (Joblin et al. 2018). This latter value is in agreement with observations of carbon radio recombination lines (CRRLs) that measure the electron density in the PDR directly (n ≳ 4−7 × 10⁵ cm⁻³) (Cuadrado et al. 2019). The thermal pressures in the EON portion of M 42, M 43, and NGC 1977 have been calculated from the parameters given in Table 5.

The magnetic field strength in the Orion Bar has been estimated from the far-IR dust polarization using the Davis-Chandrasekhar-Fermi method to be 300 μG (Chuss et al. 2019). This corresponds to a magnetic field pressure, p_magnetic = B²∕8π, of 3 × 10⁷ K cm⁻³. The magnetic field in the Veil in front of the Huygens Region is measured to be B_los ≃−50−(−75) μG from the H Iand OH Zeeman effect (Troland et al. 2016); we compute the magnetic field pressure from the lower value⁸. The turbulent pressure, , is calculated from the [C II] line width in the average spectra towards the Veil shell, M 43, and NGC 1977 (Figs. 3, 11, and 18) after correction for thermal broadening at a kinetic temperature of T_gas ~ 100 K (cf. Table 5); with a typical line width of Δv_FWHM ≃ 4 km s⁻¹ this gives σ_turb ≃ 1.7 km s⁻¹. In the Orion Bar, we use σ_turb ≃ 1.5 km s⁻¹ (Goicoechea et al. 2015).

The radiation pressure is in principle given by ; the luminosities of the central stars are given in Table 8. We assume the (projected) distances R_S ≃ 4 pc for the (southern) Veil shell, R_S ≃ 0.25 pc for M 43, and R_S ≃ 1.0 pc for NGC 1977. However, for the Orion Bar, a direct measurement of the infrared flux suggests a value that is an order of magnitude lower than the pressure calculated from the stellar luminosity and the projected distance of R_S ≃ 0.114 pc (Pellegrini et al. 2009), an incident radiation field of G₀ = 2.6 × 10⁴ (Marconi et al. 1998; Salgado et al. 2016). Hence, we choose to calculate the radiation pressure in the Orion Bar from this latter value. We note that in the other cases the radiation pressure given in Table 6 is, thus, an upper limit.

Resonant scattering of Lyα photons constitutes the major contribution to the line pressure term p_lines. We expect this term to be at most of the order of the radiation pressure and hence negligible in the cases presented here (Krumholz & Matzner 2009).

We consider now the PDR associated with the Orion Bar, which is likely characteristic for the dense PDR associated with the OMC1 core. Perusing Tables 5 and 6, we conclude that the thermal pressure derived from high-J CO lines and CRRLs exceeds the turbulent and magnetic pressure by an order of magnitude. While the [C II] line stems from the surface of the PDR, high-J CO lines and CRRLs stem from deeper within the PDR. Models of photoevaporating PDRs suggest that the pressure increases with depth into the PDR (Bron et al. 2018). The [C II]-emitting PDR surface may be well described by a lower pressure (5 × 10⁷ K cm⁻³), in which case approximate equipartition holds between the three pressure terms. Their combined pressures of 1 × 10⁸ K cm⁻³ then is balanced by the thermal pressure of the ionized gas, 8 × 10⁷ K cm⁻³, with a contribution from the radiation pressure.

For the PDR in the Veil shell, the thermal, turbulent and magnetic pressure are again in approximate equipartition. In this case, there is approximate pressure equilibrium between the combined pressures in the PDR gas and the combined pressures of the ionized gas and the hot plasma in the EON. Overall, there is a clear pressure gradient from the dense molecular cloud core behind the Trapezium stars to the Veil shell in front. This strong pressure gradient is responsible for the rapid expansion of the stellar wind bubbletowards us and sets up the ionized gas flow, which (almost freely) expands away from the ionization front at about 10 km s⁻¹ (García-Díaz et al. 2008; O’Dell et al. 2009). During the initial phase, the thermal pressure of the hot gas drives the expansion of the shell and radiation pressure provides only a minor contribution (Silich & Tenorio-Tagle 2013). Radiation pressure takes over once the plasma has cooled through energy conduction and the bubble enters the momentum-driven phase.

Turning now towards the H II regions, M 43 and NGC 1977, we note that the sum of the thermal and turbulent pressures in the PDRs of these two regions is well below the thermal pressure of the ionized gas. While, in principle, the pressure could be balanced by the magnetic field, we expect equipartition to hold as well. Rather, we take this pressure difference to imply that the overpressureof the ionized gas drives the expansion of these H II regions. We note that radiation pressure is not important for these two nebulae.

4.2 Rayleigh-Taylor instability of the Veil shell

Molecular clouds, under the influence of massive stars, can develop hydrodynamical instabilities that compress or fragment the cloud. In our [C II] data as well as in the Hα image, we see evidence that the expanding gas might suffer from Rayleigh-Taylor instability, leading to the formation of elongated structures perpendicular to the expanding surface (cf. Figs. 1a and b). These fingers are particularly clear in the EON towards the south of the Veil shell. Figure 19 shows this part of the Veil shell with its complex velocity structure. Rayleigh-Taylor instability occurs at the interface of a “heavy” fluid (here the expanding neutral [C II] gas) floating on top of a “light” fluid (here the hot plasma). This instability leads to increased turbulence and mixing of the gas. The presence of magnetic fields may alter the growth behavior of the characteristic structures, but it is yet under debate whether Rayleigh-Taylor instability will be enhanced or suppressed (Stone & Gardiner 2007; Carlyle & Hillier 2017), the details likely depending crucially on the magnetic field configuration, which is hard to determine observationally.

The scale size of the unstable region after a time t is given by: (10)

where g is the (constant) acceleration and the Atwood number is a measureof the density contrast across the contact discontinuity (Chevalier et al. 1992; Duffell 2016). In this case, , as the density contrast is large. Theparameter α is determined by experiments and numerical simulations to be ~ 0.05. Realizing that the radius of the contact discontinuity is , experiments and simulations predict that h∕R_cd ~ 0.1, which is in good agreement with the observed structures, that have a typical height of ~ 100′′ ≃ 0.2 pc in a shell with a radius of ~2.5 pc. This lends further support for the interpretation of the observed morphology as the result of the Rayleigh-Taylor instability and for the general interpretation of the dominance of the stellar wind in driving the expansion of the shell.

[C II] observations provide a unique opportunity to probe the injection of turbulence into a molecular cloud through expanding shock waves, as [C II] emission traces the surface layer that transmits the shock wave into the surrounding medium. The morphology and velocity structure of our [C II] observations seem to indicate the importance of Rayleigh-Taylor and Kelvin-Helmholtz instabilities⁹ in the transmittance of turbulence into the expanding Veil shell. The simulations of Nakamura et al. (2006) show a similar velocity structure in the synthesized spectra as our observed spectra (Fig. 20). With a typical velocity dispersion of 3 km s⁻¹, the energy injected into the Veil shell by the Rayleigh-Taylor instability can be estimated at 3² ∕13² of the kinetic energy of the Veil shell, that is 1.4 × 10⁴⁷ erg. Such macro-turbulence, associated with cloud destruction, will eventually be converted into micro-turbulence of the cloud medium (Nakamura et al. 2006). Post-processing of such numerical simulations in the [C II] emission and comparison to the data might be very illuminating.

Fig. 19 Three-color image of [C II] velocity channels of the southern Veil shell. Blue is the velocity channel vLSR =0−2 km s−1, green vLSR = 4−6 km s−1, and red vLSR = 8−10 km s−1. Arrows indicate structures that resemble fingers formed due to the Rayleigh-Taylor (RT) instability. One of these, labeled “possible pillar”, might also be formed by an overdensity encountered by the expanding shell and subsequently sculpted by the strong stellar radiation field, as it also is observed in CO (cf. Fig. 1c). The spectra shown in Fig. 20 were extracted towards the areas indicated by the numbered circles.

Fig. 20 [C II] spectra towards southern Veil shell. Each spectrum is averaged over a circle with a radius of 40′′ as indicated in Fig. 19. We note that each spectrum consists of multiple line components, which is characteristic of thermodynamic instabilities.

Fig. 21 [C II] (green), 12CO(2–1) (yellow), and H I (blue) spectra towards the center of M 43, averaged over a circle with radius r ≃40′′ and center (Δα, Δδ) =(–226′′,440′′). Plotted are also the results of Gaussian fits (dashed lines). The fit parameters are given in Table 7. The H I spectrum is multiplied by 10−2 to fit the scale.

4.3 The structure of M 43

Both [C II] and H I show three velocity components while CO shows only two (Fig. 21 and Table 7). We attribute the [C II] high-velocity component (10.7 km s⁻¹) to emission from the background molecular gas associated with OMC1, that we also recognize in CO emission (cf. Fig. 8). The 4.7 km s⁻¹ [C II] component can be attributed to a half-shell expanding towards us. As this component is also visible in H I absorption against the radio continuum emission from the H II region, this component can be placed firmly in front. We note the difference in line width between the H I and [C II] (3.2 km s⁻¹ versus 1.7 km s⁻¹), which we attribute to the importance of thermal broadening for the former. This would imply a temperature of the neutral gas of ~ 150 K.

The 8.3 km s⁻¹ [C II] component is more enigmatic. Given its large line width, we are tempted to associated this component with emission from the ionized gas. However, a 8.3 km s⁻¹ component is also present in absorption in H I and this would place it in front of the ionized gas. Perhaps, this H I component is really associated with the Veil partially overhanging M 43, whereas the [C II] component is clearly confined to the interior of M 43 (cf. Fig. 12).

Furthermore, we do not see a counterpart to the CO 9.8 km s⁻¹ component in [C II] and this gas must be associated with a structure in OMC1. Finally, the very broad, low velocity (1.7 km s⁻¹) H I emission component seems to form a coherent structure with H I emission associated with the Northeast Dark Lane (as can be seen from the comparison of the pv diagrams in Fig. 8). This H I component is also seen in dust extinction but there is no counterpart in the [C II] line, in the dust continuum at 70 μm, or in the PAH emission at 8 μm. Hence, this component is not illuminated by either θ¹ Ori C or NU Ori, placing this material firmly in front of the Veil. As has been suggested by Troland et al. (2016), this molecular gas very likely has been undisturbed by the star-formation process in OMC1.

In the [C II] line-integrated map (cf. Fig. 9), the limb-brightened shell structure of M 43 separates into two parts, with a break towards the east. This rupture is reflected in the pv diagram (Fig. 12) through the center of M 43, exhibiting a gap in the middle of the arc. In pv diagrams along the RA axis, this disruption, however, is located in the western part of M 43 (cf. Fig. C.1)

4.4 The Veil stellar-wind bubble

Güdel et al. (2008) in their influential paper on the X-ray emission from Orion revealed that the interior of the Veil shell as outlined by the [C II] line, the 8 μm PAH emission, and the 70 μm dust emission is filled with a hot (T ~ 2 × 10⁶ K), tenuous (n ~ 0.5 cm⁻³) plasma with a total mass of ~0.07 M_⊙. This hot plasma is generated by the reverse shock acting on the fast stellar wind originating from θ¹ Ori C. The observations reported in this study provide the velocity, radius and mass of the shell swept up by the bubble. Following Güdel et al. (2008), we will use the stellar wind theory of Castor et al. (1975) and Weaver et al. (1977) to analyze these data.

Based on models of Castor et al. (1975) we have obtained the following two time-dependent equations for a thin shell expanding under the influence of stellar winds, under the assumption that the stellar wind is constant; this assumption is valid during most of the lifetime of a star (Haid et al. 2018) and also in the so-called snowplow phase that is discussed by Castor et al. (1975):

where Ṁ₆ is the mass-loss rate in 10⁻⁶ M_⊙ yr⁻¹, v₂₀₀₀ is the wind velocity in 2000 km s⁻¹, and n₀ is the density of the ambient, uncompressed gas. For θ¹ Ori C, the wind parameters we adopt are Ṁ₆ ≃ 0.4 and v₂₀₀₀ ≃ 1.25 (Howarth & Prinja 1989; Stahl et al. 1996).

We thus have two equations (R_S, Ṙ_S) and two unknowns (n₀, t). We denote the solution we find for the time between the onset of expansion and the present time by t₀. Re-arranging Eqs. (11) and (12) gives

From our observations we estimate the present-day values R_S (t₀) ≃ 4 ± 0.5 pc and v_exp = Ṙ_S(t₀) ≃ 13 ± 2 km s⁻¹ at the far (south) side of the expanding shell. With these values, we obtain for the unknowns n₀ ≃ 5 ± 2 × 10¹ cm⁻³ and t₀ ≃ 0.2 ± 0.05 Myr.

With these values, the temperature of the hot plasma inside the bubble and its density are predicted by the model to be T_plasma ≃ 2.3 ± 0.2 × 10⁶ K and n_plasma ≃ 0.2 ± 0.05 cm⁻³ (Castor et al. 1975). This is in excellent agreement with the X-ray observations by Güdel et al. (2008). However, numerical models of Arthur (2012), that take into account saturation of thermal conduction at the shell interface and mass loading by proplyds in the Orion nebula, produce hotter gas with T_plasma ~ 10⁷ K. Yet, the mixing efficiency of the photo-evaporation flow from the proplyds with the stellar wind material is not well constrained and a majorsource of uncertainties.

In an alternative approach, we fit the five observables (R_S, v_S, M_S, T_plasma, n_plasma) = (2.7 pc, 13 km s⁻¹, 1500 M_⊙, 2 × 10⁶ K, 0.3 cm⁻³) with the two parameters t₀ and n₀, given the expressions by Castor et al. (1975). From a weighted least-square fit, with relative errors (0.2, 0.2, 0.5, 0.3, 0.5), we find the best fit for t₀ ≃ 0.24 ± 0.05 Myr and n₀ ≃ 1.7 ± 0.7 × 10² cm⁻³. In this fit, only the last three observables are fitted in good agreement; in contrast, the shell radius is predicted to be 3.8 pc and the expansion velocity as 9 km s⁻¹. The plasma temperature and density do not depend strongly on t₀ and n₀ (Castor et al. 1975). Adopting a lower gas mass of 600 M_⊙ results in a better fit with t₀ ≃ 0.17 ± 0.05 Myr and n₀ ≃ 1.0 ± 0.6 × 10² cm⁻³. In this scheme, however, we rely on the spherical symmetry of the problem, which is not applicable for the Orion Nebula expanding predominantly along its steepest density gradient. Considering expansion in an ambient medium with a density gradient scaling with alters the time-dependent behavior of the bubble kinematics considerably (Ostriker & McKee 1988, Sect. VII).

From the morphology of the expansion in the eastern arm of the shell (cf. Appendix B, Figs. B.1 and B.2), we obtain with R_S(t) ≃ 1.75 pc and Ṙ_S(t) ≃ 4.25 km s⁻¹ n₀ ≃ 8 × 10³ cm⁻³ and t₀ ≃ 0.25 Myr. Here, the ambient gas density is thus much higher.

The derived lifetime of the bubble, t₀ ≃ 0.2 Myr, indicates that the onset of the gas coupling to the stellar wind is rather recent, but agrees with the presumed ages of the Trapezium cluster and θ¹ Ori C, which is < 1 Myr with a median age of the stellar population of 0.3 Myr (Prosser et al. 1994; Hillenbrand 1997). In contrast, Simón-Díaz et al. (2006) give an age of θ¹ Ori C of 2.5 ± 0.5 Myr, derived from optical spectroscopy of the stars in the Trapezium cluster, which is much longer than the expansion time scale of the Veil bubble. We note, though, that the star has to leave the dense core where it was born before large-scaleexpansion can set on in a relatively dilute ambient medium. We expect this time to be of the order of 0.1 Myr, considering the peculiar velocity of θ¹ Ori C (Vitrichenko 2002; Stahl et al. 2008; O’Dell et al. 2009). From the absence of evidence for the destruction of nearby proplyds and interaction between the moving foreground layers, O’Dell et al. (2009) estimate the time scale since θ¹ Ori C has moved into lower-density material, starting to ionize the EON, is of a few 10⁴ yr. This estimate is an order of magnitude younger than the expansion time we find. However, from modelling the disk masses of proplyds close to θ¹ Ori C, Störzer & Hollenbach (1999) find that θ¹ Ori C can be significantly older than 10⁵ yr if the orbits of the proplyds are radial. Salgado et al. (2016) derive an age of 10⁵ yr for θ¹ Ori C from the mass of the Orion Bar.

For the density within the shell, the compressed gas, the simple model of Castor et al. (1975) gives (15)

If the temperature in the shell is T_s = 80 K, as assumed by Castor et al. (1975), and μ∕m_H = 1.3, assuming atomic H I gas, the density within the shell becomes . In the southern shell, this then yields n_s ≃ 400 n₀ ≃ 2 × 10⁴ cm⁻³. In the eastern bright arm of the shell, we expect a shell density of n_s ≃ 50 n₀ ≃ 5 × 10⁵ cm⁻³: here, the density is much higher but also the shell is much thinner than in the southern part. Both theoretical densities are an order of magnitude higher than what we estimated from observations.

With a mass of the expanding gas shell of M ≃ 1500 M_⊙ and an expansion velocity of v_exp ≃ 13 km s⁻¹, the kinetic energy in the shell is computed to be E_kin ≃ 2.5 × 10⁴⁸ erg. Dividing this energy by the expansion time scale provides us with a measure that we can compare to the wind luminosity of the presumed driving star θ¹ Ori C. This yields . Comparing this number with the wind luminosity, L_w ≃ 8 × 10³⁵ erg s⁻¹, we estimate a wind efficiency of about 50%. Adding the wind of θ² Ori A, which is comparable to the stellar wind of θ¹ Ori C, with L_w ≃ 2 × 10³⁵ erg s⁻¹, reduces the wind efficiency to 40%; the stellar wind parameters of θ² Ori A are Ṁ₆ ~ 0.08 and v₂₀₀₀ ~ 1.45 (Petit et al. 2012). Adding both winds in the above analysis, would increase densities by a factor 1.3, but has no effect on the expansion time.

During the initial stages, the swept up shell will be very thin and the stellar ionizing flux will be able to fully ionize the shell. Once the shell has swept up enough material, the ionization front will become trapped in the shell. The Veil shell has already reached this later stage as evidenced by the morphology of the shell in Hα and [C II] emission (cf. Fig. 1) and the narrow line width of the [C II] line. In addition, we note that bright PAH emission (cf. Fig. 1 in Pabst et al. 2019) mostly originates from the neutral gas in PDRs and is very weak if associated at all with ionized gas. For the ambient densities we estimated above, that is n₀ ≃ 2−5 × 10¹ cm⁻³, the time when the shell becomes neutral is t ≃ 0.1−0.2 Myr according to the Weaver et al. (1977) model¹⁰. This is in good agreement with our estimate of the expansion time.

As the [C II] observations reveal, the strong density gradient towards the front of the cloud has produced a rapid expansion of the bubble towards us. We note that the 13 km s⁻¹ expansion velocity exceeds the escape velocity from OMC1 (2 km s⁻¹) and eventually the swept-up mass will be lost to the environment. In addition, the bubble will break open, releasing the hot plasma and the expanding ionized gas. Inspection of the pv diagrams (Figs. C.1 and C.2) indicates that the Veil shell is quite thin in some directions and break out may be “imminent”. Eventually, the shell material, the hot plasma, and the entrained ionized gas will be mixed into the interior of the Orion-Eridanus superbubble. The next supernova that will go off, one of the Orion Belt stars in ~ 1 Myr, will sweep up this “loose” material and transport it to the walls of the superbubble. The last time this happened, about 1 Myr ago, this created Barnard’s loop.

New numerical simulations show that the effects of stellar wind on their surroundings depends strongly on the characteristics of the medium (Haid et al. 2018). Stellar winds dominate if the bubble expands into an ambient warm ionized medium (WIM, T ~10⁴ K, n ~0.1 cm⁻³). If the ambient medium is colder and neutral (CNM, T ~ 100 K, n ~ 100 cm⁻³), ionizing radiation is the dominant driving force. The warm neutral medium (WNM, T ~ 2 × 10³ K, n ~ 1 cm⁻³) is an intermediate case. However, these models are not descriptive of the situation of θ¹ Ori C. New studies need to investigate the effects of stellar winds of a star in the outskirts of a molecular cloud, an environment that is rather dilute (n ~ 10 cm⁻³) and subject to density gradients. We realize that the Veil shell provides an excellent test case for models of stellar wind expansion into a surrounding medium given the large number of observational constraints, including the mass of the shell, the temperature and density of the shell and of the plasma, the velocities of the shell and the dense ionized gas in the Huygens Region.

Fig. 22 PV diagram along position cut shown in Fig. 13 of Van der Werf et al. (2013), perpendicular to the Orion Bar through the Huygens Region.

4.5 The expanding bubble of θ² Ori A

Van der Werf et al. (2013) discuss the shell structure that seems to surround θ² Ori A, their component L, that can be identified in our [C II] data, as well. Figure 22 shows the same cut for the [C II] intensity as Fig. 13 of Van der Werf et al. (2013) for the H I optical depth. Component L is a rather small arc-like structure, r ≃50′′ ≃ 0.1 pc extending to a LSR velocity of v_LSR ≃−3.2 km s⁻¹. This would translate into an expansion velocity of v_exp ≃ 11 km s⁻¹. If this was a wind blown bubble, we would obtain from Eqs. (13) and (14), t₀ ≃ 0.006 Myr and n₀ ≃ 4 × 10⁴cm⁻³. This might be reasonable, considering that the star might only very recently have emerged from the dense material were it was born. It might also be that θ² Ori A is the ionizing source for this structure, but not its cause. Schulz et al. (2006) estimate the stellar age to be ~0.3 Myr, similar to θ¹ Ori C, but not much younger.

Also the arc-like structure D of Van der Werf et al. (2013) has been discussed as associated with θ² Ori A (García-Díaz & Henney 2007), but Van der Werf et al. (2013) attribute it to θ² Ori B, an early B-type star. However, in the H I cuts component D only appears with a radius of r ≃100′′, while we observe a larger structure, r ≃500′′ ≃ 1 pc, that has a large negative velocity (v_LSR ≃−20 km s⁻¹). The latter’s geometric center seems to lie more to the southeast, as can be seen in our pv diagrams (cf. Fig. C.1). We do not observe component D in Fig. 22, indicating that it is not illuminated by UV radiation and hence a foreground component. Assuming the larger [C II] structure to be the wind-blown bubble of θ² Ori A, comparison with models would give t₀ ≃ 0.002 Myr and n₀ ≃ 3 × 10¹cm⁻³ at the expansion velocity of v_exp ≃ 28 km s⁻¹. Possibly the wind-blown bubble of θ² Ori A has been disrupted by the stellar wind of θ¹ Ori C and merged with the second’s bubble.

4.6 The thermal bubbles of M 43 and NGC 1977

The ionizing stars of M 43 and NGC 1977 are of spectral types B0.5 and B1, respectively, and such stars have feeble stellar winds that are incapable of driving a bubble into their surroundings. We attribute the bubbles associated with these two nebulae to the thermal expansion of the ionized gas (Spitzer 1978). We will return to the stellar wind aspects later.

Simón-Díaz et al. (2011) conclude that M 43 is ionization-bounded. The Strömgren radius of the H II region environing NU Ori with the rate of ionizing photons (Simón-Díaz et al. 2011) is . Assuming the electron density Simón-Díaz et al. (2011) report, derived from optical line ratios, n_e ≃ 5 × 10² cm⁻³, , which is in excellent agreement with the radius we measure.

After the initial phase of rapid ionization, the overpressure of the ionized gas compared to its environment drives the thermal expansion of the H II region (Krumholz & Matzner 2009; Spitzer 1978). For a homogeneous region, the present-time size and velocity of expansion (R_S, Ṙ_S) are related to the expansion time and initial density (t₀, n₀):

where c_s ≃ 10 km s⁻¹ is the sound speed in the ionized gas.

For M 43, we obtain with R_S(t₀) ≃ 0.3 ± 0.05 Myr and v_exp = Ṙ_S(t₀) ≃ 6 ± 2 km s⁻¹ an expansion time of t₀ ≃ 0.02 ± 0.01 Myr and an initial density of n₀ ≃ 1.4 ± 0.7 × 10³ cm⁻³. This implies that the star NU Ori is rather young. The expected main-sequence lifetime of a ~ 20 M_⊙ star like NU Ori is ~ 10 Myr.

From the Hα emission in NGC 1977, we derive an electron density of n_e ≃ 35 cm⁻³ throughout the region, assuming a radius of r ≃ 1 pc. With the rate of ionizing photons (Kim et al. 2016; Diaz-Miller et al. 1998), the Strömgren radius of the H II region surrounding 42 Orionis is . This is smaller than the size of the shell we observe. However, Hohle et al. (2010) give an effective temperature of T_eff= 25 400 K for 42 Orionis. Also from its spectral class, B1V, we expect a much higher flux of ionizing photons, (Sternberg et al. 2003). We can estimate the rate of ionizing photons emitted by 42 Orionis (with minor assistance from other stars) from the Hα emission in the H II region. We have:

With , j_{H α} ∕j_Hβ = 2.85 (Hummer & Storey 1987), and α_B = 2.6 × 10⁻¹³ cm³ s⁻¹ (Storey & Hummer 1995) we obtain from the Hα luminsity integrated over NGC 1977, L(Hα) ≃ 45 L_⊙, a photon rate of . With this value, we derive a present-day Strömgren radius of , which is in good agreement with the observed radius of the shell (r = 1.0−1.6 pc). With the extent of the expanding shell derived from the pv diagram in Fig. 17, R_S (t₀) ≃ 1 pc ± 0.2, and the expansion velocity v_exp = Ṙ_S(t₀) ≃ 1.5 ± 0.5 km s⁻¹ the expansion time is t₀ ≃ 0.4 ± 0.2 Myr and the initial density n₀ ≃ 2.5 ± 1.2 × 10³ cm⁻³. This, as well, suggests that 42 Orionis is rather young, compared to its expected lifetime of ~ 10 Myr.

4.7 Stellar wind versus thermal expansion of bubbles

In the previous sections, we analyze the kinematics of the bubbles around M 42, M 43, and NGC 1977 in terms of stellar-wind driven and thermal-pressure driven expansions, respectively. Here, we compare these two analyses. In Tables 8 and 9, we summarize our findings for the three regions.

For M 42, the momentum and energy in the ionized gas in the champagne-flow model is about two orders of magnitude less than observed for the Veil shell: 22 M_⊙ of ionized gas moving at 17 km s⁻¹ (O’Dell 2001) versus a (half)shell of 1500 M_⊙ expanding at 13 km s⁻¹ (Table 9), the neutral Veil shell carrying a total momentum of 2 × 10⁴ M_⊙ km s⁻¹. For M 43 and NGC 1977, on the other hand, the stellar wind is too feeble to drive the expansion of a large shell and the results are in good agreement with analytical models for the thermal expansion of H II regions. Figures 23 and 24 summarize this graphically. In this depiction, we have calculated the expected kinetic energy of a stellar wind bubble using the analytical model developed by Weaver et al. (1977), (20)

and the pressure-driven expansion of ionized gas (Spitzer 1978),

where is the initial Strömgren radius and t_s = R_S,0∕c_s. Comparing the predicted shell masses, where we assume that all the swept-up material is in the shell, to the measured masses as shown in Fig. 24, we observe that the measured masses of all three shells are slightly above the predictions, by a factor of about two. The mass discrepancy in M 42 and NGC 1977 might be due to our choice of β in the SEDs and a contamination from background material that is not expanding. The mass of M 43, estimated from the gas density in the shell, is a rather rough estimate as well.

For the Veil shell, the kinetic energy is comparable to the integrated wind luminosity of a wind-blown bubble. In M 43 and NGC 1977, the observed kinetic energy of the shells much exceeds the integrated wind luminosity. For M 43 and NGC 1977, the kinetics of the shell are in good agreement with the thermal expansion of the H II region, while for M 42 this falls short by two orders of magnitude.

From models of O stars and comparison with real O stars (Pauldrach et al. 2001), we obtain an upper limit of L_w ≃ 1 × 10³⁴ erg s⁻¹. In their models of earlier B stars, however, even a B0.5V star has a significant stellar wind with L_w ≃ 1 × 10³³ erg s⁻¹ (Sternberg et al. 2003). Since no X-ray observations are available for M 43 and NGC 1977, that might reveal gas shock-heated from stellar winds in the bubble interior as in M 42, we cannot judge conclusively on this matter.

Fig. 23 Kinetic energy of expanding bubble shells versus expansion time. The lines are the predictions of wind models (blue) and models of pressure-driven expansion (red) with the parameters compiled in Tables 8 and 9.

Fig. 24 Masses of expanding bubble shells versus expansion time. The lines are the predictions of wind models (blue) and models of pressure-driven expansion (red) with the parameters compiled in Tables 8 and 9.

5 Conclusion

We have observed the [C II] emission from the Orion Nebula M 42, M 43, and NGC 1973, 1975, and 1977. Velocity-resolved observationsof the [C II] fine-structure line provide a unique tool to accurately quantify the kinematics of expanding bubbles in the ISM, in our case the large bubble associated with the Orion Nebula, the Orion Veil, and the bubbles of M 43 and NGC 1973, 1975, and 1977. Our extended map of [C II] emission can be further used to disentangle dynamic structures within the same region on larger scales than previously possible. From the consistent arc structures in the pv diagrams and the line profiles, we are able to show that the gas is indeed expanding coherently on large scales as a bubble is blown into the ISM by stellar feedback.

We have compared three regions within the Orion Nebula complex that host different central stars causing the disruption of the molecular cloud cores in which these stars are born. The Orion Nebula M 42 hosts the most massive star, an O-type star, which has a strong wind, that drives the expansion of the surrounding medium. M 43 and NGC 1977 host less massive stars of early B type with much weaker winds. Still, both are effective in blowing bubbles of significant extent.

Pabst et al. (2019) have demonstrated that the Veil shell is a coherent, half-spherical structure, enclosing the layers of ionized gas observed towards that region. The stellar wind from θ¹ Ori C seems to be very effective in sweeping up its surrounding gas into a large-scale dense shell, leaving the bubble interior filled with dilute, hot X-ray emitting gas. This precise structure, a thin, rather cold shell of swept-up gas enveloping a hot interior medium, is characteristic of wind-blown bubbles. The bubble of the Orion Nebula, though still closed, may be on the brink of entering the bursting phase, thus rejuvenating the medium of the Orion-Eridanus superbubble.

The bubbles associated with M 43 and NGC 1977 also exhibit spherical symmetry and possess a rather thin shell. In these two cases, however, the masses and kinetics of these shells are incompatible with an origin in a wind-driven shell. Instead, these shells are due to the pressure-driven expansion of the ionized gas and the observations are in good agreement with simple analytical estimates. All structures we have investigated here are rather young compared to the total lifetime of the driving (massive) stars. As stars in the Orion Nebula cluster will go off as supernovae, these structures will be disrupted as they age.

Acknowledgements

This work is based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. J.R.G. thanks the ERC and the Spanish MCIU for funding support under grants ERC-2013-Syg-610256-NANOCOSMOS and AYA2017-85111-P, respectively. Research on the interstellar medium at Leiden Observatory is supported through a Spinoza award of the Dutch Science Organisation (NWO). We are grateful for the dedicated work during the upGREAT square-degree survey of Orion of the USRA and NASA staff of the Armstrong Flight Research Center in Palmdale and of the Ames Research Center in Mountain View, and the DSI. We thank C.R. O’Dell for commenting on an earlier version of the manuscript and helpful discussions. We thank the anonymous referee for comments that helped improve the manuscript.

Appendix A SED results

Figure A.1 shows the results from an SED fit in individual points with the emissivity index β =2. The derived dust temperature decreases with increasing β, while the dust optical depth increases significantly. The integrated FIR intensity does not depend strongly on β.

Fig. A.1 SED fits to individual points throughout the Orion Nebula complex. The flux in the PACS bands is plotted as circles; fluxes of the SPIRE bands are plotted as crosses. Numbers in brackets are the RA and Dec offset, respectively, in arc seconds.

Appendix B The expanding Rim

We present here our analysis of the kinematics in the eastern bright arm, the so-called Eastern Rim, of the expanding shell. To estimate thethree-dimensional expansion velocity, we make use of the following formula: (B.1)

where v₂ − v₁ is the velocity difference of the velocity channels considered; r₁ and r₂ are the respective distances of the emission filament in these two channels. r₂ is the further outward lying filament, for an accurate calculation it must be the outermost. We can write r₁ = r₂ − Δr, where Δr is the separation between the two channel filaments. If v₂ > v₁, the shell is moving toward us and we assign a positive expansion velocity; else it is expanding away from us. Here, we use the [C II] channels at v_LSR = 0−1 km s⁻¹, v_LSR = 2−3 km s⁻¹ and v_LSR = 4−5 km s⁻¹. With v₂ − v₁ = 2 km s⁻¹, r₂ ≃ 1.75 pc and Δr ≃ 0.2 pc (cf Fig. B.1), we obtain v_exp ≃ 4.25 km s⁻¹. We note that only the distance of the furthest filament (cf. Fig. B.2) gives an accurate measure of the extent of the expanding shell, which renders the estimate from the three channels used here somewhat inaccurate; however, in higher-velocity channels the displacement between the filamentary structures is not as clear-cut as it is in these channels.

We detect the [¹³C II] F = 2−1 line in averaged spectra in this region as shown in Fig. B.3. From the peak temperatures of the [¹²C II] and [¹³C II] line (Table B.1), we compute the [C II] optical depth by Eqs. (3) and (4). The C⁺ column density is given by Eq. (6). The analysis yields τ_[C II] ≃ 3. From the column density, we estimate a gas density of the limb-brightened shell from , where [C∕H] ≃ 1.6 × 10⁻⁴ (Sofia et al. 2004) and the extent l along the line of sight is l ≃ r∕2 ≃ 1 pc. This yields n ≃ 1 × 10⁴ cm⁻³.

Fig. B.1 Color map of [C II] emission in velocity interval vLSR = 0−5 km s−1, with contours of the channel emission in vLSR = 0−1 km s−1, vLSR = 2−3 km s−1, and vLSR = 4−5 km s−1 (blue, green, and red, respectively). From the displacement of the structure of brightest emission in these channels with respect to each other, points marked by the white crosses, we can estimate the expansion velocity to be vexp ≃ 4.25 km s−1. The projected distance between the crosses is 0.2 pc, and the projected distance from the Trapezium stars of the outermost filament, traced by blue contours, is 1.75 pc.

Fig. B.2 Color map of [C II] emission in velocity interval vLSR = 0−15 km s−1, with contours of the channel emission in vLSR = 0−1 km s−1, vLSR = 2−3 km s−1, and vLSR = 4−5 km s−1 (blue, green, and red, respectively). The expanding filament grows brighter as it moves outward, probably mainly due to increasing limb-brightening, so the color map is dominated by the most outward filament.

Figure B.4 compares the CO(2–1) emission with the [C II] emission in the expanding Rim, that is observed in the velocity range v_LSR = 5−8 km s⁻¹. From the peak separation in line cut 1, d ≃70′′ ≃ 0.14 pc, we derive a density of n ≃ 9 × 10³ cm⁻³ in the Rim. From line cut 2 (d ≃110′′ ≃ 0.22 pc), we derive a density of n ≃ 6 × 10³ cm⁻³ in the southof the bright Rim.

Fig. B.3 Right: [C II] line-integrated intensity in the bright eastern Rim. White contours indicate I[C II] > 100 K km s−1.Left: averaged spectra from this region as indicated. Results of the analysis of the peak temperatures are given in Table B.1.

Fig. B.4 12CO(2–1) emission in velocity interval vLSR = 5−8 km s−1 with [C II] emission in same velocity interval in contours (dark grey: 30 K km s−1, light grey: 50 K km s−1). The inlaid panels show the line cuts along the two red lines indicated in the figure.

From the dust optical depth, we estimate a mean τ₁₆₀ ≃ 4 × 10⁻³ in the [C II]-bright regions. This translates to a gas column density of N_H ≃ 2 × 10²² cm⁻². Assuming that all carbon is ionized, this results in a C⁺ column density of , which is a factor of 2 lower than what we estimate from the [C II] spectra. Since τ₁₆₀ strongly depends on the choice of the grain emissivity index β, the agreement might be satisfactory. However, it might be a hint that we overestimate the [C II] optical depth and underestimate the excitation temperature. With these givens the physical temperature of the gas would be T ≃ 50 K, which is lower than expected from the PDR models for the [C II]-emitting (surface) layers. At these high opacities opacity broadening of the line might also become significant, which is not what we observe. The fitting procedure iscomplicated by the presence of multiple components. The main line can be fitted by two emission components or by one emission and one absorption component. However, the latter would increase the result for τ_[C II] even further. In conclusion, the values obtained here can only be a crude estimate due to the complex emission structure from this region, which might be inadmissibly averaged together in the spectra.

Appendix C PV diagrams

Figures C.1–C.4 show the pv diagrams resulting from slicing the velocity-resolved [C II] map from the Extended Orion Nebula with M 43 and NGC 1977 along the RA and Dec axis, respectively. The coherent arc structure indicative of the expansion of the Veil shell can be seen in almost all pv diagrams. The expanding structure of NGC 1977 is more difficult to recognize.

Fig. C.1 [C II] pv diagrams from the Veil shell, sliced along RA axis. Red numbers are the respective Dec offsets in arc seconds.

Fig. C.2 [C II] pv diagrams from the Veil shell, sliced along Dec axis. Red numbers are the respective RA offsets in arc seconds.

Fig. C.3 [C II] pv diagrams from NGC 1977, sliced along RA axis. Red numbers are the respective Dec offsets in arc seconds. Coordinate offsets are given with respect to the position of 42 Orionis.

Fig. C.4 [C II] pv diagrams from NGC 1977, sliced along Dec axis. Red numbers are the respective RA offsets in arc seconds. Coordinate offsets are given with respect to the position of 42 Orionis.

References

All Tables

All Figures

Fig. 1 Overview of the Orion Nebula (M 42), M 43, and NGC 1973, 1975, and 1977 in different wavelenghts. Left: DSS2 Hα emission (ESO Archive). Center: [C II] line-integrated emission (Pabst et al. 2019). Right: line-integrated 12CO(2–1) emission (Berné et al. 2014; Goicoechea et al. 2020). Hα emission stems from the ionized gas (T ~ 104 K), the [C II] lineis emitted by mostly neutral gas (T ~ 100 K), whereas CO traces the molecular gas (T ~ 30 K). M 42 and M 43 have a substantial amount of molecular gas in the background, while NGC 1977 seems to be devoid of it (the coverage of the CO map is not sufficient, however, see discussion of the expansion characteristics below). All three regions comprise ionized gas within their neutral limb-brightened shells.

In the text

Fig. 2 [C II] line-integrated intensity from the Orion Nebula (M 42), M 43, and NGC 1977. The red rectangles 1 and 2 indicate the positions of the pv diagrams shown in Figs. 5 (horizontal) and 8 (vertical). Rectangle 3 indicates the position of the pv diagram shown in Fig. 17. The red circles delineate the approximate extent of the expanding bubbles of gas that are the Veil shell, M 43, and NGC 1977, and within which the respective masses and luminosities are computed. The light-red dashed rectangle circumscribes the area of the eastern bright arm shown in Figs. B.1 and B.2. The black line indicates the position of the line cut in Fig. 4. The stars mark the position of θ1 Ori C (yellow), θ2 Ori A (orange), NU Ori (pink), and 42 Orionis (purple). Unless noted otherwise, all coordinate offsets are given with respect to the position of θ1 Ori C, (α, δ) = (5h35′ 16.46′′, − 5°23′22.85′′).

In the text

Fig. 3 Average [C II] spectrum towards the Veil shell without OMC1 and ISF. The inlaid panel shows the residual of the spectrum, in the systemic velocity of the [13C II] F = 2 − 1 line, after subtracting the [12C II] fit. The [12C II] line can be fitted by two Gaussians with TP ≃ 11.2 ± 0.3 K, vP ≃ 8.8 ± 0.1 km s−1, ΔvFWHM ≃ 5.0 ± 0.1 km s−1 and TP≃ 2.3 ± 0.1 K, vP ≃ 4.0 ± 0.4 km s−1, ΔvFWHM ≃ 9.3 ± 0.5 km s−1; the [13C II] F = 2 − 1 component is fitted by a Gaussian with TP ≃ 0.11 ± 0.02 K, vP ≃ 8.2 ± 0.3 km s−1, ΔvFWHM ≃ 2.1 ± 0.6 km s−1. The main line could also be fitted with three Gaussian components, but we compare the respective combined peak temperatures of the [12C II] and [13C II] line for a first-order estimate of the excitation conditions.

In the text

	Fig. 4 Line cut through the southwest of the Veil shell, indicated in Fig. 2. The onset of the shell is marked by a steep increase in all tracers, except CO, at Δδ ≃ −1400′′. With a line angle of 30° with respect to the vertical axis, this is at a distance of 3.25 pc from the Trapezium stars.
In the text

	Fig. 5 [C II] pv diagram along horizontal position cut 2 indicated in Fig. 2 (Δδ = −907′′ − (−831′′)). The lower panel shows the same cut with the arc structure for an expansion velocity of 13 km s−1 on a background velocity of 8 km s−1 (red dashed lines).
In the text

Fig. 6 [C II] spectra in the Veil shell taken from pv diagram in Fig. 5 with Gaussian fits. Spectra at Δα = −362′′ to Δα = −1268′′ are averaged over 151′′ × 75.5′′ in order toimprove the signal-to-noise ratio of the faint shell component to be fitted, others are averaged over 75.5′′ × 75.5′′. Coordinate offsets indicate the lower left-hand corner of the rectangle over which the spectra are averaged; Δδ = −907′′. The fit parameters are given in Table 1.

In the text

Fig. 7 Three-color image of velocity channels towards M 42 and M 43. Blue is the velocity channel vLSR =4−5 km s−1, green vLSR = 6−7 km s−1, and red vLSR = 9−10 km s−1. With increasing vLSR, the limb-brightened Veil shell filament is displaced outward, away from the bubble center. The gas of M 43 that is expanding towards us can be observed in the blue channel; the limb-brightened shell of M 43 has higher vLSR.

In the text

	Fig. 8 PV diagrams along vertical position cut 1 indicated in Fig. 2 (Δα =167′′ − 242′′) in [C II], 12CO(2–1), 13CO(2–1) (Berné et al. 2014; Goicoechea et al. 2020), and H I (Van der Werf et al. 2013).
In the text

	Fig. 9 Three-color image of M 43. Red is Hα emission from the ionized gas, green is line-integrated CO(2–1) emission from molecular gas, and blue is the [C II] line-integrated intensity, tracing the neutral gas. The line cuts in Fig. 10 are taken along lines a and b.
In the text

	Fig. 10 Line cuts along lines a and b, indicated in Fig. 9. The arrows indicate the peaks in [C II] emission (green) and CO emission (blue) that correspond to the shell.
In the text

Fig. 11 Average [C II] spectrum towards the eastern shell of M 43. The inlaid panel shows the residual of the spectrum, in the systemic velocity of the [13C II] F = 2 − 1 line, after subtracting the [12C II] fit. The [12C II] line can be fitted by a Gaussian with TP ≃ 46.2 ± 0.3 K, vP ≃ 9.6 ± 0.1 km s−1, ΔvFWHM ≃ 4.6 ± 0.1 km s−1; the [13C II] F = 2 − 1 component is fitted by a Gaussian with TP ≃ 1.1 ± 0.2 K, vP ≃ 8.6 ± 0.3 km s−1, ΔvFWHM ≃ 4.8 ± 0.8 km s−1. We note that a better fit of the [12C II] line can be obtained by fitting two components with TP ≃ 30.9 ± 1.3 K, vP ≃ 9.6 ± 0.1 km s−1, ΔvFWHM ≃ 3.3 ± 0.1 km s−1 and TP≃ 17.5 ± 1.3 K, vP ≃ 9.8 ± 0.1 km s−1, ΔvFWHM ≃ 6.8 ± 0.2 km s−1, but both fits leave significant residuals.

In the text

	Fig. 12 [C II] pv diagram along vertical position cut 1 indicated in Fig. 2 (Δα =167′′ − 242′′), zoomed in to M 43 (cf. Fig. 8). The right panel shows the same cut with the arc structure for an expansion velocity of 6 km s−1 on a background velocity of 8 km s−1 (red dashed lines).
In the text

Fig. 13 [C II] spectra along position cut 1 (cf. Fig. 12) in M 43 with Gaussian fits. Each spectrum is averaged over 75.5′′ × 25.2′′. The expanding shell has a distinguished emission peak. Additionally to the main component and the shell component, we fit a third, broad, component, that possibly stems from the ionized gas within the shell. Coordinate offsets indicate the lower left-hand corner of the rectangle over which the spectra are averaged; Δα =242′′. The fit parameters are given in Table 2.

In the text

	Fig. 14 Three-color image of NGC 1977. Red is Hα emission from the ionized gas, green is line-integrated CO(2–1) emission from molecular gas (coverage only in OMC3), and blue is the [C II] line-integrated intensity, tracing the neutral gas.
In the text

Fig. 15 Average [C II] spectrum towards [C II]-bright shell of NGC 1977 (without OMC3). The inlaid panel shows the (smoothed) residual of the spectrum, in the systemic velocity of the [13C II] F = 2 − 1 line, after subtracting the [12C II] fit. The [12C II] component can be approximated by a Gaussian with TP ≃ 16.5 ± 0.1 K, vP ≃ 12.0 ± 0.1 km s−1, ΔvFWHM ≃ 3.2 ± 0.1 km s−1; the [13C II] F = 2 − 1 component is fitted by a Gaussian with TP ≃ 0.21 ± 0.03 K, vP ≃ 11.8 ± 0.2 km s−1, ΔvFWHM ≃ 2.1 ± 0.5 km s−1.

In the text

	Fig. 16 [C II] spectrum towards NGC 1977, averaged over circular area with radius 125′′ centered at 42 Orionis ((Δα, Δδ) = (−100.4′′,1984.7′′)), with Gaussian fits. The fit parameters are given in Table 3.
In the text

	Fig. 17 [C II] pv diagram of NGC 1977 along position cut 3, indicated in Fig. 2 (Δδ =2085′′ − 2259′′). The lower panel shows the same cut with the arc structure for an expansion velocity of ± 1.5 km s−1 (red dashed lines).
In the text

Fig. 18 [C II] spectrum along pv diagram in NGC 1977 (Fig. 17) with Gaussian fits, averaged over 201′′ × 201′′. Coordinate offsets indicate the lower left-hand corner of the square over which the spectrum is averaged. The two narrow components stem from the expanding shell, while the broad component originates in the ionized gas that is contained within the shell. Since the bubble expands in both directions, the peak-velocity difference is twice the expansion velocity, vexp ≃ 1.3 km s−1. The fit parameters are given in Table 4.

In the text

Fig. 19 Three-color image of [C II] velocity channels of the southern Veil shell. Blue is the velocity channel vLSR =0−2 km s−1, green vLSR = 4−6 km s−1, and red vLSR = 8−10 km s−1. Arrows indicate structures that resemble fingers formed due to the Rayleigh-Taylor (RT) instability. One of these, labeled “possible pillar”, might also be formed by an overdensity encountered by the expanding shell and subsequently sculpted by the strong stellar radiation field, as it also is observed in CO (cf. Fig. 1c). The spectra shown in Fig. 20 were extracted towards the areas indicated by the numbered circles.

In the text

	Fig. 20 [C II] spectra towards southern Veil shell. Each spectrum is averaged over a circle with a radius of 40′′ as indicated in Fig. 19. We note that each spectrum consists of multiple line components, which is characteristic of thermodynamic instabilities.
In the text

	Fig. 21 [C II] (green), 12CO(2–1) (yellow), and H I (blue) spectra towards the center of M 43, averaged over a circle with radius r ≃40′′ and center (Δα, Δδ) =(–226′′,440′′). Plotted are also the results of Gaussian fits (dashed lines). The fit parameters are given in Table 7. The H I spectrum is multiplied by 10−2 to fit the scale.
In the text

	Fig. 22 PV diagram along position cut shown in Fig. 13 of Van der Werf et al. (2013), perpendicular to the Orion Bar through the Huygens Region.
In the text

	Fig. 23 Kinetic energy of expanding bubble shells versus expansion time. The lines are the predictions of wind models (blue) and models of pressure-driven expansion (red) with the parameters compiled in Tables 8 and 9.
In the text

	Fig. 24 Masses of expanding bubble shells versus expansion time. The lines are the predictions of wind models (blue) and models of pressure-driven expansion (red) with the parameters compiled in Tables 8 and 9.
In the text

	Fig. A.1 SED fits to individual points throughout the Orion Nebula complex. The flux in the PACS bands is plotted as circles; fluxes of the SPIRE bands are plotted as crosses. Numbers in brackets are the RA and Dec offset, respectively, in arc seconds.
In the text

Fig. B.1 Color map of [C II] emission in velocity interval vLSR = 0−5 km s−1, with contours of the channel emission in vLSR = 0−1 km s−1, vLSR = 2−3 km s−1, and vLSR = 4−5 km s−1 (blue, green, and red, respectively). From the displacement of the structure of brightest emission in these channels with respect to each other, points marked by the white crosses, we can estimate the expansion velocity to be vexp ≃ 4.25 km s−1. The projected distance between the crosses is 0.2 pc, and the projected distance from the Trapezium stars of the outermost filament, traced by blue contours, is 1.75 pc.

In the text

	Fig. B.2 Color map of [C II] emission in velocity interval vLSR = 0−15 km s−1, with contours of the channel emission in vLSR = 0−1 km s−1, vLSR = 2−3 km s−1, and vLSR = 4−5 km s−1 (blue, green, and red, respectively). The expanding filament grows brighter as it moves outward, probably mainly due to increasing limb-brightening, so the color map is dominated by the most outward filament.
In the text

	Fig. B.3 Right: [C II] line-integrated intensity in the bright eastern Rim. White contours indicate I[C II] > 100 K km s−1.Left: averaged spectra from this region as indicated. Results of the analysis of the peak temperatures are given in Table B.1.
In the text

	Fig. B.4 12CO(2–1) emission in velocity interval vLSR = 5−8 km s−1 with [C II] emission in same velocity interval in contours (dark grey: 30 K km s−1, light grey: 50 K km s−1). The inlaid panels show the line cuts along the two red lines indicated in the figure.
In the text

	Fig. C.1 [C II] pv diagrams from the Veil shell, sliced along RA axis. Red numbers are the respective Dec offsets in arc seconds.
In the text

	Fig. C.2 [C II] pv diagrams from the Veil shell, sliced along Dec axis. Red numbers are the respective RA offsets in arc seconds.
In the text

	Fig. C.3 [C II] pv diagrams from NGC 1977, sliced along RA axis. Red numbers are the respective Dec offsets in arc seconds. Coordinate offsets are given with respect to the position of 42 Orionis.
In the text

	Fig. C.4 [C II] pv diagrams from NGC 1977, sliced along Dec axis. Red numbers are the respective RA offsets in arc seconds. Coordinate offsets are given with respect to the position of 42 Orionis.
In the text