Expansion Signatures in 35 H II Regions traced by SOFIA [C II] Emission

Figure 1 presents a schematic diagram of an expanding H II region. The central ionizing source is depicted as a star at the center. Surrounding it is a region of ionized hydrogen (H II), shown as a dark blue circle, where ionized hydrogen is more prevalent than other ionized species. This H II region has an overpressure that drives a shock front into the surrounding molecular cloud, sweeping up gas into a dense shell. The inside of this shell is illuminated by far-UV (6–13.6 eV) radiation, creating a PDR boundary, shown as the outer light blue ring. As the distance from the ionizing source increases, UV photons with energies $E\geq 13.6$ eV are absorbed, but species like carbon, which has an ionization potential below 13.6 eV, can still be ionized. The small dark clouds along the outer edge of the PDR front represent sites where triggered star formation might occur as the PDR expands into the surrounding molecular cloud. A neutral shell is trapped between the ionization front and the shock front.

The time that it takes for an H II region to reach its initial Strömgren radius is given by the recombination timescale, $t_{r}=(\alpha_{B}*n_{0})^{-1}$, where $\alpha_{B}$ is the “Case B” hydrogen recombination coefficient as in Equation 1 (Spitzer Jr, 1968).

As an H II region evolves, this Strömgren sphere expands beyond its initial Strömgren radius. This expansion is driven by increased thermal pressure of the ionized gas or the pressure of the hot plasma created by the stellar wind from the ionizing source(s). Initially, this expansion is faster than the sound speed in the ambient medium, causing a shock front that propagates either in front of, with, or eventually falls behind the ionization front (Franco et al., 2000).

In theory, as the HII region expands, its internal pressure drops and as it nears the pressure of the surrounding medium, causing the shock wave to become subsonic and the expansion to stall (Franco et al., 1989).

Of the 35 H II regions studied in this paper, 12 exhibit expansion in at least one of their PV diagrams ($\sim\!34\%$ of the sample). Figure 2 shows two axes for M17 that show expansion signals, with PV diagram 2 showing only a blueshifted signal, while PV diagram 10 shows both a redshifted and a blueshifted signal. Table 2 shows the global expansion parameters for these 12 expansion candidates (ECs), along with expansion velocity values from the literature (when available). Of the 12 ECs, only G316+796 displays solely redshifted expansion with no blueshifted signature while eight display only blueshifted expansion with no redshifted signature, and three show both redshifted and blueshifted expansion. The average expansion rate for all 12 ECs is $\sim\!12.2$ $\,{\rm km\,s^{-1}}$. The average blueshifted expansion velocity of the ECs is $\sim\!10.9$ $\,{\rm km\,s^{-1}}$ while the average redshifted expansion velocity is $\sim\!13.2$ $\,{\rm km\,s^{-1}}$. The average percent of PV axes showing an expansion signature for all 12 ECs is 63%. As is expected, all four of the ECs that display expansion in 100% of their PV diagrams appear as bubbles in [C II] emission (Table 2). Of the 12 ECs, nine of them are bubbles (75%). On average, 72% of the PV cuts display expansion in the nine ECs that have bubble morphologies. On the other hand, only 34% of the PV cuts in the non-bubble ECs display expansion. Of the entire sample of 35 H II regions, 14 are bubbles (40%; see Table 1).

Table 3 shows the results of the velocity residual map analysis for eight of the 12 ECs where the emission has a bubble morphology, appearing annular in its [C II] moment-0 map, and either a radius on the sky of at least $2.5\arcmin$ or a peak [C II] moment-0 intensity of at least 400 K $\,{\rm km\,s^{-1}}$ (G083+936, G316+796, N19 and RCW 36 are too faint and/or compact). The second column labeled “$n\times m$” gives the binning parameters used for the spatial axes (see Section 4). The $v_{\rm r,map}$ and $v_{\rm b,map}$ columns give the values for the expansion velocity of the best-fit model for redshifted and blueshifted shells, respectively. The associated errors for these values can be inferred from the standard deviation plot shown in the right panel of Figure 7. Using the determined velocity of expansion and the observed S/N, the standard deviation can be determined. These standard deviations can also be used to determine the size of the error bars for the x-axis in Figure 5. The $\langle v\rangle_{\rm r,PV}$ and $\langle v\rangle_{\rm b,PV}$ columns indicate the average expansion values (across 16 PV axes per region) from the position-velocity diagram analysis outlined in Section 5.1 for comparison also for redshifted and blueshifted shells, respectively. PV axes that did not show signs of expansion are excluded from the average, which can lead to artificially inflated expansion values in regions where only a few axes exhibited expansion.

Only two of the ten expansion signatures with associated errors in Table 3 fall within the error bars of their corresponding PV analyses: the redshifted expansion in RCW 79 and the blueshifted expansion for M43; suggesting one-to-one agreement of the two methods for these signals. The remaining nine signatures fall outside their PV-derived uncertainties, but still produced well-constrained expansion models based on the [C II] moment-0 morphology. The average expansion velocity from the ten velocity residual map analyses is 12.8 $\,{\rm km\,s^{-1}}$. Two maps model redshifted shells (M17 and RCW 79) with values of 10.0 $\,{\rm km\,s^{-1}}$ and 14.5 $\,{\rm km\,s^{-1}}$, respectively. The remaining eight maps model blueshifted shells (M17, M42, M43, NESSIE-A, NGC 7538, RCW 120, RCW 79, and W40), with an average modeled blueshifted expansion velocity of 12.2 $\,{\rm km\,s^{-1}}$. To test the robustness of this analysis, we repeat it on simulated data, as described in Section 6.6. Figure 5 shows the expansion velocity derived from our PV analysis on the y axis versus the expansion velocity estimated from the velocity map on the x axis. As indicated by the one-to-one line in Figure 5, the estimated expansion velocities for almost all regions fall within $5\sigma$ for the two approaches considered.

Several previous studies have identified expansion signatures in H II regions, primarily based on velocity offsets between the ionized gas and surrounding PDRs. Our analysis builds on this body of work by providing a uniform comparison across a large sample.

Roshi et al. (2005) observed carbon recombination lines toward 18 ultracompact H II regions and identified dense PDRs in 11 sources. In nine cases, they found a consistent velocity offset of $\sim$3.3 $\,{\rm km\,s^{-1}}$ between the PDR and ionized gas, which they interpreted as evidence of expansion. Kirsanova et al. (2017) studied 14 regions in CS(2–1) and ¹³CO(1–0), detecting a 1.2 $\,{\rm km\,s^{-1}}$ offset in G183.35$-$0.58, which they interpreted as expansion of the H II region into the molecular cloud. In follow-up work, Kirsanova et al. (2020) observed S235 in [C II], interpreting the velocity structure as further evidence of expansion.

Several studies using [C II] observations with SOFIA have reported expansion in individual H II regions. Pabst et al. (2020) reported expansion velocities of $\sim$13 $\,{\rm km\,s^{-1}}$ in M42, $\sim$6 $\,{\rm km\,s^{-1}}$ in M43, and $\sim$1 $\,{\rm km\,s^{-1}}$ in NGC 1977, concluding that stellar winds likely dominate the expansion in M42, while thermal pressure may drive expansion in M43 and NGC 1977. These observations have been extended in subsequent work by Pabst et al. (2021, 2022, 2024). Luisi et al. (2021) observed an expansion velocity of $\sim$15 $\,{\rm km\,s^{-1}}$ in RCW 120 and attributed the expansion to stellar winds. Prior work by Zavagno et al. (2006) suggested that the expansion of RCW 120 may be triggering star formation along the southwestern edge of the region. Tiwari et al. (2021) identified expansion in RCW 49 with a similar velocity ($\sim$13 $\,{\rm km\,s^{-1}}$), supported by evidence of stellar wind-driven expansion from diffuse X-ray emission.

More recent SOFIA observations have reported expansion velocities of $\sim$5 $\,{\rm km\,s^{-1}}$ in RCW 36 (Bonne et al., 2022) and NGC 7538 (Beuther et al., 2022). Bonne et al. (2023b) observed simultaneous blueshifted and redshifted [C II] emission in RCW 79, identifying an expansion velocity of up to $\sim$25 $\,{\rm km\,s^{-1}}$. This represents the first detection of both blueshifted and redshifted expansion signatures in the same region using [C II]. Figueira et al. (2017) suggested that expansion in RCW 79 may be responsible for triggering star formation at the interface between the H II region and the surrounding molecular material.

Finally, Saha et al. (2024) used ALMA and the VLA to observe G24.47+0.49 in HCO⁺(1–0) and radio continuum emission. They identified a 9 $\,{\rm km\,s^{-1}}$ expanding ring of HCO⁺ surrounding the H II region and multiple collapsing dense cores along the expansion front, providing direct observational evidence of triggered star formation.

While these studies provide important context for understanding feedback-driven dynamics, they generally focus on individual regions. In contrast, our sample includes 35 H II regions analyzed using a consistent method, allowing us to examine expansion signatures statistically across a wide range of morphologies and evolutionary stages. In several cases, our results corroborate previously measured expansion velocities (e.g., M42, M43, RCW 120, RCW 79), while in others, we provide new measurements for previously unstudied regions.

Blueshifted expansion is more common than redshifted expansion in our sample of 12 ECs. Of the 15 total expansion signatures, 11 are blueshifted and only four are redshifted. We test the statistical significance of this asymmetry using a binomial test under the null hypothesis that blueshifted and redshifted signatures occur with equal probability (50%). The resulting p-value of 0.081 indicates that this asymmetry is statistically significant at the 5% level. Selection bias may influence the observed asymmetry, though the sample size of 12 ECs is too small to determine meaningful population statistics. The FEEDBACK program primarily targets well-known regions of massive star formation, which may be more likely to be “blister” H II regions like Orion, with dense background molecular clouds. Such a geometry would make them easier to detect at optical wavelengths, but the molecular material would inhibit redshifted expansion.

Eight of the 12 ECs ($\sim\!67\%$) do not display signs of expansion along every axis sampled. Estimated expansion velocities along the 16 axes differ, which we interpret as expansion asymmetry. We show this radial expansion asymmetry in Figure 6, which defines the expansion velocities along 16 different axes in the region RCW 79 (left panel).

We choose RCW 79 as an example because it has the largest expansion asymmetry of any EC in the sample. The blueshifted expansion values of RCW 79 range from $-8.5\,\,{\rm km\,s^{-1}}$ to $-14.8\,\,{\rm km\,s^{-1}}$. There are also six axes stretching from the northeast to the southwest of the region (where the material is likely denser and therefore the intensity of [C II] emission is higher, as a higher intensity of emitting radiation from ionized gas corresponds to a larger column density of the emitting species) that do not exhibit any detectable expansion. Our observations are consistent with a non-uniform density in the surrounding medium where the expansion is fastest and easiest to detect along axes pointing towards less dense surroundings. The right panel of Figure 6 shows the same plot for M42, where the expansion is found to be the most uniform of any EC in the sample ($-13.4\,\,{\rm km\,s^{-1}}$ to $-14.00\,\,{\rm km\,s^{-1}}$ with 100% of the PV diagrams showing blueshifted expansion signatures). It is possible that M42 displays a more uniform density in its surrounding medium, leading to a more uniform expansion than that seen in regions such as RCW 79.

Of the 35 H II regions observed in this study, 23 do not exhibit any detectable signs of expansion in their [C II] PV diagrams. While this might initially suggest that these regions have reached pressure equilibrium, we tested this possibility by estimating their theoretical stagnation radii, $R_{\rm stag}$.

For each of the 23 non-expanding regions, we compiled values of the ionizing photon flux, $F^{*}$, from the literature and calculated $R_{\rm stag}$ following Equation 4. In all cases where $F^{*}$ is available, we found that the radius of the H II region is smaller than its predicted stagnation radius by at least a factor of 10. This indicates that these regions have not yet reached pressure equilibrium with the surrounding ISM and are therefore not expected to have stagnated.

The lack of detected expansion signatures could be the result of several factors. First, expansion velocities may be too low to distinguish them from the systemic velocity given our spectral resolution. Second, expansion may be occurring asymmetrically or along lines of sight not probed by our PV diagrams. Finally, the surrounding environment may be too dense or clumpy to support a clearly defined expanding shell. Despite the absence of observable expansion features, the fact that these regions have not yet reached their stagnation radii suggests that they are likely still undergoing some expansion, potentially in a slow or highly non-uniform fashion.

As can be seen in the bottom panel of Figure 4, we do not always attain a best-fit model with a velocity residual map where the residuals are all zero. Although RCW 120 shows expansion in all 16 of the axes probed with PV diagrams (see Table 2), the velocity residual map still shows some subcubes with larger residual values than others. The velocity residuals are larger in the northwest of the region, where Luisi et al. (2021) suggests that the bubble is bursting open into the surrounding ISM. This could point towards velocity components in the ionized gas along the line of sight that do not follow uniform expansion in the PDR. Such components could represent the erosion of molecular clouds as suggested by Bonne et al. (2023a). The column density of [C II] is also smaller in the center of the region than in the limb-brightened edges, which results in a smaller S/N in the center. This will result in a less robust extraction of the expansion signal, which could be influencing our residual map. For this reason, when defining our best-fit model using the sum of all residuals, we remove three of the subcubes from the observed map in Figure 4 near the center of the region, which give values highly dependent on the velocity used as a cutoff threshold to search for a residual expansion signal after rest Gaussian is subtracted.

To test the robustness of the results from our velocity residual map analysis carried out in Section 5.2 we re-run the velocity residual map analysis on simulated data with a range of signal-to-noise ratios and expansion velocities. We allow the expansion velocity to vary from 0 to 18 $\,{\rm km\,s^{-1}}$ by steps of 2 $\,{\rm km\,s^{-1}}$ while we allow the S/N to range from $-2$ to 2 in log-space by steps of 0.5 using Equation 8 from Lenz & Ayres (1992):

We run the simulation 100 times for each combination of the two free parameters. We show the average difference between the simulated and fitted expansion velocities for the explored free parameter space in Figure 7. In this figure, the left panel shows the absolute difference between the simulated and fitted expansion velocities while the right panel color bar shows the absolute difference as a percent of the simulated expansion velocity. We expect that when the S/N value of the expansion signal is higher and the simulated expansion velocity is higher (the upper right corner of both graphs in Figure 7), our results will be more reliable and the pixels will have lower values in this area. This relationship is seen in both of the plots. The vertical lines shown in both plots indicate the S/N of the observations with reliable fits from this analysis. The cyan lines represent the blueshifted signals while the red lines represent the redshifted signals. The observed S/N values are estimated as the average S/N of the redshifted/blueshifted expansion signal in the brightest 50% of subcubes within the radius defined in Section 4.

Table 3 summarizes the results of the velocity residual map fitting for each region with a successful fit and also shows the S/N for the observed expansion signals. Using the range in expansion velocities from this table along with the tabulated S/N values, we can look at Figure 7 and see that most of our observations fall in the upper right-hand corner of both plots with S/N values greater than $\textrm{log}_{10}(0.5)=3$ and a sample average expansion velocity of $\sim\!12.2$ $\,{\rm km\,s^{-1}}$.

For the 23 regions without a detected expansion signal, it is likely that these non-detections could be due to a combination of three different factors. The expansion velocity could be so low (or zero, producing no expansion signal) that its expansion signal blends with the systemic velocity signal. How much this interferes with detecting the expansion signal will be a function of the line width (which we have set as a fixed parameter in the simulations) and the expansion velocity, since a smaller line width and larger expansion velocity will have a clearer expansion signal, while a wider line width and smaller expansion velocity will result in a blended rest and expansion signal. It is also possible that the expansion signal is low enough and/or the noise could be high enough that the S/N is not large enough to produce a clear expansion signal. Finally, a region with an irregular shape and non-uniform expansion will be harder to fit with the velocity residual map analysis than a region appearing as a bubble with uniform spherical expansion. It could be a combination of all three of these factors that results in no detection of expansion signals in 23 of the 35 regions in the sample ($\sim\!66\%$).

The observed PDR radii are shown in Table 4 along with stagnation radii and theoretical radii for both stellar wind and thermal pressure driven expansion. All of the observed PDR radii fall below the theoretical range for the stagnation times. It is therefore unlikely but not impossible that we observe effects of stagnation during the ionizing source’s lifetime.

One possible explanation for the observed discrepancy between the lifetime of OB stars and the dynamical ages of the regions in our sample is a negative density gradient moving away from the central ionizing source with a non-homogeneous angular density distribution, which could drive fast expansion into the surrounding, less dense medium as the H II region evolves (Zamora-Avilés et al., 2019). In addition, the dynamical age of an H II region begins when the region “bursts” out of the dense core in which the central ionizing star was born. This could be much later than when nuclear fusion is first ignited in the proto-stellar core, resulting in dynamical lifetimes much shorter than the age of the ionizing source. The low derived dynamical ages could also highlight a selection bias towards younger, brighter regions in the sample.

In addition to these factors, cloud ablation could also contribute to the large observed expansion velocities in some of the regions. As ionizing radiation erodes dense cloud surfaces, ablated material is accelerated outward, sometimes reaching velocities comparable to or exceeding the observed expansion speeds. This effect has been documented in studies such as Bonne et al. (2023a), which observed high-velocity [C II] emission in RCW 79 due to photoablation of molecular clouds and radiation-driven flows. This could explain the low dynamical time estimates, as the high-velocity gas is continuously replenished while flowing outward. Consequently, [C II] primarily traces the gas currently being expelled from the region rather than a sustained ionization front that has been expanding since the onset of the ionizing source’s lifetime. Similar mechanisms may be at play in some of the regions in our sample, particularly those with strong local density gradients that result in asymmetric erosion and acceleration of material.

Table 4 shows an estimated theoretical range of values for $R_{s,0}$, $R_{\rm therm}$, $v_{\rm therm}$, $R_{\rm wind}$, and $v_{\rm wind}$, among other parameters. These values are estimated from Equations 1,  2,  3,  6 and  7, respectively. This table also shows the observed PDR radii ($R_{\rm PDR}$) for all ECs for which values of $R_{s,0}$ could be estimated. $R_{\rm PDR}$ roughly represents the observed center of the PDR shell while $R_{\rm therm}$ represents the theoretical inner boundary, and $R_{\rm wind}$ represents the theoretical outer boundary of the PDR shell. The values are presented as pairs separated by commas that represent the minima and maxima given a range in the input parameters ($\rho_{0}$ = $3.35\times 10^{-22}$ - $3.35\times 10^{-20}$ g cm^-3 - a number density of $100-10^{4}$ cm^-3 multiplied by the mass of an $H_{2}$ molecule, $m_{H_{2}}=3.35*10^{-24}$; $t$ = $10^{4}$ yr to 1 Myr; $L_{w}=10^{34}-10^{37}$ erg s^-1). For the ECs that we are able to estimate theoretical ranges for $v_{\rm therm}$ and/or $v_{\rm wind}$ for, we compare these ranges to the PV diagram-derived expansion velocity in the region (see Section 5.1). It is worth noting that the theoretical maxima for $v_{\rm therm}$ correspond to the earliest stages of expansion, when the swept-up shell mass is still negligible, and thus is not directly applicable to the more evolved regions considered here. In this section we will compare our observed PDR radii and derived expansion velocities to the theoretical ranges for stellar wind driven expansion and thermal pressure driven expansion in an attempt to determine the dominant mechanism of expansion in the region.

We compare observed expansion velocities and PDR radii to theoretical expectations from both thermally and wind-driven models, using the ranges shown in Table 4. In general, regions with observed expansion velocities significantly exceeding the maximum theoretical $v_{\rm therm}$ are likely dominated by stellar winds, while those with velocities within the range of $v_{\rm therm}$ are more consistent with thermal pressure-driven expansion.

In M17, M42, RCW 120, NGC 7538, and RCW 79, the observed expansion velocities exceed the maximum expected from thermal pressure but fall within the range predicted for stellar wind-driven shells. This suggests that stellar winds are the dominant driver of expansion in these regions. Notably, M17 and M42 each show PDR radii and expansion velocities consistent with wind-driven expansion, aligning with prior studies (e.g., Pabst et al., 2020). RCW 120 and RCW 79 similarly display high expansion velocities and large PDR radii indicative of stellar wind influence (Luisi et al., 2021; Bonne et al., 2023a). In NGC 7538, the high expansion velocity also points to wind-driven dynamics, although the region exhibits complex internal substructure that may complicate the interpretation (Beuther et al., 2022).

By contrast, M43 shows a low expansion velocity consistent with thermal pressure-driven expansion, and no reliable $L_{w}$ estimate is available to test for wind influence. This aligns with prior work suggesting a thermally expanding shell (Pabst et al., 2020).

RCW 36 presents an intermediate case: both the observed radius and expansion velocity fall within the ranges predicted by both models. Spatially resolved studies reveal slower expansion in dense molecular structures and faster expansion in cavity regions, suggesting that the variation in expansion velocity is tied to density variations in the surrounding ISM (Bonne et al., 2022). This is consistent with the scenario proposed by Bonne et al. (2022), in which the star formed in a thin molecular sheet. Once the expanding bubble broke out of the sheet, the expansion proceeded more rapidly in the perpendicular direction, appearing stellar wind driven in all directions despite asymmetry in the velocity field. This may help explain discrepancies between the stellar age and the dynamical age inferred from expansion modeling.

Overall, most regions with high expansion velocities are consistent with stellar wind-driven shells, while a smaller subset, such as M43, are better explained by thermal pressure. RCW 36 exemplifies a case where both mechanisms appear to play a role, varying spatially across the region. When studied at X-ray wavelengths, many of the rapidly expanding bubbles (such as M17, M42, RCW 36, RCW 79, and W40) reveal the presence of hot plasma generated by stellar wind shocks (Townsley et al., 2014, 2018, 2019).

In addition to observed and theoretical PDR radii, Table 4 also shows the estimated theoretical values for $R_{\rm stag}$ (column 4, from Equation 4), and $t_{\rm stag}$ (column 7, from Equation 5) for the ECs for which information is available on their ionizing sources. Comparing the theoretical $t_{\rm stag}$ value with the observationally derived $\langle t_{\rm dyn}\rangle$ value shows that the inferred dynamical time is, on average, 10% of the expected stagnation time. For all of these calculations, it is assumed that the ambient density is $n_{0}=100$ cm^-3 and the hot and cold medium sound speeds are taken to be $c_{i}=10.74$ $\,{\rm km\,s^{-1}}$ and $c_{o}=0.34$ $\,{\rm km\,s^{-1}}$   respectively. Every EC shows signs of expansion in at least one of their PV cuts, suggesting that stagnation has not occurred in their PDRs.

Table 5 shows our ECs for which information about the ionizing source(s) is available in the literature. The spectral type and approximate main sequence lifetime of the ionizing source are shown, along with the average observed expansion velocity and dynamical time. The last column shows that $t_{\rm dyn}$ is at least one order of magnitude smaller than the lifetime of the ionizing source for all ECs analyzed. If we instead use Equation 5 to estimate $t_{\rm dyn}$ by replacing $R_{\rm stag}$ with $R_{\rm PDR}$ (see Table 4), we find that the dynamical times decrease by about a factor of ten.

In this section we discuss the results for three of our ECs that have not yet been discussed in the literature.