The SOFIA FEEDBACK Legacy Survey Dynamics and Mass Ejection in the Bipolar H ii Region RCW 36

The FEEDBACK C⁺ Legacy survey observed the [C ii] ³ P_3/2 → ³ P_1/2 transition, at 158 μm, and [O i] ³ P₁ → ³ P₂ transition, at 63 μm, simultaneously with the upGREAT instrument (Risacher et al. 2018) on board the airborne SOFIA (Young et al. 2012) during a flight on 2019 June 6 from Christchurch, New Zealand. The observation strategy was to carry out [C ii] and [O i] on-the-fly (OTF) maps and is described in detail in Schneider et al. (2020). For RCW 36, the planned map within the legacy program has a size of 144 × 144, indicated in Figure 1. Here we present the lower part of the map (144 × 72) that has been observed so far. The center of the full map is located at α_(J2000) = 08^h59^m268 and δ_(J2000) = −43° 44'141, and the emission-free OFF position is at α_(J2000) = 09^h01^m317 and δ_(J2000) = −43°22'500.

Both the reduced [C ii] (smoothed to 20'') and [O i] (smoothed to 30'') maps, used in this paper, have a spectral binning of 0.2 km s⁻¹. The typical noise rms of these smoothed data within 0.2 km s⁻¹ is 0.8–1.0 K for [C ii] and ∼0.8–1.5 K for [O i]. The calibration was done using the GREAT pipeline (Guan et al. 2012). To convert the antenna to main-beam temperatures, a forward efficiency η_f = 0.97 and main-beam efficiency η_mb = 0.65 for [C ii] and 0.69 for [O i] were applied. To improve the baseline removal of the data, a new method based on principal component analysis (PCA) was used (C. Buchbender et al. 2022, in preparation). This technique was also employed for the data published in Luisi et al. (2021), Tiwari et al. (2021), and Kabanovic et al. (2022). Specifically, PCA identifies the systematic components of baseline variation in different spectra, which allows us to accurately remove them and recover the actual spectrum. For a more general overview of PCA we refer the reader to Jolliffe (2011).

Complementary ¹²CO (3–2) and ¹³CO (3–2) data were simultaneously obtained with the LAsMA 7-pixel receiver on 2019 September 27 and 2021 July 21 in good weather conditions (precipitable water vapor (pwv) ≈1 mm) using the APEX telescope ¹³ (Güsten et al. 2006). This heterodyne spectrometer allows simultaneous observations of the two CO isotopomers in the upper (¹²CO (3–2) at 345.7959899 GHz) and lower (¹³CO (3–2), 330.5879653 GHz) sidebands of the receiver. The LAsMA back ends are fast Fourier transform spectrometers (FFTSs) with 2 × 4 GHz bandwidth (Klein et al. 2012) and a native spectral resolution of 61 kHz. The main-beam size of the APEX observations (θ_mb= 182 at 345.8 GHz) is similar to the resolution of the [C ii] data. The mapped field with a size of 20' × 15' was scanned in total power on-the-fly mode, with a spacing of 9'' between rows (while oversampling to 6'' in scanning direction, R.A.). The map data were calibrated against a sky reference position at α_(J2000) = 09^h01^m317 and δ_(J2000) = −43° 22'500. For the APEX observations, a first-order baseline was removed, and a main-beam efficiency η_mb = 0.68 was used to convert the antenna temperatures to main-beam temperatures. The baseline-subtracted spectra have been binned to a spectral resolution of 0.2 km s⁻¹. The final data cubes (with a pixel size of 91) have an rms of 0.45 K.

The Chandra X-ray Observatory ACIS camera (Garmire et al. 2003) observed a 17' × 17' field centered on RCW 36 for ∼70 ks in 2006. The MYStIX project (Kuhn et al. 2013) found 502 X-ray point sources in this field. These were excised, and the remaining data smoothed, to obtain an 0.5–7 keV diffuse X-ray image of the field (Townsley et al. 2014). ACIS is an intrinsic spectrometer, tagging every X-ray event with its energy and its position. We will use the Chandra/ACIS data to study the X-ray spectra in a variety of regions in and around RCW 36.

The RCW 36 region, being part of the Vela C molecular cloud, was observed with the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) bolometers on the Herschel telescope as part of the Herschel imaging survey of OB young stellar objects (HOBYS) program (Motte et al. 2010). In this work, we use the flux map at 70 μm with an angular resolution of ∼6'', as it clearly outlines the bipolar H ii region. Further, we use the dust temperature map at 36'' resolution and a dust column density map at 182 resolution. This dust temperature and column density map was obtained by fitting the spectral energy distribution (SED) to the 160–500 μm data, using the method presented in Palmeirim et al. (2013).

We further make use of data from the WISE telescope from 3.4 to 22 μm (Wright et al. 2010), an Hα map of RCW 36 from the Southern Hα Sky Survey Atlas (SHASSA; Gaustad et al. 2001), and a large-scale dust polarization map at 500 μm with a spatial resolution of $2\buildrel{\,\prime}\over{.} 5$ that was observed with the BLASTPol balloon-borne receiver (Fissel et al. 2016, 2019; Gandilo et al. 2016).

Figures 2 and 3 present the integrated intensity map for [C ii] and [O i], respectively. [C ii] is brightest toward the region that is surrounded by the high column density ring and shows more extended emission toward the two cavities. Outside of the cavity walls, the [C ii] intensity decreases rapidly. [O i], on the other hand, is only detected toward the high column density ring, where a high-density PDR can be expected.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The line-integrated intensity maps of ¹²CO (3–2) and ¹³CO (3–2), presented in Figure 4, show a different structure. The emission follows the ring-like Herschel dust column density maps rather closely and shows that the CO emission is brightest at the edge around the bright [C ii] emission, which would be expected for a PDR irradiated by the central cluster. This is not a resolution effect since both beams have similar spatial resolution. To the southwest of the ring there is a region (local peak in dust column density) where the ring is significantly less bright in ¹²CO (3–2) and ¹³CO (3–2) (by a factor 2–3), but the [C ii] contours show a protrusion toward the west. This may indicate that there is a puncture in the ring so that the [C ii]-emitting gas can freely expand into a lower-density region, driven by the stellar feedback from the east. On the other hand, the weaker CO emission at this location (C3 in Figure 2) can be partly caused by self-absorption. This self-absorption can be related to the protrusion, as it would give rise to warm and cold gas layers along the line of sight, while there still appears to be a significant column density based on the Herschel data. On a larger scale outside the ring in Figure 4, ¹²CO (3–2) shows filamentary emission that is perpendicular to the dense ring. Several of these filamentary structures are also detected in [C ii] and show a curvature that appears to originate in the central cluster. This suggests that preexisting filaments are being swept up.

Zoom In Zoom Out Reset image size

Recently, [C ii]-integrated intensity observations at an angular resolution of 90'' were presented in Suzuki et al. (2021) using a 100 cm balloon-borne infrared telescope. Smoothing the FEEDBACK-integrated intensity map to a resolution of 90'', presented in Appendix A, the results look extremely similar to the maps presented in Suzuki et al. (2021). The intensity maps appear to agree within 10%, which corresponds to typical calibration errors. Comparing the smoothed map in Figure 20 with the FEEDBACK map in Figure 2, it is evident that the higher spatial resolution of the SOFIA observations is necessary to resolve, e.g., the ring and filamentary structures in [C ii]. Another major advantage of the FEEDBACK observations is the high spectral resolution of the observations, which allows us to study the dynamics of the region.

The [C ii], CO, and [O i] spectra shown in Figure 2 are all very complex and cannot be described by a single Gaussian velocity component. The bulk emission of the molecular cloud/dense ring is best traced with the ¹³CO (3–2) line that is visible for all center positions (C1–C5) around 7–8 km s⁻¹ but mostly absent for all western (W1–W3) and eastern (E1–E2) positions. The ¹²CO (3–2) line profile shows significant variation depending on the position. For the spectra at the central molecular ring, we observe two apparent line components that are caused by foreground or self-absorption because the ¹³CO line peaks in the ¹²CO dip at ∼7 km s⁻¹. Note that the optically thin dense gas tracers N₂H⁺ (2–1) and NH₃(1,1), presented in Fissel et al. (2019), also peak at velocities around 7–8 km s⁻¹. The absorption thus mimics two separate velocity components that were interpreted by Sano et al. (2018) as two velocity components tracing a cloud–cloud collision. We observe that the absorption in ¹²CO continues toward the east and west because the dip in the line profile is always around 7 km s⁻¹ and the ¹³CO emission clearly peaks in the ¹²CO dip for position W1 (but is too weak for W2, E1, and E2). For all center positions, the self-absorption is dominant toward the blueshifted part of the main velocity component, which is generally associated with expansion (e.g., Castor & Lamers 1979).

The [C ii] line has a bright, prominent component (main-beam temperatures up to ∼100 K for the center positions) around 8–9 km s⁻¹. This bright component is also found in the eastern spectra but is less strong for the three western positions (W1–W3). This component is also affected by absorption, which we can prove with our [¹³C ii] line observations shown in Figure 5 and Appendix B. The intensity of the [C ii] main component drops where [¹³C ii] has its emission peak, similar to what was observed for several other bright [C ii] sources with SOFIA (e.g., Graf et al. 2012; Guevara et al. 2020; Kirsanova et al. 2020; Mookerjea et al. 2021; Kabanovic et al. 2022). There is also a separate second velocity component between 0 and 5 km s⁻¹ (hereafter the blueshifted component), which is well visible in [C ii] and even brighter than the 7–8 km s⁻¹ component in the west of the ring (W1–W3). At some locations (C3, C4, W3) this component is also detected in ¹²CO (3–2) (see Figure 2), but mostly at slightly higher velocities (3–5 km s⁻¹).

Zoom In Zoom Out Reset image size

The [O i] spectra, detected toward the dense ring, look slightly different than ¹²CO (3–2) and [C ii]. [O i] is only detected at velocities higher than 7 km s⁻¹. At lower velocities, i.e., between 3 and 7 km s⁻¹, the spectra indicate absorption since the emission in this velocity interval goes under the baseline in C1–C5. This is demonstrated in more detail in Appendix B, which shows that the integrated intensity of [O i] from 3 to 7 km s⁻¹ becomes negative toward the dense ring, in particular toward the Herschel column density peaks. This suggests significant [O i] 63 μm foreground absorption at the same velocities where [C ii] and ¹²CO (3–2) are absorbed, similar to what was recently found toward several bright [O i] 63 μm sources in Goldsmith et al. (2021). Based on detailed studies of several clouds, e.g., M17, MonR2 (Guevara et al., to be submitted), and RCW 120 (Kabanovic et al. 2022), it has been proposed that this could be due to a foreground cloud or cold atomic oxygen and C⁺ in the low-density halo of the cloud. Lastly, in several spectra of Figure 2, high-velocity blue- and redshifted wings are observed in the intervals from −8 to 0 km s⁻¹ and from 14 to 23 km s⁻¹. Even though these high-velocity wings are not always very clear on the presented figures because of the very bright emission at 5–12 km s⁻¹, associated with the central ring, they are well above the noise level (which is 0.8 K in individual spectral bins). The wings are particularly evident at positions C2, C3, and W3 in Figure 2. For the redshifted wing, it has to be noted that the [¹³C ii] F = 2–1 transition is located between 17 and 23 km s⁻¹, which provides a contribution (1.5 K km s⁻¹ on average) to the wing at these velocities.

3.3.1. The Expanding Shells

The [C ii] position–velocity (PV) diagram presented in Figure 6 is a cut through RCW 36, parallel to the central axis of the bipolar cavities. This provides a view on the large-scale kinematics in RCW 36. It displays a blueshifted expanding shell signature in the western region of the observed map, i.e., r > 0 pc in Figure 6 (indicated with "cavity 2"), where the western cavity is located. This expanding structure is similar to the ones observed toward bubble-like structures in Orion, RCW 120, and RCW 49 (Pabst et al. 2020; Luisi et al. 2021; Tiwari et al. 2021). The expansion feature is observed between 1 and 6 km s⁻¹, which suggests that this expanding shell has an expansion velocity of ∼5 km s⁻¹. Later in this section we will constrain the expansion velocity quantitatively. The emission from the shell is also clearly identified in the individual spectra of Figure 2 as a separate velocity component in the same velocity interval. This feature is most prominent in the spectra of the western cavity (W1–W3) and dominates over the component associated with the velocity of the molecular ring. The cavity to the east (r < 0 pc) has a similar, though less bright, expanding shell that covers the same velocity range; see Figure 6. As the eastern cavity is less bright, it is not as clear in the PV diagram, but when inspecting the [C ii] channel maps in Figure 7, it becomes evident that the [C ii] emission between 1 and 6 km s⁻¹ shows the same morphology as the western cavity. This blueshifted expanding shell morphology is also seen in the upper RGB image of Figure 6, where the center of both cavities shows blueshifted emission confined within the green limb-brightened cavity walls. Both the western and eastern cavities thus show a blueshifted expanding shell in [C ii] with similar expansion velocity. To quantify the expansion velocity of these shells in the cavities, we determined the velocity of peak intensity for the shells (i.e., at v_LSR < 6.5 km s⁻¹); see Figure 6. This curve was then fitted in each cavity by both a second-order polynomial and a sine function that assumes that the reference velocity for expansion is v_LSR = 6.5 km s⁻¹. To calculate the expansion velocity, we then determined the point of minimal v_LSR obtained from the fit. The difference between this minimal velocity and v_LSR = 6.5 km s⁻¹ is the expansion velocity. The results of these fits are summarized in Table 1. The expansion velocity in the cavities is then estimated by taking the average of all obtained expansion velocities, which gives 5.2 ± 0.5 ± 0.5 km s⁻¹. The first estimated error corresponds to the spread of all obtained expansion velocities. The second estimated error is a systematic error caused by the uncertainty on the rest v_LSR for the expanding shell, which should be between v_LSR = 6 and 7 km s⁻¹. The two expanding shells appear to spatially come together in projection at the central molecular ring, which is also observed in the channel maps in Figure 7. The ¹²CO (3–2) channel maps in Figure 9 show indications that the outer layer of the western expanding shell is also detected in ¹²CO (3–2) between 1.4 and 3 km s⁻¹. This is not the case for all expanding [C ii] shells (e.g., Luisi et al. 2021). In PDRs, CO generally exists in regions with A _V ≥ 2–3 (e.g., Tielens & Hollenbach 1985; Röllig et al. 2007), which implies that the expanding shell has a thickness of several A _V . This is consistent with the Herschel column density map of Vela C, which shows an A _V of 3–6 toward the cavities.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Minimal Velocity (v ${}_{\mathrm{LSR},\min }$) and Corresponding Expansion Velocity (v_exp) Obtained for the Two Fitted Functions to the Expanding Shells in the West and the East, Respectively

Function	${v}_{\mathrm{LSR},\min ,\mathrm{west}}$ (km s−1)	${v}_{\exp ,\mathrm{west}}$ (km s−1)	${v}_{\mathrm{LSR},\min ,\mathrm{east}}$ (km s−1)	${v}_{\exp ,\mathrm{east}}$ (km s−1)
sine	1.4 ± 0.1	5.1 ± 0.1 ± 0.5	0.8 ± 0.9	5.7 ± 0.9 ± 0.5
poly	1.5 ± 0.1	5.0 ± 0.1 ± 0.5	1.6 ± 0.1	4.9 ± 0.1 ± 0.5

Note. The fitted functions are a sine function (sine) and a second-order polynomial (poly). The indicated errors for the minimal v_LSR correspond to the error on the fit. The two indicated errors for the expansion velocities correspond to the fitting error and uncertainty on the assumed rest v_LSR for the shells.

Download table as:  ASCIITypeset image

3.3.2. The Kinematics of the Molecular Ring

Our observations also resolve the central molecular ring, which is best visible in the dust column density map and in the ¹³CO (3–2) emission (Figure 4). However, the ring is also associated with [C ii] emission as shown in the RGB images of Figure 6. The top panel indicates that the eastern part of the ring is more blueshifted than the western part. Since it was proposed, by Minier et al. (2013) and from Figure 1, that the eastern part of the ring is in front of the western part, this velocity structure fits with the proposed expansion of the molecular ring. A zoom into the ring kinematics by the channel maps in Figure 8 confirms this expansion and demonstrates that the expansion is also seen in ¹³CO (3–2). Unfortunately, it is very difficult to accurately estimate an expansion velocity for the ring. There are several reasons for this: self-absorption effects in all available lines limit an accurate determination of the typical velocity of the gas, stellar feedback breaking through the ring affects the velocity field, and the presence of additional velocity structure in the line of sight (see also A. Bij et al. 2022, in preparation) limits the possibility to use the moment 1 map. Therefore, we estimate the expansion velocity from inspecting the channel maps in Figure 8. The eastern area of the ring is found at 5–6 km s⁻¹, whereas the western part is at ∼8.0–8.8 km s⁻¹. This suggests a typical expansion velocity of 1.0–1.9 km s⁻¹. More generally, Figure 8 shows that [C ii] and CO have the same kinematic structure toward the ring, which implies that they trace the same expansion dynamics. The full ring is thus swept up by the expansion from the cluster with the CO emission at the outer side of the [C ii] emission as would be expected in a PDR structure with internal FUV irradiation. Lastly, in the southwest of the ring, where Figure 4 shows that the integrated [C ii] emission crosses the molecular ring, Figure 8 now shows that ¹³CO and [C ii] have the same kinematic morphology between 8.0 and 9.6 km s⁻¹.

Zoom In Zoom Out Reset image size

3.3.3. Filamentary Structures in [C ii]

The PV diagram in Figure 6 also shows significant [C ii] emission around 10 km s⁻¹, with corresponding ¹²CO (3–2) emission toward the eastern cavity. This emission is indicated with the name "filaments" in the PV diagrams. The typical elliptical signature expected in PV diagrams for expanding shells is not as evident for this [C ii] emission, which could be due to the small velocity difference with the ring (<2–3 km s⁻¹). In the channel maps in Figures 7 and 9, this 10 km s⁻¹ emission displays a filamentary structure that is not typical for expanding shells. However, Figure 9, in particular, also shows that these filaments have a curved morphology that appears to have their origin toward the central cluster. This [C ii] and CO emission might thus originate from low-density filaments, originally converging toward the ridge where RCW 36 formed, which are now being swept away at velocities <3 km s⁻¹ by the feedback from the stellar cluster. Such filamentary gas, which might originally be flowing toward a central filament or ridge, was observed for a wide variety of low- to high-mass star formation regions such as Musca, B211/3, and DR 21 (Goldsmith et al. 2008; Schneider et al. 2010; Palmeirim et al. 2013; Cox et al. 2016; Bonne et al. 2020a). The observed small expansion velocity for these filaments around RCW 36 is quite similar to the expansion velocity in the ring, which suggests that overdensities, which are clearly detected in CO, experience significantly slower expansion driven by stellar feedback. Lastly, no corresponding bright filamentary [C ii] emission at 10 km s⁻¹ is observed toward the western cavity, which might be the result of the molecular cloud morphology around RCW 36 before the onset of stellar feedback.

Zoom In Zoom Out Reset image size

The spectra in Section 3.2 showed the presence of blue- (v_LSR < 0 km s⁻¹) and redshifted (v_LSR > 14 km s⁻¹) high-velocity [C ii] wings. These high-velocity wings have a brightness up to 5–10 K in individual spectral bins (see Figure 2), which is well above the detection limit for [C ii] with a noise rms of 0.8–1.0 K. In Figure 10 the integrated intensity map of the high-velocity wings shows that the blue- and redshifted wings display a bipolar distribution, similar to what is observed for protostellar outflows in CO (e.g., Bachiller 1996; Arce et al. 2007). However, for RCW 36, there are no known driving protostars or young stellar objects (Ellerbroek et al. 2013a) at the location where the blue- and redshifted outflow wings come together in the plane of the sky. Furthermore, the high-velocity wings are only detected in [C ii]. Even at the high column density ring the wings are at no point detected in ¹²CO (see spectra C1–C3 in Figure 2), and higher angular resolution observations (A. Bij et al. 2022, in preparation) indicate that there are no protostellar sources in the ring that could drive such powerful mass ejection. The observed [C ii] wings are thus most likely not of protostellar nature.

Zoom In Zoom Out Reset image size

In the currently covered area of RCW 36, there is a significant disparity between the blue- and redshifted integrated wings. Even when integrating over the full 14–22 km s⁻¹ interval, which contains contamination from [¹³C ii], the redshifted wing remains 2 times less bright and covers a significantly smaller area than the blueshifted wing. These high-velocity wings also do not appear to be directly associated with expanding shells, as they do not show an expanding shell velocity structure in the PV diagram in Figure 6 (i.e., an elliptical velocity profile). The brightest part of the blueshifted wing can be found relatively close to the two reported HH objects in Ellerbroek et al. (2013b). However, Figure 10 also shows the position angle of the HH jets, which indicates that the wings seen in [C ii] are not associated with these HH jets. Lastly, Figure 10 shows that the brightest wings are located around the dense molecular ring and toward the cavity walls. This suggests that the wings might be the result of gas that is entrained from the ring and cavity walls by the impact of feedback processes.

Both cavities have blueshifted expanding [C ii] shells. In contrast to the shells in previously studied bubble H ii regions (e.g., Pabst et al. 2020; Luisi et al. 2021; Tiwari et al. 2021), the analysis for RCW 36 is a bit more involved, as it has two expanding cavities and an expanding molecular ring around the cluster. There are two approaches to estimate the expansion energetics. One is based on the Herschel data, and the second one is based on the [C ii] data.

First, we will work with the Herschel data of Vela C. The column density map produced from SED fitting to the 160–500 μm data might not be ideal to calculate the shell mass, as most of these wavelengths particularly trace the cold dust (<20–30 K). Furthermore, the shells in the bipolar cavities are prominent at 160 μm and shorter wavelengths, which could be expected for the strongly irradiated expanding shells if they have a relatively low A_V (≲5). In this case, the shell is mostly dominated by the presence of PDRs, which would lead to higher temperatures. A method that would be more sensitive to relatively warm shells was presented in Tiwari et al. (2021), using a column density map produced based on the Herschel 70 and 160 μm map. This makes use of the fact that the optical depth at 160 μm (τ₁₆₀) can be converted to a hydrogen nucleus column density N_H using the dust extinction cross section at this wavelength, which is given by C_ext,160/H = 1.9 ×10⁻²⁵ cm²/H (Draine 2003). The optical depth is related to the observed intensity by I_ν = B_ν (T_d )τ_ν , assuming optically thin emission and limited fore/background emission. B_ν (T_d ) is the Planck function at a specific dust temperature T_d . In order to calculate the dust temperature map, the method presented in Peretto et al. (2016) and Bonne et al. (2022), based on the intensity ratio of 70 μm (after smoothing to the resolution of the 160 μm data) and 160 μm, is used. From this dust temperature map, it is then possible to calculate the optical depth and the column density (Draine 2003).

Based on a visual inspection of the Herschel and WISE data at 160 μm and shorter wavelengths, the limb-brightened part of the shell in the cavities corresponds to a width of ∼0.3 pc, which is ∼30% of the 1.0 pc radius for both cavities; see Figure 11. This figure also defines the area of the expanding ring. The estimated masses of the ring and shells are presented in Table 2. To obtain the full mass of the expanding shells in the cavities, the calculated mass in the limb-brightened shells is multiplied by a geometric correction factor of 2.7. This correction factor was calculated using a code that determines the fraction of the shell located in the limb-brightened part of the shell in a 1D cut through the spherical half-shell and then assumes an axial symmetry. Using an expansion velocity of 1.0–1.9 km s⁻¹ for the ring and 5.2 km s⁻¹ for the two expanding shells in the cavities, we then arrive at a total kinetic energy of (2.4–2.6) × 10⁴⁷ erg; see Table 2.

Zoom In Zoom Out Reset image size

We note that the expanding central ring contains dust that is also detected at wavelengths >160 μm. When using the mass for the ring from the column density map produced by the SED fitting between 160 and 500 μm, the mass agrees roughly within a factor 2 with the method based on the 160 and 70 μm ratio (i.e., the mass of the ring obtained from SED fitting is about a factor 2 lower; see Table 2). In fact, the results for the bipolar cavities also agree within a factor 2. Identifying the best column density map to calculate the contribution from the central ring is a nontrivial issue, but these results indicate that an uncertainty of at least a factor 2 should be taken into account.

[C ii] is a direct tracer of the expanding shells at velocities between 0 and 5 km s⁻¹ in the cavities and can thus be used to calculate its mass and kinetic energy. Note, however, that currently only 50% of the planned region is observed in [C ii]. Therefore, we apply a correction factor to the calculations. The correction factor is calculated by determining the fraction of the 70 μm emission from RCW 36 covered by the current [C ii] map. Using different minimal 70 μm intensities (i.e., 500–1500 MJy sr⁻¹) to define RCW 36, we find that typically 47%–55% of RCW 36 was covered. Therefore, we apply a correction factor of 2.0 ± 0.2. In order to calculate the C⁺ column densities in the expanding shell, we need the excitation temperature. Generally, this is calculated using the [¹³C ii] line (e.g., Guevara et al. 2020). However, this is not possible for RCW 36 at the moment because there is no convincing [¹³C ii] detection corresponding to the emission between 0 and 5 km s⁻¹. Therefore, we calculate the kinetic temperature of the expanding shell from estimating the FUV field strength at the location of the expanding shells and converting this into a C⁺ kinetic temperature using the PDR Toolbox (Kaufman et al. 2006; Pound & Wolfire 2008). This method is further specified in Appendix C and results in an estimated surface temperature range of 100–500 K, which is a good proxy for the gas temperature that is traced by [C ii] emission and in agreement with work toward other H ii regions (e.g., Schneider et al. 2018; Pabst et al. 2020). Assuming optically thin [C ii] emission for the shell, which is currently justified by the nondetection of [¹³C ii] emission at velocities lower than 6 km s⁻¹, we use Equation (26) from Goldsmith et al. (2012) to calculate the C⁺ column density. We assume that [C ii] emission is thermally excited and work with a ratio C_ul/A_ul = 2 (C_ul is the collision rate and A_ul the Einstein-coefficient for radiative excitation) because H and H₂ number densities are expected to be around (3–6) × 10³ cm⁻³ in the shells. Varying this ratio between 1 and 5, the resulting column density variations remain within 30% at each temperature, which is not significantly larger than other uncertainties such as the temperature (leading to variations up to 50%) or the optically thin assumption. Decreasing the ratio of Einstein coefficients down to C_ul/A_ul = 0.3 would increase the mass by a factor 2. The column density map for the [C ii] emission between 0 and 5 km s⁻¹ is presented in Figure 12. Integrating the column density map, using the value C/H = 1.6 × 10⁻⁴ from Sofia et al. (2004), the total masses of the expanding shells in the cavities are 1.5 × 10² M_⊙, 1.1 × 10² M_⊙, and 1.0 × 10² M_⊙ at a mean kinetic temperature (T_kin) of 100, 250, and 500 K, respectively. Note that this conversion factor might give a lower limit for the total column density since the expanding shell is also weakly detected in the ¹²CO (3–2) line and the densities toward the shell of RCW 36 are far higher than the typical densities reported in Sofia et al. (2004), where C/H = 1.6 × 10⁻⁴ is mostly found to be an upper limit for the higher densities in their sample. These uncertainties might be part of the explanation why these masses, and resulting expansion energetics, are significantly lower than the values obtained with the Herschel methods. They are about a factor 4 smaller than the values obtained using the SED-fitted column density map, which gives the lowest values for the Herschel-based calculations. Another reason for this difference might be the spherical assumption when correcting the Herschel limb-brightened shells to a total mass. Recently, Kabanovic et al. (2022) demonstrated with detailed radiative transfer and self-absorption analysis that the expanding shell in RCW 120 is expanding in a flattened, and not a spherical, molecular cloud. A similar point might be valid for the cavities in RCW 36, which would reduce the mass estimated based on the Herschel data.

Zoom In Zoom Out Reset image size

At velocities around −8 to 0 km s⁻¹ and 14–22 km s⁻¹ we observe [C ii] emission features that are not associated with expanding shell(s). These blue- and redshifted high-velocity wings rather resemble what is observed for protostellar outflows in CO. However, there are no indications of classical outflow wings in ¹²CO (3–2) at any location in the map. In addition, [C ii] emission is not a typical tracer for protostellar outflows. Nevertheless, we perform a first-order approximation to calculate the mass and momentum flux of these bipolar [C ii] high-velocity wings using the method presented in Cabrit & Bertout (1992) and Beuther et al. (2002) for bipolar molecular outflows. ¹⁴ We first calculate the C⁺ column density maps associated with the high-velocity wings with the same method used in the previous section for the shells. Again, the resulting C⁺ mass associated with the wings depends on T_kin and C_ul/A_ul. T_kin, between 100 and 500 K, is again quite well constrained from the PDR models. C_ul/A_ul is less constrained, as the density and collision partners of C⁺ in these high-velocity wings are uncertain. Therefore, taking C_ul/A_ul ∼ 2 will give a reasonable view on the plausible mass range (within a factor 2). The resulting C⁺ column density maps of the velocity intervals from −8 to 0 km s⁻¹ (blueshifted wing) and from 14 to 22 km s⁻¹ (redshifted wing) are presented in Figure 13. In the redshifted wing a correction of 1.5 K km s⁻¹ is made to take into account the contribution from the [¹³C ii] line at velocities above 17 km s⁻¹. Measured with respect to the central velocity of 7 km s⁻¹, these velocity ranges give an observed velocity of 15 km s⁻¹ for both the blue- and redshifted wings, respectively. Since the driving mechanism, and thus the inclination, of the mass ejection is uncertain, we work with an average angle of 57° (Bontemps et al. 1996), which also seems reasonable for mass flows in the cavities that are close to the plane of the sky. From the integrated intensity map of the high-velocity wings, we get a radial extent of 1.5 pc. Taking again into account that currently only half of the cavity is observed, the obtained values are corrected by a factor 2.0 ± 0.2. The resulting mass, mass ejection rate, and momentum flux of the bipolar high-velocity wings, assuming a kinetic temperature of 100–500 K for [C ii], are presented in Table 3.

Zoom In Zoom Out Reset image size

Table 3. The Calculated Mass, Momentum Flux (F_m ), and Mass Ejection Rate (${\dot{M}}_{\mathrm{out}}$) of the Bipolar Wings Observed with [C ii] in the RCW 36 Cavities for Kinetic Temperatures of 100, 250, and 500 K

Tkin	Mass	Fm	${\dot{M}}_{{out}}$
100 K	4 × 102 M⊙	2 × 10−2 M⊙ km s−1 yr−1	7 × 10−4 M⊙ yr−1
250 K	3 × 102 M⊙	1 × 10−2 M⊙ km s−1 yr−1	5 × 10−4 M⊙ yr−1
500 K	3 × 102 M⊙	1 × 10−2 M⊙ km s−1 yr−1	4 × 10−4 M⊙ yr−1

Note. Based on arguments presented in the text, the errors, which are difficult to quantify for all sources of uncertainty, on these values are easily a factor 2.

Download table as:  ASCIITypeset image

The diffuse X-ray emission is presented in Figure 14. The brightest peak is centered on the two O-star candidates of the young cluster. Fainter emission is seen inside the bipolar cavities, which is explained by extinction from the cavity walls when analyzing the X-ray spectra (see Appendix D). The morphology of the diffuse X-ray emission in Figure 14 hints that the hot plasma formed in the OB cluster is leaking out of the H ii region into the ISM beyond the cavities of RCW 36.

Zoom In Zoom Out Reset image size

To constrain the excitation conditions of the diffuse X-ray emission, X-ray spectra of several regions in RCW 36 are fitted using plasma models in XSPEC (Arnaud 1996). Five diffuse extraction regions, shown in Figure 14, were defined to obtain spectra with sufficient S/N for fitting. X-ray events from hundreds of X-ray point sources were excised from these regions, to minimize their influence on the diffuse X-ray spectra. Details on our methods for fitting diffuse X-ray emission in and around massive star-forming regions are given in Townsley et al. (2019).

The extracted diffuse X-ray spectra and fitting results are presented in Appendix D and Table 8. The central region is filled by a soft, single-temperature thermal plasma plus residual emission from unresolved pre-main-sequence stars in the young cluster. The eastern and western cavities are filled with a thermal plasma of similar temperature. The intrinsic surface brightness of this emission, which corrects for the intervening absorption, is higher than that outside the cavities, confirming that the morphology of the apparent surface brightness map is influenced by intervening absorption. Hot plasma in the ISM surrounding RCW 36 shows two thermal plasma components, one similar in temperature to the RCW 36 emission and one substantially hotter. While the hotter component is likely unrelated background emission (often seen along Galactic plane sight lines), the cooler component is strong evidence that the plasma generated by massive star winds in RCW 36 is leaking out of its H ii region into the surrounding molecular cloud.

From the fitted spectra of the hot plasma, the physical properties (plasma density, pressure, etc.) can be calculated with an assumed emissivity curve and a filling factor that was taken to be unity (Townsley et al. 2003). The assumed emissivity curve from Landi & Landini (1999) assumes collisional ionization equilibrium (CIE) for the plasma. This assumption is justified since nonequilibrium effects only become notable below plasma temperatures of 10⁶ K (Vasiliev 2011, 2013), whereas the hot plasma temperature in RCW 36 ((1.7–2) × 10⁶ K) is above this value; see also Table 8. The calculation of the physical properties was not done for the regions outside the bipolar H ii region (north and southeast), as there is no evident assumption for the 3D geometry at these locations. This calculation of the physical properties also was not done for the central region directly around the OB cluster because of the nonequilibrium ionization (NEI) of the plasma and relatively low temperature. As a result, the bandpass correction factor necessary for the limited ACIS bandpass is larger than 40, which would make the results for the center highly uncertain. For the cavities, the bandpass correction factors are 4.4 and 3.4 for the east and west, respectively, which makes their estimated properties less uncertain. The resulting physical properties of the plasma are presented in Table 4. Both cavities have an estimated temperature around 2 × 10⁶ K and electron densities around 0.2 cm⁻³, which results in an estimated pressure of (0.6–0.9) × 10⁶ K cm⁻³.

The resulting estimated energy of the plasma in both cavities matches the expansion energy of the expanding shells, presented in Table 2, well (i.e., it is slightly higher than the estimated expansion energy based on [C ii] and slightly lower than the estimates based on Herschel). This suggests that under adiabatic conditions the expanding shells in the cavities might be driven by the hot plasma (e.g., Weaver et al. 1977; Lancaster et al. 2021a). This will be examined in more detail later on, as other processes such as photoionized gas pressure and direct radiation pressure could also play a role. Note that the unity filling factor is an assumption. The appropriate scaling factor (η^−1/2, with η the filling factor) to correct for a lower filling factor is also presented in Column (2) of Table 4. This indicates that the resulting plasma properties are not extremely sensitive to the filling factor.

RCW 36 consists of a bipolar cavity and a dense central ring surrounding the waist of the H ii region created by the central OB cluster. In the cavities, the expanding [C ii] shell is only observed at blueshifted velocities, which implies that it is expanding toward us. Fitting the foreground column to the Chandra data toward the cavities, in particular the eastern one, suggests that most column density, observed with Herschel, is associated with and on the near side of RCW 36 from our point of view. The observed [C ii] shells are thus expanding into the predominant nearby column density reservoir of RCW 36. This can offer an explanation why no corresponding redshifted expanding shell is detected for the cavities even though a few expanding filamentary structures are observed toward slightly more redshifted velocities. As most of the column density appears to be located in front, the gas behind the ridge might be mostly diffuse gas. As a result, the hot plasma created by the OB stars would be able to leak from RCW 36 in that direction, as sketched in Figure 15. This is in essence the same scenario observed for RCW 49 in Tiwari et al. (2021). Interestingly, in contrast to RCW 36 (two O-star candidates and ∼1 Myr), RCW 49 contains 37 OB identified stars and is slightly older with an estimated age ∼2 Myr.

Zoom In Zoom Out Reset image size

To obtain a first idea of the expansion timescale for the dense molecular ring, we use the estimated expansion velocity between 1.0 and 1.9 km s⁻¹ and an estimated radius of 0.9 pc. As it is unclear what drives the expansion of the ring, we calculated the Spitzer expansion timescale for expansion driven by ionizing radiation (Spitzer 1978), the Weaver solution for expansion driven by stellar winds (Weaver et al. 1977), and the ballistic expansion timescale that can be expected if shell expansion is driven by stellar winds in an isothermal density field with a power-law index of −2 (Geen & de Koter 2022). The Spitzer solution gives an expansion timescale of 1.2–1.7 Myr, the Weaver solution gives a timescale of 0.3–0.6 Myr, and the ballistic expansion gives a timescale of 0.5–0.9 Myr. The calculations of the Spitzer and Weaver solutions are further specified in Appendix F. Broadly speaking, the Spitzer and ballistic solutions for the dense ring appear to be in better agreement with the cluster age than the Weaver solution (Ellerbroek et al. 2013a). However, the upper value of the Spitzer solution, which might be the most representative one (see Appendix F), has an equally large deviation from the estimated cluster age as the Weaver solution. The cavities, on the other hand, have a radius of ∼1.0 pc and an expansion velocity of 5.2 ± 0.5 ± 0.5 km s⁻¹, which results in an estimated Spitzer expansion timescale of 0.2–0.4 Myr (see Appendix F), a Weaver expansion timescale of 0.1 Myr, and a ballistic expansion timescale of ∼0.2 Myr. This is roughly a factor 3–5 shorter than for the ring and suggests that the expansion in the cavities is a recent feature. Note that this assumes that the cavity expands from the center of the cavity, which might not be fully correct. Considering that the expansion starts from the center of RCW 36, we get a radius in the range of 1.5–2.0 pc. This results in an estimated ballistic expansion timescale of 0.3–0.4 Myr, which is still shorter than the expansion of the dense ring. RCW 36 thus shows inhomogeneous expansion: relatively slow expansion in the dense gas, followed by increasing expansion velocity once the expansion reaches the lower-density gas. This fits with relatively low (<3 km s⁻¹) expansion velocity of the denser filamentary structures perpendicular to the ring.

Samal et al. (2018) noted that only 16 of the 1377 H ii regions considered in their sample are visibly identified as bipolar H ii regions. A correction factor has to be taken into account based on inclination selection effects, but, as mentioned in Samal et al. (2018), this still implies that bipolar H ii regions are rather rare. From the faster expansion in the bipolar cavities of RCW 36 compared to the ring, it can be argued that the morphology of RCW 36 is changing rapidly (on ∼0.2 Myr timescales) and thus that a H ii region might only be observed as bipolar for a part of its evolution. As a short exercise, assuming a correction factor 5 for the inclination ¹⁵ and a potential correction factor 5 for the short timescale ¹⁶ would imply that up to 30% of the H ii regions might go through a bipolar phase. However, this remains speculative at the moment since this is based on a single case. More generally, the observation of expanding half-shells, instead of fully spherical shells, in [C ii] emission toward H ii regions (e.g., Anderson et al. 2019; Pabst et al. 2020; Luisi et al. 2021; Tiwari et al. 2021) and the evidence of important hot plasma leakage can fit in the view of the ISM where (massive) stars form in a complex organization of filaments and sheets. The creation of bipolar H ii regions can then occur when an O star forms close to the midplane of a sheet-like structure such that the radiation and wind break out on both sides at about the same time. The hot plasma can then start driving expanding shells in the lower-density gas outside the sheet. This is, however, speculative since it is based on a single case, while the observations also indicate that the expansion strongly depends on the density structure. Furthermore, simulations by Wareing et al. (2017) suggest that the bipolar morphology remains present for a significant part of their simulation. This is due to the apparent low expansion velocity (<1 km s⁻¹) for the cavities in their simulation, which is not in agreement with the observed velocities in the cavities for RCW 36. As the expansion dynamics might vary from region to region, a detailed analysis of age and dynamics of a sample of bipolar H ii regions will be required to fully address this question.

Analysis of the Chandra data shows that a hot plasma, created by the stellar winds, pervades RCW 36 and that it has the required energy to drive the observed expanding shells in both cavities. To discuss the role of different physical processes, the associated pressure terms were estimated and are presented in Table 5. Specific assumptions made for these calculations are specified in Appendix E. The thermal pressure from the hot plasma in both cavities is ∼(0.6–1.0) ×10⁶ K cm⁻³, which is slightly larger than the estimated thermal pressure of the photoionized gas, i.e., (4.6–5.0) × 10⁵ K cm⁻³, derived from the SHASSA Hα observations (Gaustad et al. 2001). However, as for the hot plasma, if the photoionized gas is confined in a subregion of the cavities, this will increase the estimated pressure. As a result, these terms probably dominate over the direct radiation pressure in the cavities; see Table 5. Since cavities are filled by both the hot plasma and photoionized gas, the photoionized and collisionally ionized gas can thus be in pressure equilibrium across the contact discontinuity.

The physical processes that drive expanding motion of H ii regions are still a matter of debate. Simulations generally predict that ionizing radiation is the main driver of expansion compared to stellar winds (e.g., Haid et al. 2018; Geen et al. 2021; Ali et al. 2022). Yet simulating the high dynamic range of temperatures and densities associated with the presence of stellar winds in the ISM is challenging (e.g., Dale 2015). Observations by Lopez et al. (2014) also indicate that the ionized gas pressures appear to dominate most H ii regions; however, other observational studies (e.g., Pabst et al. 2019; Luisi et al. 2021; Tiwari et al. 2021) have proposed that stellar winds drive the high-velocity expansion detected with spectrally resolved [C ii] observations. Above we already noted that thermal pressure of the ionized gas and hot plasma seem to be of the same order for RCW 36. Here, we will discuss the associated energy terms. Based on the SHASSA Hα intensity, the ionizing photon flux in each cavity is ∼8.5 ×10⁴⁶ s⁻¹; see Appendix E. This results in a total injected energy of ∼6.2 × 10⁴⁹ erg over the estimated lifetime of the cluster (i.e., 1.1 Myr) when assuming an energy of 13.6 eV for the ionizing photons. Taking into account that the coupling efficiency of ionizing radiation energy to radial expansion typically is <10⁻⁴ (Haid et al. 2018), this is about an order of magnitude too small to explain the expanding motion of the shells in the cavities. We do note that at very early evolutionary stages in Haid et al. (2018) the coupling efficiency for ionizing photons is significantly higher. However, considering a coupling efficiency of 10⁻³ over the first 0.3 Myr remains insufficient to explain the observed expansion. Therefore, unless the coupling efficiency for ionizing photons in Haid et al. (2018) is about an order of magnitude lower than in RCW 36, we propose that the hot plasma created by stellar winds drives the expansion of the blueshifted shells in the cavities.

On the other side of the shell, ram pressure, thermal pressure, and magnetic pressure from the ambient cloud confine the shell. The values in Table 5 indicate that in this low-density region the ram pressure confines the shell. Inside the compressed shell, turbulent pressure seems to be the dominant term to prevent the shell from being further compressed between the hot plasma on one side and ram pressure on the other side. Even though the estimated magnetic pressure in the shell is significantly smaller than the turbulent pressure, it might play a role in limiting the initial compression rate of the propagating shock that is driven by the expansion. Inspecting the BLASTPol magnetic field map toward RCW 36 in Figure 16, it is observed that the magnetic field morphology follows the bipolar cavities and forms a kink near the ring where the polarization fraction drops. This behavior where the magnetic field orientation is aligned with the shell morphology around a H ii region, as it is dragged along by the expansion, was recently also observed in the Keyhole Nebula (Seo et al. 2021) and NGC 6334 (Arzoumanian et al. 2021). The magnetic field morphology in the ring and whether it significantly affects the expansion remain unclear with current data. This requires higher angular resolution observations (A. Bij et al. 2022, in preparation). In the ring, it is also unclear which process (ionizing radiation or the hot plasma) drives the expansion. The number of ionizing photons and coupling efficiency is sufficient to drive the ring expansion, but the hot plasma in the ring remains unconstrained. As a result, we cannot reach a firm conclusion.

Zoom In Zoom Out Reset image size

In Figure 17, a noteworthy increase in ¹³CO (3–2) velocity dispersion is observed in the ring of RCW 36 compared to the surrounding regions of the Vela C molecular ridge, which does not seem to be explained by opacity broadening; see Appendix E. This indicates that the mechanical energy input from stellar feedback can significantly increase the turbulent energy in the compressed gas. This energy injection might thus play a role in setting the thickness of the shell, as it provides a rough pressure equilibrium with both sides of the shell. Specifically, the turbulent pressure in the ring, derived from the spectral line width, is similar to or higher than the ram pressure provided by the ambient gas; see Table 5. Otherwise, the shell would be further compressed until an increased temperature (or magnetic field strength or turbulence) due to this compression would halt it. From the information provided by the ¹³CO (2–1) line width, corrected for optically thick line broadening (Appendix E), it is possible to estimate the turbulence injection by stellar feedback in the molecular ring. The calculated total turbulent energy injection in the ring and cavity shells is given in Table 6. From Figure 17, it is reasonable to assume a line width for this gas of 0.7 km s⁻¹ before the onset of stellar feedback, which would give a total turbulent energy of 2.5 × 10² M_⊙ (km s⁻¹)² for the gas now located in the ring. Comparing this with the value for the ring in Table 6, this implies that the turbulent energy of the gas, after opacity broadening correction, has roughly increased by more than 60% compared to the onset of stellar feedback. Turbulence typically dissipates over a few crossing timescales (t_cross = d/σ_los; Ostriker et al. 2001). For the ring, with a width of ∼0.1 pc, this gives a crossing timescale of 1.1 × 10⁻¹ Myr. Using this timescale results in a turbulence dissipation rate of 2.4 × 10³³ erg s⁻¹ that has to be injected to maintain the turbulent line width. This is only a small fraction of the probable energy injection rate by the O-star candidates, which will be discussed in more detail in the next section.

Zoom In Zoom Out Reset image size

Table 6. The Turbulent Energy (E_turb), Dissipation Timescale (t_diss), Turbulent Dissipation Rate $({\dot{E}}_{\mathrm{turb}})$, and Total Injected Turbulent Energy (E_turb,tot) for the Dense Ring and the Full Expanding Shell Associated with the Bipolar Cavities

	Eturb (erg)	tdiss (Myr)	${\dot{E}}_{\mathrm{turb}}$ (erg s−1)	Eturb,tot (erg)
Ring	8.2 × 1045	1.1 × 10−1	2.4 × 1033	8.2 × 1046 (over 1.1 Myr)
Shells (from [C ii])	2.2 × 1046	2.0 × 10−1	3.2 × 1033	2.2 × 1046 (over 0.2 Myr)

Note. These calculations for the shell make use of the typical velocity dispersion of 1.5 km s⁻¹ in the expanding shell and the shell mass deduced from [C ii].

Download table as:  ASCIITypeset image

In Section 5.1 we proposed that there is leakage of the hot plasma through the diffuse medium around RCW 36, which leads to the hot plasma found outside RCW 36 in Figure 14. To quantify the leakage, we have to estimate the total energy injected in the ISM by the stellar winds of the OB cluster. This is done using several stellar atmosphere models, as there is a large uncertainty on the actual stellar wind mass-loss rates, which from observations appear to be lower than predicted by the standard models (e.g., Martins et al. 2004, 2005b; de Almeida et al. 2019). However, recently it was proposed that the observed values are underestimated because observations miss a significant part of the stellar winds, as a large fraction of the stellar wind might be in a hot phase that is not traced by the observations (Lagae et al. 2021). To cover the range of possibilities, we will use the relations that were put forward in Equation (12) of Vink et al. (2000), Equation (3) of Lucy (2010), Equation (11) of Krtička & Kubát (2017), and Equation (20) of Björklund et al. (2021), all with a terminal velocity of 2500 km s⁻¹ (e.g., Lamers et al. 1995; Krtička & Kubát 2017; Lagae et al. 2021). The formulae differ in the choice of critical parameters and their dependencies (luminosity, metallicity, effective temperature, terminal velocity, etc.), but it is beyond the scope of this paper to discuss their differences in more detail. In any case, employing a range of relations allows us to cover the standard predicted mass ejection rate from stellar models (Vink et al. 2000) to predictions that are more adapted to fit with the current observations (Björklund et al. 2021). For the 09V and 09.5V star candidates in the cluster we use the stellar parameters from Martins et al. (2005a) and Pecaut & Mamajek (2013) to take into account the effects of different stellar models. The results are summarized in Table 7, which indicates a significant change depending on the stellar model used and an even larger effect related to the stellar wind mass-loss rate prescription used. Only taking into account the two O-star candidates of RCW 36, this results in a total injected energy in the range of 4.7 × 10⁴⁷ erg to 8.3 × 10⁴⁸ erg over 1.1 Myr; see Table 7. This order-of-magnitude difference depending on the model might be overcome by the missing hot phase proposed in Lagae et al. (2021); however, future observations will have to address whether there is indeed such a large fraction of the stellar winds in the hot phase. The energy in the hot plasma of RCW 36, from Table 4, is of the order of (1–2) × 10⁴⁷ erg, implying that >50%, and potentially up to 97%, of the total injected stellar wind energy is no longer found in the hot plasma.

Table 7. The Stellar Wind Mass Ejection Rates (${\dot{M}}_{\mathrm{SW}}$) and Total Ejected Stellar Wind Energy (E_SW)

		Vink et al. (2000)	Lucy (2010)	Krtička & Kubát (2017)	Björklund et al. (2021)
Model	Star	log[${\dot{M}}_{\mathrm{SW}}$ (M⊙ yr−1)]	log[${\dot{M}}_{\mathrm{SW}}$ (M⊙ yr−1)]	log[${\dot{M}}_{\mathrm{SW}}$ (M⊙ yr−1)]	log[${\dot{M}}_{\mathrm{SW}}$ (M⊙ yr−1)]
Martins et al. (2005a)	O9V	−7.35	−8.56	−7.78	−8.31
Martins et al. (2005a)	O9.5V	−7.57	−9.00	−7.94	−8.53
Pecaut & Mamajek (2013)	O9V	−7.12	−8.56	−7.61	−8.10
Pecaut & Mamajek (2013)	O9.5V	−7.35	−8.58	−7.78	−8.31
		ESW (erg)	ESW (erg)	ESW (erg)	ESW (erg)
Martins et al. (2005a)	O9V	2.7 × 1048	1.7 × 1047	1.0 × 1048	3.0 × 1047
Martins et al. (2005a)	O9.5V	1.6 × 1048	5.7 × 1046	6.7 × 1047	1.7 × 1047
Martins et al. (2005a)	total	4.3 × 1048	2.3 × 1047	1.7 × 1048	4.7 × 1047
Pecaut & Mamajek (2013)	O9V	5.4 × 1048	1.9 × 1047	1.7 × 1048	5.4 × 1047
Pecaut & Mamajek (2013)	O9.5V	2.9 × 1048	1.7 × 1047	1.1 × 1048	3.1 × 1047
Pecaut & Mamajek (2013)	total	8.3 × 1048	3.6 × 1047	2.8 × 1048	8.5 × 1047

Note. Bottom: the resulting total ejected energy over the lifetime of the cluster is calculated assuming an age of 1.1 Myr and using a terminal wind velocity of 2500 km s⁻¹.

Download table as:  ASCIITypeset image

Stellar wind energy-loss methods were listed in Townsley et al. (2003), Harper-Clark & Murray (2009), and Rosen et al. (2014). Here we will explore the options. From Table 4, it is found that the hot plasma luminosity is a factor 3 to more than an order of magnitude lower than required to explain the missing energy in the hot plasma. The injected turbulent energy in the molecular ring, presented in Table 6, is a factor 5 to two orders of magnitude smaller than the total injected energy as well. Thermal conduction at the interface of the hot plasma and the turbulent shell layer is another potential way to lose the plasma energy. Recently, work by Lancaster et al. (2021a, 2021b) indicated that cooling can be very efficient at the interface because of turbulent substructure and that this could even explain all the observed missing hot plasma energy. However, it was also noted that thermal conduction is strongly reduced when a (weak) magnetic field is parallel to the shell (Soker 1994; Rosen et al. 2014). This is the case for RCW 36 based on Figure 16, and it should be noted that Lancaster et al. (2021a, 2021b) currently do not include the magnetic field. Simulations in Arthur et al. (2011) indicate that the magnetic field will be mostly perpendicular to the shell in the ionized gas and aligned with the shell in the neutral gas. Rosen et al. (2014) note that a parallel magnetic field with 1 μG would decrease the thermal conduction by several orders of magnitude. In RCW 36, a magnetic field strength of at least several μG is expected (Crutcher 2012), which seems to exclude thermal conduction based on this work. However, more in-depth theoretical studies of energy loss at the turbulent interface, including the magnetic field, would be invaluable to better constrain this. Lastly, wind energy transfer to dust grains through collisions was also excluded in Rosen et al. (2014) as a plausible explanation. This suggests that the majority of the stellar wind energy might have leaked through lower-density regions out of RCW 36. Therefore, we here try to estimate the plasma energy leakage from RCW 36 based on Equation (8) from Harper-Clark & Murray (2009). As we propose that the region is fully open, at one side of the region we use a shell cover fraction (C_f ) of 0.5 for RCW 36 and a radius of 1 pc. Plugging in these values with the hot plasma parameters from Table 4 gives a current energy leakage per cavity of −2.0 × 10³⁵ erg s⁻¹, which is almost two orders of magnitude larger than derived values from other processes in this work. Since leakage probably occurred over the evolution of the H ii region, we use an intermediate radius of 0.5 pc to estimate the leakage over the lifetime of RCW 36. Note that this is only an approximation, as we have no detailed information on the time evolution of RCW 36. It should also be noted that the C_f definition in Harper-Clark & Murray (2009) assumes equal leakage in all directions, whereas in RCW 36 we propose that leakage has a preferential direction. This then gives a total leaked plasma energy of −1.7 × 10⁴⁸ erg over 1.1 Myr. Compared to the values in Table 7, this suggests that leakage indeed has the potential to account for most missing hot plasma energy.

If there is indeed such strong hot plasma leakage, it could be wondered whether the leaking plasma still has an effect on the molecular cloud. Studying the large-scale molecular cloud structure around RCW 36 with sensitive dust continuum observations (i.e., Herschel and BLASTPol), it can be noted that RCW 36 is located at the center of a large, low column density region in Vela C (Hill et al. 2011; Fissel et al. 2016); see Figure 18. Furthermore, inspecting Figure 18, we observe that the leaking hot plasma shows hints that it correlates with the structure of this large, low column density region. It is thus possible that the potential large fraction of leaking hot plasma clears the low-density ambient cloud on a larger scale. When the injected feedback energy escapes into the diffuse ISM through such lower-density regions or "chimneys" in the cloud, it can thus further disrupt these regions and their mass provision to the ridge around RCW 36. From the large-scale SHASSA Hα map presented in Figure 19, we confirm that this low-density region is indeed a channel of feedback energy leaking into the larger cloud. Estimating the number of photons in these larger cavities is relatively uncertain, as their extent is not strongly constrained. Excluding the Hα emission in the cavities defined by the 70 μm emission, we obtain ionizing photon fluxes as high as (6–7) × 10⁴⁶ s⁻¹ in these larger cavities. This is fairly similar to the values obtained inside the 70 μm cavities and thus supports significant leakage from RCW 36 into the diffuse cloud.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As Vela C still is a relatively young molecular cloud where RCW 36 is the first H ii region, this suggests that the leakage is the result of the molecular cloud morphology before feedback sets in. This presence of lower-density gas (or voids) is expected if molecular clouds are indeed built out of filaments and sheets.

The bipolar wings at high velocities in [C ii] show no clear expanding shell morphology and are not driven by protostars in the ring, but rather by feedback from the OB cluster. This could explain why the wings become less clear at the high column density region of the ring in Figure 10. [C ii] high-velocity wings were also found in the bipolar regions S106 (Simon et al. 2012) and NGC 7538 IRS1 (Sandell et al. 2020). Unlike S106 (Schneider et al. 2018), [O i] is not clearly detected in the wings of RCW 36, which suggests that the density in these wings is comparatively low. In NGC 7538 IRS1, the wings are also detected in ¹²CO (3–2), which suggests that for NGC 7538 IRS1 the [C ii] wings might trace a protostellar outflow there.

To investigate whether ablation by stellar winds can drive the mass ejection observed with the [C ii] wings, we compare the energy (and momentum) rates for the stellar winds and the [C ii] wings. The energy injection rate of the most favorable stellar wind models (∼(3.9–7.5) × 10⁴² erg yr⁻¹, from Vink et al. 2000) is of the same order as the outflow energy rate seen with the [C ii] wings (∼(3–5) × 10⁴² erg yr⁻¹). As the [C ii] outflow parameters should be accurate within a factor 2–3 according to our evaluation in Section 4.3, this suggests that stellar wind ablation might be able to continuously drive the observed [C ii] wings. If these high-velocity wings are continuously ejected from RCW 36, it might account for a significant fraction of the missing plasma energy. However, it has to be taken into account that we also observe actual leakage of the hot plasma and that there are other potential stellar wind energy sinks discussed in previous sections. In addition, this conclusion is only valid for the Vink et al. (2000) models. The estimated energy injection rate by all other models would not be able to explain the observed wings. When considering the total momentum injection rate of the stellar winds, this is about two orders of magnitude lower than the momentum flux observed for the [C ii] wings.

Another potential driver would be the radiation from the OB cluster. In the simulations by Kim et al. (2018), which focus on radiation feedback in molecular clouds, they found mass ejection with velocities between ∼6 and 40 km s⁻¹ as a result of thermal and radiation pressure from ionizing radiation. Radiation might thus be able to drive the observed mass ejection from the center of RCW 36. The calculated radiation pressure (I/c) at r = 0.9 pc (corresponding to the ring) is 1.2 × 10⁶ K cm⁻³ and thus not significantly higher than the pressure from the hot plasma in the cavities. If the acceleration originates in the ionized gas closer to the OB cluster, the radiation pressure would be higher and could thus more easily drive high-velocity mass ejection into the low-density cavities. As the ionized gas does not cover the full region, because of the hot plasma, this also increases the thermal pressure of the ionized gas. However, there is very little mass inside the central ring (Minier et al. 2013), which implies that the mass ejection should happen at the surface of the ring and/or bipolar cavities. The mass ejection seen with the wings might thus be driven by a combination of stellar wind and pressure from the ionizing radiation. This combination would leave more room for the observed hot plasma leakage and other energy sinks. In general, it thus remains challenging to explain the observed wings with the constrained feedback mechanisms for RCW 36.

The mass ejection rate derived from the [C ii] wings in Table 3 is high. It suggests that (4.8–7.4) × 10² M_⊙ can be removed from the region confined by the ring in RCW 36 over the current lifetime of the OB cluster (or (0.9–1.3) × 10² M_⊙ over 0.2 Myr). This could have important implications for the star formation efficiency (SFE) in the cloud. If this mass ejection is going on for 1.1 Myr, which could fit with the missing hot plasma energy, it would have already removed about the same mass from the original star-forming region as what is now found in the molecular ring and more than an order of magnitude more than the ∼20 M_⊙ of ionized gas confined by the ring (Minier et al. 2013). We found no indications of inflow toward RCW 36 from the ambient gas, as it is actively being swept away, so the main mass inflow would happen along the ridge of Vela C. Using πR² ${v}_{\inf }\rho$, it is then possible to calculate the remaining mass inflow to the RCW 36 ring through the ridge, where R is the radius, v_inf the inflow velocity, and ρ the density. Here we assume that the 1.5 km s⁻¹ velocity gradient along RCW 36 (Fissel et al. 2019) is associated with inflow, though it is equally plausible that this is mainly the result of expansion. This gives ${\dot{M}}_{\inf }$ = 0.7 × 10⁻⁴ M_⊙ yr⁻¹ when we further assume R = 0.1 pc (Minier et al. 2013) and ${n}_{{{\rm{H}}}_{2}}$ = 2.4 × 10⁴ cm⁻³ (converted to the density assuming a mean molecular mass μ = 2.33). The inflow rate is thus more than a factor 5 lower than the mass ejected from RCW 36 and therefore has a small effect on the disruption. The observed mass ejection rate could thus stop any new star formation around RCW 36 on the timescale of 1–2 Myr.

It is also interesting that the mass ejection rate is within a factor two of the typical proposed mass inflow rates (∼10⁻³ M_⊙ yr⁻¹) in globally collapsing hubs and ridges before massive stars have formed (e.g., Schneider et al. 2010; Peretto et al. 2013). This implies that the observed mass ejection can remove the mass from the central region on roughly the same timescale that the collapsing hub/ridge provides mass inflow (i.e., a few freefall timescales of the hub), which is consistent with several simulations that emphasize the importance of stellar feedback to maintain a low SFE in molecular clouds (e.g., Dale et al. 2012; Walch et al. 2012). The observations of RCW 36 thus suggest that stellar feedback is essential in preserving the low SFE in proposed dynamic high-mass star-forming scenarios owing to effective mass removal from the originally collapsing parsec-scale region.

There are still several uncertainties related to this rapid mass removal (when the high-velocity mass removal started exactly, the currently incomplete map, etc.), but future studies using the larger statistical sample that will be provided by the FEEDBACK Legacy program will give further insight into the timescale of this rapid mass removal and how it affects the SFE in high-mass star-forming regions.