<figure><img src="image_1.png" /></figure>

HAL
open science

# Quantifying the energy balance between the turbulent ionised gas and young stars

Oleg V Egorov, Kathryn Kreckel, Simon C O Glover, Brent Groves, Francesco Belfiore, Eric Emsellem, Ralf S Klessen, Adam K Leroy, Sharon E Meidt, Sumit K Sarbadhicary, et al.

► To cite this version:

Oleg V Egorov, Kathryn Kreckel, Simon C O Glover, Brent Groves, Francesco Belfiore, et al. Quantifying the energy balance between the turbulent ionised gas and young stars. Astronomy & Astrophysics - A&A, 2023, 678, pp.A153. 10.1051/0004-6361/202346919. hal-05018654

HAL Id: hal-05018654

https://hal.science/hal-05018654v1

Submitted on 3 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Quantifying the energy balance between the turbulent ionised gas and young stars

Oleg V. Egorov¹, Kathryn Kreckel¹, Simon C. O. Glover², Brent Groves³, Francesco Belfiore⁴,
Eric Emsellem⁵,⁶, Ralf S. Klessen²,⁷, Adam K. Leroy⁸,⁹, Sharon E. Meidt¹⁰, Sumit K. Sarbadhicary⁸,⁹,
Eva Schinnerer¹¹, Elizabeth J. Watkins¹, Brad C. Whitmore¹², Ashley T. Barnes⁵,¹³, Enrico Congiu¹⁴,
Daniel A. Dale¹⁵, Kathryn Grasha¹⁶, Kirsten L. Larson¹⁷, Janice C. Lee¹⁸,¹⁹, J. Eduardo Méndez-Delgado¹,
David A. Thilker²⁰, and Thomas G. Williams²¹

*(Affiliations can be found after the references)*

Received 16 May 2023 / Accepted 17 August 2023

## ABSTRACT

**Context.** Stellar feedback is a key contributor to the morphology and dynamics of the interstellar medium in star-forming galaxies. In particular, energy and momentum input from massive stars can drive the turbulent motions in the gas, but the dominance and efficiency of this process are unclear. The study of ionised superbubbles enables quantitative constraints to be placed on the energetics of stellar feedback.

**Aims.** We directly compare the kinetic energy of expanding superbubbles and the turbulent motions in the interstellar medium with the mechanical energy deposited by massive stars in the form of winds and supernovae. With such a comparison, we aim to determine whether the stellar feedback is responsible for the observed turbulent motions and to quantify the fraction of mechanical energy retained in the superbubbles.

**Methods.** We investigated the ionised gas morphology, excitation properties, and kinematics in 19 nearby star-forming galaxies from the PHANGS-MUSE survey. Based on the distribution of the flux and velocity dispersion in the Hα line, we selected 1484 regions of locally elevated velocity dispersion ($σ(Hα) > 45$ km s⁻¹), including at least 171 expanding superbubbles. We analysed these regions and related their properties to those of the young stellar associations and star clusters identified in PHANGS-HST data.

**Results.** We find a good correlation between the kinetic energy of the ionised gas and the total mechanical energy input from supernovae and stellar winds from the stellar associations. At the same time, the contribution of mechanical energy injected by the supernovae alone is not sufficient to explain the measured kinetic energy of the ionised gas, which implies that pre-supernova feedback in the form of radiation and thermal pressure as well as winds is necessary. We find that the gas kinetic energy decreases with metallicity for our sample covering Z = 0.5–1.0Z⊙, reflecting the lower impact of stellar feedback. For the sample of well-resolved superbubbles, we find that about 40% of the young stellar associations are preferentially located in their rims. We also find a slightly higher (by ~15%) fraction of the youngest (<3 Myr) stellar associations in the rims of the superbubbles than in the centres and the opposite trend for older associations, which implies possible propagation or triggering of star formation.

**Conclusions.** Stellar feedback is the dominant source for powering the ionised gas in regions of locally (on a 50–500 pc scale) elevated velocity dispersion, with a typical coupling efficiency of 10–20%. Accounting for pre-supernovae feedback is required to set up the energy balance between gas and stars.

**Key words.** galaxies: ISM – ISM: kinematics and dynamics – ISM: bubbles – galaxies: star formation

## 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 ($\Sigma_{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 superbuble 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 ($\Sigma_g \ge 50 M_\odot pc^{-2}$). 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 $\eta$ 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<fnref n="1" /> 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.′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_{25}$). In the present paper, we use the publicly available<fnref n="2" /> data cubes and maps ('copt' dataset with homogenised PSFs) with the distribution of properties of the ionised gas.

<footnote n="1">¹ http://phangs.org</footnote>

<footnote n="2">² https://archive.eso.org/scienceportal/home?data_collection=PHANGS</footnote>
<table><caption>**Table 1.** Sample of analysed galaxies, their properties and the number of identified regions with elevated velocity dispersion.</caption><thead><tr><th rowspan="2">Galaxy</th><th rowspan="2"><i>D</i><sup>(1)</sup><br>(Mpc)</th><th rowspan="2"><i>R</i><sub>25</sub><sup>(2)</sup><br>(kpc)</th><th rowspan="2">log <i>M</i><sub>*</sub><sup>(3)</sup><br>(<i>M</i><sub>☉</sub>)</th><th rowspan="2">log(SFR)<sup>(3)</sup><br>(<i>M</i><sub>☉</sub> yr<sup>-1</sup>)</th><th rowspan="2">12 + log(O/H)<sup>(4)</sup><br>(dex at <i>R</i><sub>eff</sub>)</th><th rowspan="2">PSF<sup>(5)</sup><br>(")</th><th rowspan="2">σ(Hα)<sub>m</sub><sup>(6)</sup><br>(km s<sup>-1</sup>)</th><th rowspan="2">log Σ(Hα)<sub>DIG</sub><sup>(6)</sup><br>(erg s<sup>-1</sup> cm<sup>-2</sup> arcsec<sup>-2</sup>)</th><th colspan="3">Number of regions<sup>(6)</sup></th></tr><tr><th><i>N</i><sub>tot</sub></th><th><i>N</i><sub>HST</sub></th><th><i>N</i><sub>cls</sub></th></tr></thead><tbody><tr><td>IC5332</td><td>9.01</td><td>8.0</td><td>9.67</td><td>-0.39</td><td>8.30</td><td>0.87</td><td>23.0 ± 5.2</td><td>-16.68</td><td>44</td><td>41</td><td>16</td></tr><tr><td>NGC628</td><td>9.84</td><td>14.1</td><td>10.34</td><td>0.24</td><td>8.48</td><td>0.92</td><td>23.4 ± 5.6</td><td>-16.28</td><td>116</td><td>69</td><td>30</td></tr><tr><td>NGC1087</td><td>15.85</td><td>6.9</td><td>9.93</td><td>0.12</td><td>8.41</td><td>0.92</td><td>29.4 ± 4.9</td><td>-16.18</td><td>28</td><td>28</td><td>19</td></tr><tr><td>NGC1300</td><td>18.99</td><td>16.4</td><td>10.62</td><td>0.07</td><td>8.54</td><td>0.89</td><td>27.3 ± 5.3</td><td>-16.58</td><td>64</td><td>63</td><td>10</td></tr><tr><td>NGC1365</td><td>19.57</td><td>34.2</td><td>10.99</td><td>1.23</td><td>8.48</td><td>1.15</td><td>29.4 ± 4.8</td><td>-16.48</td><td>55</td><td>20</td><td>2</td></tr><tr><td>NGC1385</td><td>17.22</td><td>8.5</td><td>9.98</td><td>0.32</td><td>8.42</td><td>0.77</td><td>31.9 ± 4.6</td><td>-16.38</td><td>108</td><td>108</td><td>33</td></tr><tr><td>NGC1433</td><td>18.63</td><td>16.8</td><td>10.87</td><td>0.05</td><td>8.56</td><td>0.91</td><td>26.8 ± 4.7</td><td>-16.78</td><td>103</td><td>50</td><td>6</td></tr><tr><td>NGC1512</td><td>18.83</td><td>23.1</td><td>10.71</td><td>0.11</td><td>8.57</td><td>1.25</td><td>22.6 ± 5.0</td><td>-16.98</td><td>51</td><td>51</td><td>10</td></tr><tr><td>NGC1566</td><td>17.69</td><td>18.6</td><td>10.78</td><td>0.66</td><td>8.58</td><td>0.80</td><td>30.1 ± 5.5</td><td>-16.23</td><td>107</td><td>102</td><td>38</td></tr><tr><td>NGC1672</td><td>19.40</td><td>17.4</td><td>10.73</td><td>0.88</td><td>8.55</td><td>0.96</td><td>26.7 ± 4.8</td><td>-16.38</td><td>93</td><td>93</td><td>19</td></tr><tr><td>NGC2835</td><td>12.22</td><td>11.4</td><td>10.00</td><td>0.09</td><td>8.40</td><td>1.15</td><td>23.1 ± 4.6</td><td>-16.28</td><td>58</td><td>48</td><td>28</td></tr><tr><td>NGC3351</td><td>9.96</td><td>10.5</td><td>10.36</td><td>0.12</td><td>8.59</td><td>1.05</td><td>24.6 ± 5.1</td><td>-16.38</td><td>43</td><td>37</td><td>13</td></tr><tr><td>NGC3627</td><td>11.32</td><td>16.9</td><td>10.83</td><td>0.58</td><td>8.54</td><td>1.05</td><td>30.6 ± 5.2</td><td>-15.98</td><td>102</td><td>101</td><td>47</td></tr><tr><td>NGC4254</td><td>13.10</td><td>9.6</td><td>10.42</td><td>0.49</td><td>8.56</td><td>0.89</td><td>30.1 ± 4.8</td><td>-16.18</td><td>137</td><td>137</td><td>61</td></tr><tr><td>NGC4303</td><td>16.99</td><td>17.0</td><td>10.52</td><td>0.73</td><td>8.58</td><td>0.78</td><td>28.5 ± 5.5</td><td>-15.98</td><td>108</td><td>108</td><td>58</td></tr><tr><td>NGC4321</td><td>15.21</td><td>13.5</td><td>10.75</td><td>0.55</td><td>8.56</td><td>1.16</td><td>26.8 ± 5.0</td><td>-16.28</td><td>58</td><td>46</td><td>27</td></tr><tr><td>NGC4535</td><td>15.77</td><td>18.7</td><td>10.53</td><td>0.33</td><td>8.54</td><td>0.56</td><td>26.9 ± 5.7</td><td>-16.48</td><td>63</td><td>48</td><td>15</td></tr><tr><td>NGC5068</td><td>5.20</td><td>5.7</td><td>9.40</td><td>-0.56</td><td>8.32</td><td>1.04</td><td>24.3 ± 5.5</td><td>-16.18</td><td>89</td><td>89</td><td>38</td></tr><tr><td>NGC7496</td><td>18.72</td><td>9.1</td><td>10.00</td><td>0.35</td><td>8.51</td><td>0.89</td><td>28.7 ± 4.8</td><td>-16.48</td><td>57</td><td>57</td><td>16</td></tr></tbody></table>

**Notes.** $N_{tot}$, $N_{HST}$, and $N_{cls}$ denote the total number of regions with elevated velocity dispersion, the number of those with available HST data, and the number of regions containing young star clusters or associations, respectively. (1) Anand et al. (2021), includes measurements from Jacobs et al. (2009; NGC 628, NGC 1365, NGC 3351, NGC 3627), Kourkchi & Tully (2017; NGC 1087, NGC 1566, NGC 4303), Shaya et al. (2017), Kourkchi et al. (2020; NGC 1300, NGC 1385, NGC 1672, NGC 7496), Scheuermann et al. (2022; NGC 1433, NGC 1512), Nugent et al. (2006; NGC 4254) and Freedman et al. (2001; NGC 4321, NGC 4535). (2) From HyperLEDA Makarov et al. (2014); (3) Leroy et al. (2021); (4) Groves et al. (2023); (5) Emsellem et al. (2022); (6) This work.

The maps were obtained with a data analysis pipeline (DAP<fnref n="3" />; based on the GIST<fnref n="4" /> 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, $\sigma_{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(Hα) \approx \sqrt{\sigma_{obs}(Hα)^2 - \sigma_{LSF}(Hα)^2}$, where $\sigma_{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 $\sigma_{LSF}(Hα)$ is 48.9 km s⁻¹ after averaging over all sample galaxies. We do not correct $\sigma(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 $\Sigma_*$ 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 $\sigma(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 HII 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 HII 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 $\sigma(Hα) \sim 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

<footnote n="3">³ https://gitlab.com/francbelf/ifu-pipeline</footnote>

<footnote n="4">⁴ https://abittner.gitlab.io/thegistpipeline</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Statistics of the intrinsic $H\alpha$ 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\alpha$ are excluded. The blue dashed line corresponds to the mean MUSE spectral resolution (velocity dispersion) at the wavelength of $H\alpha$. The values of $\sigma(H\alpha)$ were corrected for instrumental broadening by subtracting, in quadrature, the MUSE instrumental velocity dispersion for each galaxy.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Relation between the recovered $\sigma(H\alpha)$ and the true $\sigma(H\alpha)$ 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\alpha$ is $FWHM = 2.54 \pm 0.1$ Å (according to Bacon et al. 2017, but with twice larger scatter; vertical dashed line on panel a).</figcaption></figure>

histograms in Fig. 1 show the distribution of $\sigma(H\alpha)$. Only in some environments, such as bars and galaxy centres, is $\sigma(H\alpha)$ above the MUSE resolution. Our analysis is focused, however, on areas of locally elevated intrinsic velocity dispersion (exceeding 45 km s$^{-1}$, 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 $\sigma(H\alpha)$.

To test the precision with which we can trace the relative changes of the $H\alpha$ 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 $\sigma = 10–100$ km s$^{-1}$). These profiles were then convolved with a Gaussian with the width of the MUSE instrumental profile around its mean value at $H\alpha$ (FWHM = $2.54 \pm 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<fnref n="5" />. We then added noise (both Gaussian and Poisson) according to the adopted signal-to-noise ratio (S/N) and fitted a Gaussian<fnref n="6" /> 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 $\sigma(H\alpha)$ 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\alpha$ line with MUSE with a precision better than 30% (median uncertainty is ~12%) starting from $\sigma(H\alpha) \sim 40$ km s$^{-1}$ 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 $\sigma(H\alpha)$.

<footnote n="5">5 The mosaics in PHANGS data were obtained from individual exposures obtained at 4 different orientations in order to minimise the instrumental artefacts and spatial variations of the LSF. The individual exposures in the deep field from Bacon et al. (2017) were obtained at 8 different orientations.</footnote>

<footnote n="6">6 Law et al. (2021b) compared the results of fitting the simple Gaussian and more realistic inverse gamma error distribution to the modelled line profiles and found a small (≤10%) difference in $\sigma(H\alpha)$ when both $\sigma(H\alpha)$ and S/N are low, and negligible difference at high S/N.</footnote>
Panel b shows that for the regions having $\sigma(H\alpha) > 45$ km s$^{-1}$ and $S/N > 30$, we can expect the precision of the $\sigma(H\alpha)$ measurements to be better than ~5 km s$^{-1}$, while the velocity dispersion is typically significantly overestimated with MUSE for $\sigma(H\alpha) < 25$ km s$^{-1}$ (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.

## 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<fnref n="7" /> 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<fnref n="8" /> 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 \times 10^3 M_{\odot}$ 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 $\sigma(H\alpha)$ are observed in all 19 PHANGS-MUSE galaxies. Those areas with small angular size and high $H\alpha$ 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\alpha$ flux and velocity dispersion for one galaxy from our sample: NGC 4254. The $H\alpha$ and [N II] line profiles in panel d demonstrate the reliability of the observed variations in $\sigma(H\alpha)$. 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\alpha$ velocity dispersion.

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 $\sigma(H\alpha)$ 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 - \sigma$ diagnostics

The first step in this analysis is the automated identification of regions with spatially correlated $\sigma(H\alpha)$ excess. For that, we rely on an analysis of the $I - \sigma$ diagrams representing the pixel-by-pixel variation of the velocity dispersion $\sigma(H\alpha)$ versus the logarithm of the $H\alpha$ intensity, $\log I(H\alpha)$. 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 - \sigma$ diagram as a triangle-shaped region of intermediate $I(H\alpha)$ and elevated $\sigma(H\alpha)$ 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\alpha$ 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 - \sigma$ 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 - \sigma$

<footnote n="7">⁷ NGC 628, NGC 1433, NGC 1512, NGC 1566 and part of the data for NGC 3351 and NGC 3627.</footnote>

<footnote n="8">⁸ https://archive.stsci.edu/hlsp/phangs-cat</footnote>
<figure><img src="image_5.png" /></figure>

**Fig. 3.** Distribution of the $H\alpha$ 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\alpha$ brightness. Panel c demonstrates a zoom-in $H\alpha$ 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 $\sigma(H\alpha)$ towards the bright regions. Panel d shows the examples of the $H\alpha$ and [N II] λ6548, 6584 Å line profiles extracted from three regions with low (green line) and high (orange and magenta lines) $H\alpha$ velocity dispersion. These line profiles corresponding to high $\sigma(H\alpha)$ 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.

<figure><img src="image_6.png" /></figure>

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

diagnostics in a form similar to that described in Egorov et al. (2021). The $I - \sigma$ 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.

Compared with previous studies, the application of $I - \sigma$ 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 \gtrsim 30$ and better we are able to reliably measure the relative variations of the intrinsic velocity dispersion at the level of $\sigma(H\alpha) \gtrsim 30 \text{ km s}^{-1}$ with a typical uncertainty better than 5 km s$^{-1}$. To ensure those measurement uncertainties do not bias our analysis, we masked out all pixels with $S/N < 30$ in $H\alpha$. The majority of the areas of elevated velocity dispersion that we investigate in this paper have $\sigma(H\alpha) > 45 \text{ km s}^{-1}$ and thus provide reliable measurements of $\sigma(H\alpha)$ 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 con-
structed the *I* − *σ* diagrams for each galaxy and analysed them to
select the regions with spatially correlated elevated velocity dis-
persion. As a short summary, we performed the following steps
(see Appendix A for more details):

1. Isolated the peaks of high velocity dispersion correspond-
ing 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 show-
ing 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 × σ<sub>A</sub> as size of the regions along
their major or minor axes, where σ<sub>A</sub> is the standard deviation
along these axes of the 2D distribution resulting from the statis-
tics 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 ran-
domly 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 inte-
grated 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 iso-
late 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<sub>V</sub> = 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_{\text{ion}}$, in the regions. The first one assumes a thin spherical shell geometry for all of them. Then

$$M_{\text{ion}}^{(1)} \approx 4\pi R_{\text{eff}}^2 \Delta R n_e \mu m_{\text{H}} = 4\Omega_{\text{reg}} D^2 \Delta R n_e \mu m_{\text{H}}, \quad (1)$$

where $R_{\text{eff}} = \sqrt{R_{\text{min}}R_{\text{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), $\Delta R$ is the thickness of the shell, $n_e$ is the electron density, $m_{\text{H}}$ is the mass of a hydrogen atom, the coefficient $\mu = 1.27$ accounts for a singly ionised helium contribution, $\Omega_{\text{reg}}$ is the total area of the region in the plane of the sky, and $D$ is the distance to the galaxy. To calculate $\Delta R$, we can relate the emission measure $EM$ to the surface brightness of Hα, $I(Hα)$ as (Reynolds 1977; Pengelly 1964):

$$EM = 2.75 \left( \frac{T_e}{10^4 \text{ K}} \right)^{0.9} I(H\alpha) \text{pc cm}^{-6}, \quad (2)$$

where $I(Hα)$ is in units of Rayleighs (1 Rayleigh $\approx 5.69 \times 10^{-18}$ ergs s$^{-1}$ cm$^{-2}$ arcsec$^{-2}$ 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\Delta R$ for a spherical shell geometry) as $EM = \int n_e n_{H^+} ds \approx \int n_e^2 ds \approx 2\Delta R n_e^2$. Thus, we obtain

$$\Delta R \approx 24.2 \left( \frac{T_e}{10^4 \text{ K}} \right)^{0.9} \frac{F(H\alpha)\Omega_c^{-1}}{10^{-16} \text{ erg s}^{-1} \text{ cm}^{-2} \text{ arcsec}^{-2}} \left( \frac{n_e}{\text{cm}^{-3}} \right)^{-2} \text{pc}, \quad (3)$$

where $F(Hα) = I(Hα)\Omega_c$ is the total reddening-corrected Hα flux from the central aperture used in the spectral extraction, which has an area of $\Omega_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), $\Omega_{\text{reg}}/\Omega_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

$$M_{\text{ion}}^{(1)} \approx 522 \frac{F(H\alpha)}{10^{-16} \text{ erg s}^{-1} \text{ cm}^{-2}} \left(\frac{D}{\text{Mpc}}\right)^2 \left(\frac{n_e}{\text{cm}^{-3}}\right)^{-1} M_{\odot}. \quad (4)$$

Note that this method does not depend on how well we can measure $R_{\text{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α)_{\text{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(H^0) = \alpha_B n_H^2 V \approx \alpha_B n_e^2 V$ when the volume is in ionisation equilibrium, where $\alpha_B$ is the case B recombination coefficient and $\langle n_e^2 \rangle$ is the mean squared electron density within $V^9$. If we assume, in the absence of evidence to the contrary, that the density distribution is uniform, this relation simplifies to $Q(H^0) = \alpha_B n_e^2 V$. For our adopted $T_e = 8000$ K, interpolation of the values in Storey & Hummer (1995) yields $\alpha_B \approx 3.1 \times 10^{-13} \text{ cm}^3 \text{ s}^{-1}$, and the ionising photon flux is related to the Hα luminosity by $Q(H^0) \approx 7.1 \times 10^{11} L(H\alpha)$ (Osterbrock & Ferland 2006). Given that $L(Hα) = 4\pi D^2 F(Hα)_{\text{tot}}$, we then have

$$M_{\text{ion}}^{(2)} = \mu m_{\text{H}} n_e V \approx 29.3 \frac{F(H\alpha)_{\text{tot}}}{10^{-16} \text{ erg s}^{-1} \text{ cm}^{-2}} \left(\frac{D}{\text{Mpc}}\right)^2 \left(\frac{n_e}{\text{cm}^{-3}}\right)^{-1} M_{\odot}. \quad (5)$$

As our adopted value of $M_{\text{ion}}$, we use $M_{\text{ion}}^{(1)}$ for the regions having a size along the minor axis greater than 3×PSF, and $M_{\text{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_{\text{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_{\text{ion}}$, but difficult because of several observational and physical limitations. It is typically measured using density-sensitive line ratios like $R_{\text{[SII]}} = [\text{S II}]\lambda6717\AA/[\text{S II}]\lambda6731\AA$ (e.g., Osterbrock & Ferland 2006; Kewley et al. 2019). Such diagnostics have low-density limits (e.g., $n_e \sim 30 \text{ cm}^{-3}$ for $R_{\text{[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_{\text{ion}}$ (Eqs. (4) and (5)), we estimate the electron density $n_e$ from the $R_{\text{[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_{\text{[SII]}}$ falls below the value of $1.46^{10}$ by at least $3\sigma$ (chosen to minimise the bias towards higher densities if $R_{\text{[SII]}}$ is

<footnote n="9">⁹ Accounting for helium makes only a small difference to this number.</footnote>

<footnote n="10">¹⁰ This value corresponds to the low-density limit at $T_e = 8000$ K (Barnes et al. 2021).</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 5. Mass of the ionised gas in the selected regions assuming different prescriptions for electron density: lower limit ($n_e^{(min)}$) and upper limit ($n_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 $n_{e,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.</figcaption></figure>

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_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_e^{(max)}$ derived in the latter way.

The estimates of $n_e^{(max)}$ obtained from $R_{[S II]}$ 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

$$n_e^{(\text{min})} = \sqrt{\frac{3Q(H^0)}{4\pi R_{\text{eff}}^3 \alpha_B}} \approx 1.36 \times 10^{-16} \left(\frac{L(H\alpha)}{\text{erg s}^{-1}}\right)^{0.5} \left(\frac{R_{\text{eff}}}{\text{pc}}\right)^{-1.5} \text{cm}^{-3}. \quad (6)$$

Applying this to the selected reference H II regions in the same way as above, we estimate $n_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 \times \text{PSF}$ are excluded). In the ideal case, $n_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_e^{(min)}$ corresponds to the lower limit of the electron density.

We use the estimates of $n_e^{(min)}$ and $n_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_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_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_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_e^{(min)}$) is closer to the real value, while its lower limit (based on $n_e^{(max)}$) better traces the relative changes between the regions. The median formal uncertainties of $\log(M_{ion})[n_e^{(min)}]$ and $\log(M_{ion})[n_e^{(max)}]$ are equal to 0.09 dex and 0.29 dex, respectively. Therefore our further analysis relies mainly on $M_{ion}[n_e^{(min)}]$, but we compare the results with measurements made assuming $M_{ion}[n_e^{(max)}]$.

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

Assuming that the measured elevated $H\alpha$ 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:

$$E_{\text{kin}} \approx E_{\text{turb}} \approx \frac{3}{2} M_{\text{ion}} (\sigma(H\alpha)^2 - \sigma(H\alpha)^2_{\text{m}}). \quad (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)<fnref n="11" /> and Geneva stellar evolution tracks which include rotation with two metallicities

<footnote n="11">¹¹ Note that the parameters of the star clusters have been determined using a Chabrier IMF (Chabrier 2003). This is very similar to the</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Evolution of the cumulative mechanical energy input to the ISM ($E_{\text{mech}}$) normalised by the mass of star cluster ($M_{\text{cls}}$) 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.</figcaption></figure>

($Z = 0.014 \sim Z_{\odot}$, Ekström et al. 2012, and $Z = 0.002 \sim 0.15 Z_{\odot}$, Georgy et al. 2013), for clusters with total stellar mass $\log(M_{*}/M_{\odot}) = (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_{\text{WSN}}$) and by winds only ($E_{\text{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_{\odot}$). 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_{\text{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_{\text{WSN}}$ and $E_{\text{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_{\text{WSN}}$ and $E_{\text{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.

# 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\alpha$ 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\alpha$ 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\alpha$ flux distribution and thus are classified as compact objects, with both $H\alpha$ 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\alpha$ 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).

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\text{II}]/H\alpha$ 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\text{II}]/H\alpha$ above or below 0.4, which is often used to separate supernovae remnants from photoionised nebulae (Blair et al. 1981; Blair & Long 2004).

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 $\sigma(H\alpha)$ in the rest of the regions – 73% of them have $[S\text{II}]/H\alpha > 0.4$, making them probable SNR candidates. Moreover, 74% of all regions with $[S\text{II}]/H\alpha > 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\text{II}]/H\alpha$ 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\alpha$ luminosity $L(H\alpha) < 2.5 \times 10^{37}$ erg s$^{-1}$ 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 HII regions at the ~50 pc PHANG-MUSE resolution and contribute to the ionisation balance (Scheuermann et al. 2023), but these single stars will

<footnote>Kroupa IMF, but there are small differences in their mass-to-light ratio (by a factor of ~1.08 according to Madau & Dickinson 2014). Therefore, we may expect a negligible systematic offset ($\log(\Delta M) \sim 0.03$ dex) between the measured star cluster masses and those in the STARBURST99 models.</footnote>
<figure><img src="image_9.png" /><figcaption>**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\alpha$ 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\alpha$ flux distribution, the $[S\,II]/H\alpha$ line ratio, and the $H\alpha$ 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\alpha$ flux distribution relative to the centre of each region (marked by crosses).</figcaption></figure>

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\alpha$ 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 $\sigma(H\alpha)$, 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-$\sigma$ relation observed in nearby galaxies of different types – the $H\alpha$ 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^{51}$ 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.

## 4.2. Size distribution

The distribution of the sizes (effective radii, $R_{\text{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
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

~70 pc) are described very well by a power-law: $p(R_{\text{eff}}) \propto R_{\text{eff}}^{\alpha}$. Using the POWERLAW package (Alstott et al. 2014), we got a best-fit value of $\alpha = -2.71 \pm 0.19$. The power-law breaks at radii of $R_{\text{eff}} \sim 60-80$ pc, and also at $R_{\text{eff}} \ge 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.

The measured slope of a power-law approximating our regions in Fig. 10 is very similar to $\alpha = -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 ($\alpha \approx -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 ($\alpha \approx -3.2 \pm 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 ($\alpha \approx -2.2 \pm 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_{\text{exp}}$) and kinematic age ($t_{\text{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 $\sigma(H\alpha)$ in order to estimate the corresponding $V_{\text{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

$$V_{\text{exp}} \approx 3.15 \times (\sigma(H\alpha)^2 - \sigma(H\alpha)_m^2)^{0.4} - 4.1 \text{ km s}^{-1}, \quad (8)$$

where the exact coefficients depend on the value of $\sigma(H\alpha)_m$ (unperturbed velocity dispersion in the galaxy) and are parameterised in that paper<fnref n="12" />, and $\sigma(H\alpha)$ 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_{\text{kin}} \approx 0.6R/V_{\text{exp}}$ when the bubble is in the adiabatic stage of its evolution. Figure 11 (left) shows that the expansion velocities, $V_{\text{exp}}$, of the identified regions with elevated velocity dispersion decline with their kinematic age. That is expected given that the bubbles slow down

<footnote n="12">¹² The parameterisation was obtained relying on mock spectra of toy-models of expanding bubbles in the turbulent ISM with properties representative of those found in PHANGS-MUSE galaxies. The coefficients in Eq. (8) correspond to the median $\sigma(H\alpha)_m$ in our sample.</footnote>
<figure><img src="image_11.png" /><figcaption>**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\alpha$ 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).</figcaption></figure>
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’.

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_{cls}^{SN} = t_{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):

$$V_{exp}(t) \approx 38.6 \left(\frac{L_{38}}{n_0}\right)^{1/5} t_6^{-2/5} \text{ km s}^{-1}, \quad (9)$$

where $L_{38} = L_{\text{mech}}/10^{38}$ erg s⁻¹ is the mechanical luminosity of the stellar association, $n_0$ is the ambient density, and $t_6$ is $t_{kin}$ expressed in Myr. Considering $n_e^{(\text{min})}$ as a proxy of density and $L_{\text{mech}} \approx 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} \sim 10-60 \text{ km s}^{-1}$ (mostly within $20-35 \text{ km s}^{-1}$). Thus, the large spread of $E_{\text{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_{\text{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_{\text{kin}}$ change with the metallicity of a region in Fig. 12 (left panel). Regions with lower oxygen abundance tend to have lower

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Distribution of sizes for regions of high $H\alpha$ velocity dispersion. The red line shows a result of a power-law fit to the data for $R_{\text{eff}} = 70-300$ pc. This fit agrees with a power-law approximation for the sizes of superbubbles according to hydrodynamical simulations ($\alpha = -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.</figcaption></figure>

kinetic energy, though the number of regions with $Z < 0.5 Z_{\odot}$ is rather low for drawing definitive conclusions. Note, however, that measurements of energetics of low-metallicity ($Z < 0.1 Z_{\odot}$) superbubbles available in the literature (Egorov et al. 2021; Gerasimov et al. 2022) follow the same trend as in our data.

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_{\odot}$ 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_{\text{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_{\text{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\alpha$ velocity dispersion almost independently
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Kinematic properties of the regions with elevated velocity dispersion. Left panel shows the relation between the estimated expansion velocity ($V_{exp}$) 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 $V_{exp}$ 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.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Distribution of the measured kinetic energy of the ionised gas in the identified regions of locally elevated $H\alpha$ 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_{\odot}$ 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.</figcaption></figure>

of their morphology, so they do not necessarily have the same properties as typical HII regions. In Fig. 13 we examine some of the relations between the $H\alpha$ 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 HII regions. We use a sub-sample of HII regions from the PHANGS-MUSE nebulae catalogue, com-

piled 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.

Massive stars in the young stellar associations and compact clusters produce the hydrogen-ionising photons that are responsible for the observed $H\alpha$ flux from the ionised gas surrounding
<figure><img src="image_15.png" /><figcaption>**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 $\sigma(H\alpha)$, respectively. The H$\alpha$ luminosity ($L(H\alpha)$) correlates with the total mass of the clusters ($M_{\text{cls}}$) 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. $\sigma(H\alpha)$ for the regions in this study slightly decreases with the $M_{\text{cls}}$ (middle panel) and $L(H\alpha)$ (right panel) contrary to what is observed in the reference H II regions.</figcaption></figure>

the clusters. Therefore one naturally expects to see a correlation between the mass of the star clusters and the H$\alpha$ 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$\alpha$ 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$\alpha$ 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$\alpha$ luminosities as the H II regions.

By construction, the regions considered in this paper have higher H$\alpha$ velocity dispersion than normal H II regions. The middle and right panels of Fig. 13 demonstrate that the distributions of $\sigma(H\alpha)$ 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$\alpha$ luminosity or the mass of the ionising star clusters and the H$\alpha$ 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 – $\sigma(H\alpha)$ 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 $\sigma(H\alpha)$ 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 $\sigma(H\alpha)$ between the pixels associated with the young stellar associations older than 4 Myr age and those not overlapping with any stellar associations. Surprisingly, $\sigma(H\alpha)$ 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 $\sigma(H\alpha)$), we find the separation between the two distributions is 5–11 km s$^{-1}$, and the probability that they are equal is negligible.

It is not obvious why pixels corresponding to the locations of the youngest stellar associations show lower H$\alpha$ 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<fnref n="13" />. The elevated $\sigma(H\alpha)$ 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$\alpha$ emission. Namely, we suggest that the H$\alpha$ emission around all young stellar associations is tracing the kinematics of the thin star-forming gas disc, but the diffuse H$\alpha$ outside the youngest of them is more impacted by turbulence and outflows. A higher mean $\sigma(H\alpha)$ 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

<footnote n="13">¹³ Note that the difference in absolute values of $\sigma(H\alpha)$ and $\sigma(\text{CO})$ is consistent with what measured in other objects (e.g., Girard et al. 2021).</footnote>
<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Distribution of the velocity dispersion of $H\alpha$ (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.</figcaption></figure>

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\alpha$ 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_{\text{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 $\Sigma E_{\text{kin}}$ with the mean total stellar mass surface density $\Sigma_{\star}$ (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 $\Sigma E_{\text{kin}}$ is largely independent of $\Sigma_{\star}$, implying that gravitation alone<fnref n="14" /> is unlikely to be the dominant powering source of the turbulence. In contrast, $E_{\text{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_e^{(\text{max})}$ leads to a more linear trend.

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_{\text{kin}}$. 64% of them have $[S\text{II}]/H\alpha > 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_{\text{kin}}$ are for the lower density limit $n_e^{(\text{min})}$ corresponding to an upper limit of $E_{\text{kin}}$. All regions reside below the 1-to-1 line (within the errors) if we assume an upper limit estimate of $n_e = n_e^{(\text{max})}$ instead (bottom-left panel).

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_e^{(\text{min})}$ (13% if counting only those with higher $E_{\text{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.

<footnote n="14">¹⁴ Where it was possible, we estimated the local gas mass fraction using the ALMA CO observations from Leroy et al. (2021). Accounting for molecular and ionised gas, we found the median value of ~0.23 implying that the additional contribution of gas to the total mass surface density will not affect significantly this relation.</footnote>
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Dependence of the kinetic energy of ionised gas $E_{\text{kin}}$ 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, $\Sigma_{E_{\text{kin}}}$, vs. corresponding stellar mass surface density, $\Sigma_{\star}$ (right panels). Top and bottom panels correspond to minimal and maximal values of $n_{\text{e}}$ underlying the $E_{\text{kin}}$ measurements. Colour encodes the $H\alpha$ flux surface brightness. The red line shows the median values of $E_{\text{kin}}$ or $\Sigma_{E_{\text{kin}}}$ in several bins, and the red shaded area corresponds to the 25–75 percentile interval. Spearman correlation coefficients $\rho$ are given above the plots.</figcaption></figure>

We note, however, that almost all regions reside below the 1-to-1 line if considering upper limit of the electron density ($n_{\text{e}}^{(\text{max})}$), but the correlation between $E_{\text{kin}}$ of gas and $E_{\text{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_{\text{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\alpha$ surface brightness (colour-coded) correlates with the stellar mass in the associations (see also Fig. 13). Given that the total $H\alpha$ 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\alpha$ luminosity (and thus the observed $H\alpha$ 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 HII regions (Scheuermann et al. 2023). Meanwhile, our results imply that this correlation is preserved also for the turbulent ionised gas around HII regions and for superbubbles, whose properties are not necessarily the same as those of normal HII regions (see Sect. 4.5). Given that the mass of the young stars and of the ionised gas (as well as $E_{\text{mech}}$ from the stars and $E_{\text{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: $\eta = E_{\text{kin}}(\text{gas})/E_{\text{mech}}(\text{stars})$. Our results imply that $\eta \sim 10–20\%$ (see Fig. 17), although there are regions with $\eta > 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_{\text{kin}}$ we
<figure><img src="image_18.png" /></figure>

**Fig. 16.** Dependence of the kinetic energy of ionised gas $E_{\text{kin}}$ 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_e^{(\text{min})}$ and $n_e^{(\text{max})}$ density measurements (thus representing upper and lower limits of $E_{\text{kin}}$, 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\text{II}]/H\alpha$ lines ratio tracing the dominant gas excitation mechanism ($[S\text{II}]/H\alpha > 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 $\eta = 100, 30, 5, 1\%$, respectively. Spearman correlation coefficient $\rho$ is given above the plots.

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_{\odot}$ is very low, and more data are required to confirm this finding.

### 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\alpha$ 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 – $\eta \sim 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 $\eta \sim 20–25\%$. We note that other simulations predict lower values – for example, Yadav et al. (2017) give $\eta \sim 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_e^{(\text{min})}$ generally corresponding to the upper limit of measured $\eta$. Using the upper limit on electron density $n_e^{(\text{max})}$ results in a much lower value of $\eta \sim 1\%$. We would need larger statistics and more precise measurements of the electron densities to obtain strict observational limits on the parameter $\eta$, 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
<figure><img src="image_19.png" /><figcaption>**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.</figcaption></figure>

and predicted stellar masses and ages inside the molecular gas superbubbles, that is close to our estimate.

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_{\text{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_{\text{eff}} < 3 \times \text{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<fnref n="15" /> are observed. In these regions, about 40% of the star clusters are offset from the centre of the superbubble by at least $0.5R_{\text{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\alpha$ 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

<figure><img src="image_20.png" /><figcaption>**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 $\eta = 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.</figcaption></figure>

the regions (at $R > 0.5R_{\text{eff}}$) than close to their centres ($R < 0.5R_{\text{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

<footnote n="15">¹⁵ Further results do not change if we consider higher (50 Myr) or lower (e.g., 8 Myr) age limits.</footnote>
<figure><img src="image_21.png" /><figcaption>**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.3R_{\text{eff}}$. Errors are calculated from 1000 Monte Carlo realisations with positions of each cluster randomly varying within the PSF of the MUSE observations.</figcaption></figure>

observational signatures of triggered star formation in super-
bubbles 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 morpho-
logically identified superbubbles (see, e.g., Barnes et al. 2023;
Watkins et al. 2023a, as examples of the first results).

**6. Summary**

Here we quantify the contribution of mechanical feedback from
young massive stars to the regulation of the kinematics and mor-
phology 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 observa-
tional 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 dis-
persion (σ(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 dis-
tribution 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_{\text{exp}} = 10-60 \text{ km s}^{-1}$ (mostly within $20-35 \text{ km s}^{-1}$) 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 sur-
rounding 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 exten-
sion of this analysis to dwarf galaxies will allow us to establish
a purely observational constraint on the efficiency of mechani-
cal stellar feedback in various environments, which is a critical
ingredient in cosmological simulations.

**Acknowledgements.** 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 50002206). RSK and SCOG also acknowledge computing resources provided by The Land 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

Agertz, O., Kravtsov, A. V., Leitner, S. N., & Gnedin, N. Y. 2013, ApJ, 770, 25

Allen, M. G., Groves, B. A., Dopita, M. A., Sutherland, R. S., & Kewley, L. J. 2008, ApJS, 178, 20

Alstott, J., Bullmore, E., & Plenz, D. 2014, PLoS ONE, 9, e85777

Ambrocio-Cruz, P., Le Coarer, E., Rosado, M., et al. 2016, MNRAS, 457, 2048

Anand, G. S., Lee, J. C., Van Dyk, S. D., et al. 2021, MNRAS, 501, 3621

Anderson, J. P., Gutiérrez, C. P., Dessart, L., et al. 2016, A&A, 589, A110

Bacchini, C., Fraternali, F., Iorio, G., & Pezzulli, G. 2019, A&A, 622, A64

Bacchini, C., Fraternali, F., Pezzulli, G., & Marasco, A. 2020, A&A, 644, A125

Bacon, R., Accardo, M., Adjali, L., et al. 2010, SPIE Conf. Ser., 7735, 773508

Bacon, R., Conseil, S., Mary, D., et al. 2017, A&A, 608, A1

Bagetakos, I., Brinks, E., Walter, F., et al. 2011, AJ, 141, 23

Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5

Barnes, A. T., Glover, S. C. O., Kreckel, K., et al. 2021, MNRAS, 508, 5362

Barnes, A. T., Chandar, R., Kreckel, K., et al. 2022, A&A, 662, L6

Barnes, A. T., Watkins, E. J., Meidt, S. E., et al. 2023, ApJ, 944, L22

Belfiore, F., Santoro, F., Groves, B., et al. 2022, A&A, 659, A26

Bittner, A., Falcón-Barroso, J., Nedelchev, B., et al. 2019, A&A, 628, A117

Blair, W. P., & Long, K. S. 2004, ApJS, 155, 101

Blair, W. P., Kirshner, R. P., & Chevalier, R. A. 1981, ApJ, 247, 879

Boquien, M., Burgarella, D., Roehlly, Y., et al. 2019, A&A, 622, A103

Bordalo, V., Plana, H., & Telles, E. 2009, ApJ, 696, 1668

Bresolin, F., Rizzi, L., Ho, I. T., et al. 2020, MNRAS, 495, 4347

Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000

Burkhart, B. 2021, PASP, 133, 102001

Calzetti, D., Lee, J. C., Sabbi, E., et al. 2015, AJ, 149, 51

Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245

Cedrés, B., Beckman, J. E., Bongiovanni, Á., et al. 2013, ApJ, 765, L24

Chabrier, G. 2003, ApJ, 586, L133

Chevance, M., Kruijssen, J. M. D., Krumholz, M. R., et al. 2022, MNRAS, 509, 272

Churchwell, E., Povich, M. S., Allen, D., et al. 2006, ApJ, 649, 759

Colman, T., Robitaille, J.-F., Hennebelle, P., et al. 2022, MNRAS, 514, 3670

Cosens, M., Wright, S. A., Murray, N., et al. 2022, ApJ, 929, 74

D'Agostino, J. J., Kewley, L. J., Groves, B. A., et al. 2019, MNRAS, 485, L38

Dale, J. E., Haworth, T. J., & Bressert, E. 2015, MNRAS, 450, 1199

Dawson, J. R., McClure-Griffiths, N. M., Wong, T., et al. 2013, ApJ, 763, 56

Della Bruna, L., Adamo, A., Bik, A., et al. 2020, A&A, 635, A134

Dib, S., Bell, E., & Burkert, A. 2006, ApJ, 638, 797

Egorov, O. V., Lozinskaya, T. A., Moiseev, A. V., & Smirnov-Pinchukov, G. V. 2014, MNRAS, 444, 376

Egorov, O. V., Lozinskaya, T. A., Moiseev, A. V., & Shchekinov, Y. A. 2017, MNRAS, 464, 1833

Egorov, O. V., Lozinskaya, T. A., Moiseev, A. V., & Smirnov-Pinchukov, G. V. 2018, MNRAS, 478, 3386

Egorov, O. V., Lozinskaya, T. A., Vasiliev, K. I., et al. 2021, MNRAS, 508, 2650

Egorova, E. S., Moiseev, A. V., & Egorov, O. V. 2019, MNRAS, 482, 3403

Ehlerová, S., & Palouš, J. 2005, A&A, 437, 101

Ejdetjärn, T., Agertz, O., Östlin, G., Renaud, F., & Romeo, A. B. 2022, MNRAS, 514, 480

Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146

Elmegreen, B. G., & Lada, C. J. 1977, ApJ, 214, 725

Elmegreen, B. G., & Calvo, J. 2004, ARA&A, 42, 211

Elmegreen, B. G., Martinez, Z., & Hunter, D. A. 2022, ApJ, 928, 143

Emsellem, E., Schinnerer, E., Santoro, F., et al. 2022, A&A, 659, A191

Fisher, D. B., Bolatto, A. D., White, H., et al. 2019, ApJ, 870, 46

Forbes, J. C., Emami, R., Somerville, R. S., et al. 2023, ApJ, 948, 107

Freedman, W. L., Madore, B. F., Gibson, B. K., et al. 2001, ApJ, 553, 47

Fujii, K., Mizuno, N., Dawson, J. R., et al. 2021, MNRAS, 505, 459

Georgy, C., Ekström, S., Eggenberger, P., et al. 2013, A&A, 558, A103

Gerasimov, I. S., Egorov, O. V., Lozinskaya, T. A., Moiseev, A. V., & Oparin, D. V. 2022, MNRAS, 517, 4968

Ginzburg, O., Dekel, A., Mandelker, N., & Krumholz, M. R. 2022, MNRAS, 513, 6177

Girard, M., Fisher, D. B., Bolatto, A. D., et al. 2021, ApJ, 909, 12

Girichidis, P., Offner, S. S. R., Kritsuk, A. G., et al. 2020, Space Sci. Rev., 216, 68

Green, A. W., Glazebrook, K., McGregor, P. J., et al. 2014, MNRAS, 437, 1070

Groves, B., Kreckel, K., Santoro, F., et al. 2023, MNRAS, 520, 4902

Haffner, L. M., Dettmar, R. J., Beckman, J. E., et al. 2009, Rev. Mod. Phys., 81, 969

Haid, S., Walch, S., Seifried, D., et al. 2018, MNRAS, 478, 4799

Hartmann, L., Ballesteros-Paredes, J., & Bergin, E. A. 2001, ApJ, 562, 852

Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288

Hopkins, P. F., Quataert, E., & Murray, N. 2011, MNRAS, 417, 950

Hunt, L. K., & Hirashita, H. 2009, A&A, 507, 1327

Hunter, D. A., Elmegreen, B. G., Archer, H., Simpson, C. E., & Cigan, P. 2021, AJ, 161, 175

Hunter, L. C., van Zee, L., McQuinn, K. B. W., Garner, R., & Dolphin, A. E. 2022, AJ, 163, 132

Inutsuka, S.-I., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49

Jacobs, B. A., Rizzi, L., Tully, R. B., et al. 2009, ApJ, 138, 332

Kauffmann, G., Heckman, T. M., Tremonti, C., et al. 2003, MNRAS, 346, 1055

Keller, B. W., Kruijssen, J. M. D., & Chevance, M. 2022, MNRAS, 514, 5355

Kennicutt, R. C., & J., 1984, ApJ, 287, 116

Kewley, L. J., Dopita, M. A., Sutherland, R. S., Heisler, C. A., & Trevena, J. 2001, ApJ, 556, 121

Kewley, L. J., Geller, M. J., & Barton, E. J. 2006, AJ, 131, 2004

Kewley, L. J., Nicholls, D. C., Sutherland, R., et al. 2019, ApJ, 880, 16

Klessen, R. S., & Glover, S. C. O. 2016, in Saas-Fee Advanced Course, eds. Y. Revaz, P. Jablonka, R. Teyssier, & L. Mayer, 43, 85

Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17

Koch, E. W., Rosolowsky, E. W., Lockman, F. J., et al. 2018, MNRAS, 479, 2505

Kourkchi, E., & Tully, R. B. 2017, ApJ, 843, 16

Kourkchi, E., Tully, R. B., Eftekharzadeh, S., et al. 2020, ApJ, 902, 145

Krause, M. G. H., & Diehl, R. 2014, ApJ, 794, L21

Kreckel, K., Egorov, O. V., Belfiore, F., et al. 2022, A&A, 667, A16

Kroupa, P. 2001, MNRAS, 322, 231

Krumholz, M. R., & Burkhart, B. 2016, MNRAS, 458, 1671

Krumholz, M. R., Bate, M. R., Arce, H. G., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 243

Lancaster, L., Ostriker, E. C., Kim, J.-G., & Kim, C.-G. 2021, ApJ, 914, 90

Larson, K. L., Lee, J. C., Thilker, D. A., et al. 2023, MNRAS, 523, 6061

Law, D. R., Westfall, K. B., Bershady, M. A., et al. 2021a, AJ, 161, 52

Law, D. R., Ji, X., Belfiore, F., et al. 2021b, ApJ, 915, 35

Law, D. R., Belfiore, F., Bershady, M. A., et al. 2022, ApJ, 928, 58

Lee, J. C., Whitmore, B. C., Thilker, D. A., et al. 2022, ApJS, 258, 10

Lee, J. C., Sandstrom, K. M., Leroy, A. K., et al. 2023, ApJ, 944, L17

Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3

Leitherer, C., Ekström, S., Meynet, G., et al. 2014, ApJS, 212, 14

Leroy, A. K., Schinnerer, E., Hughes, A., et al. 2021, ApJS, 257, 43

Levy, R. C., Bolatto, A. D., Sánchez, S. F., et al. 2019, ApJ, 882, 84

López-Cobá, C., Sánchez, S. F., Anderson, J. P., et al. 2020, AJ, 159, 167

Lozinskaya, T. A. 2002, Astron. Astrophys. Trans., 21, 223

Lozinskaya, T. A., Moiseev, A. V., & Podorvanyuk, N. Y. 2003, Astron. Lett., 29, 77

Luridiana, V., Morisset, C., & Shaw, R. A. 2015, A&A, 573, A42

Mac Low, M.-M., & McCray, R. 1988, ApJ, 324, 776

Mac Low, M.-M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Makarov, D., Prugniel, P., Terekhova, N., Courtois, H., & Vauglin, I. 2014, A&A, 570, A13

Martínez-Delgado, I., Tenorio-Tagle, G., Muñoz-Tuñón, C., Moiseev, A. V., & Cairós, L. M. 2007, AJ, 133, 2892

Martins, F., Schaerer, D., & Hillier, D. J. 2005, A&A, 436, 1049

Melnick, J., Tenorio-Tagle, G., & Terlevich, R. 1999, MNRAS, 302, 677

Moiseev, A. V. 2015, Astrophys. Bull., 70, 494

Moiseev, A. V., & Lozinskaya, T. A. 2012, MNRAS, 423, 1831

Moiseev, A. V., Tikhonov, A. V., & Klypin, A. 2015, MNRAS, 449, 3568

Muñoz-Tuñon, C., Tenorio-Tagle, G., Castañeda, H. O., & Terlevich, R. 1996, AJ, 112, 1636

Nath, B. B., Das, P., & Oey, M. S. 2020, MNRAS, 493, 1034

Nazé, Y., Chu, Y.-H., Points, S. D., et al. 2011, AJ, 122, 921

Nelson, D., Pillepich, A., Genel, S., et al. 2015, Astron. Comput., 13, 12

Ntormousi, E., Burkert, A., Fierlinger, K., & Heitsch, F. 2011, ApJ, 731, 13
Nugent, P., Sullivan, M., Ellis, R., et al. 2006, ApJ, 645, 841

Nusser, A., & Silk, J. 2022, MNRAS, 509, 2979

Oey, M. S., & Clarke, C. J. 1997, MNRAS, 289, 570

Oparin, D. V., & Moiseev, A. V. 2018, Astrophys. Bull., 73, 298

Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei, 2nd edn. (University Science Books)

Ostriker, E. C., Stone, J. M., & Gammie, C. F. 2001, ApJ, 546, 980

Peimbert, M. 1971, Boletin de los Observatorios Tonantzintla y Tacubaya, 6, 29

Pengelly, R. M. 1964, MNRAS, 127, 145

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Pilyugin, L. S., & Grebel, E. K. 2016, MNRAS, 457, 3678

Pilyugin, L. S., Cedrés, B., Zinchenko, I. A., et al. 2021, A&A, 653, A11

Pokhrel, N. R., Simpson, C. E., & Bagetakos, I. 2020, AJ, 160, 66

Querejeta, M., Schinnerer, E., Meidt, S., et al. 2021, A&A, 656, A133

Reynolds, R. J. 1977, ApJ, 216, 433

Robitaille, T., Rice, T., Beaumont, C., et al. 2019, Astrophysics Source Code Library [record ascl:1907.016]

Rubin, R. H. 1989, ApJS, 69, 897

Sarbadhicary, S. K., Martizzi, D., Ramirez-Ruiz, E., et al. 2022, ApJ, 928, 54

Scheuermann, F., Kreckel, K., Anand, G. S., et al. 2022, MNRAS, 511, 6087

Scheuermann, F., Kreckel, K., Barnes, A. T., et al. 2023, MNRAS, 522, 2369

Sharda, P., Ginzburg, O., Krumholz, M. R., et al. 2023, MNRAS, submitted, [arXiv:2303.15853]

Sharma, P., Roy, A., Nath, B. B., & Shchekinov, Y. 2014, MNRAS, 443, 3463

Shaya, E. J., Tully, R. B., Hoffman, Y., & Pomarède, D. 2017, ApJ, 850, 207

Sil'chenko, O. K., Moiseev, A. V., & Egorov, O. V. 2019, ApJS, 244, 6

Smirnov-Pinchukov, G. V., & Egorov, O. V. 2021, Astrophys. Bull., 76, 367

Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41

Stuber, S. K., Saito, T., Schinnerer, E., et al. 2021, A&A, 653, A172

Tamburro, D., Rix, H. W., Leroy, A. K., et al. 2009, AJ, 137, 4424

Tenorio-Tagle, G., Munoz-Tunon, C., & Cid-Fernandes, R. 1996, ApJ, 456, 264

Terlevich, R., & Melnick, J. 1981, MNRAS, 195, 839

Thilker, D. A., Braun, R., & Walterbos, R. A. M. 2000, AJ, 120, 3070

Thilker, D. A., Whitmore, B. C., Lee, J. C., et al. 2022, MNRAS, 509, 4094

Turner, J. A., Dale, D. A., Lee, J. C., et al. 2021, MNRAS, 502, 1366

Utomo, D., Blitz, L., & Falgarone, E. 2019, ApJ, 871, 17

Valdez-Gutiérrez, M., Rosado, M., Georgiev, L., Borissova, J., & Kurtev, R. 2001, A&A, 366, 35

Vasiliev, E. O., Moiseev, A. V., & Shchekinov, Y. A. 2015, Balt. Astron., 24, 213

Vasiliev, E. O., Shchekinov, Y. A., Koval', V. V., & Egorov, O. V. 2022, Astrophys. Bull., 77, 132

Viegas, S. M., & Clegg, R. E. S. 1994, MNRAS, 271, 993

Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574

Walch, S., & Naab, T. 2015, MNRAS, 451, 2757

Watkins, E. J., Barnes, A. T., Henny, K., et al. 2023a, ApJ, 944, L24

Watkins, E. J., Kreckel, K., Groves, B., et al. 2023b, A&A, 676, A67

Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377

Whitmore, B. C., Lee, J. C., Chandar, R., et al. 2021, MNRAS, 506, 5294

Williams, T. G., Kreckel, K., Belfiore, F., et al. 2022, MNRAS, 509, 1303

Xu, D., Offner, S. S. R., Guttermuth, R., & Oort, C. V. 2020, ApJ, 890, 64

Yadav, N., Mukherjee, D., Sharma, P., & Nath, B. B. 2017, MNRAS, 465, 1720

Yarovova, A. D., Egorov, O. V., Moiseev, A. V., & Maryeva, O. V. 2023, MNRAS, 518, 2256

Yu, X., Shi, Y., Chen, Y., et al. 2019, MNRAS, 486, 4463

¹ Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstr. 12-14, 69120 Heidelberg, Germany

e-mail: oleg.egorov@uni-heidelberg.de

² Institut für Theoretische Astrophysik, Zentrum für Astronomie der Universität Heidelberg, Albert-Ueberle-Strasse 2, 69120 Heidelberg, Germany

³ International Centre for Radio Astronomy Research, University of Western Australia, 7 Fairway, Crawley 6009 WA, Australia

⁴ INAF – Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50157 Firenze, Italy

⁵ European Southern Observatory, Karl-Schwarzschild Straße 2, 85748 Garching bei München, Germany

⁶ Univ. Lyon, Univ. Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France

⁷ Interdisziplinäres Zentrum für Wissenschaftliches Rechnen der Universität Heidelberg, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany

⁸ Department of Astronomy, The Ohio State University, 140 West 18th Avenue, Columbus, Ohio 43210, USA

⁹ Center for Cosmology and Astroparticle Physics, 191 West Woodruff Avenue, Columbus, OH 43210, USA

¹⁰ Sterrenkundig Observatorium, Universiteit Gent, Krijgslaan 281 S9, 9000 Gent, Belgium

¹¹ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

¹² Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹³ Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany

¹⁴ European Southern Observatory (ESO), Alonso de Córdova 3107, Casilla 19, Santiago 19001, Chile

¹⁵ Department of Physics & Astronomy, University of Wyoming, Laramie, WY 82071, USA

¹⁶ Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia

¹⁷ AURA for the European Space Agency (ESA), Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

¹⁸ Gemini Observatory/NSF's NOIRLab, 950 N. Cherry Avenue, Tucson, AZ, USA

¹⁹ Steward Observatory, University of Arizona, 933 N Cherry Ave, Tucson, AZ 85721, USA

²⁰ Department of Physics & Astronomy, Bloomberg Center for Physics and Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA

²¹ Sub-department of Astrophysics, Department of Physics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
# 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).

A large scatter in $\sigma(H\alpha)$ is usually observed at the low surface brightness end of the I – $\sigma$ 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 $\sigma(H\alpha)$ 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, $\Sigma(H\alpha)_{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, $\Sigma(H\alpha)_{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 $\Sigma(H\alpha)_{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 $\Sigma(H\alpha)_{DIG}$. However, visually those taken from Belfiore et al. (2022) are better at separating the low brightness peak of $\sigma(H\alpha)$ on I – $\sigma$ diagrams for the rest of the objects. Therefore, we used the values from Belfiore et al. (2022) to determine $\Sigma(H\alpha)_{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 $\Sigma(H\alpha)_{DIG}$ to better match the low-brightness peak in the $\sigma(H\alpha)$ 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 $\Sigma(H\alpha)_{DIG}$ for each galaxy are given in Table 1.

We estimated the mean $H\alpha$ velocity dispersion $\sigma(H\alpha)_{m0}$ for the H II regions in each galaxy as the flux-weighted mean value across all pixels brighter than $\Sigma(H\alpha)_{DIG}$. The uncertainty of the mean velocity dispersion $\Delta(\sigma_m)_{0}$ was measured as the standard deviation of $\sigma(H\alpha)$ around $\sigma(H\alpha)_{m0}$. The final adopted values of $\sigma(H\alpha)_{m}$ and $\Delta(\sigma_m)$ were obtained on the second iteration, where we excluded all pixels with $\sigma(H\alpha) > \sigma(H\alpha)_{m0} + \Delta(\sigma_m)_{0}$. The derived values of $\sigma(H\alpha)_{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 $\sigma(H\alpha)_{m}$ and SFR as reported in other works (e.g. Law et al. 2022).

Based on $\sigma(H\alpha)_{m}$ and $\Delta(\sigma_m)$, we performed a classification of all pixels brighter than $\Sigma(H\alpha)_{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 – $\sigma$ diagrams. For brighter

regions, a scatter of $\sigma(H\alpha)$ 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 $\sigma(H\alpha)$ vs $\log(I[\text{H}\alpha])$ for such regions can be parameterised by a Gaussian with the mean value and standard deviation equal to $\sigma(H\alpha)_{m}$ and $\Delta(\sigma_m)$, respectively. We refer to velocity dispersion defined by such a normal distribution for each flux value as $\sigma(H\alpha)_{\text{norm}}$.

We classify a pixel as related to Class 1 if $\sigma(H\alpha) < \sigma(H\alpha)_{\text{norm}}$ for its $H\alpha$ brightness. Note that all pixels with $\sigma(H\alpha) < \sigma(H\alpha)_{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 $\sigma(H\alpha)$ between two distributions – $\sigma(H\alpha)_{\text{norm}}$ and $\sigma(H\alpha)_{1.5\times\text{norm}}$. The latter is defined in the same way as $\sigma(H\alpha)_{\text{norm}}$, but with a standard deviation equal to $1.5\Delta(\sigma_m)$. This is an intermediate class: the elevated $\sigma(H\alpha)$ 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 $\sigma(H\alpha) > \sigma(H\alpha)_{1.5\times\text{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 – $\sigma$ diagrams.

**Top σ(Hα) and SNR/WR candidates.** This is a sub-group of the Class 3 pixels selected as having intrinsic $\sigma(H\alpha) > 100 \text{ km s}^{-1}$, or $\sigma(H\alpha) > \sigma(H\alpha)_{3\times\text{norm}}$ (where $\sigma(H\alpha)_{3\times\text{norm}}$ is the normally distributed value with a standard deviation equal to $3\Delta(\sigma_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[\text{H}\alpha])$ 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 $\sigma(H\alpha)$.

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<fnref n="16" /> 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 $\sigma(H\alpha)$ 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 $\sigma(H\alpha)$ of all pixels within the galaxy classified as Class 2 or

<footnote n="16">¹⁶ http://www.dendrograms.org/</footnote>
<table><caption>**Table A.1.** Summary of the criteria applied for classification and segmentation of the pixels in I − σ diagrams.</caption><thead><tr><th>Classification</th><th>Observed σ(Hα) range</th><th>Clipped σ(Hα) for ASTRODENDRO</th><th>Physical interpretation</th></tr></thead><tbody><tr><td>Low Σ(Hα)</td><td>Any; I(Hα) &lt; Σ(Hα)<sub>DIG</sub></td><td>σ<sub>ref</sub> − Δ(σ<sub>m</sub>)</td><td>Diffuse ionised gas</td></tr><tr><td>Class 1</td><td>σ(Hα) &lt; σ(Hα)<sub>norm</sub></td><td>σ<sub>ref</sub> − Δ(σ<sub>m</sub>)</td><td>Unperturbed H II regions</td></tr><tr><td>Class 2</td><td>[σ(Hα)<sub>norm</sub>, σ(Hα)<sub>1.5×norm</sub>]</td><td>σ<sub>ref</sub></td><td>May be illusive, but likely periphery of superbubbles; low-velocity shocks</td></tr><tr><td>Class 3</td><td>[σ(Hα)<sub>1.5×norm</sub>, σ(Hα)<sub>3×norm</sub>]</td><td>max[σ<sub>ref</sub> + 2Δ(σ<sub>m</sub>), σ(Hα)]</td><td>Expanding superbubbles; shocks; other supersonic gas motions</td></tr><tr><td>Top σ</td><td>σ(Hα) ≥ min[σ(Hα)<sub>3×norm</sub>, 100 km s<sup>-1</sup>]</td><td>max[σ<sub>ref</sub> + 3Δ(σ<sub>m</sub>), σ(Hα)]</td><td>SNR/WR or other extreme source candidates</td></tr></tbody></table>

**Notes.** − σ(Hα)norm, σ(Hα)1.5×norm and σ(Hα)3×norm are the upper boundaries of velocity dispersion defined by Gaussians in σ(Hα) vs log I(Hα) space with mean = σ(Hα)m and std = Δ(σm), 1.5Δ(σm) and 3Δ(σm), respectively. − σ(Hα)m and Δ(σm) are flux-weighted mean velocity dispersion and its standard deviation (see measured values for each galaxy in Table 1). − σref is the median σ(Hα) derived from all pixels classified as Class 2, 3 or Top σ.

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ₐ, where FWHMₐ = 2√2 ln 2σₐ, and σₐ 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 ΔV ∝ √(1 − r²/R²), 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) σₐ ≈ 0.32 × R for various expansion velocities and spatial resolutions<fnref n="17" />. 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.

<footnote n="17">¹⁷ We used SDSS5-LVMDATASIMULATOR package (https://github.com/sdss/lvmdatasimulator) for doing this analysis as it allows one to easily construct such toy models with various parameters and numbers of spatial resolution elements across the nebula and obtain the simulated IFU spectra for them.</footnote>
<figure><img src="image_22.png" /><figcaption>**Fig. A.1.** *l* − *σ* 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.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. A.1. Continued.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. A.1. Continued.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. A.1. Continued.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. A.1. Continued.</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>Fig. A.1. Continued.</figcaption></figure>