Comparing the pre-SNe feedback and environmental pressures for 6000 H ii regions across 19 nearby spiral galaxies

ABSTRACT

1 INTRODUCTION

High-mass stars (>8 M_⊙) are fundamental for driving the evolution of galaxies across cosmic time, due to the large amount of energy and momentum – stellar feedback – that they inject into the interstellar medium (ISM) during their relatively short lifetimes (e.g. Krumholz et al. 2014). This is crucial, as in the absence of any stellar feedback, the ISM would rapidly cool and form stars at a high efficiency, consuming most available gas in galaxies on a short time-scale incompatible with observations (e.g. White & Rees 1978). Recent simulations (e.g. Dale, Ercolano & Bonnell 2012, 2013; Raskutti, Ostriker & Skinner 2016; Gatto et al. 2017; Rahner et al. 2017, 2019; Kim, Kim & Ostriker 2018; Kannan et al. 2020; Kim, Ostriker & Filippova 2021b; Jeffreson et al. 2021) and observational evidence (e.g. Grasha et al. 2018, 2019; Kruijssen et al. 2019; Chevance et al. 2020b,c; Kim et al. 2021a; Barrera-Ballesteros et al. 2021a,b) suggest that feedback in the early (pre-supernova) stages of high-mass stars plays a critical role in destroying molecular clouds, and hence producing the low star formation efficiencies inferred for giant molecular clouds (GMCs) in the Milky Way (e.g. Zuckerman & Evans 1974; Krumholz & Tan 2007; Evans et al. 2009; Murray 2011; Longmore et al. 2013; Evans, Heiderman & Vutisalchavakul 2014; Lee, Miville-Deschênes & Murray 2016; Barnes et al. 2017) and in many other nearby galaxies (e.g. Leroy et al. 2008, 2017; Utomo et al. 2018; Schruba, Kruijssen & Leroy 2019; Sánchez 2020; Sánchez et al. 2021).

Stellar feedback from young stars and stellar clusters is heavily associated with H ii regions. In the idealized picture, Strömgren (1939) described H ii regions as static, uniform density, spherical regions of ionized gas with a radius set by the balance of ionization and recombination rates. However, over the following decades our understanding of several physical effects has led to departures from this simple static model: the dynamical expansion of an H ii region, if the pressure in the surrounding neutral medium cannot confine its ionized gas, deviates from sphericity due to nonuniform density, injection of energy and momentum by a stellar wind, absorption of hydrogen-ionizing photons by dust grains and radiation pressure acting on gas and dust (see e.g. Kahn 1954; Savedoff & Greene 1955; Mathews 1967, 1969; Gail & Sedlmayr 1979). More recently, many works have focused on observationally quantifying the impact of the various feedback mechanisms on driving the expansion of feedback-driven bubbles by detailed studies of their feedback mechanisms, ionization structures, morphologies, dynamics, and the stellar content across the Milky Way (e.g. Rugel et al. 2019; Watkins et al. 2019; Barnes et al. 2020; Olivier et al. 2021), the Small and Large Magellanic Clouds (Oey 1996a,b; Pellegrini, Baldwin & Ferland 2010; Lopez et al. 2011, 2014; Pellegrini et al. 2012; Chevance et al. 2016; McLeod et al. 2019), and in nearby galaxies (e.g. McLeod et al. 2020, 2021).

The dynamics and expansion of H ii regions may be driven by several possible sources of internal energy and momentum injection. By definition, H ii regions are filled with warm (∼10⁴ K) ionized hydrogen, which imparts an outward gas pressure (e.g. Spitzer 1978). Yet, in addition, several other forms of stellar feedback can drive the dynamics of H ii regions and deposit energy and momentum in the surrounding ISM: the direct radiation of stars (e.g. Dopita et al. 2005, 2006; Krumholz & Matzner 2009; Fall, Krumholz & Matzner 2010; Murray, Quataert & Thompson 2010; Peters et al. 2010; Commerçon, Hennebelle & Henning 2011; Hopkins, Quataert & Murray 2011; Rathjen et al. 2021), the dust-processed infrared radiation (e.g. Thompson, Quataert & Murray 2005; Murray et al. 2010; Andrews & Thompson 2011; Skinner & Ostriker 2015; Reissl et al. 2018; Tsang & Milosavljević 2018), stellar winds and supernovae (SNe; e.g. Yorke et al. 1989; Harper-Clark & Murray 2009; Rogers & Pittard 2013), and protostellar outflows/jets (e.g. Quillen et al. 2005; Cunningham et al. 2006; Li & Nakamura 2006; Nakamura & Li 2008; Wang et al. 2010; Rosen et al. 2020). While we have a good understanding of how individual stars or massive stellar populations produce each of these effects, the field still lacks a substantial number of quantitative observations for a diverse sample of H ii regions and their environments.

In this work, we investigate the role of early stellar feedback within a sample of H ii regions identified across the discs of 19 nearby spiral galaxies (see Table 1). To do so, this work exploits optical integral field unit spectroscopy (IFU; see Sánchez 2020 for a recent review of nearby galaxy IFU studies) from the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) instrument mounted on the Very Large Telescope (VLT) obtained as part of a VLT large programme (PI: Schinnerer). The spatial sampling and large field of view of MUSE allow us to analyse the properties of the ionized gas at high resolution (∼50–100 pc) in systems as far as 19 Mpc (see Fig. 1; Emsellem et al. 2021). We use these observations to place limits on the sizes, luminosities, and ultimately the feedback-related pressure terms (i.e. the direct radiation pressure, the pressure from stellar winds, and the ionized gas pressure) for each of the H ii regions. Contrasting these with the local environmental pressure, we can capture a snapshot of the physical and dynamical state of the ionized gas at the later evolutionary stages of H ii regions (around a few Myr).

Open in new tabDownload slide

Open in new tabDownload slide

Figure 1.

Overview of the H ii region samples for two of galaxies studied in this work (NGC 1672 and NGC 1300). Upper left: Three colour image composed of 814 nm (red), 555 nm (green) and 435 nm (blue) wideband emission from the PHANGS-HST survey (Lee et al. 2021), and the HST continuum-subtracted 658 nm or H |$\alpha$| (pink) narrowband emission. Upper right: MUSE H |$\alpha$| emission map obtained as part of the PHANGS-MUSE survey (Emsellem et al. 2021). Also shown are the beam sizes of the MUSE H |$\alpha$| observations in the lower left corner, and a scale bar of 5 kpc in the lower right corner. Lower left: Full H ii region sample identified within each galaxy (Santoro et al. 2021; see also Section 2.3). Circle sizes represent the measured sizes of the H ii regions (r_eff; see Section 3). Lower right: Sample of H ii regions that have measurements of both the electron density (n_e) and effective radius (r_eff) shown by blue circles. Also shown are the samples that are resolved and below the low-density limit as orange circles, and that are not resolved and above the low-density limit as green circles (see Section 3.4).

This paper is organized as follows. In Section 2, we introduce the sample of 19 nearby galaxies and the MUSE observations that are used to identify and study their H ii region populations. In Section 3, we outline assumptions for the unresolved density distribution within each of the identified H ii regions and use these to estimate their physical properties. In Section 4, we place limits on the internal pressures within each H ii region. We compare how the pressure components vary across the galaxies, how the total internal pressure compares to the external pressure, and how these results compare to samples within the literature in Section 5. Finally, the main results of this paper are summarized in Section 6.

2 SAMPLE SELECTION

To study the young stellar feedback mechanisms, we require a large sample of H ii regions that have accurate measurements of ionizing photon flux, electron density, and/or their size (see Section 4). In this section, we outline the sample of galaxies studied in this work, introduce the MUSE/VLT observations taken as part of the PHANGS (Physics at High Angular Resolution in Nearby GalaxieS) survey (see Emsellem et al. 2021; Lee et al. 2021; Leroy et al. 2021a), and outline how these observations are used to identify a catalogue of around 23 699 H ii regions (see Fig. 1).

2.1 Galaxy sample of PHANGS-MUSE

The parent galaxy sample of the overall PHANGS programme was constructed according to the criteria outlined in Leroy et al. (2021a) and Emsellem et al. (2021). Briefly, the PHANGS galaxies were selected to be observable by both ALMA (Leroy et al. 2021a,b) and MUSE (Emsellem et al. 2021; −75^○ ≤ δ ≤ +25^○), nearby (5 Mpc ≤ D ≤ 17 Mpc),¹ to allow star-forming regions and molecular clouds to be resolved at high spatial resolution (∼100 pc), at low to moderate inclination to limit the effects of extinction and line-of-sight confusion (i < 65^○), and to be massive star-forming galaxies with log (M_⋆/M_⊙) ≳ 9.75, and log (sSFR/yr⁻¹) ≳ −11. In this work, we use a subset of 19 PHANGS galaxies that have been observed with the MUSE spectrograph on the VLT (see Emsellem et al. 2021, for full details). The sample of PHANGS-MUSE galaxies is given in Table 1 along with their key properties. For each galaxy, we tabulate the central right ascension (RA) and declination (Dec.) from Salo et al. (2015), and the inclination (i) and position angle (PA) based on PHANGS CO(2–1) kinematics from Lang et al. (2020). Insufficient CO emission was detected in IC 5332 to allow the kinematics to be constrained, and so for this galaxy we use values from i and PA from Querejeta et al. (2015). We also show the source distances that are taken from the compilation of Anand et al. (2021), which along with their Tip of the Red Giant Branch method (TRGB) estimates also include distances taken from Freedman et al. (2001), Nugent et al. (2006), Jacobs et al. (2009), Kourkchi & Tully (2017), Shaya et al. (2017), and Kourkchi et al. (2020). The deprojected galactocentric radii (in parsec) quoted in this work use these central positions, orientations, and distance estimates. In Table 1, we list the average metallicity within the region of each galaxy mapped by MUSE, computed using the Scal method of Pilyugin & Grebel (2016), as discussed in more detail in Kreckel et al. (2019). We also show (in logarithmic units) mass estimates of the atomic gas (⁠|$M_{\text{H}\, \small {I}}$|⁠; Makarov et al. 2014), molecular gas (|$M_{\mathrm{ H}_2}$|; Leroy et al. 2021a), and stars (M_⋆; Leroy et al. 2021a), and the average star formation rate (SFR; Leroy et al. 2021a).

2.2 MUSE observations

The MUSE Integral Field Unit (IFU) provides a 1 arcmin × 1 arcmin field of view, 0.2 arcsec pixels, and a typical spectral resolution (FWHM) of ∼2.5 Å (or ∼100 km s⁻¹) covering the spectral range 4800–9300 Å. Observations of the 19 galaxies are reduced using the pymusepipe package. pymusepipe was developed specifically for these observations by the PHANGS team,² and is a python wrapper around the main processing steps of the data reduction conducted by the MUSE pipeline (MUSE DRS; Weilbacher et al. 2020) accessed via EsoRex command-line recipes. A complete discussion of the processing and reduction of the MUSE observations is presented in Emsellem et al. (2021). The final reduced cubes have an angular resolution ranging between ∼0.6 arcsec (for the subset observed using ground-layer correction adaptive optics) and ∼1.2 arcsec (see Table 2), which at the distances of our sample corresponds to physical scales of 25–70 pc at the corresponding galaxy distance (see Table 1).

To identify and determine the properties of the H ii regions across each galaxy using their optical spectroscopic features, Emsellem et al. (2021) first produce emission line maps covering the full field of view. To do so, all data cubes are processed using the penalized pixel fitting python package (ppxf; Cappellari & Emsellem 2004; Cappellari 2017). These fits include E-MILES simple stellar population models (Vazdekis et al. 2016) and a set of emission lines that are treated as additional Gaussian templates. A detailed description of the spectral fitting process is presented in Emsellem et al. (2021).

2.3 Sample of H ii regions

We make use of an ionized nebulae catalogue derived from point spread function (PSF)-homogenized H |$\alpha$| line maps (the copt data products described in Emsellem et al. 2021), as described in detail in Santoro et al. (2021). Briefly, Santoro et al. (2021) first run an implementation of hiiphot (Thilker, Braun & Walterbos 2000), which has been adapted for use with integral field data. The final catalogue contains not only H ii regions, but also planetary nebulae and supernova remnants. The ionized nebulae spatial masks are then applied to the original data cube to extract integrated spectra for each object. The emission lines within these spectra are subsequently fitted using the same procedure as the one described in Section 2.2. This fitting procedure includes both strong lines (e.g. H |$\beta$|⁠, [O iii] |$\beta 5007$|⁠, H |$\alpha$|⁠, [N ii] |$\beta 6583$|⁠, [S ii] |$\beta 6716$|⁠, [S ii] |$\beta 6731$|⁠) as well as fainter, temperature-sensitive auroral lines (e.g. [N ii]λ5755, [S iii]λ6312). We then employ the following selection criteria to identify H ii regions within the nebular catalogue:

• We remove all nebulae that are flagged as being smaller than a single PSF (see Table 2) separated from the edge of the field of view in order to ensure that we do not include regions that have artificially small sizes because they lie partially outside of the field of view.

• We require all of the strong lines used to compute line ratio diagnostics (see below) to be detected with signal-to-noise greater than three.³

• We apply three emission line ratio diagnostic diagrams (BPT diagrams; Baldwin, Phillips & Terlevich 1981) to separate the nebulae photoionized by high-mass stars from those ionized by other sources (e.g. Active Galactic Nuclei). Nebulae are classified as H ii regions if they fall below the Kauffmann et al. (2003) line in the |$[\rm{O}\, \small{III}]\,\beta 5007/\rm{H\,\beta }$| versus |$[\rm{N}\, \small {II}]\beta 6583/\rm{H\,\alpha }$| diagram and below the Kewley et al. (2001) line in the |$[\rm{O}\, \small{III}]\beta 5007/\rm{H\,\beta }$| versus |$[\rm{S}\, \small {II}]\beta 6716{+}\beta 6731/\rm{H\,\alpha }$| and |$[\rm{O}\, \small{III}]\beta 5007/\rm{H\,\beta }$| versus |$[\rm{O}\, \small {I}]\beta 6300/\rm{H\,\alpha }$| diagrams.⁴

• We discard regions with H |$\alpha$| velocity dispersions that exceed 100 km s⁻¹, which are likely compact supernova remnants.⁵

We note that Santoro et al. (2021) conducted the source identification on the H |$\alpha$| emission maps from MUSE that have not been extinction corrected. The extinction correction is then determined for the H |$\alpha$| flux within the source masks using the Balmer decrement (see Section 3.3). Hence, there is no correction for the surface brightness dimming within extincted regions during the source identification stage of our analysis. Due to the high sensitivity of our MUSE observations, however, we can typically recover the H |$\alpha$| emission out to the point where the H ii regions merge with the diffuse ionized gas (DIG). In addition, we do not directly account for the contribution of the DIG in our determination of the H ii region line fluxes, as this was minimized by the selection of the source identification parameters in the HIIphot package. That said, the contribution of the DIG to the H ii regions studied in this work is expected to be low, given that here we have explicitly chosen to study the brightest regions due to our high flux selection thresholds (signal-to-noise greater than three). Doing so is also not trivial, and can introduce large uncertainties into the remaining fluxes. Lastly, the physical interpretation of removing the DIG in our pressures analysis is not clear (see Section 4). The leakage of ionized gas from the H ii regions is the main contributor to the DIG (Belfiore et al. 2021), yet this gas could still provide a contribution to the pressure terms we measure.

The above selection criteria remove 7798 sources from the initial catalogue of 31 ×  497 objects, hence leaving ∼75 per cent of the sources (23 699) as H ii regions. The number of H ii regions within each galaxy is summarized in Table 2. To achieve a final sample of ∼6000 H ii regions, which is used throughout this work, we imposed further selection criteria on size and density measurements (outlined in Section 3). The total number of identified H ii regions within each galaxy is presented in Table 1, and ranges from a few hundred to a few thousand. In Fig. 1, we show the distribution of the H|$\alpha$| emission line and H ii region catalogue compared to an optical HST three colour composite (Lee et al. 2021) for two galaxies in our sample (NGC 1300 and NGC 1672).

3 PHYSICAL PROPERTIES

We use the wealth of information provided by the MUSE observations to estimate several fundamental physical properties for the H ii regions in our catalogue: the H ii region sizes (r_eff), electron densities (n_e), and ionization rates (Q), as well as the mass (M_cl), bolometric luminosity (L_bol), total mass loss rate (⁠|$\dot{M}$|⁠), and mechanical luminosity (L_mech) of the cluster or association powering each region.

With regard to the H ii region sizes, we face the complication that the size that we ideally want to measure is the characteristic radius at which the majority of the mass of the ionized gas is located, since this is more relevant for understanding the dynamics of the H ii region and the interplay between the different pressure terms than the maximum physical extent of the ionized region. For an H ii region that is well described by the classical Strömgren sphere solution (Strömgren 1939), or one with a shell-like morphology, this is comparable to the extent of the H ii region, but for a partially embedded or blister-type H ii region, this is not necessarily the case.

At the resolution of our MUSE observations, we cannot easily distinguish these different H ii region morphologies, and so we instead consider two limiting cases for the distribution of ionized gas within the observed H ii regions. In one limit, we assume that the ionized gas is smoothly distributed throughout the measured volume of the H ii region, and the properties of the H ii regions are hence determined for the maximum radius (i.e. the measured r_eff; Section 3.1). In the other limit, we assume that the ionized gas is clumpy, with most of the mass located in dense clumps that lie close to the centre of the H ii region. In this case, the properties are determined using the minimum volume within which these clumps can be accommodated while remaining consistent with the measured electron density (Section 3.2) and H |$\alpha$| flux. The radial length-scale associated with this minimum volume is hereafter denoted as r_min. These are, of course, not the only possibilities – for instance, a shell-like H ii region may have most of its gas in dense clumps that are located far from the ionizing source at r ∼ r_eff – but for our purposes we restrict our attention to these two limiting cases as they will later allow us to put upper and lower limits on the various pressure terms. These assumptions are illustrated in Fig. 2 and are discussed further in Section 3.4.

Figure 2.

Schematic diagram representative of an H ii region in our sample. As the linear resolution of our observations is of the order a few tens of parsec at the distance of the galaxies in our sample, in this work we consider two limiting cases for the unresolved density distribution within an H ii region (see Section 3). Left-hand panel: On the one extreme, we posit that the density distribution is smooth, and the properties of the H ii regions are hence determined for the effective spherical radius determined from the resolved H ii regions size (i.e. the measured r_eff). Right-hand panel: On the other extreme, we posit that the H ii regions have a clumpy density structure, with the clumps located close to the centre of the region, as might be expected, e.g. in an H ii region in the process of breaking out from a molecular cloud. In this case, the properties are determined using the minimum volume within which these dense clumps can be distributed while remaining consistent with the measured H |$\alpha$| flux and the density inferred from the [S ii] doublet. The corresponding size scale in this case is r_min. Our calculated pressure terms therefore each have two limits: a maximum (P_max) at r_min and a minimum (P_min) at r_eff (Section 4).

Open in new tabDownload slide

3.1 Measured effective radii – r_eff

We estimate the effective angular radii for the H ii regions in the catalogue by circularizing the area contained within each H ii region; |$r_\mathrm{eff}= \sqrt{a/\pi }$|⁠,⁶ where a is the area enclosed by the boundary identified using the hiiphot routine (i.e. above some intensity threshold, not a fitted ellipse). An inherent problem with many automated source identification algorithms, such as with HIIphot, is their tendency to separate compact emission into distinct sources that have sizes comparable to the point spread function (or resolution) of the input observations. The result of this is the identification of a large sample of unresolved or only marginally resolved (point) sources, the sizes for which are either unconstrained or highly uncertain. As we require an accurate measure of H ii region sizes for our pressure analysis (see Section 4.1.1), we consider regions to be resolved only if their effective radius satisfies the resolution criterion r_eff ≥ FWHM_PSF, i.e. the effective diameter must be at least two times the resolution limit of the observations for each galaxy (see Table 2). This threshold was chosen to include the H ii regions that are significantly more extended than the observational limits, yet without substantially limiting our sample for the most distant galaxies. We consider all regions smaller than the r_eff < FWHM_PSF limit to be unresolved, and do not make use of the values of r_eff derived for these regions in our later analysis. For the subset of these unresolved regions for which we can determine the electron density (see Section 3.2), we can place a lower limit on their sizes, as discussed later in Section 3.4. Unresolved regions without a well-determined electron density cannot be assigned meaningful values of either r_eff or r_min and are not considered further in our analysis.

The physical effective radius of each H ii region in units of parsec is determined using the source distance given in Table 1. We find that the r_eff size range across the whole sample of H ii regions (including both resolved and unresolved sources) is 16.6–388.4 pc (median: 54.6 pc), while for the resolved sub-sample of H ii regions it is 26.2–388.4 pc (median: 80.64 pc). In Fig. 3, we show two distributions of r_eff: one for regions that are resolved and that have a measured electron density (blue histogram with solid outline) and one for regions that are resolved but that do not have a measured electron density (orange histogram). In the figure, we also show the distribution of r_min for those regions in which it can be calculated (see Section 3.4 below).

Figure 3.

Histograms of the size distribution for the H ii region sample for assumptions of a (a) smooth and (b) clumpy unresolved density distribution (see Fig. 2). Upper panel: We show r_eff for the resolved (r_eff ≥ FWHM_PSF; see Section 3.1) sources as both blue (solid outline) and orange histograms. The blue and orange histograms differentiate the samples with and without n_e,max electron density measurements, respectively (see Section 3.2). Lower panel: Where possible, we also determine a size (r_min) using the n_e,max determined from |$R_{[\text{S}\, \small {II}]}$| and the H |$\alpha$| flux (see Section 3.4). These r_min distributions are shown in green and blue (dashed line). Note that the H ii regions that are both resolved and are below the low-density limit have both r_min and r_eff estimates, and, hence, the blue histograms have the same sample size.

Open in new tabDownload slide

3.2 Measured electron densities – n_e,max

To calculate the electron density of the H ii regions in our sample (n_e,max), we use the pyneb package (Luridiana, Morisset & Shaw 2015). pyneb is a python module for the analysis of emission lines. It solves the equilibrium equations and determines level populations for one or several user-selected model atoms and ions. We use pyneb to solve for the electron density within each H ii region given the flux ratio |$R_{[\rm{S}\, \small {II}]}= F_{[\rm{S}\, \small {II}]\beta 6716} / F_{[\rm{S}\, \small {II}]\beta 6731}$|⁠, and a value for the electron temperature (T_e; Belfiore et al. in preparation).⁷In the lower panel of Fig. 4, we show the PyNeb solutions for n_e as a function of |$R_{[\rm{S}\, \small {II}]}$| for two values of the electron temperature that are representative of the extremes of the temperature distribution for the H ii region sample (6000–16 000 K; Belfiore et al. in prep). |$R_{[\rm{S}\, \small {II}]}$| has two limiting values: a high-density limit at |$R_{[\rm{S}\, \small {II}]}\sim 0.4$| that is reached at number densities above a few 1000 cm⁻³ and a low-density limit at |$R_{[\rm{S}\, \small {II}]}\sim 1.45$| that is reached at number densities below a few 100 to 10s cm⁻³. For H ii regions with densities between these two limits, measuring the value of |$R_{[\rm{S}\, \small {II}]}$| allows us to infer the [S ii] emission-weighted mean density of the ionized gas. In the case of a clumpy H ii region, this density will primarily reflect that of the gas in the clumps and may significantly exceed the mean density of the ionized gas in the H ii region as a whole. For that reason, we refer to this estimate loosely as the ‘maximum’ electron density for the H ii region, which we will later compare with a ‘minimum’ electron density estimate derived using a different technique. Note also that the inferred density depends on the electron temperature, but as Fig. 4 shows, this dependence is relatively weak.

Figure 4.

Upper panel: Histogram of the |$R_{[\rm{S}\, \small {II}]}= F_{[\rm{S}\, \small {II}]\beta 6716} / F_{[\rm{S}\, \small {II}]\beta 6731}$| line ratio across the H ii region sample. We show the distribution of |$R_{[\rm{S}\, \small {II}]}$| for the full sample of H ii regions as a grey histogram. The vertical dashed red line and shaded region shows the low-density limit (also shown in the lower panel). The blue and green histograms show the distributions of H ii regions that are statistically distinguishable from the low-density limit (see Fig. 6), and have resolved and unresolved sizes (see Section 3.1). The orange histogram shows the H ii regions that are indistinguishable from the low-density limit, yet are resolved (see Fig. 6). Lower panel: The dashed curve shows the conversion between |$R_{[\rm{S}\, \small {II}]}$| and electron density (n_e, see right y-axis) as determined from pyneb assuming two electron temperatures that span the observed range of approximately 6000–16 000 K.

Open in new tabDownload slide

In the upper panel of Fig. 4, we show the distribution of values for |$R_{[\rm{S}\, \small {II}]}$| that we measure for the whole sample of H ii regions (grey). We see that the distribution peaks just above |$R_{[\rm{S}\, \small {II}]}\gt 1.4$|⁠, i.e. at a value comparable to the one we expect to recover on the low-density limit. However, we also see that many of the values we measure for |$R_{[\rm{S}\, \small {II}]}$| lie above this limiting value (indicated by the red shaded region in the figure). These values are unphysical and so we assume that they are due to the statistical errors in our measurements of the fluxes of the [S ii] lines, which introduce an error into the calculated line ratio. To check this, we calculate for each region the absolute difference of |$R_{[\rm{S}\, \small {II}]}$| from the low-density limiting value, normalized by the uncertainty in the value of |$R_{[\rm{S}\, \small {II}]}$| for that region (⁠|$\sigma _{R_{[\rm{S}\, \small {II}]}}$|⁠). This uncertainty is calculated using the formal errors in the fluxes of the two [S ii] lines, adjusted upwards by a factor of 1.43 to account for the fact that these formal errors are still somewhat underestimated in the latest version of the MUSE data reduction, likely due to imperfect sky substraction.⁸ We use a temperature-dependent low-density limit of |$R_{[\rm{S}\, \small {II}]}= 1.49\!-\!3.94\times 10^{-6}\, T_\mathrm{e}$|⁠. For regions that do not have reliable estimates of the electron temperature, as can happen if the |$[\rm{N}\, \small {II}]\beta 5755$| auroral line is not detected, we adopt a representative electron temperature of 8000 K, which yields |$R_{[\rm{S}\, \small {II}]}= 1.46$|⁠. Finally, we account for the uncertainty in |$R_{[\rm{S}\, \small {II}]}$| that arises due to the statistical error in the T_e measurement by combining this in quadrature with the line flux uncertainties when computing |$\sigma _{R_{[\rm{S}\, \small {II}]}}$|⁠.

We show the distribution of the normalized absolute differences in Fig. 5. The blue histogram corresponds to H ii regions with values of |$R_{[\rm{S}\, \small {II}]}$| above the low-density limiting value, while the orange histogram shows the H ii regions with values of |$R_{[\rm{S}\, \small {II}]}$| below this limit. We see that the distribution of |$R_{[\rm{S}\, \small {II}]}$| values in the unphysical region above the low-density limit is Gaussian, with a standard deviation of 1, consistent with what we would expect if all of these regions have a true value of |$R_{[\rm{S}\, \small {II}]}$| at or very close to the low-density limiting value. For H ii regions with measured |$R_{[\rm{S}\, \small {II}]}$| below this limit, we recover a Gaussian distribution for low values of the normalized deviation and a clear non-Gaussian tail for higher deviations. In order to exclude regions which are consistent with Gaussian noise around the low-density limit, we select H ii regions that are at least |$3\sigma _{R_{[\rm{S}\, \small {II}]}}$| away from the low-density limit, where the Gaussian distribution becomes sub-dominant (see dashed vertical grey line in Fig. 5). The |$R_{[\rm{S}\, \small {II}]}$| distributions for the samples of H ii regions, that are significantly below the low-density limit, are shown as blue and green histograms in Fig. 4. Note that these fall below the histogram for the whole sample (shown in grey) as, even below the low-density limit, the uncertainties of |$R_{[\rm{S}\, \small {II}]}$| can be large, causing some values to be indistinguishable from the low-density limit. We also show the distribution of H ii regions that are indistinguishable from the low-density limit and have resolved sizes, for which we will calculate lower limits for the electron density (see Section 3.4).

Figure 5.

Histogram of the absolute difference of the |$R_{[\rm{S}\, \small {II}]}= F_{[\rm{S}\, \small {II}]\beta 6716} / F_{[\rm{S}\, \small {II}]\beta 6731}$| line ratio from the low-density limit, normalized by the error in the line ratio for each H ii region. The orange and blue histograms show the distributions of |$R_{[\rm{S}\, \small {II}]}$| above and below an electron temperature-dependent low-density limit (see Section 3.2). The vertical dotted lines show a difference of 1 and 3 |$\sigma _{R_{[\rm{S}\, \small {II}]}}$| from the low-density limit. The dashed curve shows the Gaussian function fit to the histogram distribution of H ii regions above the |$R_{[\rm{S}\, \small {II}]}$| limit, which has a standard deviation of ∼1 |$\sigma _{R_{[\rm{S}\, \small {II}]}}$|⁠. We highlight the non-Gaussian tail for those sources below the low-density limit.

Open in new tabDownload slide

The histogram distribution of allowed n_e,max measurements below the low-density limit are shown in Fig. 6. Here, we differentiate those that are resolved and have both n_e,max and r_eff measurements (in blue), and those that are unresolved and have only n_e,max measurements (see Section 3.1). We find that the n_e range across these samples of H ii regions is 13–537 cm⁻³ (median: 40.5 cm⁻³), respectively. In Fig. 6, we also show the distribution of the (minimum) electron densities (n_e,min), which are determined using a different method for the sample of resolved sources (see Section 3.4). These values of the electron density are similar to those determined from other IFU studies of resolve H ii regions within nearby galaxies (e.g. NGC 628: Rousseau-Nepton et al. 2018; NGC 300: McLeod et al. 2020). We note, however, that they sit at the lower end of the values estimated from lower spectral and angular resolution studies (e.g. Sánchez et al. 2012, 2015; Espinosa-Ponce et al. 2020). The comparison between resolved and unresolved studies presents an interesting avenue for future investigations.

Figure 6.

Histograms of the electron density (n_e) distribution for the H ii region sample for assumptions of a (a) smooth and (b) clumpy unresolved density distribution (see Fig. 2). Upper panel: The blue (solid line) and green histograms show the n_e,max distributions for H ii regions that are statistically distinguishable from the low-density limit (see Fig. 6). The blue and green histograms differentiate the samples with resolved and unresolved sizes (see Section 3.1). Lower panel: We derive an n_e,min for the resolved H ii regions using the measured effective radius (r_eff) and extinction-corrected H |$\alpha$| flux (⁠|$F_\mathrm{H\,\alpha }$|⁠; see Section 3.4). We show the distribution of n_e,min as orange and blue (dashed line) histograms. Note that the H ii regions that are both resolved and are below the low-density limit have both n_e,max and n_e,min estimates, and, hence, the blue histograms have the same sample size.

Open in new tabDownload slide

Overall, out of our initial sample of 23 699 H ii regions, a total of 5810 (∼25 per cent) have measurements of either r_eff or n_e,max. 5669 (98 per cent) of these regions have sizes above our resolution limit and hence have valid r_eff measurements, whilst 2379 (41 per cent) have densities large enough to allow us to distinguish |$R_{[\rm{S}\, \small {II}]}$| from its low-density limiting value, thereby allowing us to determine n_e,max for these regions. Finally, we have both measurements for a total of 2238 H ii regions (see Table 2 for a summary). We remind the reader that these measurements correspond to our two limiting cases (see Fig. 2): the ionized gas is smoothly distributed throughout the measured volume of the H ii region (i.e. r_eff), or the ionized gas is clumpy and fills only some fraction of the H ii region close to the source (i.e. n_e,max). Note that the requirement of a resolved size or accurate electron density measurement, biases these samples to the brightest and largest H ii regions within each galaxy. In Section 3.4, we use the estimates of maximum sizes, r_eff (or electron densities, n_e,max), of the H ii regions to place lower limits on the electron densities, n_e,min (or sizes, r_min); or in other words, we also determine the electron density and sizes for both our assumptions of a smooth and clumpy ionized gas density distribution. To do so, however, we must first outline how the ionising photon rate within each H ii region is calculated.

3.3 Ionization rate – Q

To calculate the ionization rate (Q), we first calculate the H |$\alpha$| luminosity from |$L_{\rm{H\,\alpha }} = 4 \pi D^2 F_{\rm{H\,\alpha }}$|⁠, where |$F_{\rm{H\,\alpha }}$| is the extinction-corrected H |$\alpha$| flux computed by Santoro et al. (2021) and D is the distance to each galaxy given in Table 1. The extinction correction is computed by measuring the reddening using the H α/H β ratio measured for each H ii region and applying a correction assuming the O’Donnell (1994) reddening law with R_V = 3.1 and a theoretical H α/H β = 2.86. For optically thick nebulae (case B recombination; Osterbrock & Ferland 2006) at T_e = 10 000 K, the ionization rate is given as |$Q \approx L_{\rm{H\,\alpha }} (\alpha _\mathrm{B}/\alpha ^\mathrm{eff}_{\rm{H\,\alpha }}) / (h\nu _{\rm{H\,\alpha }}) = L_{\rm{H\,\alpha }} / (0.45 h\nu _{\rm{H\,\alpha }})$|⁠, where the total recombination coefficient of hydrogen is |$\alpha _\mathrm{B} = 2.59\times 10^{-13}$| cm³ s⁻¹, the effective recombination coefficient (i.e. the rate coefficient for recombinations resulting in the emission of an H |$\alpha$| photon) is |$\alpha ^\mathrm{eff}_{\rm{H\,\alpha }} \approx 1.17\times 10^{-13}$| cm³ s⁻¹ (Osterbrock & Ferland 2006), |$\nu _{\rm{H\alpha }}$| is the frequency of the H |$\alpha$| emission line and h is the Planck constant.

We find that Q ranges from 10^49.0 to 10^52.8 s⁻¹ across the sample of H ii regions with both n_e,max and r_eff measurements. We find that Q ranges from 10^49.1 to 10^51.6 s⁻¹ across the n_e,max only sample, and 10^47.6 to 10^51.0 s⁻¹ across the r_eff only sample. The fact that we recover systematically lower Q values for the r_eff only sample is easily understood: the H ii regions that are weaker in H |$\alpha$| emission (⁠|$Q \propto F_{\rm{H\alpha }}$|⁠) are also weaker in the [S ii]|$\beta 6716$| and [S ii]|$\beta 6731$| emission lines, causing larger errors on the |$R_{[\rm{S}\, \small {II}]}$| ratio, making it harder to distinguish |$R_{[\rm{S}\, \small {II}]}$| in these regions from the low-density limit. Finally, note that the value of Q we derive for each H ii region does not depend on the escape fraction of ionising photons from that H ii region (f_esc), since Q here refers only to the ionization rate of gas within the H ii region.

3.4 Minimum radii and electron densities – r_min and n_e,min

For the H ii regions with sizes greater than the resolution limit, the measured r_eff represents an estimate of their maximum extent. However, as mentioned previously, in cases where the H ii region is clumpy and the clumps are close to the ionising source, the average distance of the clumps from the source is a more appropriate measure of the H ii region size from the point of view of understanding its dynamics. Our observations do not have sufficient resolution to allow us to measure this distance directly. However, for regions where we have a measure of the electron density from the [S ii]doublet, we can put a lower limit on this size, which we hereafter denote as r_min. We can do this because the extinction-corrected H|$\alpha$| luminosity of an H ii region is determined by three quantities: the electron temperature (measured as explained in the previous section), the root-mean-squared density of the gas and the volume of the H ii region. If we assume that the root-mean-squared density is the same as the density we measure from [S ii], then we can straightforwardly solve for the volume, and hence the size of the region if we approximate it as a sphere. The rms density could of course be lower than our [S ii]-derived density if the [S ii]-bright clumps fill only a small fraction of the volume, but is unlikely to be larger than this value. The estimate of the H ii region size that we get from this argument is therefore a lower limit on the true size, complementing the upper limit we get from r_eff. Note also that we can also derive a value of r_min for unresolved H ii regions for which we cannot measure an accurate r_eff, so long as we have a measure of n_e,max for these regions. Finally, we can also apply the same logic to derive an estimate of the minimum density of our resolved H ii regions, n_e,min, by fixing the volume and solving for the density.

The values of r_min we derive from this approach range from a few parsecs to a few tens of parsecs, as illustrated in Fig. 3. For the resolved regions, they are typically around a factor of ten smaller than r_eff. We also see that the values of r_min that we derive for resolved H ii regions are generally larger than those we derive for unresolved regions. This is a consequence of the H ii region size–luminosity relationship: larger H ii regions tend to also be brighter, and hence the minimum volume of dense gas required to produce their observed H|$\alpha$| luminosities is larger.

In Table 2, we list the number of H ii regions in each galaxy for which we can derive both r_min and n_e,min (2238 regions in total), only n_e,min (3431 regions) or only r_min (141 regions).

3.5 Stellar population models

Lastly, we require a final set of properties for our H ii region sample before determining their internal pressure terms, which we obtain from synthetic stellar population modelling. Namely, in this section we estimate their bolometric luminosity (L_bol), cluster mass (M_cl), mass-loss rate (⁠|$\dot{M}$|⁠), and mechanical luminosity (L_mech). We employ the starburst99 model (Leitherer et al. 1999),⁹ adopting the default parameter set and varying the cluster mass between 10⁴ and 10⁶ M_⊙. Of note within the default parameter set, we use the Evolution wind model (Leitherer, Robert & Drissen 1992) for the calculation of the wind power, the mass-loss rates from the Geneva models with no rotation, an instantaneous star formation burst populating a with Kroupa initial mass function (IMF; Kroupa 2001). We assess the evolution of the ionization, (wind) feedback power and mass-loss across a time range of ∼0–10 Myr (i.e. requiring the quanta, snr, power, yield, spectrum, ewidth outputs from starburst99).

3.5.1 Bolometric luminosity – L_bol

We wish to estimate the bolometric luminosity of the cluster(s) responsible for the H ii regions in order to determine the radiation pressure (see Section 4.1.1). In Fig. 7, we show the bolometric luminosity and H |$\alpha$| luminosity as a function of time for a range of cluster masses (also see Agertz et al. 2013), computed with the assumption that no ionising photons escape from the H ii region (i.e. that f_esc = 0). We see that the bolometric luminosity remains relatively constant across the first 10 Myr (varying by only ∼1 dex), whereas the H|$\alpha$| luminosity drops significantly after ∼2 Myr (∼3 dex). Since our H ii region sample is constructed from observed H |$\alpha$| emission, by definition they must have lifetimes less than the lifetime of ionising radiation. We then consider the maximum age (t_max) as the time when |$L_{\rm{H\,\alpha }}$| has dropped significantly from the zero-age value (⁠|$L_{\rm{H\,\alpha }}(t_0)$|⁠). We set t_max at 4 Myr, where |$L_{\rm{H\,\alpha }}$| has decreased by half an order of magnitude (a factor of around 3). This t_max includes the zero-age main sequence (ZAMS) and Wolf–Rayet phases of the high-mass stars (as labelled in Fig. 7). We note that the choice of t_max is somewhat arbitrary; however, we do not see any appreciative change in our results by increasing it to, e.g. 5 Myr, where |$L_{\rm{H\,\alpha }}$| has dropped by an order of magnitude, or decreasing it to, e.g. 0 Myr, to only include the ZAMS.

Figure 7.

Early time evolution of the luminosity from the starburst99 model (Leitherer et al. 1999). Left-hand and centre panels: We show the H |$\alpha$| luminosity (⁠|$L_{\rm{H\alpha }}$|⁠) and the bolometric luminosity (L_bol) for cluster with masses of M_⋆ = 10⁴–10⁶ M_⊙, in steps of 0.5 dex. Right-hand panel: We show the ratio of |$L_\mathrm{bol}/L_{\rm{H\alpha }}$|⁠, which is the same for all cluster masses. We label the average |$L_\mathrm{bol}/L_{\rm{H\alpha }} \sim 88$| between t = 0 Myr and t_max = 4 Myr, where t_max is the time for |$L_{\rm{H\alpha }}$| to drop by half an order of magnitude from the zero-age main sequence (see Section 3.5.1). We also label the peak in the mass loss rate seen at ∼3.5 Myr that corresponds to the time at which winds are the most effective (Leitherer et al. 1999; Rahner et al. 2017, 2019), as the most massive O stars are in their Wolf–Rayet phase but have not yet exploded as supernovae. The estimate of L_bol from |$L_{\rm{H\alpha }}$| for each H ii region is used in the calculation of the direct radiation pressure (P_rad; see Section 4.1.1).

Open in new tabDownload slide

In the right-hand panel of Fig. 7, we show the ratio of the bolometric to H|$\alpha$| luminosity. We find an average value of |$L_\mathrm{bol}/L_{\rm{H\alpha }} \approx 88$| between t₀ = 0 Myr and t_max ∼ 4 Myr. To check if this is reasonable, we compare to other estimates of |$L_\mathrm{bol}/L_{\rm{H\alpha }}$|⁠. First, we can derive a robust lower limit on L_bol by assuming that the only contribution to it is from the ionising radiation of the high-mass stars. In that case, L_bol = 〈hν〉Q, where 〈hν〉 ∼ 15 eV is a reasonable estimate for the mean energy of an ionising photon, and Q is the ionization rate derived in Section 3.3 above. This gives a lower limit to the conversion factor of |$L_\mathrm{bol}/L_{\rm{H\alpha }} \approx 18$|⁠. Secondly, there is the conversion presented by Kennicutt & Evans (2012), which accounts for a stellar population that fully samples the IMF and the stellar age distribution. This conversion is given as |$L_\mathrm{bol}/L_{\rm{H\alpha }} \approx 138$|⁠. As the H ii regions in this work are assumed to be relatively young (given their bright H |$\alpha$| emission), we expect the correct conversion to be somewhere between these two estimates. Hence, |$L_\mathrm{bol}/L_{\rm{H\alpha }} \approx 88$| is reasonable for our sample, and is used throughout this work to estimate the L_bol for each H ii region.

Finally, we note that although we assume here that f_esc = 0 for simplicity, we know that in reality some ionizing photons will escape from the H ii regions into the diffuse ISM. Estimates of the average value of f_esc for a population of H ii regions vary (see e.g. the discussion in Chevance et al. 2020a) but are typically in the range f_esc = 0.3–0.6, and so accounting for this would increase our estimates of L_bol by around a factor of 2. In practice, the impact on the radiation pressure will be smaller than this, since the photons that escape from the H ii region obviously do not contribute to the radiation pressure, and so we feel justified in neglecting this complication in our current study.

3.5.2 Cluster mass, mass-loss rate, and mechanical luminosity – |$M_\mathrm{cl}, \dot{M}, \rm{ and } L_\mathrm{mech}$|

We investigate how changing the cluster mass (M_cl) of the starburst99 model varies the ionization rate (Q). We can then use the measured ionization rate of each H ii region in the catalogue to estimate its cluster mass, which we then also use to estimate its mass loss rate (⁠|$\dot{M}$|⁠) and mechanical luminosity (L_mech; also see e.g. Dopita et al. 2005, 2006). The upper left panel of Fig. 8 shows Q as a function of time (t) for a range of M_cl. Similarly to |$L_{\rm{H\alpha }}$|⁠, we see that Q is higher for higher M_cl, and suffers a strong decrease after ∼4 Myr. As before, we then average Q within a time of 0–4 Myr. This time averaged Q is then plotted as a function of the M_cl,6 = M_cl/10⁶ M_⊙ in the lower left panel of Fig. 8. Plotted in log–log space, we see the relation is linear, with a constant of Q/M_cl,6 = 10^52.5 (s⁻¹/M_⊙). We use this conversion factor with the estimate of Q (Section 3.3) to determine M_cl for each of the H ii regions in our sample (see Fig. B1). As Q is directly estimated from the observed L_H α emission, we also outline that L_H α/M_cl,6 = 10^40.5 (erg s⁻¹/M_⊙) and, for completeness, L_bol/M_cl,6 = 10^42.5 (erg s⁻¹/M_⊙).

Figure 8.

Upper panels (left to right): Early time evolution of the ionization rate (Q), mass-loss rate (⁠|$\dot{M}$|⁠), and mechanical luminosity (L_mech) from the starburst99 model (Leitherer et al. 1999) for cluster with masses of M_cl = 10⁴–10⁶ M_⊙. The shaded region represents t < t_max, where t_max is the time for |$L_{\rm{H\alpha }}$| to drop by half an order of magnitude from the zero-age main sequence (see Section 3.5.1). We label the peak in the mass loss rate seen at ∼3.5 Myr that corresponds to the time at which winds are the most effective (Leitherer et al. 1999; Rahner et al. 2017, 2019), as the most massive O stars are in their Wolf–Rayet phase but have not yet exploded as supernovae. Lower panels (left to right): The average Q, |$\dot{M}$| and L_mech between t = 0 Myr and t_max = 4 Myr plotted as a function of the cluster mass (M_cl,6 = M_cl/10⁶ M_⊙). We use the measured Q for each H ii region to infer M_cl. We then use this M_cl to estimate L_mech and |$\dot{M}$|⁠, which are used in the calculation of the wind ram pressure (P_wind; see Section 4.1.2).

Open in new tabDownload slide

In the upper central and right-hand panels of Fig. 8, we show the time evolution of the mass-loss rate (⁠|$\dot{M}$|⁠) and mechanical luminosity (L_mech), respectively. We see that |$\dot{M}$| has an overall increase relative to its zero-age main sequence when averaged over the shown time-scale of 10 Myr, which is in contrast to the sharp declines seen in the |$L_{\rm{H\alpha }}$|⁠, M_cl and L_mech. The peak seen at ∼3.5 Myr corresponds to the time at which winds are the most energetic (Leitherer et al. 1999; Rahner et al. 2017, 2019), as the most massive O stars are in their Wolf–Rayet phase but have not yet exploded as supernovae. The time-averaged |$\dot{M}$| and L_mech are shown as a function of M_cl in the lower centre and right-hand panels, respectively. We see that |$\dot{M}/M_\mathrm{cl,6} = 10^{-2.3}$| (M_⊙ yr⁻¹/M_⊙) and L_mech/M_cl,6 = 10^40.0 (erg s⁻¹/M_⊙). We use these conversion factors to estimate L_mech and |$\dot{M}$| for each H ii region within the sample, which are used in the following section to estimate the wind ram pressure (see Section 4.1.2).

4 PRESSURE CALCULATION

4.1 Internal pressure components

In this section, we will place quantitative observational constraints on the main feedback mechanisms driving the expansion of our large sample of H ii regions. We will use these constraints to then examine if the feedback mechanisms differ with evolutionary time-scale. In this section, we will also identify the local environmental conditions surrounding the H ii regions. By contrasting the internal and external properties of the H ii region, we will investigate different dependences on initial and current environmental conditions. To do so, we first determine the components of the internal pressure within an H ii region (see also Lopez et al. 2011, 2014; Pellegrini, Baldwin & Ferland 2011; McLeod et al. 2019; Barnes et al. 2020; Olivier et al. 2021).

In this work, we consider three pressure terms that can be determined from our catalogue of H ii regions:

• thermal gas pressure (P_therm),

• direct radiation pressure (P_rad),¹⁰

• wind ram pressure (P_wind).¹¹

We then assume that the total internal pressure of an H ii region is equal to the sum of these three components P_tot = P_therm + P_rad + P_wind; i.e. assuming all components act independently and combine constructively to create a net positive internal pressure. The calculation of these various internal pressure components is outlined in this section. Note that throughout this work we will refer to the pressure terms in units of K cm⁻³ or e.g. P/k_B (where |$P/k_\mathrm{B}\, [\rm{{\rm K\,cm^{-3}}}] = P/1.38 \times 10^{-16}$| [dyn cm⁻²]).

The following calculations are simplistic in the sense that they do not account for the leaking of radiation or material into the diffuse ionized gas (Kim, Kim & Ostriker 2019; Belfiore et al. in preparation), and the cancellation of radiation forces from distributed sources (e.g. Kim et al. 2018), which may act to reduce our calculated pressures (Chevance et al. 2020a). In addition, we cannot constrain the unresolved density distribution within the ionized gas (e.g. Kado-Fong et al. 2020), yet, in the previous section, we have placed limits on the various physical properties for the H ii regions assuming that they have a smooth or clumpy density profile (see Fig. 2). Throughout this next section, we continue to use these two simple assumptions when calculating the various pressure terms. For each H ii region, we define a maximum pressure (P_max) calculated for the smallest volume (i.e. using n_e,max and r_min), and a minimum pressure (P_min) calculated for the largest volume (i.e. using n_e,min and r_eff).

4.1.1 Direct radiation pressure – P_rad

4.1.2 Wind (ram) pressure – P_wind

4.1.3 Thermal gas pressure – P_therm

4.2 External (dynamical) pressure components

In this section, we outline the method used to calculate the external pressure components acting against the internal pressures (outlined above), to confine the H ii regions and limit their expansion. To do so, we use the dynamical equilibrium pressure (P_de), which is an indirect measurement of the ambient pressure consisting of the sum of thermal, turbulent, magnetic pressure, and the ambient radiation and cosmic rays (see e.g. Kim, Ostriker & Kim 2013; Kim & Ostriker 2015). The most simplistic ‘classic’ form of P_de includes the gas self-gravity and the weight of the gas in the potential well of the stars, and is commonly adopted within the literature (e.g. Spitzer 1942; Elmegreen 1989; Elmegreen & Parravano 1994; Gallagher et al. 2018; Schruba et al. 2019; Barrera-Ballesteros et al. 2021a).

The kpc-scale molecular gas surface density, Σ_mol,1 kpc, is calculated from the CO (2–1) intensity I_CO,1 kpc from the PHANGS-ALMA survey (Leroy et al. 2021a), assuming a constant CO (2–1)/(1–0) ratio of 0.7 (den Brok et al. 2021), and the metallicity-dependent CO-to-H₂ conversion factor (α_CO) described in Sun et al. (2020). Radial metallicity measurements were estimated using the galaxy mass–metallicity relation reported by Sánchez et al. (2019), and a universal radial metallicity gradient (Sánchez et al. 2014). The kpc-scale atomic gas surface density, Σ_atom,1 kpc, is calculated from the H i 21 cm line intensity I_{H i,1 kpc} using data from the PHANGS-VLA project (PI: D. Utomo), the EveryTHINGS project (PI: K. Sandstrom), as well as existing data from VIVA (Chung et al. 2009), THINGS (Walter et al. 2008), and VLA observations associated with HERACLES (Leroy et al. 2013). The kpc-scale stellar mass surface density, Σ_*,1 kpc, is calculated from the (dust-corrected) 3.6 |$\mu$|m specific surface brightness I_3.6,1 kpc from Spitzer, assuming a mass-to-light ratio of 0.47 (McGaugh & Schombert 2014). All surface densities were corrected for the projection effect from the galaxy inclination.

Sun et al. (2020) provide the values of P_de averaged within kpc-sized apertures for 12 of the 19 galaxies studied in this work. Of the 7 galaxies without P_de measurements, NGC 1365, NGC 1433, NGC 1512, NGC 166, NGC 1672, and NGC 7496 have no available high-resolution H i observations, and IC5332 lacks any significant CO (2–1) emission in the PHANGS-ALMA data (Leroy et al. 2021a).

In this work, we want to compare the estimated internal pressure in each H ii region to this kpc-scale estimate of P_de. We note, however, that multiple H ii regions could be located in the same kpc-sized aperture, in which case the P_de values used for such comparison are identical. Moreover, we note that these kpc-scale estimates do not account for the smaller scale density fluctuations on the scales of the H ii regions. Sun et al. (2020) did introduced a modified, cloud-scale dynamical equilibrium pressure, 〈P_de,θpc〉_1 kpc (their equation 15), which treats the clumpy molecular ISM and diffuse atomic ISM separately, allowing them to have a different geometry (also see e.g. Ostriker et al. 2010; Schruba et al. 2019). Sun et al. (2020) find that the 〈P_de,θpc〉_1 kpc range between factors of 2 and 10 higher than P_de. Either P_de or 〈P_de,θpc〉_1 kpc could be an appropriate estimate of the ambient pressure of an H ii region depending on its location; if embedded inside a cloud, then 〈P_de,θpc〉_1 kpc may be relevant (i.e. akin to the initial conditions of the H ii region), yet if outside a cloud (i.e. a more evolved state), P_de would be more appropriate. As the H ii regions studied here have been identified from H |$\alpha$| emission, they are not highly obscured, and therefore are most likely not embedded within molecular clouds. Whilst this is true for the majority of cases, there is a known crossover between the H |$\alpha$| emitting phase and the embedded phase, which typically corresponds to around a third of the total H |$\alpha$| emitting lifetime (e.g. Kim et al. 2021a). The effect of the local environment and the initial conditions of the H ii regions will be assessed in detail in future work, and here we adopt P_de for the external dynamical pressure.

5 PRESSURES COMPARISON

5.1 Global variations in the pressure components

In this section, we compare the global variations of the pressure components across and between the galaxies in our sample. In Fig. 9, we show the total distribution of each pressure component for all galaxies. In the violin plot (left-hand panel), we use a kernel density estimation (KDE) to compute the smoothed distribution for both the minimum (P_min) and maximum (P_max) pressure limits.¹² In the right two panels of Fig. 9, we compare the histogram distributions for P_min and P_max separately. We list the mean and standard deviation of each P_min and P_max pressure component across the whole galaxy sample in Table 3.

Figure 9.

Distribution of the various pressure components determined across the H ii region sample for all galaxies for assumptions of a (a) smooth and (b) clumpy unresolved density distribution (see Fig. 2). We show the internal pressure components of the direct radiation pressure (P_rad) in orange, thermal pressure from the ionized gas (P_therm) in blue, and ram pressure from the wind (P_wind) in purple. We show the distribution of the external confining pressure or dynamical pressure (P_de) in green. Left-hand panel: The KDE distributions are represented as violin plots, where the mean value of each distribution is highlighted by a horizontal bar and the width of each distribution corresponds to the logarithmic scale histogram distribution. The distributions for both the maximum pressure limit (P_max) calculated for the smallest volume, and minimum pressure limit (P_min) calculated for the largest volume are shown with dashed and solid outlines, respectively. The dashed, faded arrows highlight that P_max and P_min represent upper and lower limits on the internal pressure. Right-hand panels: We show the histogram distributions for P_min and P_max separately.

Open in new tabDownload slide

The first thing to note in Fig. 9 is the one to two orders of magnitude difference between the mean values of P_min and P_max. To understand this difference, we outline here how the P_max/P_min are related. For radiation and winds, P_max/P_min = (r_eff/r_min)² = (n_e,max/n_e,min)^4/3, and from Fig. 6 we see that there is around 1.5 dex difference between the centres of the n_e,min and n_e,max distributions. This would then translate to the two orders of magnitude ratio of P_max/P_min for direct radiation and winds seen in Fig. 9 (1.8 dex given in Table 3). For the thermal pressure, P_therm,max/P_therm,min = n_e,max/n_e,min, which then would also be consistent with the 1.2 dex difference observed in Fig. 9 (also see Table 3).

We now compare the relative difference between the various pressure components considering either their maximal (P_max) or minimal (P_min) values. In Fig. 9, we see that the maximum internal pressures are all relatively similar, with mean values of around P_max/k_B ∼ 10⁶ K cm⁻³. On the other hand, the minimum values appear relatively different, with the direct radiation and wind pressures having values around P_min/k_B ∼ 10⁴ K cm⁻³ and the thermal pressures being around a factor of 4 higher (P_therm,min/k_B ∼ 10^4.6 K cm⁻³). Comparing to the external dynamical pressure, we see that that typically P_min < P_de < P_max. Interestingly, we see that values of P_de determined towards those H ii regions with P_max measurements are slightly (0.3 dex) higher than towards those with P_min measurements.

Lastly, in Fig. 9, we compare the external pressure (P_de) associated with each H ii region. Note that as each H ii region may not have both P_min and P_max measurements (see Table 2), the differences in the P_de distributions are caused by these different samples rather than a difference in the P_de measurement method (Section 4.2). In addition, several galaxies were not included in the sample from Sun et al. (2020) due to no available H i observations or due to the lack of CO significant emission (e.g. IC 5332), and hence have no P_de measurement. We see that the mean of our P_de distribution is higher by ∼0.5 dex than that shown in fig. 1 of Sun et al. (2020). This is due to the fact that the majority of H ii regions within our samples are identified towards the spiral arms and centres of the galaxies, which have systematically higher P_de values than the galaxy averages. This can be seen in the upper right panel of Fig. C1, which shows the P_de apertures taken from Sun et al. (2020) overlaid with the H ii region sample within NGC 4321.

In Fig. 10, we show the pressure components determined within the H ii regions for each galaxy. The upper and lower panels show the separate distributions for P_max and P_min, respectively. Here, we again see that all the pressure terms are similar for P_max, yet for P_min, P_therm are consistently larger than P_rad and P_wind. Moreover, we see that in general P_min < P_de < P_max for individual galaxies. We do not see any significant deviations from these trends within the individual galaxies, and are careful to compare the distributions between galaxies given the systematic biases of our H ii sample. We list the mean and standard deviation of each pressure component for individual galaxies in Table 3.

Figure 10.

Distribution of the various pressure components determined across the H ii region sample within individual galaxies for assumptions of a (a) smooth and (b) clumpy unresolved density distribution (see Fig. 2). We show the internal pressure components of the direct radiation pressure in orange (P_rad), thermal pressure from the ionized gas in blue (P_therm), and ram pressure from the wind in purple (P_wind). Note that several galaxies were not included in the sample from Sun et al. (2020) due to no available H i observations or due to the lack of significant CO emission (e.g. IC 5332), and hence have no P_de measurement. We show the distribution of the external confining pressure, or dynamical pressure (P_de), in green. The KDE distributions are represented by violin plots, where the width of each distribution corresponds to the logarithmic scale histogram distribution. Upper panel: The distributions for the maximum pressure limit (P_max) are shown in colour, and minimum pressure limit (P_min) in faded grey. Lower panel: The distributions for the P_min are shown in colour, and P_max in faded grey.

Open in new tabDownload slide

5.2 Pressure components as a function of size and position

In this section, we assess how the various pressure components vary as a function of the sizes of the H ii regions and their position within the host galaxies. In Fig. 11, we show the minimum and maximum pressure limits for P_rad, P_therm, and P_wind as a function of the radius of the H ii regions, where P_max is plotted at r_min and P_min is at r_eff. Due to the high density of individual measurements on this plot, we show the KDE distribution as contours that increase to include 99, 90, 75, 50, and 25 per cent of the data points of each P_max or P_min.

Figure 11.

Pressure components as a function of the radius (log–log axes). We show the internal pressure components of the thermal pressure from the direct radiation pressure (P_rad), the ionized gas (P_therm), and wind ram pressure (P_wind) as orange, blue, and green contours, respectively (left-hand to right-hand panels). The contours show the Gaussian KDE distribution in levels that include 99, 90, 75, 50, and 25 per cent of the points. We show the distributions for both the maximum pressure limit (P_max) calculated for the smallest volume, and minimum pressure limit (P_min) calculated for the largest volume (as labelled).

Open in new tabDownload slide

The first thing to note in Fig. 11 is that the distributions for P_rad and P_wind are very similar; albeit P_rad ≳ P_wind. This is due to the fact that they are both calculated using the H |$\alpha$| emission and have the same radial dependence of r⁻² (Section 4). On the contrary, P_therm uses the n_e calculated from the [S ii] line ratio, and the pressure calculation has no radial dependence, hence it is independent of P_rad and P_wind. The distribution of P_therm in Fig. 11 is therefore different to both P_rad and P_wind.

It is worth quickly reviewing the biases within our H ii region samples before continuing the discussion of Fig. 11 further. First, we identified our initial sample of H ii regions using an automated algorithm on H |$\alpha$| emission observations, which have a finite resolution and sensitivity. Hence, this means we could be missing detections of unresolved and weak H ii regions, or multiple compact and clustered H ii regions within complex environments (e.g. galaxy centres). Moreover, we may be missing lower surface brightness, larger H ii regions. In other words, our samples of H ii regions are biased to the brightest H ii regions across the galaxies (see Section 2.3). Secondly, in determining P_max or P_min, we have split this initial sample into two sub-samples: those with r_min and/or r_eff measurements, respectively (see Table 2 for these sample sizes within each galaxy). In the case of P_max, for example, this requires determination of n_e,max, which is calculated from the |$R_{[\rm{S}\, \small {II}]}$| ratio (Section 3.2). Therefore, the distribution of P_max does not include sources for which this ratio is statistically indistinguishable from the low-density limit, and hence measurements of P_max for larger, lower density H ii regions will be missing. We are then cautious in drawing conclusions for the radial size dependence of the pressure terms due to the number of biases affecting the size distribution of our H ii region samples.

The above discussion may then explain the different dependences we observe between P_max or P_min for each pressure component. We see that P_max suffers a moderate decline with increasing H ii region radius, which one may expect as e.g. P_rad ∝ r⁻² by definition (Section 4.1.1). P_min, on the other hand, shows an increase with increasing H ii region radius. The reason for this increase is not clear, but we speculate that this could be a result of the lack of either large, diffuse H ii regions or compact, clustered H ii regions in our sample.

In Fig. 12, we show how the pressure components for each galaxy vary as a function of the galactocentric radius normalized to R_eff (see Table 1). In this figure, the size of each point is proportional to the r_eff or r_min of each H ii region. We see that our H ii region sample spans log (R_gal/R_eff) of −2 to 1, and hence covers a large range of galactic environments; from central molecular zones to outer edges of discs (see also Fig. 1).

Open in new tabDownload slide

Open in new tabDownload slide

Figure 12.

Pressure components as a function of the galactocentric radius (normalized to the R_eff radius; see Table 1) for each galaxy in the sample. We show the internal pressure components of the direct radiation pressure (P_rad), the thermal pressure from the ionized gas (P_therm), and wind ram pressure (P_wind) as orange, blue and green points, respectively (left-hand to right-hand panels). We show the external confining pressure, or dynamical pressure (P_de), as the green points. The size of each point has been scaled to the radius of the H ii region, and corresponds to the scale shown in the lower left of each panel. We highlight the points that correspond to the maximum pressure limit (P_max) and minimum pressure limit (P_min; Section 4.1) as circles and squares, respectively.

In the majority of galaxies, we see that the external pressure (P_de) shows a systematic increase by several orders of magnitude towards the centres. A systematic increase in P_de is expected as the gas and stellar surface densities increase towards galaxy centres. However, notable exceptions are NGC 2835 and NGC 5068 that appear to have a relatively constant P_de across the whole disc (within a 1 dex scatter), which could be a result of these having lower than average atomic, molecular and stellar masses for the sample. Moreover, we see that NGC 3627 has a large scatter P_de (around 4 dex) within the disc, which could be a result of the strong bar (e.g. Bešlić et al. 2021), or the strong ongoing tidal interaction in the Leo triplet (e.g. Zhang, Wright & Alexander 1993). Again, however, we caution any interpretation of the galaxy-to-galaxy variations seen here, given the systematic biases affecting the H ii region sample (Section 2.3 and 2).

Interestingly, we also see that both the internal P_max or P_min pressures generally show systematic increases towards galaxy centres (e.g. see NGC 1365 and NGC 4535), albeit with some significant scatter within discs (e.g. NGC 1566). The increase in P_rad, P_therm, and P_wind towards centres is however smaller than the relative increase in P_de. For example, in the case of NGC 4321, within the disc at R_gal/R_eff = 1, P_de ∼ 10^4.5 K cm⁻³ and P_rad,max ∼ 10⁶ K cm⁻³, while near the centre at R_gal/R_eff = 0.1, P_de ∼ 10^7.0 K cm⁻³ and P_rad,max ∼ 10^6.5 K cm⁻³ (see also NGC 3351 and NGC 4254). This is then a factor of ∼300 increase in P_de towards the centre, yet only a factor of ∼3 increase in P_rad,max (similarly minor relative increases are observed for P_therm and P_wind, and the measurements for P_min).

It is not entirely clear why H ii regions should be more highly internally pressured within galaxy centres. For example, galaxy centres typically have higher metallicities, and hence cooling is more efficient within H ii regions and electron temperatures are lower. On the other hand, the more highly pressured environment causes higher local gas densities (see for electron density measurements e.g. Herrera-Camus et al. 2016). The latter could plausibly be an order of magnitude or more (i.e. effecting both P_therm and P_rad), whilst the former is at most a factor of two (P_therm ∝ T_e). From an observational side, we have the most trouble getting consistent boundaries (and hence sizes) for the H ii regions within the centres, as they’re often quite clustered and sitting on a high diffuse ionized gas background (Santoro et al. 2021). If larger objects are preferentially identified in the centres (i.e. smaller regions merged into one larger region), then this would act to increase the relative difference in the pressure compared to the discs. One thing to note is that r_min never gets below ∼10 pc, but in the Galactic Centre H ii region sizes are at most a ∼few 1pc in size (e.g. Barnes et al. 2020). In Section 5.5, we return how this can be addressed in future.

5.3 Total internal pressure as a function of external pressure

Figure 13.

Total internal pressure (P_tot = P_therm + P_wind + P_rad) of each H ii region as a function of its external or dynamical pressure (P_de) for assumptions of a (a) smooth and (b) clumpy unresolved density distribution (see Fig. 2). The contours show the Gaussian KDE distribution in levels that include 99, 90, 75, 50, and 25 per cent of the points. Above the dashed line the H ii regions would be overpressured (blue shaded), and below the H ii regions are underpressured (orange shaded). Upper panel: We show the distributions for the maximum pressure limit (P_max) calculated for the smallest volume. Lower panel: We show the minimum pressure limit (P_min) calculated for the largest volume.

Open in new tabDownload slide

We see that P_de spans four orders of magnitude in Fig. 13, whereas P_tot covers only one and two orders of magnitude for P_tot,min and P_tot,max, respectively. In addition, here we see a very gradual increase of both P_tot limits as a function of P_de. This significantly larger range of P_de compared to P_tot, and the tentative correlation between the two pressures, is suggestive that the ambient environmental pressure potentially has only a minor effect in regulating the internal pressures of H ii regions. For example, we posit a scenario where a high ambient environmental pressure could confine an H ii region, and therefore cause the H ii region to become more highly pressured for a given size. If this is indeed the case, this effect would be relatively minor for the larger H ii regions we observe. There is a potential caveat to discuss here, however, that the P_de the H ii regions feel is different from what we consider in this work (see Section 4.2). Here we use the kpc-scale average P_de from Sun et al. (2020) (i.e. P_de,1 kpc; which would be relevant if H ii regions are located randomly within the ISM disc), and not the cloud-scale average 〈P_de,θpc〉_1 kpc (which would be relevant if most H ii regions are within/near ISM overdensities). 〈P_de,θpc〉_1 kpc would have less radial variation, as at larger radii ISM overdensities are less common and, therefore, more impactful when calculating 〈P_de,θpc〉_1 kpc (i.e. because of the luminosity weighting) than at small radii where the ISM is more densely packed at a fixed measurement scale of θ ∼ 100 pc (Sun et al. 2020). As the H ii regions studied here are large and, therefore, evolved, we expect this effect to be significant for a small sample that will be investigated further in a future work.

Focusing on the P_tot,min limits in Fig. 13 (lower panel), we see that just under half of the H ii regions (1927 in total) are overpressured relative to their environment. These H ii regions would therefore still be expanding, despite their large measured sizes of several tens to a few 100 pc (r_eff). These large and overpressured H ii regions could then be expanding into superbubble-like structures, which are interesting targets to study the effect of large-scale expansion in the (ionized, atomic, and molecular) gas spatial and kinematic distributions. The remaining 2206 H ii regions appear to be underpressured relative to their environment, highlighting that these have most likely stopped expanding. If they have not already, these H ii regions will begin to dissipate without further energy and moment injection from young stars. The low-density cavities of these large H ii regions present the perfect environments into which future SNe can quickly expand.

Comparing to the P_tot,max limits in Fig. 13 (upper panel), we see that the majority of H ii regions are now overpressured. This is expected as P_tot,max is estimated for H ii region size scales of a few to a few tens of parsecs, and hence H ii regions that are still relatively young and expanding. Interestingly, there is a small number of H ii regions (15) that have P_tot,max < P_de and are therefore underpressured relative to their environment. To determine where these H ii regions reside within each galaxy, in Fig. 14, we show the ratio of the internal pressure over the external pressure, log (P_tot/P_de), as a function of the galactocentric radius (also see Table 3). The horizontal dashed line shows where both P_de and P_tot are equal, log (P_tot/P_de) = 0, and the region where P_de < P_tot is shaded blue and P_de > P_tot is shaded in orange. Here, the size of each point has been scaled to the effective radius (r_eff) of the H ii region, and corresponds to the size scale shown in the lower left of each panel. Note that the galaxies IC 5532, NGC 1365, NGC 1433, NGC 1512, NGC 1566, NGC 1672, and NGC 7496 have been omitted from this analysis due to their lack of available P_de measurements.

Open in new tabDownload slide

Open in new tabDownload slide

Figure 14.

Overpressure, P_tot/P_de = (P_therm + P_wind + P_rad)/P_de, of each H ii region as a function of the galactocentric radius for each galaxy for assumptions of a (a) smooth and (b) clumpy unresolved density distribution (see Fig. 2). We plot the galactocentric radius normalized to the R_eff radius (see Table 1). The size of each point has been scaled to the radius of the H ii region, and corresponds to the scale shown in the lower left of each panel. The horizontal dashed black line shows where the external pressure is equal to the internal pressure. Above the dashed line the H ii regions would be overpressured (blue shaded), and below the H ii regions are underpressured (orange shaded).

Fig. 14 shows that log (P_tot/P_de) systematically increases with increasing galactocentric radius across the sample. With the exception of NGC 3627, we see that in the six galaxies with P_tot,max < P_de, this occurs at a radius around R_gal/R_eff < 0.1, approximately corresponding to the central <1 kpc (see Table 1). This then highlights that centres are interesting high-pressured regions in which to assess the effects of stellar feedback (see e.g. Barnes et al. 2020). In the case of the strongly barred galaxy NGC 3627, we previously mentioned the large scatter in the P_de measurements at a galactic radius are coincident with the prominent bar-end features (see e.g. Beuther et al. 2017; Bešlić et al. 2021). The build-up of gas at the bar-end regions causes an increase in the gas density, and hence an increase in the dynamical pressure similar to that within the galaxy centres. It is interesting to then assess if, more generally, the increase in gas density towards the galaxy centres and bar-end regions has a significant effect on the over- (under-) pressure of a H ii region.

Fig. 15 shows the overpressure of each H ii region as a function of the total gas mass surface density (⁠|$\Sigma _\mathrm{gas} = \Sigma _\mathrm{H_2} + \Sigma _{\rm{H}\, \small {I}}$|⁠). Where the molecular (⁠|$\Sigma _\mathrm{H_2}$|⁠) and atomic (⁠|$\Sigma _{\rm{H}\, \small {I}}$|⁠) mass surface densities are taken from Sun et al. (2020), and have been measured over the same 1 kpc hexagonal grid as the P_de measurements. Here, we see that both the P_max and P_min distributions show a decreasing P_tot/P_de with increasing Σ_gas (modulo the alternative case described above that 〈P_de,θpc〉_1 kpc is larger than P_de). This then shows that the radial trends shown in Fig. 14 also apply between galaxies, and galaxies (or environments in general) with higher global gas surface densities are less overpressured (more underpressured). The simple interpretation of this is that how quickly/easily H ii regions can expand depends on the global gas surface density; i.e. more dense environments may inhibit rapid expansion (e.g. also see Dopita et al. 2005, 2006; Watkins et al. 2019). Studying the impact of stellar feedback (e.g. Grudić et al. 2018; Kim et al. 2018; Fujimoto et al. 2019; Li et al. 2019; Keller & Kruijssen 2020) and its effect on the molecular cloud lifecycle (Chevance et al. 2020b), in setting the initial conditions for star formation (e.g. Faesi, Lada & Forbrich 2018; Sun et al. 2018, 2020; Schruba et al. 2019; Jeffreson et al. 2020) and the subsequent star formation efficiency (e.g. Krumholz & McKee 2005; Blitz & Rosolowsky 2006; Federrath & Klessen 2012), within dense regions is particularly important, because ISM pressures observed within starburst systems, and at the peak of the cosmic star formation history, are several orders of magnitude higher than those observed in disc galaxies today (e.g. Genzel et al. 2011; Swinbank et al. 2011, 2012; Tacconi et al. 2013).

Figure 15.

Overpressure, P_tot/P_de = (P_therm + P_wind + P_rad)/P_de, of each H ii region as a function of the total gas mass surface density, |$\Sigma _\mathrm{gas} = \Sigma _\mathrm{H_2} + \Sigma _{\rm{H}\, \small {I}}$| (measured on kpc-sizescales; Sun et al. 2020). The contours show the Gaussian KDE distribution in levels that include 99, 90, 75, 50, and 25 per cent of the points. We show both the maximum pressure limits (P_max) and minimum pressure limits (P_min) as separate distributions (see Section 4.1). Above the dashed line the H ii regions would be overpressured (blue shaded), and below the H ii regions are underpressured (orange shaded).

Open in new tabDownload slide

5.4 Pressure components within the literature

We now compare the internal pressure components to those determined in previously studied samples of H ii regions taken from the literature. Although the wind pressure, and additional internal pressure components such as that from the heated dust, have been determined within the literature, the methodologies for calculating these differ; e.g. P_wind has been inferred from the shocked X-ray-emitting gas (e.g. Lopez et al. 2011, 2014). Therefore, here we focus on only the P_therm and P_rad pressure components from the literature, as these have been determined using a methodology consistent to that used in this work.

We make use of measurements of P_therm and P_rad for a sample of extragalactic H ii regions within the Small and Large Magellanic Clouds from Lopez et al. (2014), the Large Magellanic Cloud from McLeod et al. (2019),¹³ and NGC 300 from McLeod et al. (2020). We also compare to Galactic measurements focusing on H ii regions within in the central regions (R_gal < 100 pc) of the Milky Way from Barnes et al. (2020),¹⁴ and the disc of the Milky Way from Olivier et al. (2021).

In total, the literature sample comprises 293 H ii regions. This represents the entire sample of H ii regions with consistent internal pressure estimates currently available within the literature. The addition of the ∼6000 H ii regions studied in this work represent a significant advancement in the number of measurements available, and for the quantitative assessment of (large-scale) H ii region dynamics.

Fig. 16 shows relative contributions of P_therm and P_rad as a function of the radius for all the literature measurements mentioned above, and including the measurements determined in this work. In the left-hand panel, we show the sum of P_therm and P_rad, and colour the points by the reference. In the centre panel, we show only the distribution of P_therm (in blue), whereas in the right-hand panel we only show the distribution of the P_rad (in orange). For reference, we also show the distribution of the P_therm+P_rad in these panels in faded grey.

Figure 16.

The relative contributions of the thermal (P_therm) and direct radiation (P_rad) internal pressure components as a function of H ii region size (r) for a sample of galactic and extragalactic H ii regions taken from the literature. We show as points the Small and Large Magellanic Clouds (Lopez et al. 2014), the Large Magellanic Cloud (McLeod et al. 2019), NGC 300 (McLeod et al. 2020), the centre (R_gal < 100 pc) of Milky Way (Barnes et al. 2020), and disc of the Milky Way (Olivier et al. 2021). The contours show the Gaussian KDE of the distribution of pressure terms determined in this work, where the levels are as shown in Fig. 11. We show the distributions for both the maximum pressure limit (P_max) calculated for the smallest volume, and minimum pressure limit (P_min) calculated for the largest volume as labelled. Left-hand panel: We show the sum P_therm + P_rad for each of the literature samples. Centre and right-hand panels: We show the separate P_therm and P_rad distributions in blue and orange, respectively. The shaded grey distributions are of the P_therm + P_rad for reference. Shown as diagonal black dashed and dotted lines are power-law dependences of y ∝ x⁻¹ and y ∝ x⁻².

Open in new tabDownload slide

The first thing to note in the left-hand panel of Fig. 16 is that the Galactic sources have a much smaller scale and more highly pressured measurements than the extragalactic sources (r_eff ≲ 1 pc). This is likely due to two reasons. First, this could be a resolution effect. As it is possible to more easily achieve a higher spatial resolution with observations within the Milky Way due to its relative proximity compared to other galaxies. The observed extragalactic H ii regions could then fragment on smaller scales, and these H ii regions could be more highly pressured. Secondly, the H ii regions within the sample we observe could be at a later stage in their evolution compared to the Galactic samples, and hence larger and less pressured. Differentiating between these two possibilities is ultimately the aim of the future investigations discussed in the next section.

As previously noted, radii larger than ∼1 pc have somewhat larger thermal pressures compared to their direct radiation pressure. Although, this is only true when considering their P_min estimates. Whereas, their P_max are similar (see Fig. 14). On scales of 0.01–1 pc, however, P_therm and P_rad are comparable, and on the smallest scale of <0.01 pc, P_therm < P_rad (also see Barnes et al. 2020; Olivier et al. 2021). Shown in Fig. 16 are diagonal black dashed and dotted lines highlighting power-law dependences of y ∝ x⁻¹ and y ∝ x⁻² for reference (note that these are not fits to the data). We see that P_therm approximately follows a r⁻¹ relation. The P_rad follows r⁻², albeit with a significant scatter for the larger H ii regions (i.e. for the extragalactic observations). These are both in agreement with their expected radial dependences (e.g. by construction from equation 4; see Barnes et al. 2020). Moreover, these are in broad agreement with Olivier et al. (2021), who compared their Milky Way measurements to the LMC and SMC measurements from Lopez et al. (2014). These authors found radial power-law dependences for P_rad of r^−1.36 and P_therm of r^−0.74 (also thermal dust pressure is found to scale as r^−1.43), showing that the thermal pressure is typically sub-dominant on the smallest scales, yet does not decrease as rapidly with increasing size, and hence becomes mildly dominant on the largest scales. The transition in pressure terms shows that the impact of the different feedback mechanisms evolves and that, for example, direct radiation pressure is more significant at early times (e.g. Arthur & Hoare 2006; Tremblin et al. 2014). However, we note that we do not consider trapped radiation pressure in this work (e.g. in the form of dust heating), which could contribute significantly at later times (e.g. Lopez et al. 2014; Olivier et al. 2021).

5.5 Open questions and future avenues

This work represents a significant milestone in observationally quantifying the feedback properties in a large sample of evolved H ii regions. However, there are limitations to our study that leave several questions unanswered. We end this section by outlining these questions, and noting the possible future avenues for building on our analysis.

First, we had to make several simplistic assumptions about the unresolved sub-structure of the H ii regions within our sample. With these current data, it is difficult to address, for example, which of the smooth (i.e. P_min) versus clumpy (i.e. P_max) models may be favoured? Comparison to the literature trends in Fig. 16, shows that P_max is more similar to the P_therm measurements, whereas P_min could favour P_rad. In addition, the literature shows that Galactic centre H ii regions have maximum sizes of few |${\sim}1$| pc (Barnes et al. 2020), yet our r_min within these central regions is still ≳10 pc. Could we then be overestimating the sizes, particularly in this environment? This could be due to confusion from the DIG, which is particularly extended within centres, and clustering, which could cause multiple smaller H ii regions to be merged in our observations. To address these questions, higher spatial resolution data sets are required to constrain the true sizes and separations of the H ii regions, and hence allow us to place tighter constraints on the internal pressures. Such observations could be obtained from the Hubble Space Telescope, the James Webb Space Telescope (JWST), or integral field spectroscopy observations of more nearby targets (e.g. NGC 300; McLeod et al. 2020, SIGNALS; Rousseau-Nepton et al. 2019, and the SDSS-V Local Volume Mapper survey; Kollmeier et al. 2017).

Secondly, we have inferred several key properties of the stellar populations within the H ii regions using the starburst99 models. When doing so, we used a representative age range for the sample of 0–4 Myr, defined by the time for |$L_{\rm{H\alpha }}$| to drop by half an order of magnitude from the zero-age main sequence (see Section 3.5.1). Could our results differ when accounting for the H ii regions at various evolutionary stages? This could be addressed by using age estimates of the cluster associations within the H ii regions from the, e.g. PHANGS-HST survey (e.g. Lee et al. 2021). In addition, by confirming the presence of the young stellar population driving the H ii regions, we will be able to mitigate contamination from other sources of ionization in our sample (e.g. shocks; e.g. see Espinosa-Ponce et al. 2020 and references therein). Efforts to link the ionized gas properties from PHANGS-MUSE and cluster properties from PHANGS-HST are currently underway.

Thirdly, in this work, we have focused on the internal pressure components of P_therm, P_rad, and P_wind, yet could the additional contribution from, e.g. trapped radiation from heated dust be important in driving the large-scale expansion of the H ii regions (e.g. Krumholz & Matzner 2009; Draine 2011)? The inclusion of heated dust pressure is, however, difficult for distant extragalactic sources. As, for temperatures of ∼100 K the blackbody function peaks within the mid-infrared (∼30 |$\mu$|m). Hence, modelling the heated dust emission requires observations within the infrared regime. Currently, available data sets are limited in resolution (e.g. from the Spitzer and Herschel space observatories), yet this may be possible in the near future using the scheduled JWST observations for this sample of galaxies (PI: Lee).

Lastly, it would be interesting to see how the balance of internal and external pressures of the H ii regions varies with both local and global galactic environments. We could ask: does the pressure balance within H ii regions differ for arm, inter-arm and central regions? Are the H ii regions still embedded within molecular clouds, and how does this affect their expansion? To assess the effect of the local environment, one can compare to the measurements of the 〈P_de,θpc〉_1 kpc characterized at scales of θ ∼ 100 pc, which accounts for the clumpy molecular and diffuse atomic ISM (Sun et al. 2020; also see Barrera-Ballesteros et al. 2021a,b). We could then assess the pressure balance relative to 〈P_de,θpc〉_1 kpc for those H ii regions that are still potentially embedded, which could be identified as having a higher extinction or associated with CO(2–1) emission from the PHANGS-ALMA survey (Leroy et al. 2021a). To assess the effect of the global galactic environment on the H ii region properties, we can compare the pressure balance to the environmental masks produced by Querejeta et al. (2021). These masks were produced using the Spitzer 3.6 |$\mu$|m images, and differentiate stellar structures that form the galaxy centres, bars, spiral arms, and inter-arms regions. In this work, we inferred that a higher fraction of H ii regions may be underpressured within the galaxy centres, yet it would be interesting to assess if this could be found within these various other galaxy environments.

6 SUMMARY

In this paper, we compare the internal and external pressures acting on a sample of 5810 H ii regions across 19 nearby spiral galaxies. The H ii region sample is identified using Multi Unit Spectroscopic Explorer (MUSE) data taken as part of the Physics at High Angular resolution in Nearby GalaxieS survey (PHANGS-MUSE; Emsellem et al. 2021). We constrain the internal pressure components of the thermal pressure from the warm ionized gas (P_therm; Section 4.1.3), direct radiation pressure (P_rad; Section 4.1.1), and mechanical wind pressure (P_wind; Section 4.1.2), which we compare to the confining external pressure of their host environment, or their dynamical pressure (P_de; Section 4.2). With the MUSE observations, we cannot constrain the unresolved density distribution within the ionized gas, and hence we place upper and lower limits on each of the internal pressure components. The lower limit (P_min) corresponds to the assumption of a smooth density profile, where the measured radius (r_eff) is assumed to be representative of the H ii volume over which the pressure is acting (Section 3.1). The upper limit (P_max) corresponds to the assumption of a more clumpy density profile, where the minimum radius (r_min) is derived from the electron density measurement (Section 3.4). Of the sample of 5810 H ii regions studied in this work, 2238 H ii regions have both P_max and P_min measurements, whereas 3431 have only P_max, and 141 have only P_min (see Table 2). Due to our observational selection criteria (Section 3.4), these samples are biased towards the brightest and largest H ii regions within the galaxies. The main conclusions from the analyses of these samples are summarized below.

We assess the relative differences of the P_max or P_min measurements for each pressure term. We see that the maximum internal pressures are all relatively similar, with mean values of around P_max/k_B ∼ 10⁶ K cm⁻³. On the other hand, the minimum values appear relatively different, with the direct radiation and wind pressures having values around P_min/k_B ∼ 10⁴ K cm⁻³ and the thermal pressures being around a factor of 4 higher (P_therm,min/k_B ∼ 10^4.6 K cm⁻³). This shows that at best the pressure terms are comparable if they have a compact density distribution (i.e. at P_max). However, it is likely that the vast majority of H ii regions have at least some extended structure (i.e. tending to P_min) that would then cause P_therm to become dominant.

Comparison to a sample of H ii region pressure measurements available within the literature shows that on the scales of several tens to a couple of hundred parsecs P_therm is expected to be the highest internal pressure (e.g. Lopez et al. 2011, 2014; McLeod et al. 2019, 2020). In addition, we compare to H ii regions within the Milky Way and more nearby galaxies such as the LMC and SMC; combined with the presented measurements of this work, the sample covers spatial scales that span a total of six orders of magnitude (0.001–300 pc). Indeed, above scales of around 0.1 to 1 pc the thermal pressure is marginally dominant, yet below 0.1 pc the direct radiation pressure is dominant (Lopez et al. 2014; McLeod et al. 2019; Barnes et al. 2020; Olivier et al. 2021). We note that due to inconsistencies within the literature, this comparison does not include the indirect (trapped) radiation pressure from heated dust or the contribution of winds.

We compare our total internal pressures (P_tot = P_therm + P_rad + P_wind) within each H ii region to the external pressure (P_de) of their host environment, which we take directly from Sun et al. (2020, but also see Barrera-Ballesteros et al. 2021a,b for similar calculations of P_de). We see that for the P_tot,min limits (see Fig. 13) just under half of the H ii regions (1927 in total) are overpressured relative to their environment (i.e. P_tot,min > P_de). These H ii regions would still be expanding despite their large measured sizes of several tens to a few 100 pc (r_eff), and would, therefore, represent interesting targets to study the effect of large-scale expansion in the (ionized, atomic, and molecular) gas spatial and kinematic distributions. The remaining 2206 H ii regions appear to be underpressured relative to their environment, highlighting that these have most likely stopped expanding.

We find that for the P_tot,max limits, the majority of H ii regions are now over-pressured. This is expected as P_tot,max is estimated assuming an H ii region size scale typically of the order a few to a few tens of parsecs. In this case, H ii regions would be still relatively young and expanding. Interestingly, however, there is a small number of compact H ii regions (15) that are underpressured relative to their environment. Plotting the ratio of the internal pressure over the external pressure, log(P_tot/P_de), as a function of galactocentric radius (see Fig. 14), we see that the majority of these compact underpressured H ii regions reside within galaxy centres. This then highlights that centres are interesting high-pressured regions in which to assess the effects of stellar feedback (see e.g. Barnes et al. 2020). To assess the effect of environment more generally, we investigate if the increase in gas density has a significant effect on the over- (under-) pressured nature of an H ii region (see Fig. 15). We see that regions of galaxies (or environments in general) with higher gas surface densities have fewer overpressured H ii regions (and more underpressured H ii regions). The simple interpretation of this is that a more dense environment may inhibit rapid expansion, and thus limit the effect of stellar feedback. This is of particular importance not only for current-day star formation, but also has implications for cosmic time-scales, given that ISM pressures and densities observed at the peak of the cosmic star formation history are several orders of magnitude higher than those observed in disc galaxies today (e.g. Genzel et al. 2011; Swinbank et al. 2011, 2012; Tacconi et al. 2013).