© The Authors 2023

1. Introduction

Massive stars (M > 8 M_⊙) inject a large amount of energy and momentum into the interstellar medium (ISM) through several channels – ionising radiation, stellar winds, and supernova (SN) explosions (see, e.g., Mac Low & Klessen 2004; Krumholz et al. 2014; Klessen et al. 2016; Girichidis et al. 2020). This stellar output leads to the formation of expanding bubbles and superbubbles in the surrounding ISM with sizes from a few parsecs to 1–2 kpc (e.g., Oey & Clarke 1997; Nath et al. 2020; Watkins et al. 2023a). These bubbles are seen in many star-forming galaxies, including our own, and have been observed in multiple ISM phases: atomic hydrogen, as traced by the H I 21 cm line (e.g., Ehlerová & Palouš 2005; Bagetakos et al. 2011; Pokhrel et al. 2020), ionised gas (e.g., Valdez-Gutiérrez et al. 2001; Lozinskaya et al. 2003; Ambrocio-Cruz et al. 2016; Egorov et al. 2014, 2017, 2018; Gerasimov et al. 2022; Barnes et al. 2022), warm dust (e.g., Churchwell et al. 2006; Watkins et al. 2023a; Barnes et al. 2023), and even in molecular gas (Xu et al. 2020; Watkins et al. 2023b). Observationally, expanding superbubbles can be identified by their shell-like morphology and their emission-line profiles showing multiple components, corresponding to the approaching and receding sides of the bubble. When observed at low spectral and spatial resolution, the shell morphology can be lost, and the line profile broadened due to the overlapping of the individual unresolved components (e.g., Tenorio-Tagle et al. 1996; Smirnov-Pinchukov & Egorov 2021). In such circumstances, the superbubbles can be observed as regions of locally elevated line-of-sight velocity dispersion, including individual supernova remnants (SNRs, see e.g., Vasiliev et al. 2015).

Outside of the expansion of the superbubbles, the ISM itself has broad kinematics, likely turbulent in nature, though this is not easy to quantify from observations (see Elmegreen & Scalo 2004; Burkhart 2021, for reviews). The turbulent motions in nearby galaxies usually manifest themselves in the supersonic velocity dispersion of different tracers (e.g., for reviews Melnick et al. 1999; Green et al. 2014; Moiseev et al. 2015; Cosens et al. 2022; Law et al. 2022). It is still not clear what mechanisms are responsible for the observed supersonic velocity dispersion of the ISM in galaxies, although stellar feedback and the gravitational potential of the galaxy are usually considered the main probable drivers. On the one hand, the mean velocity dispersion in the discs of the galaxies correlates with the total star formation rate (SFR) and its surface density (Σ_SFR) (see, e.g., Yu et al. 2019; Law et al. 2022). It has also been shown for several samples of star-forming galaxies that at the current level of the SFR, SNe can produce sufficient energy to drive the turbulent motions in the atomic gas (Dib et al. 2006; Tamburro et al. 2009; Bacchini et al. 2019, 2020). At a high SFR, SNe produce even more energy and momentum than is required to explain the observed velocity dispersion of H I (Utomo et al. 2019; Sarbadhicary et al. 2022), with the excess of the energy likely deposited in the ionised gas or carried away by superbubble outflows. Girard et al. (2021) show that SN feedback alone can explain the observed velocity dispersion of molecular gas in their sample of galaxies with a high specific SFR, although it is not sufficient to explain ionised gas kinematics there.

In some cases, however, SNe alone do not explain the observed mean gas velocity dispersion, for example, in several gas-rich galaxies (e.g., Fisher et al. 2019 for ionised gas, Hunter et al. 2021; Elmegreen et al. 2022 for H I) and in the outer regions of some galaxies (Tamburro et al. 2009; Koch et al. 2018). In their analysis of four dwarf galaxies, Hunter et al. (2022) found that the H I velocity dispersion on scales of ∼400 pc correlates with the SFR 100–200 Myr ago, but no correlation was present between the previous star formation activity and the Hα velocity dispersion. Colman et al. (2022) conclude that stellar feedback alone is not enough to explain the ISM structure in the Large Magellanic Cloud (LMC) on scales larger than 60 pc – if the ISM structure is generated by turbulence, then another large-scale driving mechanism is needed. Krumholz & Burkhart (2016) tested models of feedback and gravity-driven turbulence and conclude that the gravity-driven model better explains the underlying observational measurements of the velocity dispersion and SFR, at least for the rapidly star-forming and high-velocity dispersion galaxies.

Several theoretical models and simulations suggest that both stellar feedback and gravitational instabilities are responsible for the turbulence of the ISM, but at different scales. According to Nusser & Silk (2022), turbulence is driven mostly by stellar feedback in normal star-forming disc galaxies, while gravitational instabilities become important at high gas surface densities (Σ_g ≳ 50 M_⊙ pc⁻²). Ejdetjärn et al. (2022) find that a difference in the velocity dispersion in warm (ionised) and cold (atomic and molecular) gas can be explained if stellar feedback is responsible for the high mean velocity dispersion in Hα, while the gravitational potential is a probable driver of that for cold gas. Gas accretion is also a possible driver of the turbulence in galaxies (e.g., Klessen & Hennebelle 2010; Forbes et al. 2023; Ginzburg et al. 2022; Sharda et al. 2023). In fact, an increase in the gas velocity dispersion and the non-circular motions, caused by the recent accretion of external gas from the intergalactic medium or nearby companions, has been observed at small (< 1 kpc) scales (e.g., Sil’chenko et al. 2019; Egorova et al. 2019; Pilyugin et al. 2021).

In modern cosmological simulations, stellar feedback is often introduced as a sub-grid model (e.g., Hopkins et al. 2011; Nelson et al. 2015; Pillepich et al. 2018), and a proper understanding of what fraction of the energy injected by massive stars into the ISM is radiated away and what fraction is retained in the ISM in the form of kinetic energy is important for calibrating these models. Different models predict different values of this mechanical energy coupling efficiency – from $\eta =E_{\text{kin}}^{(\text{gas})}/E_{\text{mech}}^{(\text{stars})}\sim 1\%$ to ∼40% (e.g., Krause & Diehl 2014; Sharma et al. 2014; Yadav et al. 2017; Lancaster et al. 2021), where 100% means that all mechanical energy injected by stellar sources is converted to kinetic and turbulent energy of the gas. In particular, Sharma et al. (2014) demonstrate that about 40% of the energy can be retained in superbubbles during a few tens of million years, while single SNe lose almost all their energy in ∼1 Myr. The particular value of η inferred from the models changes with the ISM gas density, age of the cluster, and number of OB stars, but also with the resolution of simulations (Yadav et al. 2017). It is therefore especially important to quantify this parameter from observations.

The main goal of the present paper is to establish the connection between the small-scale (∼30 − 500 pc) supersonic motions of the ionised gas (in the form of turbulence or expanding superbubbles) and the stellar population for a sample of nearby (D < 20 Mpc) star-forming spiral galaxies. We investigate 19 galaxies from the PHANGS survey¹ for which both integral-field spectroscopy (IFS) with the Multi Unit Spectroscopic Explorer (MUSE) at the Very Large Telescope (VLT; Emsellem et al. 2022) and multi-band imaging with the Hubble Space Telescope (HST; Lee et al. 2022) are available. With this dataset in hand, we can simultaneously analyse the properties of the ionised gas and the young stellar associations, and link their properties to each other. In particular, after identifying the regions dominated by non-circular supersonic motions of the ionised gas via their elevated velocity dispersion, we can compare their kinetic energy with the energy input produced by massive stars to judge the importance of stellar feedback in powering the regions of high velocity dispersion and estimate the efficiency of the mechanical stellar feedback.

The paper is organised as follows. Section 2 briefly describes the PHANGS-MUSE and PHANGS-HST observations and their data reduction. Section 3 presents the data analysis, including the procedure for identifying the regions of locally increased velocity dispersion in the ionised gas, isolating the young stellar association linked to them, and measuring the parameters used in the present paper. In Sect. 4 we consider the derived properties of the selected regions of high velocity dispersion. In Sect. 5 we discuss the results, and Sect. 6 summarises our main findings.

2. Observational data

2.1. PHANGS-MUSE

This paper analyses IFS data from the MUSE/VLT (Bacon et al. 2010) instrument as part of the PHANGS-MUSE survey (PI: Schinnerer). For details of the survey strategy, log of observations, and description of the data reduction process, we refer the reader to Emsellem et al. (2022), while here we briefly describe those aspects important for the present analysis.

2.1.1. Description of the data cubes and their analysis

The important advantage of MUSE is the simultaneous combination of high angular resolution (PSF is given in Table 1 for each galaxy), small pixel size (0$\stackrel{″}{.}$2), and large field of view (1′×1′). PHANGS-MUSE data cubes represent the mosaics of three to 15 MUSE fields for 19 nearby star-forming galaxies and cover a large part of their discs (extending to at least half of their optical radii, R₂₅). In the present paper, we use the publicly available² data cubes and maps (‘copt’ dataset with homogenised PSFs) with the distribution of properties of the ionised gas. The maps were obtained with a data analysis pipeline (DAP³; based on the GIST⁴ code by Bittner et al. 2019), which we also apply further to analyse the extracted spectra of the regions of interest. We derive the reddening-corrected fluxes (based on the observed Balmer decrement and a Cardelli et al. (1989) extinction curve), line-of-sight velocity, and velocity dispersion of the emission lines by fitting a single-component Gaussian model to the stellar population subtracted spectra. Details of the data reduction, mosaicking, and DAP are given in Emsellem et al. (2022). The typical spatial resolution of the PHANGS-MUSE data is 50 pc, and the 3σ flux sensitivity in Hα is 0.7 − 1.7 × 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻² (4 − 7 × 10³⁷ erg s⁻¹ kpc⁻²).

In the initial steps of the present analysis, we used the spatial distribution of the Hα flux, F(Hα), and the measured velocity dispersion, σ_obs(Hα), provided in the PHANGS-MUSE public data release DR1 (internal data release DR2.2). To remove the contribution of instrumental broadening, we derived the intrinsic velocity dispersion as $\sigma (\mathrm{H}\alpha )\approx \sqrt{{\sigma }_{\mathrm{obs}}{(\mathrm{H}\alpha )}^2-{\sigma }_{\mathrm{LSF}}{(\mathrm{H}\alpha )}^2}$, where σ_LSF(Hα) corresponds to the width of the MUSE instrumental profile at the observed wavelength of the Hα line using the parametrisation of the line-spread function (LSF) as a function of wavelength from Bacon et al. (2017). The mean value of σ_LSF(Hα) is 48.9 km s⁻¹ after averaging over all sample galaxies. We do not correct σ(Hα) for thermal and fine structure broadening, as their contribution is negligible compared to the values measured in this paper. The maps of the stellar mass surface density Σ_⋆ measured from the stellar continuum spectra and provided in the internal data release DR2.2 are used for the analysis in Sect. 5.2. We also analysed the integrated spectrum of each region with locally elevated σ(Hα). For this purpose, we ran the DAP on the extracted integrated spectra for each of these regions with the same settings as for the main PHANGS-MUSE products. As a result, we obtained information on the fluxes of the main diagnostic emission lines used in our later analysis (Hβ, [O III] 5007 Å, [S II] 6717, 6731 Å, Hα, [N II] 6584 Å), and their velocity dispersion (corrected for instrumental broadening in the same way as described above).

Our analysis further relies on the catalogue of nebulae detected within PHANGS-MUSE galaxies (hereinafter – the PHANGS-MUSE nebular catalogue; Groves et al. 2023). It was created by identifying H II regions using the HIIPHOT code (Thilker et al. 2000), and includes derived physical properties of these regions. In the present paper, we use the information on the electron density (derived from the ratio of the [S II] 6717/6731 Å lines using PYNEB; Luridiana et al. 2015), size, total reddening-corrected Hα flux of the region, and gas-phase metallicity (traced by oxygen abundance 12 + log(O/H) measured using the S-calibration from Pilyugin & Grebel 2016). We consider only those nebulae classified as photoionised H II regions (based on the diagnostic BPT diagrams; Baldwin et al. 1981). Details of the PHANGS-MUSE nebular catalogue compilation and the parameters are provided in Groves et al. (2023).

2.1.2. On the reliability of the MUSE velocity dispersion measurements

Studies of the velocity dispersion of ionised gas with MUSE are limited due to its relatively modest spectral resolution – the MUSE instrumental broadening is significantly larger than the typical values of σ(Hα) ∼ 15 − 30 km s⁻¹ observed in the ISM of nearby star-forming galaxies (e.g., Terlevich & Melnick 1981; Moiseev et al. 2015; Green et al. 2014; Yu et al. 2019; Law et al. 2022). This is also true for the PHANGS-MUSE galaxies – the histograms in Fig. 1 show the distribution of σ(Hα). Only in some environments, such as bars and galaxy centres, is σ(Hα) above the MUSE resolution. Our analysis is focused, however, on areas of locally elevated intrinsic velocity dispersion (exceeding 45 km s⁻¹, see Sects. 3.1 and 5.2), and thus the limited MUSE spectral resolution does not strongly impact our results. Nevertheless, it is important to estimate how precisely we can measure the relative variations of σ(Hα).

Fig. 1. Statistics of the intrinsic Hα velocity dispersion in PHANGS-MUSE galaxies (sorted by their stellar mass, see Table 1) for different environments (the blue colour is for the disc, and orange is for the central regions and bars, according to the classification from Querejeta et al. 2021). Pixels with a S/N < 30 in Hα are excluded. The blue dashed line corresponds to the mean MUSE spectral resolution (velocity dispersion) at the wavelength of Hα. The values of σ(Hα) were corrected for instrumental broadening by subtracting, in quadrature, the MUSE instrumental velocity dispersion for each galaxy.

To test the precision with which we can trace the relative changes of the Hα velocity dispersion with MUSE, we performed a Monte Carlo (MC) simulation. We generated 50 000 synthetic Gaussian line profiles having the same sampling as the MUSE data (1.25 Å per pixel) and a uniformly distributed line width (corresponding to σ = 10 − 100 km s⁻¹). These profiles were then convolved with a Gaussian with the width of the MUSE instrumental profile around its mean value at Hα (FWHM = 2.54 ± 0.05 Å) according to the parameterisation from Bacon et al. (2017). We consider a two times larger scatter in this profile (0.1 Å) given the higher spatial variation of the LSF in our data, likely due to the smaller number of rotations per pointing⁵. We then added noise (both Gaussian and Poisson) according to the adopted signal-to-noise ratio (S/N) and fitted a Gaussian⁶ to the obtained line profile. We repeated this procedure five times for different S/N and compared the measured LSF-subtracted velocity dispersion from the MC-sampled profiles to the true values. Figure 2 shows how the relative error of the measured σ(Hα) changes with velocity dispersion and with S/N (and also that for the absolute error for S/N = 30). Qualitatively, these distributions agree with what was obtained by Law et al. (2021b) from a similar test for MaNGA data. As can be seen, we can trace the velocity dispersion in the Hα line with MUSE with a precision better than 30% (median uncertainty is ∼12%) starting from σ(Hα) ∼ 40 km s⁻¹ for a S/N = 30 or higher, and increasing the S/N above 30 does not improve the precision significantly even for much higher values of σ(Hα). Panel b shows that for the regions having σ(Hα) > 45 km s⁻¹ and S/N > 30, we can expect the precision of the σ(Hα) measurements to be better than ∼5 km s⁻¹, while the velocity dispersion is typically significantly overestimated with MUSE for σ(Hα) < 25 km s⁻¹ (typical values for H II regions). We also note that while the MUSE instrumental profile is quite well described by a Gaussian, the true instrumental LSF is slightly more square in shape according to Bacon et al. (2017). This effect is not considered here but in principle, can produce slight additional systematic offsets of the measurements.

Fig. 2. Relation between the recovered σ(Hα) and the true σ(Hα) depending on the S/N. Panel a: relative errors for various S/N. Coloured areas encompass the 5th–95th percentile interval and dotted lines trace the median values. Panel b: probability density of absolute errors for S/N = 30. The adopted MUSE resolution at the wavelength of Hα is FWHM = 2.54 ± 0.1 Å (according to Bacon et al. 2017, but with twice larger scatter; vertical dashed line on panel a).

2.2. PHANGS-HST

All 19 galaxies considered in this paper have publicly available HST observations obtained as part of the PHANGS-HST survey (PI: Lee). Some of the data⁷ were obtained within the LEGUS (Legacy ExtraGalactic Ultraviolet Survey) programme (Calzetti et al. 2015). Images in 5 broad-band filters were obtained with Wide Field Camera 3 (WFC3) or Advanced Camera for Surveys (ACS) in bands roughly corresponding to the NUV, U, B, V, and I bandpasses. A complete description of the PHANGS-HST survey and details on the data reduction are provided in Lee et al. (2022).

At the HST resolution (PSF ∼ 0.08″ = 2.0–7.6 pc at the distance of our galaxies), stellar associations can be identified by morphologically associating neighbouring young (NUV-bright) sources. Larson et al. (2023) have constructed a catalogue of stellar associations, which were defined as local peaks in the NUV band, together with location masks derived by applying a watershed algorithm. Multi-scale catalogues were produced for fixed scales of 8, 16, 32 and 64 pc. Here we use the catalogue corresponding to the 32 pc spatial scale, which is typically a factor of two better than the MUSE spatial resolution, and thus allows us to cleanly associate the ionised gas with individual young stellar associations (the choice of spatial scale does not affect the results presented below).

In addition to the young stellar associations, we also use the catalogue⁸ of compact star clusters described in Whitmore et al. (2021) and Thilker et al. (2022). From that catalogue, we selected only those human-classified clusters of Class 1 (symmetric, centrally concentrated), Class 2 (asymmetric, centrally concentrated) and Class 3 (compact stellar associations) that are not part of any of the young stellar associations from Larson et al. (2023).

The catalogues of stellar associations and compact star clusters provide information on their mass and age, computed from fitting the observed spectral energy distribution (SED) in all 5 bands with CIGALE (Boquien et al. 2019) assuming theoretical models appropriate for a single stellar population from Bruzual & Charlot (2003). The details of the SED fitting of the PHANGS-HST data are given in Turner et al. (2021). In the present analysis, we use the measurements of the total stellar mass and the luminosity-weighted age of each association (or star cluster), and the corresponding stellar association masks. We note that both catalogues (stellar associations and compact star clusters) are used for the analysis, although for brevity we sometimes refer to stellar associations or star clusters only. The completeness limit for the young clusters and associations catalogues is about 3 − 5 × 10³ M_⊙ at the typical distance of PHANGS galaxies (Turner et al. 2021).

3. Data analysis

We observe that the velocity dispersion of the ionised gas is not homogeneously distributed across the discs of our galaxies. Multiple areas of spatially coherent and locally elevated σ(Hα) are observed in all 19 PHANGS-MUSE galaxies. Those areas with small angular size and high Hα surface brightness can be indicative of the presence of SNRs, Wolf-Rayet (WR) or Luminous Blue Variable (LBV) stars, or other high energy stellar sources (e.g., Moiseev & Lozinskaya 2012; Yarovova et al. 2023). More extended areas of elevated velocity dispersion are often observed in the inter-arm regions and are associated with the diffuse ionised gas (DIG, Belfiore et al. 2022). However, some of these regions are clearly connected to nearby H II regions and identify the locations of expanding ionised superbubbles. Figure 3 shows the distribution of Hα flux and velocity dispersion for one galaxy from our sample: NGC 4254. The Hα and [N II] line profiles in panel d demonstrate the reliability of the observed variations in σ(Hα). The presence of a broad underlying component or double-peaked emission lines (from the approaching and receding sides of a superbubble) is evident in the emission line profiles extracted from regions of high Hα velocity dispersion.

Fig. 3. istribution of the Hα surface brightness (panel a) and velocity dispersion (panel b) for NGC 4254 shown as an example. Contours on these panels are the same and correspond to the lines of a constant Hα brightness. Panel c demonstrates a zoom-in Hα map of the area shown by dashed rectangle in panels a, b with the same brightness limits as in panel a. As follows from panel b, velocity dispersion tends to be higher in the low surface brightness regions (dominated by DIG), but there are also areas of elevated σ(Hα) towards the bright regions. Panel d shows the examples of the Hα and [N II] λ6548,  6584 Å line profiles extracted from three regions with low (green line) and high (orange and magenta lines) Hα velocity dispersion. These line profiles corresponding to high σ(Hα) are either broadened (magenta line) or can be decomposed by two components (orange line). Exact areas of the spectra extraction (3 × 3 pixels) are shown on panel c by squares of the same colour as the spectra.

In this paper, we aim to isolate the regions of correlated supersonic motions in the ionised gas (caused by either turbulence or expanding superbubbles) within the PHANGS-MUSE data and to compare them with the presence and properties of the massive stars identified within the PHANGS-HST data. The analysis of both datasets can be briefly summarised in the following three major steps. First, we identified high σ(Hα) regions by their local excess of ionised gas velocity dispersion. Then we searched for stellar associations residing within each identified region. Finally, we derived the physical properties of each region and of the stellar associations and calculated the kinetic energy of the ionised gas and the total mechanical energy input produced by the stellar associations. These steps are described in detail in the rest of this Section.

3.1. Identification of expanding superbubbles and turbulent motions in the ionised ISM with I − σ diagnostics

The first step in this analysis is the automated identification of regions with spatially correlated σ(Hα) excess. For that, we rely on an analysis of the I − σ diagrams representing the pixel-by-pixel variation of the velocity dispersion σ(Hα) versus the logarithm of the Hα intensity, log I(Hα). This approach was introduced by Muñoz-Tuñon et al. (1996) and further developed by Moiseev & Lozinskaya (2012). As was shown in these papers, an expanding superbubble appears on the I − σ diagram as a triangle-shaped region of intermediate I(Hα) and elevated σ(Hα) in comparison to H II regions, which tend to have a relatively uniform distribution in their velocity dispersion across almost the whole range of Hα fluxes. Moiseev & Lozinskaya (2012, see also Yarovova et al. 2023) demonstrated that stellar-like objects (e.g., SNRs, WR, and LBV stars) that are strongly impacting the ISM can also be easily identified on these diagrams as elongated diagonal stripes. The I − σ diagnostics were successfully implemented for the detection and analysis of supersonic ionised gas motions in individual star-forming complexes of nearby galaxies and in irregular and blue compact dwarf galaxies based on high spectral resolution data from Fabry-Perot interferometer (FPI) and IFU (GMOS, KCWI) observations (e.g., Martínez-Delgado et al. 2007; Bordalo et al. 2009; Egorov et al. 2018, 2021; Bresolin et al. 2020; Gerasimov et al. 2022). In this paper, we apply the I − σ diagnostics in a form similar to that described in Egorov et al. (2021). The I − σ diagram constructed for NGC 4254 is shown in Fig. 4 as an example (see Appendix A for other PHANGS-MUSE galaxies). The interpretation of the different zones and the analysis of these diagrams are described in detail in Appendix A, while here we outline the main steps in the analysis.

Fig. 4. Localisation of the regions of locally elevated Hα velocity dispersion (cyan ellipses) in NGC 4254 galaxy, identified based on the ‘intensity – velocity dispersion’ (I − σ) diagnostics, overlaid on the Hα surface brightness map (left panel) and the classification map (central panel; see text). The I − σ diagram is shown in the right panel. The black solid line shows the mean value of σ(Hα) in the galaxy, and the dashed grey lines show its 1σ uncertainty. See Fig. A.1 for the rest of the PHANGS-MUSE galaxies.

Compared with previous studies, the application of I − σ diagnostics to the PHANGS-MUSE data has one obvious limitation – the spectral resolution, which can be insufficient to resolve the complex kinematics of the ionised gas in star-forming complexes. As we demonstrated in Sect. 2.1, with S/N ≳ 30 and better we are able to reliably measure the relative variations of the intrinsic velocity dispersion at the level of σ(Hα)≳30 km s⁻¹ with a typical uncertainty better than 5 km s⁻¹. To ensure those measurement uncertainties do not bias our analysis, we masked out all pixels with S/N < 30 in Hα. The majority of the areas of elevated velocity dispersion that we investigate in this paper have σ(Hα) > 45 km s⁻¹ and thus provide reliable measurements of σ(Hα) and even allow for the identification of kinematically distinct components for some of the regions (Fig. 3). Moreover, thanks to the high angular resolution of the MUSE data, the local peaks in the velocity dispersion are well distinguished, making it easier to identify expanding structures on the I − σ plot.

In contrast to the analysis of individual star-forming complexes and dwarf galaxies, massive galaxies provide a more diverse set of environments, with large variations in the properties of the ISM (i.e. temperature, density, and metallicity) that affect the gas pressure and potentially the observed mean velocity dispersion of the H II regions (e.g., Barnes et al. 2021). In galaxies with strong bars, massive bulges, or inner discs, each morphological component produces its own unique distribution of the data on the I − σ plot, making the uniform application of this diagnostic across the whole galaxy very difficult. In general, the measured σ(Hα) tends to be higher in the galactic centres (including nuclear rings, lenses etc.) and bars of the PHANGS-MUSE galaxies than in their discs (Fig. 1). To remove this bias, we exclude from further analysis all the regions related to these morphological components according to the Spitzer 3.6 μm based environmental masks defined by Querejeta et al. (2021). It is possible to extend the current study to these environments by analysing them separately, but for this work we focus only on regions in the discs. We skip the masking procedure for NGC 628 and for the two lowest mass galaxies (IC 5332 and NGC 5068) where the σ(Hα) distribution looks similar in all environments. Additional, manually defined masks were applied to NGC 1365, NGC 1672, and NGC 7496 to remove very extended regions of high σ(Hα) associated with AGN outflows (Stuber et al. 2021).

After applying these S/N and environmental masks, we constructed the I − σ diagrams for each galaxy and analysed them to select the regions with spatially correlated elevated velocity dispersion. As a short summary, we performed the following steps (see Appendix A for more details):

1. Isolated the peaks of high velocity dispersion corresponding to the DIG with surface brightness below Σ(Hα)_DIG (given in Table 1), which were excluded from the analysis;

2. Derived the mean velocity dispersion (σ(Hα)_m) and its standard deviation and defined the normal distribution of σ(Hα) around σ(Hα)_m as a function of Hα intensity for each galaxy;

3. Divided individual pixels into three major classes based on their Hα intensity and σ(Hα) relative to the defined normal distribution around σ(Hα)_m (these classes are shown by different colours in Fig. 4, and particular criteria are summarised in Table A.1);

4. Applied ASTRODENDRO (Robitaille et al. 2019) to the obtained classification map, isolating and tracing the variations of σ(Hα) for pixels with elevated velocity dispersion (Class 2 and Class 3 and SNR/WR candidates in Fig. 4). Defined the ellipses encircling the selected regions;

5. Refined the estimates of the sizes for those regions showing clear peaks (either central or associated with the swept-up shells) in the radial distribution of their Hα flux. If no such peak is found, we adopt 6.2 × σ_A as size of the regions along their major or minor axes, where σ_A is the standard deviation along these axes of the 2D distribution resulting from the statistics computed by ASTRODENDRO. Significantly overlapping or elongated regions are excluded.

As a result, we selected ∼30 − 130 regions of locally elevated σ(Hα) in each of the 19 PHANGS-MUSE galaxies (exact numbers are given in Table 1), and 1484 regions in total. This sample is not complete (some candidates were excluded as they do not pass all our criteria, and also better resolution data would allow us to detect more regions); however, the number of false-detections is negligible (< 3 regions per galaxy or ∼1% in total, see Appendix A).

3.2. Connecting the identified regions with the stellar associations

As a next step, we searched for young stellar associations and compact star clusters which might be connected to the identified regions of locally elevated ionised gas velocity dispersion. For each region, we identified the stellar associations where the masks (or position of the star cluster) overlap with the elliptical borders calculated in the previous section. We also calculated the overlapping fraction as the ratio of the area of the stellar mask within the region to the total area of the association. Those associations residing mostly outside the elliptical borders (overlapping by less than 20%; considering a higher or lower percentage within ±10% does not change our results) were excluded from further analysis as they are unlikely to be physically associated, but rather just chance projections (this criterion was skipped, however, for the centrally concentrated regions). We also assume that the turbulent motions decay within a given timescale, which we set here equal to 10 Myr. The exact timescale depends on the spatial scale, and velocity of the gas motions: a typical value is 1 − 10 Myr, where the longest value corresponds to the characteristic crossing timescale for galaxy discs in the vertical direction (Agertz et al. 2013; Ostriker et al. 2001). Therefore we considered only the mechanical energy injection from the star clusters during the last 10 Myr (see Sect. 3.3), and excluded from the calculations any stellar associations and star clusters older than 50 Myr (∼10 Myr higher than the lifetime of an 8 M_⊙ star which still can explode as a Type II SN).

For 188 regions with high σ(Hα) we are unable to establish a link with a star cluster because they are outside the area covered by HST observations. For the rest of the regions, we find that 511 (∼39% of the sample covered by HST) are associated with at least one stellar association or star cluster. This is almost twice what would be expected if the same stellar associations were randomly distributed across the discs. About 70% of these regions (364) contain only relatively young stellar associations (younger than 10 Myr), and 264 contain stellar associations younger than 4 Myr.

3.3. Estimation of physical parameters for the identified regions and stellar associations

For each identified region with high σ(Hα) we extracted an integrated spectrum across an elliptical aperture covering a third of the region’s effective radius. In choosing the central aperture we tried to minimise any contamination from nearby H II regions or the denser rims in cases where the region is a superbubble (to isolate the emission from its approaching and receding sides). The integrated spectra were analysed with the DAP in the same way as was done for the individual pixels in the PHANGS-MUSE datacubes (see Sect. 2.1 and Emsellem et al. 2022). From these spectra, we extracted the Hα surface brightness and emission line fluxes to characterise each region. We corrected their values for interstellar extinction using the observed Balmer decrement and the (Cardelli et al. 1989) reddening curve with R_V = 3.1 as for our Galaxy, and assuming an intrinsic ratio of I(Hα)/I(Hβ) = 2.86 for case B recombination (Storey & Hummer 1995).

The oxygen abundance 12 + log(O/H) is taken from the smoothed two-dimensional distributions published by Williams et al. (2022). These were derived by applying Gaussian Process Regression to the metallicity measurements made for individual H II regions for all 19 PHANGS-MUSE galaxies. These values rely on the S calibration from Pilyugin & Grebel (2016), requiring measurements of the [S II] λ6717, 6731 Å, [N II] λ6584 Å, [O III] λ5007 Å, Hα and Hβ emission lines. We use these estimates instead of obtaining metallicity measurements from the integrated region spectra, as the H II regions are less biased by the DIG, although both methods gave us consistent results.

3.3.1. Mass of the ionised gas

We consider two ways to estimate the total mass of the ionised gas, M_ion, in the regions. The first one assumes a thin spherical shell geometry for all of them. Then

$$\begin{array}{c}\hfill M_{\mathrm{ion}}^{(1)}\simeq 4\pi R_{\mathrm{eff}}^2\mathrm{\Delta }R{n}_{\mathrm{e}}\mu m_{\mathrm{H}}=4{\mathrm{\Omega }}_{\mathrm{reg}}{D}^2\mathrm{\Delta }R{n}_{\mathrm{e}}\mu m_{\mathrm{H}},\end{array}$$(1)

where $R_{\mathrm{eff}}=\sqrt{R_{\min }{R}_{\mathrm{maj}}}$ is the effective radius of the region, computed as the geometric mean of the minor and major semi-axes of the ellipses encircling the region (see Sect. 3.1), ΔR is the thickness of the shell, n_e is the electron density, m_H is the mass of a hydrogen atom, the coefficient μ = 1.27 accounts for a singly ionised helium contribution, Ω/italic>_reg is the total area of the region in the plane of the sky, and D is the distance to the galaxy. To calculate ΔR, we can relate the emission measure EM to the surface brightness of Hα, I(Hα) as (Reynolds 1977; Pengelly 1964):

$$\begin{array}{c}\hfill EM=2.75{\left(\frac{T_{\mathrm{e}}}{{10}^4\mathrm{K}}\right)}^{0.9}I(\mathrm{H}\alpha ){\mathrm{pc}\mathrm{cm}}^{-6},\end{array}$$(2)

where I(Hα) is in units of Rayleighs (1 Rayleigh ≃ 5.69 × 10⁻¹⁸ ergs s⁻¹ cm⁻² arcsec⁻² for the Hα line), and T_e is the electron temperature. Assuming a homogeneous density distribution (that represents the average conditions), the emission measure is related to the thickness of the nebula (equal to 2ΔR for a spherical shell geometry) as $EM=\int n_{\text{e}}{n}_{{\text{H}}^{+}}\text{d}s\simeq \int n_{\text{e}}^2\text{d}s\simeq 2\mathrm{\Delta }R{n}_{\text{e}}^2$. Thus, we obtain

$$\begin{array}{c}\hfill \mathrm{\Delta }R\simeq 24.2{\left(\frac{T_{\mathrm{e}}}{{10}^4\mathrm{K}}\right)}^{0.9}\frac{F(\mathrm{H}\alpha ){\mathrm{\Omega }}_{\mathrm{c}}^{-1}}{{10}^{-16}\mathrm{erg}{\mathrm{s}}^{-1}c{m}^{-2}arcse{c}^{-2}}{\left(\frac{n_{\mathrm{e}}}{{\mathrm{cm}}^{-3}}\right)}^{-2}\mathrm{pc},\end{array}$$(3)

where F(Hα) = I(Hα)Ω/italic>_c is the total reddening-corrected Hα flux from the central aperture used in the spectral extraction, which has an area of Ω/italic>_c in the plane of the sky in squared arcseconds. Given that in our analysis we choose the central aperture for spectra extraction to have a size of 1/3 of the full size of the region (see Sect. 3.1), Ω/italic>_reg/Ω/italic>_c = 9 for every region. Assuming that T_e = 8000 K (typical for our sample of mostly slightly sub-solar metallicity H II regions, but can be higher in shocked or in lower metallicity regions; Kreckel et al. 2022), we can combine Eqs. (1) and (3) to obtain

$$\begin{array}{c}\hfill M_{\mathrm{ion}}^{(1)}\simeq 522\frac{F(\mathrm{H}\alpha )}{{10}^{-16}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}}{\left(\frac{D}{\mathrm{Mpc}}\right)}^2{\left(\frac{n_{\mathrm{e}}}{{\mathrm{cm}}^{-3}}\right)}^{-1}{M}_{\odot }.\end{array}$$(4)

Note that this method does not depend on how well we can measure R_eff for our regions, which can be especially uncertain for small, poorly resolved regions (see Barnes et al. 2021). However, it assumes a spherical shell geometry, with a uniform distribution of density (and brightness) within the shell. As previously mentioned in Sect. 3.1, at least 171 of the regions demonstrate a shell-like morphology and thus are likely superbubbles, for which this approximation is arguably valid. This is not necessarily true for all regions.

Given the spatial resolution of our data, we cannot provide better constraints on the shape of the high velocity dispersion regions. Also, the thin-shell approximation might not work for the small regions in our sample (i.e. where the size of the aperture for spectra extraction is smaller than the PSF). Therefore, we also consider a second method that does not depend on the geometry of the region but instead requires us to measure the total Hα flux from the whole region, F(Hα)_tot (measured from the spectra extracted from the elliptical apertures corresponding to the borders of a region). From atomic physics, we know that the number of hydrogen-ionising photons consumed within a region of volume V is given by $Q({\text{H}}^0)={\alpha }_{\text{B}}\langle n_{\text{H}}^2\rangle V\simeq {\alpha }_{\text{B}}\langle n_{\text{e}}^2\rangle V$ when the volume is in ionisation equilibrium, where α_B is the case B recombination coefficient and $⟨n_{\mathrm{e}}^2⟩$ is the mean squared electron density within V⁹. If we assume, in the absence of evidence to the contrary, that the density distribution is uniform, this relation simplifies to $Q({\text{H}}^0)={\alpha }_{\text{B}}{n}_{\text{e}}^{2}V$. For our adopted T_e = 8000 K, interpolation of the values in Storey & Hummer (1995) yields α_B ≃ 3.1 × 10⁻¹³ cm³ s⁻¹, and the ionising photon flux is related to the Hα luminosity by Q(H⁰)≃7.1 × 10¹¹L(Hα) (Osterbrock & Ferland 2006). Given that L(Hα) = 4πD²F(Hα)_tot, we then have

$$\begin{array}{c}\hfill M_{\mathrm{ion}}^{(2)}=\mu m_{\mathrm{H}}{n}_{\mathrm{e}}V\simeq 29.3\frac{F{(\mathrm{H}\alpha )}_{\mathit{tot}}}{{10}^{-16}\mathrm{erg}{\mathrm{s}}^{-1}{\mathrm{cm}}^{-2}}{\left(\frac{D}{\mathrm{Mpc}}\right)}^2{\left(\frac{n_{\mathrm{e}}}{{\mathrm{cm}}^{-3}}\right)}^{-1}{M}_{\odot }.\end{array}$$(5)

As our adopted value of M_ion, we use $M_{\mathrm{ion}}^{(1)}$ for the regions having a size along the minor axis greater than 3 × PSF, and $M_{\mathrm{ion}}^{(2)}$ for smaller regions. This allows us to minimise contamination by overlapping H II regions in the case of large regions, and to correctly measure M_ion for small regions where the thin-shell approximation does not work due to insufficient resolution.

3.3.2. Estimates of the electron density

Precise measurements of the electron density are crucial for deriving M_ion, but difficult because of several observational and physical limitations. It is typically measured using density-sensitive line ratios like R_[SII] = [S II]λ6717Å/[S II]λ6731Å (e.g., Osterbrock & Ferland 2006; Kewley et al. 2019). Such diagnostics have low-density limits (e.g., n_e ∼ 30 cm⁻³ for R_[SII]), below which changes in n_e have no effect on the line ratios. The densities typically measured in extragalactic H II regions are close to this low-density limit (e.g., Kennicutt 1984; Hunt & Hirashita 2009) making the calculations uncertain. In addition, density inhomogeneities in the nebulae can introduce systematic biases in the measurements of n_e and other physical properties (e.g., Peimbert 1971; Rubin 1989; Viegas & Clegg 1994; Cedrés et al. 2013). Here we use two different density diagnostics that we consider as lower and upper limits.

For both calculations of M_ion (Eqs. (4) and (5)), we estimate the electron density n_e, SII from the R_[SII] using the PYNEB package (Luridiana et al. 2015) and assuming T_e = 8000 K. The electron density of the ionised gas in our regions is close to the low-density limit. Because of this, we are able to reliably measure it only for a limited sample of ∼20% of regions in which the observed R_[SII] falls below the value of 1.46¹⁰ by at least 3σ (chosen to minimise the bias towards higher densities if R_[SII] is underestimated because of the observational uncertainties). To use all the identified regions in the remainder of the analysis, we assume they have in general the same n_e as their nearest-neighbour H II regions. We selected reference H II regions from the nebular catalogue located within two effective radii of each region centre (this distance was increased up to four radii if none were found at the smaller distance), calculated n_e, SII for each of these H II regions in the same way as above, and adopted our final $n_{\mathrm{e}}^{(\max )}$ as the median value among them. We are able to measure n_e, SII of at least one reference H II region for all but 16 of our regions, and of at least 3 for 36% of our sample. Figure 5 compares the estimates of M_ion made assuming n_e calculated in different ways. For the regions where direct measurements of n_e, SII from their spectra are possible, we obtain very similar M_ion as from the density adopted from the reference H II regions. Therefore, for consistency, we consider hereafter the electron density $n_{\mathrm{e}}^{(\max )}$ derived in the latter way.

Fig. 5. Mass of the ionised gas in the selected regions assuming different prescriptions for electron density: lower limit ($n_{\mathrm{e}}^{(\min )}$) and upper limit ($n_{\mathrm{e}}^{(\max )}$) derived from the reference H II regions using Eq. (6) and from R[S II] = [S II]λ6717Å/[S II]λ6731Å flux ratio, respectively, and ne, SII derived from R[S II] of the observed spectra of the region itself. Different morphological types introduced in Sect. 4.1 are shown by different colours. Diagonal lines show 1-to-1 relation. Spearman correlation coefficients ρ are given on the plots.

The estimates of $n_{\mathrm{e}}^{(\max )}$ obtained from R_[SII] are related to dense clumps in the nebulae where the emission of the [S II] lines originates. Thus, the derived values are not necessarily equal to the mean volume electron density n_e, but represent its upper limit because of clumping of the ionised gas. Assuming a spherical shape for the reference H II regions, we can measure the volume-averaged electron density as

$$\begin{array}{c}\hfill n_{\mathrm{e}}^{(\min )}=\sqrt{\frac{3Q({\mathrm{H}}^0)}{4\pi R_{\mathrm{eff}}^3{\alpha }_{\mathrm{B}}}}\simeq 1.36\times {10}^{-16}{\left(\frac{L(\mathrm{H}\alpha )}{\mathrm{erg}{\mathrm{s}}^{-1}}\right)}^{0.5}{\left(\frac{R_{\mathrm{eff}}}{\mathrm{pc}}\right)}^{-1.5}{\mathrm{cm}}^{-3}.\end{array}$$(6)

Applying this to the selected reference H II regions in the same way as above, we estimate $n_{\mathrm{e}}^{(\min )}$ for all regions in our sample. 53% have at least 3 reference H II regions with reliable estimates (H II regions with an effective radius R_eff < 0.5 × PSF are excluded). In the ideal case, $n_{\mathrm{e}}^{(\min )}$ yields the true value of the mean volume density, however, only about 25% of the H II regions in PHANGS-MUSE galaxies are well-resolved and have reliable measurements of R_eff (Barnes et al. 2021). Thus many of the reference H II regions in our sample are under-resolved, in which case R_eff is likely overestimated (Barnes et al. 2022). Moreover, the validity of Eq. (6) depends on the actual distribution of the ionised gas in H II regions and can underestimate the result if the ionised gas does not occupy the entire volume of the nebula. Thus the value of $n_{\mathrm{e}}^{(\min )}$ corresponds to the lower limit of the electron density.

We use the estimates of $n_{\mathrm{e}}^{(\min )}$ and $n_{\mathrm{e}}^{(\max )}$ for each region to measure the upper and lower limits of M_ion, respectively. Figure 5 shows the comparison between different methods. While the measurements are correlated, we find an order of magnitude difference between M_ion calculated with these two prescriptions for n_e. In the case of the well-resolved reference H II regions, $n_{\mathrm{e}}^{(\min )}$ should yield the volume-averaged density, less biased by the presence of dense clumps. We thus consider it as a better proxy of the true density of the ionised gas in our selected regions of high velocity dispersion. Meanwhile, $n_{\mathrm{e}}^{(\max )}$ demands a knowledge of the filling factor to be converted to a volume-averaged density, which cannot be properly derived for our sample. At the same time, $n_{\mathrm{e}}^{(\max )}$ correlates extremely well with n_e, SII derived directly from the observed spectra and unaffected by the limited angular resolution. Thus, we expect that the upper limit of M_ion (based on $n_{\mathrm{e}}^{(\min )}$) is closer to the real value, while its lower limit (based on $n_{\mathrm{e}}^{(\max )}$) better traces the relative changes between the regions. The median formal uncertainties of $log(M_{\mathrm{ion}})[n_{\mathrm{e}}^{(\min )}]$ and $log(M_{\mathrm{ion}})[n_{\mathrm{e}}^{(\max )}]$ are equal to 0.09 dex and 0.29 dex, respectively. Therefore our further analysis relies mainly on $M_{\mathrm{ion}}[n_{\mathrm{e}}^{(\min )}]$, but we compare the results with measurements made assuming $M_{\mathrm{ion}}[n_{\mathrm{e}}^{(\max )}]$.

3.3.3. Kinetic energy of the ionised gas and mechanical energy input from stars

Assuming that the measured elevated Hα velocity dispersion relative to the mean value in a galaxy (or in the immediate surroundings in the more general case) is a proxy for turbulent motions or is associated with gas expansion, we can estimate the kinetic energy of these supersonic motions in the ionised gas as:

$$\begin{array}{c}\hfill E_{\mathrm{kin}}\simeq E_{\mathrm{turb}}\simeq \frac{3}{2}{M}_{\mathrm{ion}}(\sigma {(\mathrm{H}\alpha )}^2-\sigma {(\mathrm{H}\alpha )}_{\mathrm{m}}^2).\end{array}$$(7)

To estimate the energy input from massive stars, we rely on the parameters of the stellar associations, described in Sect. 3.2. The age and total mass of the stars were taken from the PHANGS-HST catalogue of stellar associations (see Sect. 2.2). Then we used models produced by STARBURST99 (v.7.0.1, Leitherer et al. 1999, 2014) to calculate the cumulative mechanical energy produced by the winds and supernovae from a star cluster with parameters (mass, age, metallicity) matched to each of the selected stellar associations. Namely, we calculated a grid of models assuming instantaneous star formation, a Kroupa (Kroupa 2001) initial mass function (IMF)¹¹ and Geneva stellar evolution tracks which include rotation with two metallicities (Z = 0.014 ∼ Z_⊙, Ekström et al. 2012, and Z = 0.002 ∼ 0.15 Z_⊙, Georgy et al. 2013), for clusters with total stellar mass log(M_*/M_⊙) = (3.0…7.0) with steps of 0.25 dex and age varying from 0.2 to 60 Myr with steps of 0.2 Myr. The output parameters from this grid are the total mechanical energies produced by both stellar winds and supernovae (E_WSN) and by winds only (E_W) during the whole lifetime of a cluster. The evolution of these cumulative energies is dependent on the input parameters and is shown in Fig. 6 for several models as examples. All the regions analysed in this paper have metallicity significantly higher than the value used for the low-metallicity grid, and only 86 regions (4%) have 12 + log(O/H) < 8.35 (Z < 0.45 Z_⊙). Therefore we use a grid of models computed for solar metallicity to evaluate the energy input for all the clusters. Given the differences in the energy input between the high- and low-metallicity models (Fig. 6), we can expect an overestimation of the cumulative E_mech for some of the young clusters by a factor of few. The results do not change significantly if using the interpolated values between these two available grids. The exact values of E_WSN and E_W were obtained by interpolating the model values to the measured mass and age of each association. As we explained in the previous Section, for further analysis we take into account only values of E_WSN and E_W injected by star clusters during the last 10 Myr of their evolution. Accounting for the total energy does not significantly change the results as the majority (77%) of the star clusters in our analysis are younger than 10 Myr.

Fig. 6. Evolution of the cumulative mechanical energy input to the ISM (Emech) normalised by the mass of star cluster (Mcls) for different metallicity (traced by different colours) according to STARBURST99 (Leitherer et al. 1999, 2014) models. The contribution produced by stellar winds alone is shown by the solid line, while the dashed line corresponds to the impact of both stellar winds and supernovae.

4. Properties of the ionised gas in regions with high velocity dispersion

Following the analysis described in the previous section, we selected 1484 regions of locally elevated velocity dispersion from all 19 galaxies. Here we consider the properties of these regions and classify them according to their most likely powering mechanism.

4.1. Morphology, ionisation state, and links with stellar associations

Based on the results of our automated technique for refining the size of each region (by searching for the peak in the radial Hα flux distribution, see Appendix. A) we can divide the whole sample into three categories. In Fig. 7 we show three regions exemplifying each of these classes. The first are classified as probable expanding superbubbles, as their radial Hα flux distribution shows weak flux in the centre increasing to a peak at the edges, corresponding to a shell-like morphology. We identified 171 such regions. Many regions (572) show a central peak in the Hα flux distribution and thus are classified as compact objects, with both Hα flux and velocity dispersion elevated in their centres. This can be due to the presence of a small unresolved expanding bubble, and in particular of SNRs. Finally, those regions without identified peaks in their radial Hα flux distribution are classified as ambiguous. They still could be superbubbles where the rim is not clearly identified, or they could also represent regions of more complex morphology with elevated velocity dispersion due to locally induced turbulence by stellar feedback or gravity. Note that sometimes superbubbles are more easily identifiable by their kinematics than by their morphology, especially in an inhomogeneous environment (e.g., Nazé et al. 2001).

Fig. 7. Examples of regions classified as shells (top row), ambiguous (middle row) and compact (bottom row). Images from left to right correspond to a three-colour multiband image of the region (Hα image from HST (PI: R. Chandar, GO-10402) is in red, FUV from AstroSat (Hassani et al., in prep.) in green, CO from ALMA (Leroy et al. 2021, in blue), the Hα flux distribution, the [S II]/Hα line ratio, and the Hα velocity dispersion. Cyan ellipses show the adopted borders of the region and overlaid contours correspond to the stellar associations’ masks. The HST stellar associations shown in yellow are selected as associated with the region. The right-hand plots show the radial Hα flux distribution relative to the centre of each region (marked by crosses).

In Fig. 8 we show the statistical distribution of these three classes of regions, considering their link with the young stellar associations and the [S II]/Hα line ratio, which is a common tracer for the presence of shocks (e.g., Allen et al. 2008). We distinguish here the number of regions having [S II]/Hα above or below 0.4, which is often used to separate supernovae remnants from photoionised nebulae (Blair et al. 1981; Blair & Long 2004).

Fig. 8. Relative distribution of the number of regions with different morphological types, ionisation states and associated young stellar associations. Three groups are shown, corresponding to the three morphological types (see text in Sect. 4.1), and their sizes correlate with the relative number of corresponding regions. Different colours correspond to select properties: red and orange colours correspond to regions with [S II]/Hα > 0.4, while blue and green – to those with [S II]/Hα < 0.4; regions that have at least one stellar association within their borders are shown by blue and orange colours, while green and red correspond to those without identified young stars). Violet colour corresponds to the regions outside the HST footprint (and thus without information about star clusters there). Grey lines mark each 10% level.

We identified 511 regions (less than half of the whole sample) with locally elevated velocity dispersion that are linked with at least one young stellar association or compact star cluster. However, this does not rule out that stellar feedback is driving the high σ(Hα) in the rest of the regions – 73% of them have [S II]/Hα > 0.4, making them probable SNR candidates. Moreover, 74% of all regions with [S II]/Hα > 0.4 and covered by HST observations do not host stellar associations within their borders.

Compact regions make up 39% of our sample and most are not linked to any stellar associations. Meanwhile, about 43% of the compact regions without stellar associations have [S II]/Hα line ratios consistent with photoionisation. We suggest that at least some of these regions could be poorly resolved bubbles around individual massive O stars. Indeed, ∼86% of them have an Hα luminosity L(Ha) < 2.5 × 10³⁷ erg s⁻¹ that can be produced by a single OV5 star (adopting the rate of production of hydrogen-ionising photons from Martins et al. 2005). The nebulae surrounding such isolated stars would likely be observed as isolated compact H II regions at the ∼50 pc PHANG-MUSE resolution and contribute to the ionisation balance (Scheuermann et al. 2023), but these single stars will naturally not be detected as stellar associations. The brightest compact regions could originate from multiple stars with a single SNR – in such systems the line flux ratios might be still dominated by photoionisation. Planetary nebulae (PNe) could also be identified as compact regions with low [S II]/Hα and would not be linked with the young stellar associations, although we find that only two such regions overlap with PNe identified in PHANGS-MUSE galaxies (Scheuermann et al. 2022).

In contrast to the compact and ambiguous group, regions with shell-like morphologies are generally associated with at least one stellar association or star cluster (84% of the regions covered by HST observations). We consider these regions as our best candidates for expanding superbubbles, and discuss them in Sect. 5.3.

To determine the excitation state of regions with high σ(Hα), we plot them on diagnostic BPT diagrams (Fig. 9) together with other types of nebulae (including H II regions, SNRs and PNe) from the PHANGS-MUSE nebular catalogue (also applying the same environmental masks). Our regions tend to be shifted towards the higher ionising radiation domain (right-hand side of the diagrams). This is in agreement with the BPT-σ relation observed in nearby galaxies of different types – the Hα velocity dispersion usually increases towards the top-right side of the diagrams, which is typically explained by a growing contribution of shocks to the emission line excitation (Oparin & Moiseev 2018; D’Agostino et al. 2019; López-Cobá et al. 2020; Law et al. 2021a). Interestingly, as seen in the central row of Fig. 9, regions without corresponding star clusters preferentially occupy the upper-right part of the BPT diagrams, while regions linked to star clusters are more similar to H II regions. Furthermore, regions without star clusters have (in general) lower kinetic energy, with values close to 10⁵¹ erg, typical of the energy ejected during a SN explosion. Altogether, these findings support our claim that most of the regions not linked to young stellar associations are SNRs, and thus their high velocity dispersion results from the influence of stellar feedback.

Fig. 9. Diagnostic BPT diagrams showing the optical line ratios for all identified regions of locally elevated velocity dispersion (colour-coded circles). Colour denotes either Hα velocity dispersion (top panels), number of young stellar associations and star clusters within the region (middle panels) or estimated kinetic energy of the gas (bottom panels). For reference, we show as grey circles nebulae from the PHANGS-MUSE catalogue (including H II regions, supernovae remnants, planetary nebulae; see Groves et al. 2023) excluding those falling in the masked areas (see Sect. 3.1). Demarcation lines separate areas on the diagrams typically occupied by regions with different excitation mechanisms: the black solid curve from Kewley et al. (2001) separates H II regions from those with a dominant contribution from other excitation mechanisms (e.g., DIG, shocks); composite mechanism of excitation is probable in the area between this and the grey curve from Kauffmann et al. (2003); AGNs usually occupy the area above the black straight line from Kewley et al. (2006).

4.2. Size distribution

The distribution of the sizes (effective radii, R_eff) of the identified regions of locally elevated velocity dispersion is shown in Fig. 10. As can be seen, the sizes of the large structures (above ∼70 pc) are described very well by a power-law: $p(R_{\mathrm{eff}})\propto R_{\mathrm{eff}}^{\alpha }$. Using the POWERLAW package (Alstott et al. 2014), we got a best-fit value of α = −2.71 ± 0.19. The power-law breaks at radii of R_eff ∼ 60 − 80 pc, and also at R_eff ≳ 300 pc. We cannot properly resolve smaller structures because of the limited angular resolution. Most (∼86%) of the smallest regions (the first bin in Fig. 10) are detected in the four closest galaxies with D < 10 Mpc. We note that a log-normal distribution describes the size distribution better than a power-law at its small and large ends.

Fig. 10. Distribution of sizes for regions of high Hα velocity dispersion. The red line shows a result of a power-law fit to the data for Reff = 70 − 300 pc. This fit agrees with a power-law approximation for the sizes of superbubbles according to hydrodynamical simulations (α = −2.7 obtained by Nath et al. 2020). The distribution of larger and smaller regions is better approximated by log-normal distribution (the blue line). The dashed black line and the grey area show the median spatial resolution and its standard deviation across our sample.

The measured slope of a power-law approximating our regions in Fig. 10 is very similar to α = −2.7, which is derived from hydrodynamical simulations based on the evolution of superbubbles by Nath et al. (2020), as well as from the size distribution of the H I holes in nearby galaxies (α ≃ −2.9; Bagetakos et al. 2011; Oey & Clarke 1997). The observed size distribution is also in good agreement with a power-law (with a slope ∼ − 3) defining the ISM density inhomogeneities due to supersonic turbulence (see Burkhart 2021, for a review). Watkins et al. (2023b) analysed the ALMA CO data for the same sample of galaxies and found a slightly steeper distribution (α ≃ −3.2 ± 0.4) for the molecular gas superbubbles. Molecular gas is hard to detect in the rims of the larger bubbles, which naturally explains their steeper slope. A shallower distribution (α ≃ −2.2 ± 0.1) was measured by Watkins et al. (2023a) for superbubbles in one galaxy from our sample (NGC 628) visually identified by the morphology of the dust in PHANGS-JWST data (Lee et al. 2023). Future analysis of superbubbles identified in JWST images for other galaxies from our sample will clarify the origin of the differences (if real) in their size distribution in different tracers.

4.3. Expansion velocities and kinematic age

Assuming that the elevated velocity dispersion has resulted from an expanding bubble, we can estimate its expansion velocity (V_exp) and kinematic age (t_kin). At the MUSE spectral resolution, we generally cannot resolve the Hα line profile into kinematically distinct components corresponding to the approaching and receding sides (see, however, one example in Fig. 3), and thus direct measurements of the expansion velocity are not possible. However, we can use the information about spatial variation of σ(Hα) in order to estimate the corresponding V_exp following a recipe from Smirnov-Pinchukov & Egorov (2021; their Eqs. (7)–(10)). For a Gaussian line profile and the MUSE spectral resolution, we obtain

$$\begin{array}{c}\hfill V_{\exp }\simeq 3.15\times {(\sigma {(\mathrm{H}\alpha )}^2-\sigma {(\mathrm{H}\alpha )}_{\mathrm{m}}^2)}^{0.4}-4.1\mathrm{km}{\mathrm{s}}^{-1},\end{array}$$(8)

where the exact coefficients depend on the value of σ(Hα)_m (unperturbed velocity dispersion in the galaxy) and are parameterised in that paper¹², and σ(Hα) is measured from the integrated spectrum from the entire region.

In the classical analytical solution by Weaver et al. (1977), the kinematic age of a bubble with radius R expanding in a homogeneous medium can be derived as t_kin ≃ 0.6R/V_exp when the bubble is in the adiabatic stage of its evolution. Figure 11 (left) shows that the expansion velocities, V_exp, of the identified regions with elevated velocity dispersion decline with their kinematic age. That is expected given that the bubbles slow down once they have swept up sufficient mass during their evolution. Equation (8) and estimates of t_kin are valid under the assumption of a spherical thin bubble expanding in a homogeneous medium, which is not necessarily true for all our regions (see Sect. 4.1). We note, however, that in Fig. 11, all regions follow the same trends as those classified as ‘Shells’.

Fig. 11. Kinematic properties of the regions with elevated velocity dispersion. Left panel shows the relation between the estimated expansion velocity (Vexp) and kinematic age of each region (grey points). Blue points correspond to the regions classified as ‘Shells’, and light-blue crosses – to the values of Vexp for that sub-sample predicted from the analytical model (Eq. (9)) based on their estimated kinematic age and the properties of the star clusters. Top histograms in right panel shows the distribution of the kinematic age (blue colour) and the mean age of the clusters (weighted by their total mechanical energy input during the last 10 Myr; orange colour). Orange colour on the bottom histograms considers only clusters older than 4 Myr and corresponds to the time passed since the first SN (assuming it is happening at 4 Myr according to Fig. 6). Only the regions with at least one stellar association or star cluster are shown.

The right panel of Fig. 11 compares the distributions of the derived kinematic age (blue histogram) and mean age of the young star clusters associated with each region, weighted by their total mechanical energy input over the last 10 Myr (orange histogram on the top plot). This plot demonstrates that the typical kinematic age of the observed supersonic ionised gas motions is generally lower than the mean age of the star clusters and stellar associations residing within the considered regions. This discrepancy can be understood if the development of these ionised gas motions (either turbulent or superbubble expansion) starts due to the influence of pre-SN feedback, but their main driver is energy and momentum from SN explosions. Therefore we expect that the kinematic age should correlate better with the time passed since the first SN explosion for older regions, or with the mean age of the clusters within the region if SNe have not exploded yet and the supersonic gas motions are mostly driven by pre-SN feedback. Thus, the bottom plot of Fig. 11 (right) considers only clusters older than 4 Myr and shows their modified mean age $t_{\text{cls}}^{\text{SN}}=t_{\text{cls}}-4$ Myr (the threshold of 4 Myr approximately corresponds to the time when the first SN occurs according to the STARBURST99 models shown in Fig. 6). We find that the estimated kinematic ages are in better agreement with such modified ages of the clusters tracing the time passed from the first SN event.

To verify that the measured values of the expansion velocities are consistent with the properties of the young stellar population, we derive V_exp at a given t_kin from an analytical solution of a superbubble driven by multiple winds and SNe (Mac Low & McCray 1988):

$$\begin{array}{c}\hfill V_{\exp }(t)\simeq 38.6{\left(\frac{L_{38}}{n_0}\right)}^{1/5}{t}_6^{-2/5}\mathrm{km}{\mathrm{s}}^{-1},\end{array}$$(9)

where L₃₈ = L_mech/10³⁸ erg s⁻¹ is the mechanical luminosity of the stellar association, n₀ is the ambient density, and t₆ is t_kin expressed in Myr. Considering $n_{\text{e}}^{(\text{min})}$ as a proxy of density and $L_{\text{mech}}\simeq E_{\text{mech}}^{\text{SN}}/t_{\text{cls}}^{\text{SN}}$, we predict V_exp for Shells similar to what is observed.

While the velocity dispersion (and hence the derived expansion velocities) of the regions in this study is higher than it is in normal H II regions (see Sect. 4.5), it does not vary significantly between the regions and results in V_exp ∼ 10 − 60 km s⁻¹ (mostly within 20 − 35 km s⁻¹). Thus, the large spread of E_kin for these regions is driven mostly by the differences in their gas mass, not in their kinematics. We also note that the derived kinematic age is typically below 10 Myr, which supports our assumption made in Sect. 3.2 that only energy input from the star clusters during the last 10 Myr is relevant for establishing the energy balance between ionised gas and stars.

4.4. Kinetic energy of the ionised gas and its dependence on metallicity

The estimated kinetic energy E_kin of the ionised gas in the identified regions spans about three orders of magnitude – from ∼5 × 10⁵⁰ erg to ∼5 × 10⁵³ erg. We consider how the measured values of E_kin change with the metallicity of a region in Fig. 12 (left panel). Regions with lower oxygen abundance tend to have lower kinetic energy, though the number of regions with Z < 0.5 Z_⊙ is rather low for drawing definitive conclusions. Note, however, that measurements of energetics of low-metallicity (Z < 0.1 Z_⊙) superbubbles available in the literature (Egorov et al. 2021; Gerasimov et al. 2022) follow the same trend as in our data.

Fig. 12. Distribution of the measured kinetic energy of the ionised gas in the identified regions of locally elevated Hα velocity dispersion (left panel) and the masses of the star clusters in these regions (right panel) shown for different metallicity cuts. These regions tend to have lower kinetic energy in the low-metallicity environment, but also the typical masses of clusters decline. The bin for Z < 0.1 Z⊙ corresponds to the estimates of kinetic energy for 16 ionised gas superbubbles identified in the dwarf galaxies Sextans A (Gerasimov et al. 2022) and DDO 53 (Egorov et al. 2021) and is shown here for comparison.

Such a metallicity trend can be a consequence of either a lower mechanical energy input injected by young stars into the ISM, or higher energy losses, or it can result from differences in the distribution of cluster masses in low- and high-metallicity regions. The right panel of Fig. 12 shows that the latter is a likely explanation of what we see in the PHANGS-MUSE data. Indeed, clusters at Z < 0.5 Z_⊙ tend to be less massive than at higher metallicities, which naturally leads to lower mechanical energy input from these clusters. We do not see any noticeable differences in the efficiency of stellar feedback (considering the ratio of E_kin and cumulative energy input produced by clusters, see Sect. 5.2), thus there is no evidence of differences in energy losses. We cannot, however, completely rule out the probable effect of a declining mechanical energy input per star cluster. A decrease in the mechanical energy input from the star clusters at low-metallicity is present in the results of STARBURST99 simulations in Fig. 6. Theoretical models suggest that the energy contributed by stellar winds (Vink et al. 2001) declines in the low-metallicity regime. The observed number of SNe type II in low metallicity environments is also low (Anderson et al. 2016), however the mass range and fraction of stars ending as SNe do not vary significantly with metallicity (Heger et al. 2003). Note that the observed trend of E_kin with metallicity in Fig. 12 is derived from directly measurable values in the MUSE spectra only, it does not depend on any model assumptions (except the assumption of spherical geometry), and thus provides a purely observational indication for variations in mechanical stellar feedback as a function of metallicity.

4.5. Ionised gas properties compared to normal H _II regions

The regions analysed in the current study are selected by their locally elevated Hα velocity dispersion almost independently of their morphology, so they do not necessarily have the same properties as typical H II regions. In Fig. 13 we examine some of the relations between the Hα luminosity, velocity dispersion and mass of the associated star clusters. We compare these relations in the regions of elevated velocity dispersion with what is found for H II regions. We use a sub-sample of H II regions from the PHANGS-MUSE nebulae catalogue, compiled by Scheuermann et al. (2023), representing regions associated with a single star cluster. We can thus easily establish the link between the ionised gas and young stars for this reference sample.

Fig. 13. Properties of the regions of elevated velocity dispersion in comparison to those of normal H II regions. Blue and orange contours show the statistical density distribution of the H II regions (associated with a single star cluster, from the catalogue compiled by Scheuermann et al. 2023) and the regions with locally elevated σ(Hα), respectively. The Hα luminosity (L(Hα)) correlates with the total mass of the clusters (Mcls) for both samples (left panel), but the slope is different implying differences in the escape fraction of the ionising photons and/or contribution of the additional ionising sources. σ(Hα) for the regions in this study slightly decreases with the Mcls (middle panel) and L(Hα) (right panel) contrary to what is observed in the reference H II regions.

Massive stars in the young stellar associations and compact clusters produce the hydrogen-ionising photons that are responsible for the observed Hα flux from the ionised gas surrounding the clusters. Therefore one naturally expects to see a correlation between the mass of the star clusters and the Hα luminosity of the H II regions. This is however not necessarily true for our regions of elevated velocity dispersion as they are not classical H II regions. Given their kinematically disturbed nature, shock heating could contribute to their observed Hα emission. Also given that these regions are likely superbubbles or turbulent volumes of gas, leakage of ionising photons into and out of the region can also play a role. The left panel of Fig. 13 shows that the Hα luminosity of the regions with high-velocity ionised gas motions still correlates with the total mass of the young stellar associations within the regions. The slope of this correlation is, however, shallower than for bona-fide H II regions. At the same time, these regions span the same range of Hα luminosities as the H II regions.

By construction, the regions considered in this paper have higher Hα velocity dispersion than normal H II regions. The middle and right panels of Fig. 13 demonstrate that the distributions of σ(Hα) for these objects have almost no overlap with the H II region sample. The MUSE spectral resolution is not sufficient to measure the well-known positive correlation between the Hα luminosity or the mass of the ionising star clusters and the Hα velocity dispersion for H II regions (e.g., Terlevich & Melnick 1981; Melnick et al. 1999; Moiseev et al. 2015). However, we see a negative shallow trend for the regions with elevated velocity dispersion – σ(Hα) is, in general, lower for the brighter regions and/or regions associated with more massive star clusters. Therefore, from the distributions shown in Fig. 13 we conclude that the observed properties and, in particular, the correlation between properties of ionised gas and massive stars slightly differ in these dynamically active regions from what is typically observed in H II regions.

5. Discussion

5.1. Velocity dispersion at the location of the young stellar associations

If the dynamics of the ionised gas is strongly affected by stellar feedback, we expect to see the velocity dispersion correlate with the presence of young stars. Using the positions of the young stellar associations (a dominant source of stellar feedback) identified in the HST images, we check if an elevated velocity dispersion towards these stellar associations is observed. Fig. 14 (left-hand panel) shows the distribution of σ(Hα) in the individual pixels residing within the borders of young stellar associations (younger and older than 4 Myr shown by orange and blue colour, respectively) and outside the stellar associations. There is no measurable difference in σ(Hα) between the pixels associated with the young stellar associations older than 4 Myr age and those not overlapping with any stellar associations. Surprisingly, σ(Hα) is systematically lower for stellar associations younger than 4 Myr. Applying a t-test to the bootstrapped sample of 1000 MC realisations of these distributions (accounting for the uncertainties of the measurements of σ(Hα)), we find the separation between the two distributions is 5 − 11 km s⁻¹, and the probability that they are equal is negligible.

Fig. 14. Distribution of the velocity dispersion of Hα (left) and CO (right) emission within the borders of young stellar associations of different ages (orange and blue colours for those younger and older than 4 Myr, respectively) and outside them (green colour). Only pixels with S/N > 30 in the corresponding emission line are considered. Note that increasing the S/N limit does not affect the result.

It is not obvious why pixels corresponding to the locations of the youngest stellar associations show lower Hα velocity dispersion. If this results from lower turbulence towards the young stellar associations, then one would expect to see the same trend in the cold gas too. However, a similar comparison with the PHANGS-ALMA data (Leroy et al. 2021), using the CO line kinematics for the same sample of galaxies (Fig. 14, right-hand panel), shows no obvious difference¹³. The elevated σ(Hα) towards stellar associations older than 4 Myr (and the regions without any young stars) could be related to an increased contribution of the DIG to the total Hα emission. Namely, we suggest that the Hα emission around all young stellar associations is tracing the kinematics of the thin star-forming gas disc, but the diffuse Hα outside the youngest of them is more impacted by turbulence and outflows. A higher mean σ(Hα) relative to H II regions is often reported for the DIG (e.g., Moiseev 2015; Levy et al. 2019; Della Bruna et al. 2020; Law et al. 2022), although it remains an open question whether stellar feedback or large-scale processes are responsible for this. The same effect, in principle, could be a consequence of superbubbles beeing powered by off-centre stellar associations (see Sect. 5.3), and thus a large kinematic separation between the approaching and receding walls is observed outside the young stellar associations. This implies that the influence of stellar feedback on the ionised gas kinematics can be observed better at some separation from the source of feedback.

5.2. The driving source of turbulent motions in star-forming galaxies

As discussed in Sect. 1, several processes have been proposed as drivers of high velocity dispersion: stellar feedback, gravitational instabilities, and external processes. In this paper we identified 1484 regions with locally enhanced Hα velocity dispersion, and linked 511 of them to young stellar associations or compact star clusters. Comparing the measured properties of the ionised gas with the modelled energy input from the stellar population (see Sect. 3.3), we can determine which of these mechanisms (winds and SNe or gravity) is driving the observed supersonic gas motions.

We compare the measured kinetic energy E_kin of the ionised gas with the total mass of stars in the young stellar associations within the regions (Fig. 15, left panels) and the surface density of the kinetic energy ΣE_kin with the mean total stellar mass surface density Σ_⋆ (right-hand panels). The latter is taken from the PHANGS-MUSE DAP maps, and measured from the local stellar continuum spectrum (see Sect. 2.1) and averaged over the area corresponding to each region. The measured ΣE_kin is largely independent of Σ_⋆, implying that gravitation alone¹⁴ is unlikely to be the dominant powering source of the turbulence. In contrast, E_kin strongly correlates with the mass of the young stellar associations and clusters, pointing to stellar feedback playing a critical role. We emphasise, however, that this is valid for discs of galaxies as bars and centres are excluded from our analysis. The result is independent of the adopted method for the n_e measurements except that $n_{\mathrm{e}}^{(\max )}$ leads to a more linear trend.

Fig. 15. Dependence of the kinetic energy of ionised gas Ekin in the regions of locally elevated velocity dispersion on the mass of the stellar associations within them (left panels), and their gas kinetic energy per area, ΣEkin, vs. corresponding stellar mass surface density, Σ⋆ (right panels). Top and bottom panels correspond to minimal and maximal values of ne underlying the Ekin measurements. Colour encodes the Hα flux surface brightness. The red line shows the median values of Ekin or ΣEkin in several bins, and the red shaded area corresponds to the 25–75 percentile interval. Spearman correlation coefficients ρ are given above the plots.

In Fig. 16 we compare the kinetic energy of the ionised gas with the total mechanical energy input produced by massive stars during the lifetime of the association (or the last 10 Myr if the star cluster is older than that) in the form of winds and supernovae (left panels), and supernovae alone (right panels). One can see that in most regions, the cumulative mechanical energy input from young stars is sufficient to explain the observed kinetic energy of the ionised gas. There are, however, ∼16% of the regions residing in the top-left panel above the 1-to-1 line (∼4% if counting only those farther than the quoted uncertainties) – in these cases young star clusters and associations do not provide sufficient mechanical energy input to explain the measured value of E_kin. 64% of them have [S II]/Hα > 0.4, indicative of shock-heated gas due to possible recent explosion of SNe, so we expect to underestimate the number of star clusters for these regions. Another explanation for these outliers is that we do not account for the pressure-driven expansion due to the presence of hot gas, which might be significant for these small regions (Barnes et al. 2021). Finally, we note that our measurements of E_kin are for the lower density limit $n_{\mathrm{e}}^{(\min )}$ corresponding to an upper limit of E_kin. All regions reside below the 1-to-1 line (within the errors) if we assume an upper limit estimate of $n_{\mathrm{e}}=n_{\mathrm{e}}^{(\max )}$ instead (bottom-left panel).

Fig. 16. Dependence of the kinetic energy of ionised gas Ekin in the regions of locally elevated velocity dispersion on the total mechanical energy input from the stellar associations in the form of supernovae and stellar winds (left panels), and supernovae alone (right panels). Top and bottom panels correspond to the $n_{\mathrm{e}}^{(min)}$ and $n_{\mathrm{e}}^{(max)}$ density measurements (thus representing upper and lower limits of Ekin, respectively). Median value of logarithmic errors is shown in bottom-right corner of each panel. All regions with errors 2 times larger than presented are excluded from the plots. Colour encodes the [S II]/Hα lines ratio tracing the dominant gas excitation mechanism ([S II]/Hα > 0.4 is likely produced by shocks; left panels), or the relative contribution of the SNe to the total mechanical energy input according to the STARBURST99 models for the corresponding mass and age of each cluster right panels. Blue, orange, green and red lines correspond to mechanical stellar feedback energy efficiency η = 100, 30, 5, 1%, respectively. Spearman correlation coefficient ρ is given above the plots.

The situation for the energy balance changes if we account only for the SNe contribution (right panels of Fig. 16). In that case, a third of the considered regions reside above the one-to-one relation assuming $n_{\mathrm{e}}^{(\min )}$ (13% if counting only those with higher E_kin than the quoted uncertainties) – the star clusters do not produce enough mechanical energy via SNe alone to explain the observed kinetic energy of the ionised gas in these regions. The colours encode the relative contribution of supernovae to the total mechanical energy input from the star cluster during its lifetime, calculated based on STARBURST99 models given the measured ages and masses. The mechanical energy input is mostly dominated by stellar winds for the outlier regions. This means that the contribution of pre-supernovae mechanical feedback (in the form of stellar winds) must be accounted for in order to understand what is driving the small-scale turbulence in the ISM. We note, however, that almost all regions reside below the 1-to-1 line if considering upper limit of the electron density ($n_{\mathrm{e}}^{(\max )}$), but the correlation between E_kin of gas and E_mech from SNe is much weaker in both cases than that for the combined contribution of SNe and winds.

Significant contribution of pre-SNe feedback in total energy balance does not necessarily mean that it is also critical for reproducing the gas kinematics in simulations. As we demonstrated in Sect. 4.3, our estimates of the kinematic ages for our regions better agree with what is obtained from the SED fitting of star clusters when we consider only time passed since the first SN explosion. We suggest that SNe override the previous impacts of the star clusters onto the gas kinematics and that the latter is mostly regulated by SNe feedback (in agreement with simulations, e.g., Walch & Naab 2015). However, pre-SNe feedback disperses molecular clouds ionising them and blowing-out gas away from the clusters (Chevance et al. 2022), and therefore the SNe explode in the pre-cleared medium. Also, the continuous impact of pre-SNe feedback during a few Myr can produce bigger superbubbles with more gas in their rims, which explains large differences in resulting E_kin in Fig. 16, in comparison with SNe contribution only. Therefore, accounting in simulations for the impact of pre-SNe feedback on the ISM can be essential for reproducing the realistic gas morphology and variation of SFR in the galaxies (in agreement with simulations by, e.g., Haid et al. 2018; Keller et al. 2022).

In general, there must be an energy balance between the young stellar associations and the turbulent ionised ISM. Indeed, in Fig. 15 (left panel) the Hα surface brightness (colour-coded) correlates with the stellar mass in the associations (see also Fig. 13). Given that the total Hα flux is also a parameter in the calculations of the ionised gas mass (see Eq. (4)), the observed correlations in Figs. 15 and 16 are likely driven by the fact that more massive associations produce more ionised gas around them. As a consequence of this, the total Hα luminosity (and thus the observed Hα flux) increases as the stellar association mass increases. This is often predicted by photoionisation models (e.g., in STARBURST99 among many others), and this correlation is also observed on the scale of individual H II regions (Scheuermann et al. 2023). Meanwhile, our results imply that this correlation is preserved also for the turbulent ionised gas around H II regions and for superbubbles, whose properties are not necessarily the same as those of normal H II regions (see Sect. 4.5). Given that the mass of the young stars and of the ionised gas (as well as E_mech from the stars and E_kin of the gas, respectively) are measured from unrelated data in a way independent of each other, we are confident that the correlation between these parameters is real and not a statistical effect. Thus, we can use our measurements to estimate observationally what fraction of the mechanical energy produced by stars is converted into gas kinetic energy in order to support the supersonic gas motions: η = E_kin(gas)/E_mech(stars). Our results imply that η ∼ 10 − 20% (see Fig. 17), although there are regions with η > 100% where the measured mechanical feedback alone cannot explain the gas kinetic energy (assuming the lower limit of electron density as discussed above). Note that in our measurement of E_kin we consider the energy excess in each region with respect to that in the less perturbed ISM, so the energy coupling should be closer to ∼10% if we also account for the surrounding gas. This is in good agreement with simulations by Ejdetjärn et al. (2022), who found that ∼10% of the energy from stellar feedback is necessary to drive the supersonic velocity dispersion of the ionised gas. We do not find any obvious change in the energy coupling with metallicity, which suggests that the lower measured kinetic energy of the ionised gas in the low-metallicity turbulent regions (see Sect. 4.4 and Fig. 12) results from a change in the cumulative mechanical energy input by stellar feedback and not the energy losses. However, the number of regions with metallicity Z < 0.5 Z_⊙ is very low, and more data are required to confirm this finding.

Fig. 17. Distribution of the estimated mechanical stellar feedback energy efficiency for the regions shown in top-left panel of Fig. 16 grouped by metallicity.

5.3. Ionised superbubbles: Energy balance and possible signs of star formation triggering

As mentioned in Sect. 3.1, we identified a peak in the Hα flux radial distribution corresponding to a swept-up shell for 171 of our regions. These regions are our best candidates for expanding superbubbles, and it is interesting to compare their properties with the remaining regions.

In Fig. 18 we show the same comparison of the kinetic energy of ionised gas and the mechanical energy input from stellar associations as in Fig. 16, but in this figure, we include only those regions with shell-like morphologies. They follow the same trend as the rest of the regions. We measure similar mechanical energy efficiencies to those found for the whole sample – η ∼ 10 − 20%. This value is consistent with predictions from numerical simulations, where Sharma et al. (2014) followed the evolution of superbubbles powered by multiple (from tens to a few hundred) OB stars and obtained an asymptotic value of η ∼ 20 − 25%. We note that other simulations predict lower values – for example, Yadav et al. (2017) give η ∼ 6 − 10%, which is also close to the peak in Fig. 18. Our conclusion is valid for our adopted estimate of the electron density $n_{\mathrm{e}}^{(min)}$ generally corresponding to the upper limit of measured η. Using the upper limit on electron density $n_{\mathrm{e}}^{(max)}$ results in a much lower value of η ∼ 1%. We would need larger statistics and more precise measurements of the electron densities to obtain strict observational limits on the parameter η, or to identify the main source of its variations. We note also that our measurements are related to ionised gas only. Based on observations of the same sample of galaxies as in this work, Watkins et al. (2023b) found a coupling efficiency of ∼5 − 12% required to put into agreement the derived and predicted stellar masses and ages inside the molecular gas superbubbles, that is close to our estimate.

Fig. 18. Comparison of the kinetic energy of ionised gas and the mechanical energy input by young stars for our best candidate superbubbles (as exhibiting shell-like morphology). Top plot: same comparison of energies as in Fig. 16, superbubbles follow the same trend as the other regions. Colour encodes the age of the youngest association in the superbubble. Bottom plot: distribution of the estimated efficiency of the energy consumption in superbubbles. Its peak is slightly below η = 0.25 that was obtained from the modelling by Sharma et al. (2014) as an asymptotic value for the expanding superbubbles driven by multiple (up to a few hundred) OB stars in normal star-forming galaxies.

For the sample of 171 likely superbubbles, we can resolve their morphology and thus check if the stellar associations driving their expansion are located close to their centres. For that, we calculated the separation as the angular distance between the centres of the regions and the position of the stellar associations (the geometric centre of the corresponding stellar association mask), deprojected them based on the position angle and axes ratio of the regions and translated to physical distance normalised to R_eff of each region. For separations less than the PSF, we assume a separation equal to zero. We excluded from our analysis any region where R_eff < 3 × PSF to ensure the reliability of our results. As a result, we are left with 40 superbubble-like regions where stellar associations or star clusters younger than 40 Myr¹⁵ are observed. In these regions, about 40% of the star clusters are offset from the centre of the superbubble by at least 0.5R_eff. This implies that the stellar sources powering the expanding superbubbles are often located in their rims (see also Egorov et al. 2017; Gerasimov et al. 2022; Barnes et al. 2023). We speculate that in such cases we observe the superbubbles (or outflows) driven by young stars at the edges of molecular clouds (‘blister’-like regions). If such offsets are also common in other objects from our entire sample, this can partially explain the difference between the Hα velocity dispersion towards the youngest stellar associations and the regions without any young stars (see Sect. 4 and Fig. 14). Namely, more dynamically perturbed gas is observed in the lower density environment outside the powering stellar associations, which reside in the denser clouds.

We also find a weak trend with the age of the stellar association, shown in Fig. 19: fraction of the associations younger than 3 Myr is higher (by ∼15%) at the periphery of the regions (at R > 0.5R_eff) than close to their centres (R < 0.5R_eff). A higher fraction of the younger associations on the edge of superbubbles points towards sequential star formation (Elmegreen & Lada 1977), or triggering of star formation as a result of the evolution of the superbubble (e.g., Hartmann et al. 2001), though it is usually hard to distinguish between these processes (Dale et al. 2015). The triggering of star formation in the periphery of expanding supershells as a result of their interaction with the surrounding ISM, or collision with each other, is predicted in simulations (e.g., Ntormousi et al. 2011; Inutsuka et al. 2015; Vasiliev et al. 2022), and indications of this have been observed in very nearby dwarf galaxies (below 4 Mpc, see, e.g., Lozinskaya 2002; Egorov et al. 2017, 2018; Fujii et al. 2021). Dawson et al. (2013) demonstrated that about 12 − 25% of the molecular gas in the supergiant shells in the LMC was formed as a result of their interaction with the ISM. Our identification of systematically larger separations from the centre of superbubbles for the youngest stellar associations extends the observational signatures of triggered star formation in superbubbles to more massive spiral galaxies. However, the statistics are limited, in principle by the limited angular resolution of the underlying data. Investigation of superbubbles in the delivered and scheduled data for nearby galaxies from JWST imaging will allow for a more detailed investigation into the age distribution of star clusters and stellar associations, in relation to morphologically identified superbubbles (see, e.g., Barnes et al. 2023; Watkins et al. 2023a, as examples of the first results).

Fig. 19. Relative fraction of the stellar associations younger than 3 Myr at different separations between clusters/associations and centres of the regions. Only a sub-sample of 40 well-resolved superbubbles and star clusters younger than 40 Myr are considered. Each separation bin corresponds to 0.3Reff. Errors are calculated from 1000 Monte Carlo realisations with positions of each cluster randomly varying within the PSF of the MUSE observations.

6. Summary

Here we quantify the contribution of mechanical feedback from young massive stars to the regulation of the kinematics and morphology of the ionised ISM. Namely, we identify areas of locally elevated Hα velocity dispersion (likely a result of the presence of unresolved expanding superbubbles, SNRs, or turbulent ionised gas motions), connect them to young stellar associations, and quantify the energy balance between the mechanical energy in the ionised gas and what is contributed by the observed stellar populations in these regions. Our analysis combines observational data from MUSE/VLT IFS with HST multi-wavelength imaging for 19 nearby spiral galaxies obtained as a part of the PHANGS survey. Our main results are summarised below:

• We identify 1484 regions of locally elevated Hα velocity dispersion (σ(Hα) > 45 km s⁻¹), with tens in each galaxy. About half of them contain at least one young stellar association or compact star cluster within their borders. We argue that most of those regions not connected with star clusters are SNRs – they typically have a centrally concentrated Hα flux distribution and emission line flux ratios consistent with shock excitation.

• We find that the kinetic energy of the ionised gas in the regions of locally elevated Hα velocity dispersion strongly correlates with the mass of the young stellar associations within these regions and with their total mechanical energy input into the surrounding ISM (in the form of stellar winds and SNe based on STARBURST99 models). At the same time, the correlation with mechanical energy produced by young stars is significantly weaker if we only consider the impact of SNe. This implies that pre-SNe feedback is an essential contributor to the energy balance.

• About ∼10 − 20% of the total mechanical energy input is retained in the regions in the form of kinetic turbulent energy in the ionised gas. This result is in agreement with the coupling efficiencies obtained from numerical simulations for feedback-driven superbubbles (e.g., as examples of the first results Sharma et al. 2014) or ISM turbulence (e.g., Ejdetjärn et al. 2022).

• The measured kinetic energy of the ionised gas changes with metallicity, which can be a consequence of both the lower masses of the clusters and weaker stellar feedback in the low metallicity environments.

• Among the high velocity dispersion regions, only 171 (12%) exhibit a clearly distinguishable shell-like morphology in their Hα flux distribution. These regions are probably expanding ionised superbubbles. In general, the size distribution of the identified regions is well described by a power law with a slope ∼ − 2.7, in agreement with observations (e.g., Oey & Clarke 1997; Bagetakos et al. 2011; Watkins et al. 2023a) and simulations (Nath et al. 2020) of superbubbles in star-forming galaxies. Thus, we cannot exclude that most of the identified regions of high velocity dispersion appear to be expanding superbubbles, in many cases the shell-like morphology is not resolved or overlaps with other structures.

• The estimated expansion velocity V_exp = 10 − 60 km s⁻¹ (mostly within 20 − 35 km s⁻¹) for the superbubbles in our sample and corresponding kinematic ages are typically less than 10 Myr (mostly within 1 − 6 Myr). The derived ages agree better with the time passed since the first SN explosion than with the total lifetime of the star clusters. This implies that once SNe occur, their impact becomes dominant in the regulation of the current ionised gas kinematics in the analysed superbubbles.

• About 40% of the well-resolved superbubbles in our sample have young stellar associations located on their rims, not in their centres. The youngest associations are more often observed at the periphery of the superbubbles, while older stellar associations are mostly observed closer to their centres. This is an indication of sequential or triggered star formation occurring at the edges of superbubbles in our sample.

Overall, we establish a direct observational connection between the mechanical feedback from young massive stars onto the surrounding ionised gas for a large and homogeneous sample of regions, including expanding superbubbles. We demonstrate that accounting for pre-SNe mechanical feedback in form of stellar winds is important, as SNe alone fail to explain the measured kinetic energy of the ionised gas.

We also show that the mechanical stellar feedback energy efficiency does not change with metallicity. However, our data are very limited for metallicities Z < 0.5 Z_⊙. A further extension of this analysis to dwarf galaxies will allow us to establish a purely observational constraint on the efficiency of mechanical stellar feedback in various environments, which is a critical ingredient in cosmological simulations.

Acknowledgments

This work was carried out as part of the PHANGS collaboration. Based on observations collected at the European Southern Observatory under ESO programmes 094.C-0623 (PI: Kreckel), 095.C-0473, 098.C-0484 (PI: Blanc), 1100.B-0651 (PHANGS-MUSE; PI: Schinnerer), as well as 094.B-0321 (MAGNUM; PI: Marconi), 099.B-0242, 0100.B-0116, 098.B-0551 (MAD; PI: Carollo) and 097.B-0640 (TIMER; PI: Gadotti). This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. These observations are associated with programme 15654. This research made use of astrodendro, a Python package to compute dendrograms of Astronomical data (http://www.dendrograms.org/). KK, OVE, EJW, JEM-D gratefully acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in the form of an Emmy Noether Research Group (grant number KR4598/2-1, PI Kreckel). RSK and SCOG acknowledge funding from the European Research Council via the ERC Synergy Grant “ECOGAL” (project ID 855130), from the Heidelberg Cluster of Excellence (EXC 2181 - 390900948) “STRUCTURES”, funded by the German Excellence Strategy, and from the German Ministry for Economic Affairs and Climate Action in project “MAINN” (funding ID 50OO2206). RSK and SCOG also acknowledge computing resources provided by The Länd through bwHPC and DFG through grant INST 35/1134-1 FUGG and for data storage at SDS@hd through grant INST 35/1314-1 FUGG.

References

Appendix A: I − σ diagrams for all PHANGS-MUSE galaxies

The regions analysed in this paper were selected based on the I − σ diagrams introduced and described in Section 3.1. Here we describe the details of the performed analysis and how the selection criteria were determined. The I − σ diagrams and classification maps for all PHANGS-MUSE galaxies are shown in Fig. A.1 (NGC 4254 is shown in Fig. 4).

Fig. A.1. I − σ diagram (right), classification map (centre) and the Hα map (left) for all PHANGS-MUSE galaxies except NGC 4254 (shown in Fig. 4). Cyan ellipses encircle the identified regions of the locally elevated Hα velocity dispersion.

Fig. A.1. Continued.

Fig. A.1. Continued.

Fig. A.1. Continued.

Fig. A.1. Continued.

Fig. A.1. Continued.

A large scatter in σ(Hα) is usually observed at the low surface brightness end of the I − σ diagrams (e.g. Moiseev & Lozinskaya 2012), and this is also the case for most of the PHANGS-MUSE galaxies. This scatter can result from both physical (e.g. the impact of shocks) and observational (low signal-to-noise ratio) reasons, though we remind that we consider only the pixels with S/N > 30. Here we exclude σ(Hα) peaks associated with the low surface brightness DIG. We tested two quantitative criteria for isolating the low surface brightness DIG. Both criteria suggest a constant level of surface brightness, Σ(Hα)_DIG, is suitable for each galaxy in order to separate the DIG from the rest of the emission structures. First, we adopted a median surface brightness for the DIG, Σ(Hα)_DIG, as estimated for the PHANGS-MUSE galaxies in Belfiore et al. (2022). In that paper, the DIG was isolated from H II regions by a morphological analysis. Another criterion is based on the spectral properties of the DIG, as it exhibits elevated relative intensities in low-ionisation ions compared to H II regions, and usually resides above the maximum starburst line from Kewley et al. (2001) on the [O III]/Hβ vs [S II]/Hα diagnostic (e.g. Haffner et al. 2009; Belfiore et al. 2022). We derived Σ(Hα)_DIG as the value corresponding to the 75th percentile of all pixels residing above the maximum starburst line. For half of the galaxies in our sample, both methods gave very similar values of Σ(Hα)_DIG. However, visually those taken from Belfiore et al. (2022) are better at separating the low brightness peak of σ(Hα) on I − σ diagrams for the rest of the objects. Therefore, we used the values from Belfiore et al. (2022) to determine Σ(Hα)_DIG. For five galaxies (NGC 1087, NGC 1300, NGC 1385, NGC 3351, and NGC 7496) we additionally increased (by 0.1 − 0.5 dex) the value of Σ(Hα)_DIG to better match the low-brightness peak in the σ(Hα) distribution. According to Belfiore et al. (2022), the DIG emission in three of these galaxies has a very wide and shallow surface brightness distribution. The final adopted values of Σ(Hα)_DIG for each galaxy are given in Table 1.

We estimated the mean Hα velocity dispersion σ(Hα)_m0 for the H II regions in each galaxy as the flux-weighted mean value across all pixels brighter than Σ(Hα)_DIG. The uncertainty of the mean velocity dispersion Δ(σ_m)₀ was measured as the standard deviation of σ(Hα) around σ(Hα)_m0. The final adopted values of σ(Hα)_m and Δ(σ_m) were obtained on the second iteration, where we excluded all pixels with σ(Hα) > σ(Hα)_m0 + Δ(σ_m)₀. The derived values of σ(Hα)_m characterise the mean velocity dispersion of the ionised ISM in each galaxy and are given in Table 1. These values agree well with the empirical relations between σ(Hα)_m and SFR as reported in other works (e.g. Law et al. 2022).

Based on σ(Hα)_m and Δ(σ_m), we performed a classification of all pixels brighter than Σ(Hα)_DIG into four classes (shown by different colours in Figs. 4 and A.1) in a way similar to that in Egorov et al. (2021). The particular criteria and physical meanings are summarised in Table A.1 and described in detail below.

Class 1: Unperturbed σ(Hα). Normal H II regions tend to occupy the elongated area of relatively high surface brightness and low velocity dispersion on I − σ diagrams. For brighter regions, a scatter of σ(Hα) measurements is usually smaller than for faint regions, which can also be seen for most galaxies in Fig. 4. The upper boundary of the distribution of σ(Hα) vs log(I[Hα]) for such regions can be parameterised by a Gaussian with the mean value and standard deviation equal to σ(Hα)_m and Δ(σ_m), respectively. We refer to velocity dispersion defined by such a normal distribution for each flux value as σ(Hα)_norm.

We classify a pixel as related to Class 1 if σ(Hα) < σ(Hα)_norm for its Hα brightness. Note that all pixels with σ(Hα) < σ(Hα)_m are also considered as belonging to this class. This class contains all bright H II regions and surrounding ionised ISM that do not exhibit signs of perturbation in their kinematics.

Class 2: Intermediate σ(Hα). This includes all pixels with values of σ(Hα) between two distributions – σ(Hα)_norm and σ(Hα)_{1.5 × norm}. The latter is defined in the same way as σ(Hα)_norm, but with a standard deviation equal to 1.5Δ(σ_m). This is an intermediate class: the elevated σ(Hα) may still be illusive (due to the contamination by pixels in the wings of the distribution defining Class 1 objects), but can reflect real physical changes in the velocity dispersion due to turbulence, shocks, the presence of slowly expanding superbubbles, or is observed at the edges of the more dynamically active Class 3 regions.

Class 3: Elevated σ(Hα). This class is assigned to pixels with σ(Hα) > σ(Hα)_{1.5 × norm} (thus the contamination by wrongly classified ‘unperturbed’ pixels should not exceed 7% assuming their normal distribution). The corresponding regions clearly show remarkable local peaks of velocity dispersion, that can be driven by turbulence, shocks or moderately expanding superbubbles (due to the presence of emission from both approaching and receding sides). In agreement with the toy models in Moiseev & Lozinskaya (2012), Muñoz-Tuñon et al. (1996), and with the simulations in Vasiliev et al. (2015), the regions corresponding to this class usually show triangular shapes on the I − σ diagrams.

Top σ(Hα) and SNR/WR candidates. This is a sub-group of the Class 3 pixels selected as having intrinsic σ(Hα) > 100 km s⁻¹, or σ(Hα) > σ(Hα)_{3 × norm} (where σ(Hα)_{3 × norm} is the normally distributed value with a standard deviation equal to 3Δ(σ_m)). The cores with the most extreme peaks in velocity dispersion are selected in this subclass. As was demonstrated by Moiseev & Lozinskaya (2012), these pixels might correspond to the presence of energetic stellar objects like SNRs, WR, LBV etc. Those pixels that pass the criteria for Classes 2 or 3 and have log(I[Hα]) above that for any Class 1 pixel are also considered as part of this subclass, as they are definitely outliers from the normal distribution of σ(Hα).

The regions of correlated high velocity dispersion are relatively easy to identify by eye in the constructed classification maps (central panel of Figs. 4 and A.1). To automatically select these regions and derive their boundaries, we used the ASTRODENDRO¹⁶ package. However, the application of ASTRODENDRO directly to the classification map resulted in a very rough selection of the regions (many of the individual Class 3 regions were merged into a single region), and only the strongest peaks of velocity dispersion were often identified when we applied ASTRODENDRO to the σ(Hα) map directly. To overcome these issues, we created a synthetic map based on both the classification and velocity dispersion distribution by clipping the velocity dispersion range for each of the classes (as summarised in Table A.1). Namely, for each galaxy, we derived the median σ(Hα) of all pixels within the galaxy classified as Class 2 or above (hereinafter σ_ref), and set σ(Hα) = σ_ref − Δ(σ_m) for all Class 1 and DIG regions, and σ_ref for Class 2 regions. For Class 3 and SNR/WR candidates, we left their observed σ(Hα), but with the lower limit equal to σ_ref + 2Δ(σ_m) and σ_ref + 3Δ(σ_m), respectively (approximately equal to the mean thresholds used for the classification in the previous step).

Applying ASTRODENDRO to the constructed synthetic maps, we required the minimal area of adjoined pixels to be no less than π(PSF/2)² (where the full width at half maximum (FWHM) size of the PSF – point spread function – characterises the angular resolution of the MUSE data and is given in Table 1), and the minimal deviation of the identified region from the background should be 0.9Δ(σ_m). Each of the identified regions was then approximated by an ellipse having minor and major diameters equal to 2.65 × FWHM_A, where ${\mathrm{FWHM}}_A=2\sqrt{2ln2}{\sigma }_A$, and σ_A is the standard deviation along the minor and major axes of the 2D distribution resulting from the statistics computed by ASTRODENDRO. The ellipses defined in such a way correspond to the baselines of the σ(Hα) distribution in the region and encircle the entire area of the elevated velocity dispersion. In the case of a bubble, this corresponds to its rim. To prove this, we used a toy model of a thin sphere expanding in a homogeneous medium. The velocity separation between the components from the approaching and receding sides of such a simulated bubble is $\mathrm{\Delta }V\propto \sqrt{1-r^2/R^2}$, where r is the distance in the sky plane from the centre of the bubble with radius R. We analysed the radial distribution of the velocity dispersion measured from the simulated spectra integrated along the line-of-sight and found that the weighted standard deviation of the σ(Hα) (calculated in the same way as by ASTRODENDRO) σ_A ≃ 0.32 × R for various expansion velocities and spatial resolutions¹⁷. We excluded all regions having an axis ratio less than 0.35 (chosen based on visual inspection of the results obtained with the different parameters). Even if real, such elongated regions may add significant scatter to the further analysis because of the uncertainties of their size measurements.

Where possible, we refined the sizes of the elliptical regions identified by ASTRODENDRO. We considered the oversampled radial Hα flux distribution for each region. For those regions representing the expanding superbubbles we expect to see a peak in the Hα flux radial distribution associated with the rims of the swept-up shell, and thus we define the radial distance to this peak as the radius of the superbubble (see example in Fig. 7). We were able to clearly identify peaks in the radial flux distribution for 171 regions, which we use as our best candidates for expanding ionised superbubbles (see Section 5.3). For most of these regions, the change in their size after this refinement was not significant (within 15% for half of the regions, and > 30% for only 20 regions). On the other hand, we also identified many centrally concentrated regions, that could be SNR or other objects with unresolved morphology (see Section 4.1). In such a case, the analysis of the classification maps described above can significantly overestimate the size of the regions, and their real radius is better approximated by the standard deviation of the Hα radial flux distribution (Barnes et al. 2022). We left the size of a region unchanged in the cases when we found neither a central peak nor an associated shell in the radial Hα flux distribution.

Finally, we checked for overlapping ellipses and excluded the smaller regions when the overlap was more than 20%. As a result, we identified 1484 regions across all 19 galaxies. The cyan ellipses in Figs. 4 and A.1 show the locations of the identified regions of locally elevated velocity dispersion, and their edges correspond to the adopted sizes for each region.

The number of regions with ‘normal’ velocity dispersion erroneously selected with our criteria is negligible. This is estimated by simulating the synthetic maps of velocity dispersion randomly distributed around σ(Hα) with standard deviation equal to Δ(σ_m) (changing with the Hα flux in the same way as observed, see Fig. 4) and applying the same criteria to these data. For each of the 100 MC realisations, the number of erroneously selected regions per galaxy does not exceed 3, and the size of such regions is close to the PSF.

All Tables

All Figures

Fig. 1. Statistics of the intrinsic Hα velocity dispersion in PHANGS-MUSE galaxies (sorted by their stellar mass, see Table 1) for different environments (the blue colour is for the disc, and orange is for the central regions and bars, according to the classification from Querejeta et al. 2021). Pixels with a S/N < 30 in Hα are excluded. The blue dashed line corresponds to the mean MUSE spectral resolution (velocity dispersion) at the wavelength of Hα. The values of σ(Hα) were corrected for instrumental broadening by subtracting, in quadrature, the MUSE instrumental velocity dispersion for each galaxy.

In the text

Fig. 2. Relation between the recovered σ(Hα) and the true σ(Hα) depending on the S/N. Panel a: relative errors for various S/N. Coloured areas encompass the 5th–95th percentile interval and dotted lines trace the median values. Panel b: probability density of absolute errors for S/N = 30. The adopted MUSE resolution at the wavelength of Hα is FWHM = 2.54 ± 0.1 Å (according to Bacon et al. 2017, but with twice larger scatter; vertical dashed line on panel a).

In the text

Fig. 3. istribution of the Hα surface brightness (panel a) and velocity dispersion (panel b) for NGC 4254 shown as an example. Contours on these panels are the same and correspond to the lines of a constant Hα brightness. Panel c demonstrates a zoom-in Hα map of the area shown by dashed rectangle in panels a, b with the same brightness limits as in panel a. As follows from panel b, velocity dispersion tends to be higher in the low surface brightness regions (dominated by DIG), but there are also areas of elevated σ(Hα) towards the bright regions. Panel d shows the examples of the Hα and [N II] λ6548,  6584 Å line profiles extracted from three regions with low (green line) and high (orange and magenta lines) Hα velocity dispersion. These line profiles corresponding to high σ(Hα) are either broadened (magenta line) or can be decomposed by two components (orange line). Exact areas of the spectra extraction (3 × 3 pixels) are shown on panel c by squares of the same colour as the spectra.

In the text

Fig. 4. Localisation of the regions of locally elevated Hα velocity dispersion (cyan ellipses) in NGC 4254 galaxy, identified based on the ‘intensity – velocity dispersion’ (I − σ) diagnostics, overlaid on the Hα surface brightness map (left panel) and the classification map (central panel; see text). The I − σ diagram is shown in the right panel. The black solid line shows the mean value of σ(Hα) in the galaxy, and the dashed grey lines show its 1σ uncertainty. See Fig. A.1 for the rest of the PHANGS-MUSE galaxies.

In the text

Fig. 5. Mass of the ionised gas in the selected regions assuming different prescriptions for electron density: lower limit ($n_{\mathrm{e}}^{(\min )}$) and upper limit ($n_{\mathrm{e}}^{(\max )}$) derived from the reference H II regions using Eq. (6) and from R[S II] = [S II]λ6717Å/[S II]λ6731Å flux ratio, respectively, and ne, SII derived from R[S II] of the observed spectra of the region itself. Different morphological types introduced in Sect. 4.1 are shown by different colours. Diagonal lines show 1-to-1 relation. Spearman correlation coefficients ρ are given on the plots.

In the text

Fig. 6. Evolution of the cumulative mechanical energy input to the ISM (Emech) normalised by the mass of star cluster (Mcls) for different metallicity (traced by different colours) according to STARBURST99 (Leitherer et al. 1999, 2014) models. The contribution produced by stellar winds alone is shown by the solid line, while the dashed line corresponds to the impact of both stellar winds and supernovae.

In the text

Fig. 7. Examples of regions classified as shells (top row), ambiguous (middle row) and compact (bottom row). Images from left to right correspond to a three-colour multiband image of the region (Hα image from HST (PI: R. Chandar, GO-10402) is in red, FUV from AstroSat (Hassani et al., in prep.) in green, CO from ALMA (Leroy et al. 2021, in blue), the Hα flux distribution, the [S II]/Hα line ratio, and the Hα velocity dispersion. Cyan ellipses show the adopted borders of the region and overlaid contours correspond to the stellar associations’ masks. The HST stellar associations shown in yellow are selected as associated with the region. The right-hand plots show the radial Hα flux distribution relative to the centre of each region (marked by crosses).

In the text

Fig. 8. Relative distribution of the number of regions with different morphological types, ionisation states and associated young stellar associations. Three groups are shown, corresponding to the three morphological types (see text in Sect. 4.1), and their sizes correlate with the relative number of corresponding regions. Different colours correspond to select properties: red and orange colours correspond to regions with [S II]/Hα > 0.4, while blue and green – to those with [S II]/Hα < 0.4; regions that have at least one stellar association within their borders are shown by blue and orange colours, while green and red correspond to those without identified young stars). Violet colour corresponds to the regions outside the HST footprint (and thus without information about star clusters there). Grey lines mark each 10% level.

In the text

Fig. 9. Diagnostic BPT diagrams showing the optical line ratios for all identified regions of locally elevated velocity dispersion (colour-coded circles). Colour denotes either Hα velocity dispersion (top panels), number of young stellar associations and star clusters within the region (middle panels) or estimated kinetic energy of the gas (bottom panels). For reference, we show as grey circles nebulae from the PHANGS-MUSE catalogue (including H II regions, supernovae remnants, planetary nebulae; see Groves et al. 2023) excluding those falling in the masked areas (see Sect. 3.1). Demarcation lines separate areas on the diagrams typically occupied by regions with different excitation mechanisms: the black solid curve from Kewley et al. (2001) separates H II regions from those with a dominant contribution from other excitation mechanisms (e.g., DIG, shocks); composite mechanism of excitation is probable in the area between this and the grey curve from Kauffmann et al. (2003); AGNs usually occupy the area above the black straight line from Kewley et al. (2006).

In the text

Fig. 10. Distribution of sizes for regions of high Hα velocity dispersion. The red line shows a result of a power-law fit to the data for Reff = 70 − 300 pc. This fit agrees with a power-law approximation for the sizes of superbubbles according to hydrodynamical simulations (α = −2.7 obtained by Nath et al. 2020). The distribution of larger and smaller regions is better approximated by log-normal distribution (the blue line). The dashed black line and the grey area show the median spatial resolution and its standard deviation across our sample.

In the text

Fig. 11. Kinematic properties of the regions with elevated velocity dispersion. Left panel shows the relation between the estimated expansion velocity (Vexp) and kinematic age of each region (grey points). Blue points correspond to the regions classified as ‘Shells’, and light-blue crosses – to the values of Vexp for that sub-sample predicted from the analytical model (Eq. (9)) based on their estimated kinematic age and the properties of the star clusters. Top histograms in right panel shows the distribution of the kinematic age (blue colour) and the mean age of the clusters (weighted by their total mechanical energy input during the last 10 Myr; orange colour). Orange colour on the bottom histograms considers only clusters older than 4 Myr and corresponds to the time passed since the first SN (assuming it is happening at 4 Myr according to Fig. 6). Only the regions with at least one stellar association or star cluster are shown.

In the text

Fig. 12. Distribution of the measured kinetic energy of the ionised gas in the identified regions of locally elevated Hα velocity dispersion (left panel) and the masses of the star clusters in these regions (right panel) shown for different metallicity cuts. These regions tend to have lower kinetic energy in the low-metallicity environment, but also the typical masses of clusters decline. The bin for Z < 0.1 Z⊙ corresponds to the estimates of kinetic energy for 16 ionised gas superbubbles identified in the dwarf galaxies Sextans A (Gerasimov et al. 2022) and DDO 53 (Egorov et al. 2021) and is shown here for comparison.

In the text

Fig. 13. Properties of the regions of elevated velocity dispersion in comparison to those of normal H II regions. Blue and orange contours show the statistical density distribution of the H II regions (associated with a single star cluster, from the catalogue compiled by Scheuermann et al. 2023) and the regions with locally elevated σ(Hα), respectively. The Hα luminosity (L(Hα)) correlates with the total mass of the clusters (Mcls) for both samples (left panel), but the slope is different implying differences in the escape fraction of the ionising photons and/or contribution of the additional ionising sources. σ(Hα) for the regions in this study slightly decreases with the Mcls (middle panel) and L(Hα) (right panel) contrary to what is observed in the reference H II regions.

In the text

Fig. 14. Distribution of the velocity dispersion of Hα (left) and CO (right) emission within the borders of young stellar associations of different ages (orange and blue colours for those younger and older than 4 Myr, respectively) and outside them (green colour). Only pixels with S/N > 30 in the corresponding emission line are considered. Note that increasing the S/N limit does not affect the result.

In the text

Fig. 15. Dependence of the kinetic energy of ionised gas Ekin in the regions of locally elevated velocity dispersion on the mass of the stellar associations within them (left panels), and their gas kinetic energy per area, ΣEkin, vs. corresponding stellar mass surface density, Σ⋆ (right panels). Top and bottom panels correspond to minimal and maximal values of ne underlying the Ekin measurements. Colour encodes the Hα flux surface brightness. The red line shows the median values of Ekin or ΣEkin in several bins, and the red shaded area corresponds to the 25–75 percentile interval. Spearman correlation coefficients ρ are given above the plots.

In the text

Fig. 16. Dependence of the kinetic energy of ionised gas Ekin in the regions of locally elevated velocity dispersion on the total mechanical energy input from the stellar associations in the form of supernovae and stellar winds (left panels), and supernovae alone (right panels). Top and bottom panels correspond to the $n_{\mathrm{e}}^{(min)}$ and $n_{\mathrm{e}}^{(max)}$ density measurements (thus representing upper and lower limits of Ekin, respectively). Median value of logarithmic errors is shown in bottom-right corner of each panel. All regions with errors 2 times larger than presented are excluded from the plots. Colour encodes the [S II]/Hα lines ratio tracing the dominant gas excitation mechanism ([S II]/Hα > 0.4 is likely produced by shocks; left panels), or the relative contribution of the SNe to the total mechanical energy input according to the STARBURST99 models for the corresponding mass and age of each cluster right panels. Blue, orange, green and red lines correspond to mechanical stellar feedback energy efficiency η = 100, 30, 5, 1%, respectively. Spearman correlation coefficient ρ is given above the plots.

In the text

	Fig. 17. Distribution of the estimated mechanical stellar feedback energy efficiency for the regions shown in top-left panel of Fig. 16 grouped by metallicity.
In the text

Fig. 18. Comparison of the kinetic energy of ionised gas and the mechanical energy input by young stars for our best candidate superbubbles (as exhibiting shell-like morphology). Top plot: same comparison of energies as in Fig. 16, superbubbles follow the same trend as the other regions. Colour encodes the age of the youngest association in the superbubble. Bottom plot: distribution of the estimated efficiency of the energy consumption in superbubbles. Its peak is slightly below η = 0.25 that was obtained from the modelling by Sharma et al. (2014) as an asymptotic value for the expanding superbubbles driven by multiple (up to a few hundred) OB stars in normal star-forming galaxies.

In the text

Fig. 19. Relative fraction of the stellar associations younger than 3 Myr at different separations between clusters/associations and centres of the regions. Only a sub-sample of 40 well-resolved superbubbles and star clusters younger than 40 Myr are considered. Each separation bin corresponds to 0.3Reff. Errors are calculated from 1000 Monte Carlo realisations with positions of each cluster randomly varying within the PSF of the MUSE observations.

In the text

	Fig. A.1. I − σ diagram (right), classification map (centre) and the Hα map (left) for all PHANGS-MUSE galaxies except NGC 4254 (shown in Fig. 4). Cyan ellipses encircle the identified regions of the locally elevated Hα velocity dispersion.
In the text

	Fig. A.1. Continued.
In the text

	Fig. A.1. Continued.
In the text

	Fig. A.1. Continued.
In the text

	Fig. A.1. Continued.
In the text

	Fig. A.1. Continued.
In the text