© The Authors 2025

1. Introduction

In this paper our aim is to study two aspects of galaxy structure. These are the gas inhomogeneities that cause dynamical friction in gaseous discs of galaxies and the corresponding secular evolution of these gas discs. Specifically, in this paper we focus on superbubbles (SBs) created by supernova (SN) explosions in gas discs as the source of dynamical friction.

Superbubbles are cavities in the gas distribution; they span from hundreds to thousands of parsecs. Nath et al. (2020) compared the size distributions of observed HI holes from THINGS (Walter et al. 2008) with those derived from detailed hydrodynamical simulations, and concluded that SBs are caused by multiple SNe in a large OB association or in multiple associations. Many observations document the phenomenon as dark regions in the gas discs of galaxies (e.g. Lara-López et al. 2023). Recent studies using JWST observations give excellent samples of SBs in NGC 628 (Watkins et al. 2023). In addition, studies such as Brinks et al. (2007), Walter et al. (2008) provide good observations of holes in the HI distribution in NGC 6946 (see Fig. 1). The association of HI holes and SN explosions has been studied in detail (e.g. Lara-López et al. 2023; Sarbadhicary et al. 2023; Barnes et al. 2023).

Fig. 1. Illustration of SBs in galaxy NGC 6946 from the intensity map of the THINGS survey (Walter et al. 2008). The inset in the top right corner is a zoomed-in image of a fraction of the galaxy where the locations of three SBs are shown (pink circles). A more detailed map of all the bubbles in the galaxy is shown in Brinks et al. (2007).

From a physical perspective, SBs are underdense regions within the gaseous medium, filled by hot pressurised gas. This is created by multiple SN explosions that inflate a cavity within the ambient medium, driving a shock wave into the latter and accumulating gas at its edges (Mac Low & McCray 1988). When the SN energy input ceases, the SB starts to dissolve. The hot interior exchanges the heat with edges via thermal conduction (Yadav et al. 2017, hereafter Y17), causing heat to diffuse towards the cold shell and the shell material to evaporate towards the interior of the bubble (Cowie & McKee 1977). During the two phases, the shell of the SB is either expanding or shrinking, causing disturbances in density and gravitational potential. An overview of these processes is provided by Keller (2017). Physically, the situation is rather complex, as the heat conduction may depend on magnetic fields, and the initial gas distribution may be clumpy instead of homogeneous. These aspects may be relevant and were studied in Keller et al. (2014).

Our aim was to expand the study of SBs and their gravitational interaction with the galaxy disc, as the current status of this research is not extensive. García-Conde et al. (2024), while not aiming to solve the SBs influence to the structure of galaxies, notes that in low-density environments, and in connection to bending waves, SNe can redistribute gas and cause acceleration field changes. In the context of dwarf and high-redshift clumpy galaxies, the strong scattering of giant clouds and holes may help to form exponential-like discs (Struck & Elmegreen 2017; Wu et al. 2020). However, these works concentrated on turbulent rather than quiet disc galaxies.

As a rule, the gas disc and the stellar disc of a galaxy rotate at rather different velocities: the stellar disc has a slower rotation. Therefore, there is, on average, a systematic mismatch between the velocities of the stars and the SBs, which can give rise to a net momentum and/or energy transfer between the stars traversing the SB and the gaseous environment. In the case of a massive perturber in a stellar environment, this is called dynamical friction (Chandrasekhar 1943). In Kipper et al. (2020, 2023) we developed a method to calculate dynamical friction effects in a non-homogeneous gravitational potential. For the present paper we used the method for the case of SBs. For this we needed to know the density contrast (positive or negative) of a SB relative to the host environment.

Dynamical friction acts in two ways: it changes the velocities of stars and the velocities of SBs, both of which result in the evolution of the disc. These changes are slow (secular). There are various other processes causing the secular evolution of galaxies (Kuzmin 1961, 1963; Binney 2013; Kipper et al. 2021; Huško et al. 2023) (also see the review by Kormendy 2013). Getting a grasp on these processes is necessary in order to reach a level of understanding to reproduce all galaxy aspects in simulations. A good understanding of the baryonic effects also allows dark matter (DM) studies to be improved (Maxwell et al. 2015; Benito et al. 2020). The list of internal processes that contribute to secular evolution includes the scattering of stars, characterised by the Fokker-Planck equation (see Binney & Tremaine 2008); churning, which is the change in angular momentum by spiral, bar, or other non-axisymmetry (Sellwood & Binney 2002; Minchev & Famaey 2010; Vera-Ciro et al. 2014; Patil et al. 2023); blurring, where the breadth of an orbit grows (Schönrich & Binney 2009) out to very large galactic radii (Lian et al. 2022); and local disc instabilities (Huško et al. 2023). These processes provide a significant redistribution of matter in galaxies. In the present paper we study an addition to these effects that might be the reason behind a slow redistribution of stars and gas. Similarly to giant molecular clouds acting as scatterers of stars, we investigate how much the SBs alter the gas distribution.

The purpose of the current paper is to investigate whether the time-dependent gravitational potential disturbances help the transfer of angular momentum between the gas disc and stellar disc in a smooth and relatively unperturbed (regular) galaxy. Thereafter, we consider whether these contributions are sufficient to cause a significant migration of the gas disc.

In Sect. 2 we characterise our description of the SB and its evolution. In Sect. 3 we show how passage through the SB influences the single-star kinematics numerically, and we provide analytical approximations for it. In Sect. 4 we show numerically the transfer of momentum or angular momentum in various environments. We discuss the results in Sect. 5, and summarise in Sect. 6.

2. Description of a superbubble

Although SBs may have a complex morphology (Kim et al. 2017), in essence they are underdense regions embedded within a background gaseous environment (the host). In our present treatment we consider this background to be homogeneous, and describe the SB as a spherical region with (formally) negative interior density relative to the background, surrounded by a slightly overdense shell of material displaced from the interior. We subsequently consider the effect of this configuration on the local acceleration field relative to that of the unperturbed host.

2.1. The density profile

We approximated the SB density profile at any moment in time using two parameters: its inner and outer radius, R_i and R_o, respectively. The inner radius is the radius closer to the centre of the cavity before the overdensity region, and the outer radius the one closer to the interstellar medium (ISM) at the edges of the SB. For convenience, we also denote their midpoint R_m = (R_i + R_o)/2. The approximation of the SB mimics the profile from Y17, where they simulated the evolution of a SB in a homogeneous gas environment with initial density ρ_env = 0.015 M_⊙ pc⁻³. The profile itself is shown in Fig. A.1 of Y17. Our approximation assumes that the SNe have pushed out all the gas inside R_i. Hence, the physical density inside it is zero, but in the form of density-correction ρ^* ≡ ρ − ρ_env its value is −ρ_env. In our description the gas from inside is pushed between R_i and R_o, where we describe the density distribution as a quadratic polynomial. Outside of R_o, the SB has had no influence and the density-correction is zero. Combining all the regions together, we have the formula

$$\begin{array}{c}\hfill {\rho }^{\ast }\equiv \rho -{\rho }_{\mathrm{env}}=\{\begin{array}{cc}-{\rho }_{\mathrm{env}},\hfill & R<R_{\mathrm{i}}\hfill \\ B-A{(R-R_{\mathrm{m}})}^2,\hfill & R_{\mathrm{i}}\le R\le R_{\mathrm{o}}\hfill \\ 0,\hfill & R>R_{\mathrm{o}}\hfill \end{array}.\end{array}$$(1)

The coefficients A and B are found from mass conservation: the mass shifted from within radius R_i equals that added between R_i and R_o and at ρ^*(R_i) = ρ^*(R_o) = 0. Figure 2 shows this profile visually.

Fig. 2. Density contrast profile approximation of the SB at t = 30 Myr inspired by the simulated SB behaviour in Y17 and described in Sect. 2.1.

2.2. The evolution of superbubbles

Superbubbles are short-lived compared to galaxies; they vanish soon after the SNe of the star-forming region cease. We split the description of the evolution into two phases: the expanding phase, where the SNe are pushing the gas out, and the vanishing phase, where the SNe support ceases, and the SBs start to shrink. We call this phase the collapsing phase to emphasise our assumption. We denote the time of the split between the phases as t_max.

The precise evolution of the SBs is uncertain. The formation of the SBs is studied in many magnetohydrodynamic simulations (Kim & Ostriker 2017; Girichidis et al. 2016), but their complete evolution depends not only on the SN explosions, but also on the physics responsible for the collapse of the SB cavity. In the literature, the first stages of SB evolution have been worked through thoroughly. In the case of a single SN explosion, the main physical processes have been analysed (see e.g. Ryden & Pogge 2021), but in the case where many SNe explode in spatial and temporal vicinity, it becomes more complicated, and these solutions have only been simulated (e.g. Kim et al. 2017). In the collapsing phase, turbulences and shearing motions redistribute former shell material back into the SB volume (Kim & Ostriker 2017). Turbulences are chaotic phenomena that do not follow strict descriptions, which makes analytical approximation somewhat difficult. From the analytical side, the collapsing stage in atmospheric conditions was studied by Zel’dovich & Raizer (1967), among others, who showed that after the shockwave travels through the atmosphere, the outer edge of the shell stalls and inner edge turns into a cooling wave that moves the interiors of the cavity. However, this does not necessarily align with what is happening in the case of SBs, as the cooling is different for atmospheric chemistry and the ISM.

The Y17 simulation describes only the first 30 Myr of the SB evolution as this is the approximate age when the least massive core-collapse SNe explode. Thus, we denote this as the expansion phase. We approximated the evolution of both R_i and R_o as a power law with two and three free parameters, respectively. For the vanishing of the SB phase, which was not covered by the Y17 simulation, we assumed a scenario describable with a characteristic speed, which we treated as constant. It mimics a collapse occurring for both R_i and R_o with a constant speed v_col. We denote the lifetime of the SB as t_end, defined as the time when the R_i shrinks back to 0 kpc. The analytical form for R_i is

$$\begin{array}{c}\hfill R_{\mathrm{i}}^{\mathrm{fid}}=\{\begin{array}{ccc}{R}_{\mathrm{u}}{(t/t_{\mathrm{i},\mathrm{e}})}^{{\alpha }_{\mathrm{i},\mathrm{e}}},\hfill & \hfill & t\le t_{\max }\hfill \\ R_{\max ,\mathrm{i}}-v_{\mathrm{col}}(t-t_{\max ,\mathrm{i}})\hfill & t_{\max }<\hfill & t\le t_{\mathrm{end}}\hfill \\ 0,\hfill & \hfill & t>t_{\mathrm{end}}\hfill \end{array}\end{array}$$(2)

and for R_o

$$\begin{array}{c}\hfill R_{\mathrm{o}}^{\mathrm{fid}}=\{\begin{array}{cc}{R}_{\mathrm{o}0}+R_{\mathrm{u}}{(t/t_{\mathrm{o},\mathrm{s}})}^{{\alpha }_{\mathrm{o}}},\hfill & t\le t_{\max }\hfill \\ R_{\max ,\mathrm{o}}-v_{\mathrm{col}}(t-t_{\max ,\mathrm{i}})\hfill & t>t_{\max }\hfill \end{array}.\end{array}$$(3)

We approximated the parameters in this formula to match the evolution provided in Y17, and derived the following values: t_i, e = 290 Myr, α_i, e = 0.47, R_o0 = 60 pc, t_o, s = 56.7 Myr, α_o = 0.815. In addition, we denoted R_max, i and R_max, o as the SB inner and outer size at t_max; v_col as the collapsing speed; and the superscript ‘fid’ as the fiducial run of the simulation. R_u is a constant with unit-length for dimensional considerations. The evolution of $R_{\mathrm{i}}^{\mathrm{fid}}$ and $R_{\mathrm{o}}^{\mathrm{fid}}$ is shown in Fig. 3 as dashed and solid red lines, respectively.

Fig. 3. Acceleration field of SB. The solid and dashed red lines are the outer and inner radius of the fiducial model from Y17, and are described with Eqs. (8) and (9). The colour scale shows the radial acceleration values described in Sect. 3.

Our currently used approximation of the collapsing phase is when the inner and outer parts of the SB collapse with the same speed. This is not necessarily the case. In order to see how different assumptions change, we analysed our results (see Appendix A.2) and we tested the cases where (i) R_i and R_o collapse at different speeds; (ii) the SB does not collapse, but mixes with turbulence characterised by modulation of the amplitude of acceleration; and (iii) SB disappears due to shearing motions caused by differential rotation of the galaxy.

Although Y17 concentrates on their fiducial model, they provided some scaling relations for different SBs. In the current case, the relevant relation is the scaling of the SB radii depending on the number of progenitor SNe. The scaling behaves as the modifier of radii: R_i and R_o are found from fiducial values by multiplying them with the term s_r, or equivalently

$$\begin{array}{c}\hfill s_{\mathrm{r}}=\frac{R_{\mathrm{i}}}{R_{\mathrm{i}}^{\mathrm{fid}}}=\frac{R_{\mathrm{o}}}{R_{\mathrm{o}}^{\mathrm{fid}}}·\end{array}$$(4)

Y17 shows that s_r depends on the number progenitor SNe as $N_{\mathrm{SN}}^{1/5}$. The fiducial model corresponded to 100 SNe, and hence

$$\begin{array}{c}\hfill s_{\mathrm{r}}={\left(\frac{N_{\mathrm{SN}}}{100}\right)}^{1/5}.\end{array}$$(5)

To match the number of SNe with the mass of their progenitor open clusters we use the initial mass function (IMF) Ψ dm. We match the IMF normalisation to the mass of the SB progenitor open cluster (M_OC) by introducing a normalising constant n_IMF so that

$$\begin{array}{c}\hfill {\displaystyle \underset{0}{\overset{\infty }{\int }}}m{n}_{\mathrm{IMF}}\mathrm{\Psi }(m)\mathrm{d}m=M_{\mathrm{OC}}.\end{array}$$(6)

Secondly, we are interested in the number of stars that can produce core-collapse SNe. In terms of IMF it is an integral between the suitable mass range (m_SN, min…m_SN, max):

$$\begin{array}{c}\hfill N_{\mathrm{SN}}={\displaystyle \underset{m_{\mathrm{SN},\min }}{\overset{m_{\mathrm{SN},\max }}{\int }}}{n}_{\mathrm{IMF}}\mathrm{\Psi }(m)\mathrm{d}m=M_{\mathrm{OC}}\frac{{\displaystyle \underset{m_{\mathrm{SN},\min }}{\overset{m_{\mathrm{SN},\max }}{\int }}}\mathrm{\Psi }(m)\mathrm{d}m}{{\displaystyle \underset{0}{\overset{\infty }{\int }}}m\mathrm{\Psi }(m)\mathrm{d}m}·\end{array}$$(7)

The division of integrals provides the number of SNe per unit cluster mass, which we denote as ϑ so that N_SN = ϑM_OC.

For the numerical evaluation of ϑ we need to specify Ψ, m_SN, min, and m_SN, max. For the Ψ we selected the Kroupa IMF (Kroupa 2001), the m_SN, min = 8 M_⊙, but the m_SN, max is not well determined. There have been dimming periods observed of massive stars, and hints of faint or absent SNe with masses above m_SN, max = 17 M_⊙ (Reynolds et al. 2015; Byrne & Fraser 2022), which shows that there is a lack of SNe with masses above m_SN, max = 17 M_⊙, although this has not yet been confirmed. Therefore, it is possible that the highest mass stars do not explode as SNe and we should adopt m_SN, max = 17 M_⊙, but it is unclear if this affects all massive stars or only some fraction of them. The conventional text-book value is still m_SN, max = ∞ M_⊙ for integrations. In these cases, the values are ϑ₁₇ = 0.0054 M_⊙⁻¹ and ϑ_∞ = 0.0084 M_⊙⁻¹. In the current paper we adopt the more classical ϑ = ϑ_∞.

Combining all of this together, we see that for the cluster mass of M_OC the inner and outer radii are determined from

$$\begin{array}{cc}\hfill R_{\mathrm{i}}& =R_{\mathrm{i}}^{\mathrm{fid}}{(\vartheta M_{\mathrm{OC}}/100)}^{1/5},\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill R_{\mathrm{o}}& =R_{\mathrm{o}}^{\mathrm{fid}}{(\vartheta M_{\mathrm{OC}}/100)}^{1/5}.\hfill \end{array}$$(9)

2.3. Selecting open clusters able to produce SBs

The inflation of a SB is not guaranteed in all open clusters as not all them are able to produce a sufficient number of SNe for that. According to Y17, the number of SNe necessary for the production of SBs is N_SN, lim ≈ 100. The mass of the open cluster corresponding to this is

$$\begin{array}{c}\hfill M_{\mathrm{OC},\lim }\equiv N_{\mathrm{SN},\lim }/\vartheta \approx 11800{\mathrm{M}}_{\odot }.\end{array}$$(10)

We now turn our attention to the fraction of open clusters able to produce a SB and what their average mass is and the corresponding parameter s_r defined in Eq. (5). We take the open cluster IMF (≡Ψ_OC) from the Milky Way (MW) observations. The MW is in a calm non-merger state that produces a basis to study the secular evolution aspects in galaxies. The data for Ψ_OC is taken from Bhattacharya et al. (2022) and is shown in Fig. 4.

Fig. 4. Cumulative distribution of open cluster IMF from Bhattacharya et al. (2022) used to determine the open clusters able to produce SNe in Sect. 2.3. The vertical dashed line shows MOC, lim defined with Eq. (10).

Only one open cluster in the Ψ_OC is above the M_OC, lim indicating that very few open clusters are able to form a stable SB. The mass of the cluster is 14 108 M_⊙, providing the effective value of

$$\begin{array}{c}\hfill s_{\mathrm{r}}=1.19.\end{array}$$(11)

In a subsequent analysis, we will need the number of SBs produced per formation of a unit stellar mass (≡κ) that we can derive from the open cluster IMF. For our case, the total mass of the Bhattacharya et al. (2022) catalogue is 137 514 M_⊙, indicating that κ = 1/137 514 ≈ 7.3 × 10⁻⁶ M_⊙⁻¹ SBs are formed per unit stellar mass formation.

3. Dynamical coupling of a SB with its environment

An evacuated cavity within the interstellar gas introduces a perturbation into the local galactic potential, which alters the trajectories of particles that are making up the ambient medium (e.g. stars and DM). In the case of a non-zero relative velocity between the SB and the ambient medium, this effect resembles the more familiar case of dynamical friction due to a massive object (perturber), by giving rise to an energy–momentum exchange between the perturber and its surroundings.

For the dynamical friction perspective, the main difference between a massive perturber and the SB is that in the latter case, mass is merely displaced from its interior to the SB edges, hence inflating the SB and introducing zero net mass into the environment. For the idealised spherically symmetric SB and Newton’s shell theorem, the SB has no dynamical effect on particle trajectories outside its own boundaries. However, within the SB, the same theorem implies that the perturbation to a particle’s trajectory arises due to the missing mass within the spherical region interior to the particle’s radial coordinate relative to the SB centre, which means that the particle experiences an acceleration

$$\begin{array}{c}\hfill \mathit{a}=-\frac{GM(<r)}{r^3}\mathit{r},\end{array}$$(12)

where

$$\begin{array}{c}\hfill M(<r)={\displaystyle \underset{0}{\overset{r}{\int }}}4\pi r^{\prime 2}{\rho }^{\ast }\mathrm{d}r\prime =-\frac{4\pi r^3{\rho }_{\mathrm{env}}}{3}·\end{array}$$(13)

The last part of the equation holds only for the interiors of the SBs.

Compared with the case of a massive object propagating through a medium, there are a few significant differences when considering dynamical friction with a SB as the perturber:

1. Since the effective perturbing mass is negative, the particle trajectories tend to diverge rather than converge when traversing the SB. As a result, the wake behind the perturber is underdense rather than overdense relative to the unperturbed state.

2. The spatial extent of the SB’s dynamical effect is limited to its own size, in contrast to a massive perturber in which case its influence formally extends to infinity. Therefore, no divergence arises when computing the energy–momentum exchange within an extended medium and no factor resembling the Coulomb logarithm appears in the final expression.

3. In contrast with most massive perturbers, the SB can evolve on timescales of interest for the dynamical friction problem (i.e. the size of the SB at the time of entry and exit of a given particle is generally different). As a result, the encounter is not conservative even in the reference frame of the SB. This effect can dominate the energy exchange compared with interaction with an unevolving perturber, where only momentum is exchanged in the perturber frame (in the limit m ≪ M).

The acceleration field for a particle at different times and locations within the evolving SB is shown in Fig. 3. In keeping with the above discussion, the acceleration peaks near the inner edge of the SB wall where the interior (missing) mass is highest. However, from the perspective of dynamical coupling between the SB and the environment, a more relevant quantity is the net energy gain or loss of a particle during the entire time spent within the SB. We show in Appendix A.2 that in the frame of the SB, this is determined by the relative radii of the SB at the time of entry and exit of the particle; if R_exit > R_enter, then net energy gain is positive and vice versa (see Eq. A.14). Since the SB evolves through both expansion and contraction phases, one can see that the dynamical friction force can have either sign, even in a uniform medium, in contrast with a massive perturber.

The change in a star’s velocity parallel to its initial direction as it traverses the SB is shown in Fig. 5, for different impact parameters and phases in the SB evolution at the time of entry. At early entry times and for a sufficiently rapidly moving star (upper panel), we have R_exit > R_enter, and the velocity gain is positive regardless of the impact parameter b. In contrast, stars that enter the SB in later phases will be able to traverse it only when the latter has already evolved well into the contraction phase. In this case, whether the ratio R_exit/R_enter is above or below unity depends on the impact parameter, with lower values of b resulting in longer residence times inside the SB, lower values of R_exit/R_enter, and a greater energy loss.

Fig. 5. Velocity changes for stars passing through the SB. The top part of the panel shows stars that enter the SB with v0 = (15, 0) km s−1, while in the bottom half v0 = (5, 0) km s−1. The colour depicts the time the star enters the SB (see colour scale at right).

Stars with lower velocities than the characteristic expansion or contraction speed of the SB will in most cases experience a net loss of energy and momentum¹, since the SB will have had time to complete the expansion phase and contract to a size R_exit < R_enter by the time the star exits it (Fig. 5, lower panel).

The significance of the energy gain or loss of a given star compared to its initial kinetic energy is seen from Eq. (A.14), which can be written as

$$\begin{array}{c}\hfill \frac{v_{\mathrm{out}}^2}{v_{\mathrm{in}}^2}=1+\frac{[v_{\mathrm{ch}}^2(R_{\mathrm{exit}})-v_{\mathrm{ch}}^2(R_{\mathrm{enter}})]}{v_{\mathrm{in}}^2},\end{array}$$(14)

where v_in and v_out refer to the stellar velocities at the time of entry and exit, respectively, and $v_{\mathrm{ch}}(r)=\sqrt{8\pi G{\rho }_{\mathrm{env}}{r}^2/3}$. Thus, stars with v_in ≫ max[v_ch(R_enter),v_ch(R_exit)] only obtain a minor correction to their initial energy and vice versa (unless, by chance, R_enter ≈ R_exit). The same applies to their trajectories, and hence v_ch delineates the boundary between the regime where the phase space density is only slightly deformed by the SB and the collision-like regime, whereby stars are displaced by a considerable distance within phase space.

3.1. A single SB in an environment

The total momentum transfer between a SB and a stellar population is found by double integration. First, over the phase space of stars and second, the SB lifetime. The outcome can be expressed as

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathit{P}={\displaystyle \int \mathrm{\Delta }\mathit{p}}\mathrm{d}{F}_{\mathrm{in}}\mathrm{d}{S}_{\mathrm{b}}\mathrm{d}{t}_{\mathrm{in}},\end{array}$$(15)

where Δp is the momentum transfer from a single encounter, dF_in is the flux of stars through the unit area of the SB surface, dS_b is the surface element (see Eqs. (A.16) and (A.19)), and t_in is the entry time.

We note that in general the stellar flux F_in through the SB surface depends on the location on this surface; obtaining a non-zero value for ΔP relies on this being the case (otherwise the integral in Eq. (15) vanishes by symmetry). Typically, this reflects the relative motion between the SB and the stellar component, whereby more stars are encountered on one hemisphere of the SB. If there exists a frame in which the stellar distribution is approximately homogeneous and isotropic, then ΔP is non-zero only if the SB is moving relative to this frame. From symmetry considerations one finds that the momentum transfer ΔP is parallel to the direction of the SB’s motion in this frame. In contrast with the classical massive perturber, however, ΔP along this direction can take either sign depending on the specific parameters.

To illustrate the momentum exchange between the SB and the stellar population, we consider an idealised configuration in which a spherical SB with a well-defined (thin) edge moves through a homogeneous stellar field with a Maxwellian velocity distribution. This set-up allows a semi-analytic treatment outlined in Appendix A.2.1, which can serve as a validation of the full numerical model.

Figure 6 shows the momentum gain experienced by stars as a function of their entry time, multiplied by the rate at which they are entering the SB at that time; the total gain or loss ΔP over the SB lifetime is found by integration over t_enter. In broad terms, the stellar population experiences a net momentum gain (along the direction of the SB propagation) in the early phases of the SB evolution and vice versa. This is consistent with the discussion in the previous section, whereby a star entering the SB in its growing phase will more frequently have R_exit/R_enter > 1 and gain energy by the encounter. As expected, most of the momentum exchange is done by stars that traverse the SB when the latter is not far from its maximum size, while stars entering either very early or very late in the SB life cycle have R_enter, R_exit ≪ R_max and their energy–momentum transfer is correspondingly small (see Eqs. A.13 and A.14).

Fig. 6. Momentum exchange rate between the stellar population and an idealised thin-edged SB with constant expansion and contraction velocities vexp = 20 km s−1 and vcol = 10 km s−1, respectively, and expansion phase duration of 30 Myr. The ambient gas density is 0.015 M⊙ pc−3. The stellar density is 0.1 M⊙ pc−3, with a Maxwellian velocity distribution with σ = 25 km s−1. The relative velocity between the SB and the stellar distribution is vb = 10 km s−1. The solid line shows the semi-analytical solution given by Eq. (A.21); the symbols correspond to the full numerical solution.

Finally, the angular momentum can be acquired in the same way as momentum. The integration is done over the angular momenta changes of individual particles (ΔL_z′):

$$\begin{array}{c}\hfill \mathrm{\Delta }{L}_z={\displaystyle \int \mathrm{\Delta }}{L}_z\prime \mathrm{d}{F}_{\mathrm{in}}\mathrm{d}{S}_{\mathrm{b}}\mathrm{d}{t}_{\mathrm{in}}.\end{array}$$(16)

3.2. Including all the SBs

Superbubbles do not usually occur as single isolated entities; instead, in most star-forming galaxies, they are generated at a certain rate throughout the gas disc. The resulting dynamical coupling between the stellar and gas components can be described by some average momentum exchange rate between the two components per unit disc area. Specifically for the case of circular orbits within the gas disc, the exchange can be characterised in terms of either linear or angular momentum, collectively denoted as ΔΘ = {ΔP_x, ΔL_z}.

3.2.1. Transfer of linear momentum

First, we consider the exchange of linear momentum, taking ΔΘ = ΔP_x. Within a disc-like environment, we consider a gas slab with area A and height h, and denote the rate of SB formation per unit area as f_b. As the star formation activity is concentrated towards the mid-plane of the gas disc (the region of coldest gas), we assume that the SB centres are likewise in the mid-plane, disregarding any vertical variation.

The rate of SB formation in the above region is therefore f_bA, each SB contributing an exchange of momentum ΔP_x. This yields a momentum gain rate of gas within the region −f_bAΔP_x and a corresponding rate for the stellar population (with opposite sign), or equivalently, an effective force per unit area, dF_x/dA = −f_bΔP_x.

The gas mass per unit area (i.e. surface density) is Σ_gas = ρ_gash. Hence the gas deceleration rate is

$$\begin{array}{c}\hfill s\equiv \frac{\mathrm{d}{v}_{\mathrm{gas}}}{\mathrm{d}t}=-\frac{1}{{\rho }_{\mathrm{gas}}h}\frac{\mathrm{d}{F}_x}{\mathrm{d}A}=-\frac{f_{\mathrm{b}}\mathrm{\Delta }{P}_x}{{\rho }_{\mathrm{gas}}h}·\end{array}$$(17)

3.2.2. Transfer of angular momentum

The angular momentum gain or loss of the gas disc can be evaluated by taking ΔΘ = ΔL_z; the effect can be characterised by the relative rate of change of specific angular momentum of a gas element, or equivalently, the rate of radial migration of a disc annulus.

The former essentially yields a timescale (more precisely, its inverse) in which a significant change in angular momentum can take place,

$$\begin{array}{c}\hfill \zeta \equiv \frac{{\dot{L}}_z}{L_z}=-\frac{\mathrm{\Delta }{L}_z{f}_{\mathrm{b}}}{{\mathrm{\Sigma }}_{\mathrm{gas}}R{v}_{\mathrm{c}}},\end{array}$$(18)

where R is the galactocentric distance, and we have assumed that the rotational velocity is very close to the circular speed v_c² = −R∂Φ/∂R determined by the galactic potential.

To determine the resulting radial migration speed, one can use the continuity and angular momentum conservation equations for the gas disc in cylindrical coordinates,

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}{\mathrm{\Sigma }}_{\mathrm{gas}}+\frac{1}{R}\frac{\partial }{\partial R}\left(R{\mathrm{\Sigma }}_{\mathrm{gas}}{v}_{\mathrm{R}}\right)=0\end{array}$$(19)

and

$$\begin{array}{c}\hfill \frac{\partial }{\partial t}\left({\mathrm{\Sigma }}_{\mathrm{gas}}{R}^2\mathrm{\Omega }\right)+\frac{1}{R}\frac{\partial }{\partial R}\left(R{\mathrm{\Sigma }}_{\mathrm{gas}}{v}_{\mathrm{R}}{R}^2\mathrm{\Omega }\right)=T_{\mathrm{A}},\end{array}$$(20)

respectively, where v_R is the radial velocity and the right-hand side is the imposed torque per unit surface area. The parameter Ω denotes the angular velocity. Using Eq. (19), the angular momentum equation can be rewritten in the form

$$\begin{array}{c}\hfill v_{\mathrm{R}}\frac{\partial }{\partial R}\left(R^2\mathrm{\Omega }\right)=\frac{T_{\mathrm{S}}}{{\mathrm{\Sigma }}_{\mathrm{gas}}}\equiv T,\end{array}$$(21)

where T is the torque per unit mass, and we assume ∂(R²Ω)/∂t = 0, which is justified if the orbital velocity only depends on R.

In the present context, the torque per unit mass from dynamical friction is T = ΔL_zf_b/Σ_gas, whereby Eq. (21) yields

$$\begin{array}{c}\hfill v_{\mathrm{R}}=\frac{\mathrm{\Delta }{L}_z{f}_{\mathrm{b}}}{{\mathrm{\Sigma }}_{\mathrm{gas}}\partial (R^2\mathrm{\Omega })/\partial R}=\zeta R{[1+\frac{\mathrm{d}ln{v}_{\mathrm{c}}}{\mathrm{d}lnR}]}^{-1},\end{array}$$(22)

where we have used v_c = RΩ.

Although the inward speed is a very intuitive quantity, a more frequently used measure of gas migration is the flux of mass through some radius. The mass transfer inwards is included by converting v_R to mass flux by multiplying it with surface density and the area shifted over the time interval:

$$\begin{array}{c}\hfill \mathrm{d}{M}_{\mathrm{gas}}={\mathrm{\Sigma }}_{\mathrm{gas}}2\pi R{v}_{\mathrm{R}}\mathrm{d}t.\end{array}$$(23)

3.2.3. The SB formation rate

Both SB influence estimators, (17) and (18), depend on the SB rate, to which we provide three estimates: from direct observations, from the Kennicutt–Schmidt relation and from the SB saturation limit.

Ehlerová & Palouš (2013) measured the surface density of SBs with radius over 100 pc in the solar neighbourhood to be about 4 per kpc². To convert this into the SB rate, we use the lifetime of SB t_end = 71.6 Myr (see Sect. 2.2), from which during 31.6 Myr it has an inner radius exceeding 100 pc. To have four superbubbles per square kiloparsec, the formation rate of SBs should be f_b = 127 kpc⁻² Gyr⁻¹.

The Kennicutt–Schmidt relation can be used for SB rate estimates as well. In Sect. 2.3, we defined the conversion rate (κ), showing how many SBs form per formed stellar mass. Hence, we can convert the surface density of the star formation rate to the SBs formation rate by simple multiplication. Using the Kennicutt–Schmidt relation for non-starburst galaxies from de los Reyes & Kennicutt (2019) with the surface density of gas 13.6 M_⊙ pc⁻² (McKee et al. 2015) we obtain the rate of SBs f_b = 42 kpc⁻² Gyr⁻¹.

If the formation rate of SBs is high and the lifetime of SBs long, they start to overlap, and new ones form on the edges of older SBs (Barnes et al. 2022). This does not allow the SB to complete its full evolution cycle, and our approximation of a single SB influence breaks down. To estimate the SB rate when the saturation of independent SBs occurs, we assume that a higher SB rate does not produce a larger Θ transfer. We derive the saturation limit by simulating the formation of SBs at random positions and times, and f_b measures the number of random points inside of SBs (defined as the distance from a random point to the SB centre is less than R_m) on average. At a SB formation rate of f_b = 25 kpc⁻² Gyr⁻¹ the number of overlapping SBs will be over unity on average, and we use this rate as the saturation limit. Figure 7 shows the saturation level of our simulation.

Fig. 7. Average number of SBs that are in a random point based on simple numerical simulation.

Considering that all these estimates are diverse and almost independent, we choose the most conservative one for further analysis: f_b = 25 kpc⁻² Gyr⁻¹. Changing the value of f_b changes the secular evolution speed proportionally.

4. Examples of SBs in various environments

4.1. SBs in a homogeneous isotropic environment

This section shows the amount of de-acceleration (s) SBs induce in gas. We assume a homogeneous and isotropic environment as it is an approximation to a more realistic environment. First, we solve the SB influences numerically using (17), and then construct an approximation formula to simplify its future evaluations.

The set-up of the environment is a slab of gas, which we complement with homogeneous stellar distribution with isotropic velocities. The parameter v denotes the average velocity between the gas and stars, an analogue to asymmetric drift; σ denotes the stellar velocity dispersion. We described the calculation prescription in Sect. 3.2.1, and the implementation parameters are in Table 1.

Figure 8 shows the resulting de-acceleration of the gas slab, both in kinematic and SB property aspects. For the context, MW-like galaxy that formed its discs around redshift 1.5 − 2 (van der Wel et al. 2014) giving them the time-span to modify their discs for about 10 Gyr once they evolve to redshift 0. If we assume a flat rotation curve with a velocity plateau of 230 km s⁻¹, then during that time the gas can migrate from due to SB ≈ 10s/230, or 20% if s = 5 km² s⁻² kpc⁻¹. The evaluation assumes a homogeneous host environment, which we discard in the next section. Figure 8 shows that the highest impact comes from the parameter space where the stellar disc is cold, its velocity dispersion is below 10 km s⁻¹, but with some dependence on the relative velocity of the gas disc and the stellar disc. The bottom panel show the dependence on the SB properties, with slow-collapsing SBs being more efficient in transferring momentum: a longer age extends the time of momentum transfer. Analogously, the larger size (described by s_r) increases the cross-section of the SB, allowing more stars to transfer their momenta. Larger SBs also have larger amounts of displaced gas mass and they provide higher accelerations (see Sect. 2).

Fig. 8. De-acceleration s in a homogeneous environment. The long-term effects can be estimated using the relation km s−1 per Gyr as it approximates km2 s−2 kpc−1. The top panel shows the dependence of s on the kinematical parameters of the homogeneous environment: the velocity and velocity dispersion. The magenta point shows the difference between MW gas and stellar disc values. The bottom panel shows the dependence on the SB parameters of collapse speed and scaling, with kinematic parameters v and σ fixed at the magenta point in the top panel. The de-acceleration field can be approximated using Eq. (24).

To evaluate the de-acceleration with more ease, we fitted an empirical approximation formula to the top panel in Fig. 8. We made an analytic approximation for the simple case using the SymbolicRegression package (Cranmer 2023) in Julia (Bezanson et al. 2017). The mean and standard deviation of the residuals of the fit is −0.0013 km² s⁻² kpc⁻¹, 0.257 km² s⁻² kpc⁻¹. The best model has the form

$$\begin{array}{c}\hfill s_{\mathrm{approx}}=\frac{v_{\mathrm{bub}}+\sigma +a_1}{a_3(v_{\mathrm{bub}}-a_2)[a_5+exp\left(\frac{\sigma +a_4}{v_{\mathrm{bub}}}\right)]},\end{array}$$(24)

where a₁ = −112 km s⁻¹, a₂ = −1.64 km s⁻¹, a₃ = 0.650 km² s⁻² kpc⁻¹, a₄ = 5.94 km s⁻¹, and a₅ = −0.217. This equation assumes that the velocities are in units of km s⁻¹, distances in kpc, and the resulting de-acceleration in km² s⁻² kpc⁻¹.

4.2. SBs in a test galaxy

A homogeneous environment can approximate regions without a steep density gradient and host acceleration field. However, the orbits in a galaxy can have a circular nature, and they can produce mismatches with the approximation. In this section we characterise a situation where Θ = L_z, describing the SB effects by changes in angular momentum.

4.2.1. Underlying galaxy

The numerical estimations of the order of magnitude effects require specifying a galaxy. We use a simple analytical potential for a dynamical model evaluation that mimics a disc galaxy.

The galaxy has two density components: DM is described with a Navarro-Frenk-White (NFW) profile (Navarro et al. 1996), and stellar distribution is described with Miyamoto-Nagai (MN) profile (Miyamoto & Nagai 1975). The gas disc is assumed to produce an insignificant contribution to total density, and hence, its contribution to gravitational potential is neglected. The MN profile has the potential

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{MN}}=\frac{-G{M}_{\mathrm{MN}}}{\sqrt{R^2+{[a_{\mathrm{MN}}+\sqrt{z^2+b_{\mathrm{MN}}^2}]}^2}}·\end{array}$$(25)

The chosen values for the MN profile were a_MN = 3.0 kpc, b_MN = 0.4 kpc, and M_MN = 7.11 × 10¹⁰ M_⊙. The NFW profile has density in the form of

$$\begin{array}{c}\hfill \rho =\frac{{\rho }_{\mathrm{NFW}}}{\frac{r}{R_{\mathrm{NFW}}}{(1+\frac{r}{R_{\mathrm{NFW}}})}^2}\end{array}$$(26)

with the parameter values ρ_NFW = 0.001 × 10¹⁰ M_⊙ kpc⁻³ and R_NFW = 15 kpc.

To evaluate the SB influence on the galaxy, we need to know the stellar and DM particle kinematics. They are found from the Jeans equations (Jeans 1915) by assuming stationarity and cylindrical or spherical symmetry for the disc and DM halo, respectively.

The first equation for solving disc kinematics is the radial Jeans equation

$$\begin{array}{c}\hfill \frac{\partial (\rho {\sigma }_{\mathrm{R}}^2)}{\partial R}=-(1-{\beta }^2)\rho \frac{\partial \mathrm{\Phi }}{\partial R}-\frac{1-{\alpha }^2}{R}(\rho {\sigma }_{\mathrm{R}}^2),\end{array}$$(27)

with the simplifying assumptions that tangential dispersion is proportional to the radial dispersion and rotational velocity to circular velocity:

$$\begin{array}{cc}\hfill {\sigma }_{\theta }& =\alpha {\sigma }_{\mathrm{R}},\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill v_{\mathrm{rot}}& =\beta v_{\mathrm{circ}}.\hfill \end{array}$$(29)

In our applications the velocity ellipsoid is slightly elongated in the radial direction, α = 0.8, and the asymmetric drift (the relative speed between the circular and rotational velocities) is about ≈10 km s⁻¹ as we chose β = 0.95. The vertical motions originate from solving the vertical Jeans equations by assuming that the tilt term in Jeans equations is zero:

$$\begin{array}{c}\hfill \frac{\partial (\rho {\sigma }_z^2)}{\partial z}=-\rho \frac{\partial \mathrm{\Phi }}{\partial z}·\end{array}$$(30)

We solved the equations using the Runge-Kutta-Fehlberg method by integrating ρσ_r² or ρσ_z² along the coordinate axis² and then solving for the kinematical property.

By assuming non-rotating isotropic velocity distribution, the solution for the spherical Jeans equations for the DM kinematics was provided in Kipper et al. (2019). Our aim is also to test if the rotation of DM influences the angular momentum transfer by gas. We wanted to achieve a situation where as many other parameters as possible (aside from the rotation) are shared with the non-rotating solution. The rotation can be induced by modifying the spherical solution by lowering one side of the v_rot distribution (p_rot,  DM) on one side and increasing it by the same amount on the other side, similar to what was done in Joshi & Widrow (2023). Effectively, it means that the overall potential and density remain the same if the rotational direction of a fraction of stars is switched³. The factor p_rot, DM describes the modification of the distribution function

$$\begin{array}{c}\hfill f_{\mathrm{rot}}=\frac{1}{1+exp(-kx)}-\frac{1}{2},\end{array}$$(31)

with x ≡ |v_rot|/σ_DM. The value of p_rot, DM is multiplied by 1 ± f_rot, depending on whether v_rot is positive or not. In the implementation, k equals ±5, which transforms the mean value to ±0.37σ. The solution of the kinematics is in Fig. 9.

Fig. 9. Kinematic solution of the test-galaxy from solving the Jeans equations. The rotation curve of the disc is β = 0.95 times the circular velocity curve.

4.2.2. The gas disc changes

Figure 10 shows the influence of the SBs on the gas disc in two representations: relative loss of angular momentum and radial velocity. We note two important things here. Firstly, there is a strong radial dependence in the inner region of the galaxy. The innermost region was not calculated as the SB would not fit the region without breaking down the shape assumptions. Secondly, the contributions by the DM and stellar components have opposite signs. This is specific only to the region where particles pass through the SB very quickly and the gravitational potential is cylindrical and does not occur in the homogeneous and isotropic medium. There are multiple effects in play simultaneously, so pinpointing the exact culprit for this effect is not feasible in the current study.

Fig. 10. Dependence of gas disc expansion due to SBs in terms of relative angular momentum change (see Eq. (18)) and of mass transfer (see Eq. (23)) in the top and bottom panel, respectively.

The rotation of the DM profile might produce some changes in the stellar disc (e.g. Herpich et al. 2015; Kataria & Shen 2022; Li et al. 2023). One physical reason is that the dynamical friction is sensitive to the kinematics of the particles, not only density (Chandrasekhar 1943). Therefore, we investigated the role of DM rotation on the gas disc via the dynamical friction. Figure 11 shows the difference between different DM rotation models. We conclude that the contribution from DM kinematics exists, but that it is weak.

Fig. 11. DM rotation influence to the transfer of angular momentum between DM and the gas disc (ζ). Different k values correspond to different rotations, with value ⟨v⟩ = 0.37σ, the dispersion of the DM halo is shown in Fig. 9.

5. Discussion

5.1. Observational verification of the SB profile

In Sect. 2.1 we describe the density profile of a SB as derived in Y17. In this section we verify the density’s consistency with the observations. Watkins et al. (2023) analysed the galaxy NGC 628 in terms of its SB content from the observations from the PHANGS46–JWST using MIRI’s F770W filter. We took advantage of this observation and focused on the SB situated at approximately 3 kpc, south-east from the centre of this galaxy. This SB is an ideal example of these processes as it has an almost circular shape.

We measured the surface brightness in circular regions from the centre of the SB, matching our coordinate centre. The Y17 results assume initially homogeneous gas distribution, but in our case, this is not fully satisfied if we want to investigate the ISM regions further than the extent of the SB. To test if this was an issue, we measured the surface brightness only in a well-defined shell-like structure as seen in the JWST image, also done in the work of Mayya et al. (2023). This approach accounts for the influence of various external factors, including the collapse of edges, gas diffusion, and other effects.

We generated the intensity of the profile from the intensity map of the image by averaging intensities at different radii from the shell centre. Thereafter, we binned these radii so that we could get a smoother profile. For the binned data, we calculated the mean and the error of the mean by bootstrapping Efron (1979). This gave us the observed surface brightness profile within errors, which is proportional to the density (ρ) profile given by Y17. The theoretical density profile is approximated in Sect. 2.1 following the Y17 simulation. This simulation assumes a homogeneous environment for all the initial conditions. The observations are more complex as the ISM and gas outside the SB are not homogeneous as well, and there are other effects contributing, such as the higher gas density on the spiral where one side of the SB resides or some higher rates of star formation at the edges of the SB from the pressured gas or other specific gas properties. We denote this conversion factor between the surface density of the Y17 fiducial profile and the surface brightness of the observations by C. To exclude ISM effects in the observations, we take into account radii up to 0.7 kpc, a bit further than the well-defined luminosity shell of the profile. Additionally, the observations have a background intensity that we denote as B.

The conversion between the line of sight and SB-centric coordinate in a galaxy is $l=\pm \sqrt{r^2-R^2}$, where r is the 3D distance from the centre of the SB and R is the distance from the centre of the circle as measured in the photometric image. Thus, the line-of-sight integral is defined as

$$\begin{array}{c}\hfill I(R)=C{\displaystyle {\int }_{-l_{\lim }}^{+l_{\lim }}}\rho [r(l)]\mathrm{d}l+B.\end{array}$$(32)

This gives us five free parameters, three that are the result of the observational image, C, B, l, and two that describe the SB, its expansion time t and the scaling factor s_r (analytically described in Sect. 2.2). These parameters are modelled within large and uniform priors. Our aim was to match the observed data with the theoretical model by following a Bayesian approach (Laplace 1774). We considered errors as coming from the standard deviation derived from bootstrapping plus a systematic uncertainty to account for the fact that the SB is not isolated or in a homogeneous environment.

The following equation gives us the likelihood of the posterior. We used 10⁵ iterations in an adaptive Monte Carlo scheme provided by Julia, as implemented by Vihola (2020):

$$\begin{array}{c}\hfill log\mathcal{L}\propto -0.5{\displaystyle \sum _{i=1}^N}\left[{\left(\frac{I_{\mathrm{obs},i}-I_{\mathrm{model},i}}{\mathrm{\Delta }{I}_i}\right)}^2\right].\end{array}$$(33)

Here I_obs denotes the surface brightness coming from the JWST image in the radii bins and the Markov chain Monte Carlo inference estimates of the free parameters for the I_model, and ΔI_i represent the associated errors.

The constraints of the time of the expansion of the SB to t_obs = 0.030 ± 0.002 Gyr indicate that the time of the observation of the specific SB’s evolutionary stage is at its biggest peak before starting to shrink. We also constrained the scaling factor to s_r = 1.18 ± 0.04, translating to ≈230 progenitor SN explosions. For the rest of the parameters, we obtained the line-of-sight distance l_lim = 0.37 ± 0.02 kpc and the luminosity background information of the JWST image as B = 10.79 ± 0.03 Jy/sr. The latter is within the measurements consistent with regions outside the galaxy.

As demonstrated in Fig. 12, the theoretical model closely matches the observational data, providing an agreement between the approach of the Y17 profile and the observed data. This verifies that the theoretical density profile from Y17 describes accurately the SB phenomena.

Fig. 12. Comparison between the surface brightness from the JWST image and the surface brightness from the Y17 best fit of the model. The observational flux of the SB (gold line) is plotted between three sigma (light gold dashed lines) and the standard deviation plus the systematic error of the bootstrap analysis for the binned data. The best fit of the model incorporating the five free parameters is superimposed on them (blue dashed line).

5.2. Robustness of the assumptions and approximations

Our analysis relies on a few assumptions about the SB itself, its environment and the stars interacting with it. The SB description is based on the Y17 simulation and their assumptions about the SB propagate here. Their dominant assumptions are homogeneous initial gas distribution with gas density analogous to the solar neighbourhood and the stars exploding near the centre being distributed uniformly in time.

5.2.1. Environment assumptions

The gas distribution in real galaxies is not homogeneous; the radial and vertical gradients of the gas density contradict our assumptions. The Sérsic profile describes the HI surface brightness distribution, whereas the Sérsic index characterises the steepness or gradient of the profile. Low values indicate nearly constant core with truncation, and large values (∼4) indicate cuspy profiles. The Sérsic index of averaged matter distribution of spirals from the THINGS survey is between 0.18 − 0.36 (Brinks & Portas 2017), showing the gradient of the gas disc to be much shallower than a typical stellar disc (Sérsic index ∼1), and thus the homogeneity assumption in the radial direction is a good approximation. The vertical gradient of gas discs is much steeper. The scale heights of different types of gas vary between 75 pc for molecular gas to 3000 pc for a hot ionised medium (Ryden & Pogge 2021). The highest contribution comes from a warm neutral medium with a scale height of ∼300 pc, which corresponds to the fiducial SB inner radius at the maximum extent. To see if the disc-height limited profile still matches the theoretical profile, we modelled observations with our SB model. In Sect. 5.1 we show that the profiles match best when the integration extent over the line of sight is 370 pc, which slightly exceeds the scale height of the gas disc, indicating that the SB almost fits the gas discs.

We intrinsically assume that the shape of the SB remains spherical. One reason for deviations from sphericity is shearing motions from galactic rotation. We tested the extent of shear distortions by assuming a flat rotation curve with v₀, the SB is at galactocentric radius R₀, the radius of SB is R_SB, and the time of shearing is Δt. In that case, the far edge of the SB can shift by v₀R_SBΔt/R₀ ≈ 225 pc from the centre when assuming values Δt = 30 Myr, R₀ = 8 kpc, R_SB = 300 pc, and v₀ = 200 km s⁻¹⁴. When trying to see whether the elongation and tilt are present in observations (e.g. the Barnes et al. 2023 images), we can see it for some SBs, which is a plausible deviation from our assumption. However, Binney & Tremaine (2008) showed that the gravitational potential ellipticity is of the order of one-third to one-half of the density ellipticity. Our current example has an ellipticity of 0.43, which translates to a potential ellipticity of 0.14 − 0.23 or an axis ratio of ≈0.82. Considering other uncertainties of the SB rate, evolution, and structure precision, we do not consider this deviation to be of first-order importance.

The secular evolution of the gas disc due to SBs is proportional to the rate of SB formation in the galaxy (see Eqs. (17) and (18)). Although our estimates of the SB formation rate have a high scatter (see Sect. 3.2.3), we would like to point out the possibility that there are spatial variations of the rate within a galaxy and between the galaxies. Finn et al. (2024) show that most massive molecular clouds, producing a greater number of likely open clusters able to develop into SBs, occur more likely in prominent spirals, contrary to inter-arm regions or regions in flocculent galaxies. In the latter cases, we suspect that our de-accelerations are overestimated. However, although smaller SBs also have smaller effects on the galaxy evolution (see Fig. 8), there are more of them and the cumulative effect can be comparable.

5.2.2. SB assumptions

In the case of the SB physics, we assume that the shape of the SB is the same as Y17 described, and that it has the same scaling relation they provided. We describe the suitability test of the profile in Sect. 5.1.

Although the observed and theoretical profiles match, the link between the SB and the SNe/SFR needed to create them needs a discussion. In Sect. 2.3, we adopted the approach that 100 SNe produce a stable SB and evaluated the results accordingly. Fewer SNe are also likely to make a SB, but it would not be very stable. The smaller number of SNe can be significant; for example, Barnes et al. (2022) showed many SBs in a galaxy when observed at high resolution. We suspect the smaller ones can also transfer angular momentum, but less effectively (see Fig. 8). This paper does not provide the analysis for that.

Figure 8 shows that in addition to the size of the SB, the collapsing speed (v_col) strongly influences how much the SB affects the transfer of momentum. If the corresponding time of the SB is twice as fast, then the transfer efficiency reduces, but the SB rate (f_b) increases⁵. Kinugawa et al. (2024) demonstrated that due to mass transfer between binaries, the time when the SNe can occur in an open cluster gets extended, possibly delaying the collapse of the SBs.

Superbubbles are destroyed in simulations due to shearing motions within the disc and/or are mixed into the ISM by hydrodynamical instabilities and turbulence, for example. In our analytical approach we cannot mimic it well, but we included an analysis to see what would happen to the transfer rate when the edges diffuse instead of collapse (i.e. when the inner and outer edge of the profile drift apart until the inner radius reaches zero). This analysis is provided in Appendix B.

5.3. Implications

Della Croce et al. (2024) reported that the outskirts of open clusters are expanding, the outer parts faster than the inner. The trend lasts for 30 Myr. Based on our description of the SB, the first 30 Myr is the expanding phase of the SB, causing an additional acceleration field (see Fig. 3). We tested the significance of the SB acceleration on a fictive open cluster by generating an open cluster with mass 10⁴ M_⊙ and radius 5 pc, and evaluated the relative acceleration originating from the SB: |a_SB|/(|a_SB|+|a_OC|) = 0.18 ± 0.09. The ± describes different stars in the cluster. We conclude that the SB can provide an additional acceleration field to alter the kinematics of open clusters. However, the further acceleration does not necessarily translate to expansion, but to extended times in the outer parts of the orbits or possibly easier evaporation.

There have been indications that the rotation of a DM halo affects the baryonic parts of the galaxy (Herpich et al. 2015; Kataria & Shen 2022; Li et al. 2023), or the opposite (Nuñez-Castiñeyra et al. 2023). In Fig. 11 we evaluate the extent to which SBs can cause the angular momentum transfer, which is about 5% per 10 Gyr. By measuring extremely precisely how much the gas discs shrink in a galaxy absent of spirals and bars (e.g. using chemical imprints), it is, in theory, possible to infer the rotation of the DM halo.

6. Summary

Galaxies evolve via stochastic and secular processes, for example two-body relaxation, virialisation, momentum transfer by spirals and bars, major mergers, minor mergers. The last item in the list can provide torques via dynamical friction, changing the angular momentum of the merger with the disc. We studied additional effects related to dynamical friction that are caused not by external satellites but by gas inhomogeneities, namely SBs introduced by multiple SN explosions originating from open clusters.

We characterised the SB influence by acceleration field correction originating from density field correction. We applied this correction to stellar and DM particle trajectories that pass through the SB and evaluated the velocity or angular momentum changes. By stacking together all the particles, including the time-dependence of the SB, and counting the constant formation-destruction cycle of the SBs, we measured the extent of the dynamical friction from the SBs.

Our summary of the SB influences on the galaxy is as follows:

• The acceleration field of the SB correction is directed outwards: if the gas mass is missing from the SB centre, the correction to mass is negative, and the radial acceleration is also negative (see Sect. 2 for the derivation of the effect).

• Stars entering the SB slowly can reverse their movement direction, indicating that stars are almost co-moving with the SB (see Sect. 3).

• The SB effects can be approximated empirically (see Eq. (24)). We provided an analytical approximation for a gas slab in stellar distribution that is homogeneous and isotropic.

• For the test galaxy, we estimated the change to the gas disc when both the SB and the galactic potential influence the orbits of stars. The relative rate of angular momentum changes is up to 4% Gyr⁻¹, which can account for significant changes over the galaxy’s lifetime.

• Although SBs are able to produce changes in galaxy disc over extended timescales, there are processes, such as an instability that drives both turbulence and inward mass migration in galaxies (Goldbaum et al. 2015, 2016), that do it faster. Instabilities drive gas inwards at the rate Ṁ = 0.1 − 1 M_⊙ yr⁻¹, compared to SBs Ṁ ≲ 0.02 M_⊙ yr⁻¹ (see the bottom panel of Fig. 8).

• Dynamical friction is a process that depends on the velocity of the host environment (Chandrasekhar 1943). As we are studying processes similar to dynamical friction, we estimated the possible contribution of the rotation of a DM halo to the transfer, and whether this can be used to infer the kinematics of DM. The differences are of the order of ∼0.3% Gyr⁻¹, or ∼3% over the time an average disc has existed. Measuring these rates is extremely difficult, but in theory the DM rotation influences the gas disc via SBs.

Data availability

The codes and scripts used to produce these results are available in GitLab https://gitlab.ut.ee/rain.kipper/onlybubbles

Acknowledgments

We thank Punyakoti Ganeshaiah Veena for helpful discussion. The present study was supported by the ETAG projects PSG700, PRG1006 and PRG2159. This work was partially supported by the Estonian Ministry of Education and Research CoE grant “Foundations of the Universe” TK202 and the European Union’s Horizon Europe research and innovation programme (EXCOSM, grant No. 101159513). AT was also supported by the COST Action CA21126 – Carbon molecular nanostructures in space (NanoSpace)-COST (European Cooperation in Science and Technology).

References

Appendix A: Effect of dynamical friction on a spherical cavity propagating through a uniform medium

Consider an idealised situation, in which a spherical cavity within the ISM is propagating through a uniform ambient medium. This medium experiences no significant direct interaction with the ISM (except via gravity). For example, this cavity may result from energy deposition from multiple SNe exploding within a stellar cluster (i.e. a SB; see e.g. Yadav et al. 2017), which can inflate a quasi-spherical region with a density several orders of magnitude below the ambient ISM. The ambient medium under consideration may be associated with the DM halo of the galaxy as well as field stars unassociated with the cluster and/or bubble. If the relative velocity between the bubble and the ambient medium is non-zero on average (this notion will be made more precise below), the particles comprising this medium that propagate through the bubble exert a non-zero net force on the bubble (more precisely, the bubble walls) due to dynamical friction.

First, let us consider a single particle of mass m propagating through a spherical cavity of radius R with impact parameter b. According to Newton’s shell theorem, the gravitational influence of a spherical density distribution can be reduced to the case where the mass internal to the test particle radius r from the symmetry centre is concentrated into said centre. The theorem is also applicable to a spherically symmetric “hole” within an otherwise uniform medium (ISM); in this case, if r < R, the effective internal mass is M(< r) =  − 4πr³ρ_ISM/3, where ρ_ISM is the ISM mass density, and vanishes if r > R (here we have assumed that the material evacuated from the cavity is concentrated in a narrow region near its edge and that the total mass is conserved). Hence, while the particle is at r > R, it propagates within the large-scale galactic gravitational potential as if the cavity was not present, while at r < R it experiences an effective extra contribution from the ”missing" mass internal to radius r.

A.1. Cavity of fixed radius

Let us first consider the case in which the radius of the cavity does not change appreciably over timescales of interest. If the mass of the particle entering the cavity satisfies m ≪ M, then the total energy of the particle before and after the encounter remains constant in the frame of the perturber (cavity; i.e. Δv = 0). If the particle velocity further satisfies the condition v² ≫ Gρ_ISMR², its trajectory is altered only slightly upon traversing the cavity (compared with the case where the bubble is absent). Hence, the change Δv in particle velocity is approximately perpendicular to v after the encounter. The perpendicular component of acceleration at r < R can be written as

$$\begin{array}{c}\hfill {\dot{v}}_{\perp }\approx \frac{GM(<r)}{r^2}cos\varphi =\frac{4\pi G{\rho }_{\mathrm{ISM}}b}{3},\end{array}$$(A.1)

where ϕ is the angle between the pericentre and the current position of the particle (measured from the centre of the cavity), and b is the impact parameter. The total change Δv_⊥ is obtained by integrating over the time interval that the test particle spends within the cavity. Writing dt = d(btanϕ/v)≈bdϕ/cos²ϕ and observing that r < R requires $-\sqrt{R^2-b^2}/b<tan\varphi <\sqrt{R^2-b^2}/b$, one obtains

$$\begin{array}{cc}\hfill \mathrm{\Delta }{v}_{\perp }& ={\displaystyle {\int }_{{\varphi }_{\min }}^{{\varphi }_{\max }}}\frac{4\pi G{\rho }_{\mathrm{ISM}}b}{3}\frac{bd\varphi }{v{cos}^2\varphi }={\frac{4\pi G{\rho }_{\mathrm{ISM}}{b}^2}{3v}tan\varphi |}_{-\sqrt{R^2-b^2}/b}^{\sqrt{R^2-b^2}/b}\hfill \\ \hfill & =\frac{8\pi G{\rho }_{\mathrm{ISM}}}{3v}b\sqrt{R^2-b^2}.\hfill \end{array}$$(A.2)

Since the scalar v remains unchanged by the encounter, the parallel component of the velocity change can be approximated as

$$\begin{array}{c}\hfill \mathrm{\Delta }{v}_{\Vert }=\sqrt{v^2-v_{\perp }^2}-v\approx -\frac{v}{2}{\left(\frac{\mathrm{\Delta }{v}_{\perp }}{v}\right)}^2=-\frac{32{\pi }^2{G}^2{\rho }_{\mathrm{ISM}}^2}{9v^3}{b}^2(R^2-b^2).\end{array}$$(A.3)

The corresponding transfer of parallel momentum is Δp_∥ = mΔv_∥.

The rate of momentum transfer (i.e. the effective drag force) for a collimated beam of particles of density n is obtained by multiplying (A.3) by the particle flux nv and integrating over the impact parameter:

$$\begin{array}{cc}\hfill F_{\Vert ,\mathrm{coll}}& ={\displaystyle \int \mathrm{\Delta }}{p}_{\Vert }nv2\pi bdb=-\frac{64{\pi }^3{G}^2{\rho }_{\mathrm{ISM}}^2}{9v^2}\rho {{\displaystyle \int 0}}_R^{}{b}^3(R^2-b^2)db\hfill \\ \hfill & =-\frac{16{\pi }^3{G}^2{\rho }_{\mathrm{ISM}}^2\rho R^6}{27v^2}.\hfill \end{array}$$(A.4)

Here ρ = mn is the mass density of incident particles and 2πbdb is a surface element perpendicular to the particle flux.

For an arbitrary distribution of particles, one should substitute

$$\begin{array}{c}\hfill \rho \to m{\displaystyle \int f}(v)d^{3}v,\end{array}$$(A.5)

where f(v) is the phase space density, to obtain

$$\begin{array}{c}\hfill F=-\frac{16{\pi }^3{G}^2{\rho }_{\mathrm{ISM}}^2{R}^6}{27}m{\displaystyle \int \frac{f(v)}{v^2}}\omega d^{3}v,\end{array}$$(A.6)

where ω = v/v is a unit vector in the direction of the incident particle velocity.

If the incident particle distribution is isotropic, then Eq. (A.6) shows that the net momentum exchange vanishes, as is expected from symmetry considerations. Next, let us consider the simplest non-trivial case, in which the phase space density is isotropic in a frame that propagates with a velocity v₀ relative to the perturber (cavity). By symmetry, the dynamical friction force must be parallel to v₀, hence we consider only the F‖v₀ component when taking the angular integral. Making use of the fact that f(v) is invariant upon Galilean transformations (or indeed Lorentz transformations if written in terms of p) and so is d³v, one obtains

$$\begin{array}{c}\hfill F=-\frac{16{\pi }^3{G}^2{\rho }_{\mathrm{ISM}}^2{R}^6}{27}m{\displaystyle \int \frac{f^{\prime }(v^{\prime })}{v^2}}cos\theta d^3{v}^{\prime },\end{array}$$(A.7)

where v′ = v−v₀ is the particle velocity in the isotropy frame and cos θ = ω · v₀/v₀ = v · v₀/(vv₀). Using the cosine rule to express cos θ = (v² + v₀² − v^′ 2)/2vv₀ and v² = v^′ 2 + v₀² + 2v′v₀cos θ′, where cos θ′ = v′ · v₀/(v′v₀), one can write Eq. (A.7) as

$$\begin{array}{cc}\hfill F& =-\frac{16{\pi }^3{G}^2{\rho }_{\mathrm{ISM}}^2{R}^6}{27}m{\displaystyle \int f^{\prime }}(v^{\prime })v^{\prime 2}d{v}^{\prime }\hfill \\ \hfill & \times {\displaystyle \int \frac{1+qcos{\theta }^{\prime }}{v_0^2{(1+q^2+2qcos{\theta }^{\prime })}^{3/2}}}d{\varphi }^{\prime }dcos{\theta }^{\prime }.\hfill \end{array}$$(A.8)

Here we have expressed the velocity space element as d³v′=v^′ 2dv′dϕ′d cos θ′, where dϕ′d cos θ′ is the solid angle in the isotropy frame, and defined q ≡ v′/v₀. The integral over the azimuthal angle ϕ′ is trivial and contributes a factor 2π. The integral over cos θ is also elementary, yielding 2/v₀², hence the final result is

$$\begin{array}{cc}\hfill F& =-\frac{64{\pi }^3{G}^2{\rho }_{\mathrm{ISM}}^2{R}^{6}m}{27v_0^2}{\displaystyle \int f^{\prime }(v^{\prime })v^{\prime 2}}d{v}^{\prime }.\hfill \end{array}$$(A.9)

A.2. Cavity of time-dependent radius

In the case of an expanding or contracting cavity the entry and exit radii of a particle traversing it will be different and the kinetic energy of the latter is no longer conserved by the encounter even if m ≪ M.

In the frame of the cavity, the radial component of a particle’s equation of motion within the cavity reads

$$\begin{array}{c}\hfill \frac{d{v}_r}{\mathit{dt}}=\frac{d}{\mathit{dt}}(v·\frac{r}{r})=F_r+\frac{L^2}{r^3},\end{array}$$(A.10)

where $L=r^2\dot{\theta }$ is the (conserved) angular momentum of the particle. Using the relation dv_r/dt = (1/2)dv_r²/dr and writing F_r = −dΦ(r)/dr, one obtains

$$\begin{array}{c}\hfill \frac{d}{\mathit{dr}}(\frac{v_r^2}{2}+\mathrm{\Phi }+\frac{L^2}{2r^2})=0,\end{array}$$(A.11)

or simply v²/2 + Φ = E = constant, where E is the conserved total energy.

For the cavity we can write

$$\begin{array}{c}\hfill \mathrm{\Phi }(r)=-\frac{2\pi G{\rho }_{\mathrm{ISM}}}{3}{r}^2=-\frac{1}{4}{v}_{\mathrm{ch}}^2\frac{r^2}{R_{\mathrm{enter}}^2},\end{array}$$(A.12)

where R_enter is the entry radius and we have introduced a characteristic velocity

$$\begin{array}{c}\hfill v_{\mathrm{ch}}=\sqrt{\frac{8\pi G{\rho }_{\mathrm{ISM}}{R}_{\mathrm{enter}}^2}{3}}.\end{array}$$(A.13)

Applying the energy conservation condition, one obtains the kinetic energy gain as a function of entry and exit radii into/from the bubble,

$$\begin{array}{c}\hfill v_{\mathrm{out}}^2=v_{\mathrm{in}}^2+\frac{v_{\mathrm{ch}}^2}{2}[{\left(\frac{R_{\mathrm{exit}}}{R_{\mathrm{enter}}}\right)}^2-1].\end{array}$$(A.14)

For a given R_enter, the exit radius R_exit can be determined by solving Eq. (A.11) for r(t) and subsequently eliminating t from r(t) = r_exp(t), where r_exp(t) characterises the evolution of the bubble radius. In general, the latter step has to be performed numerically. For a spherically symmetric bubble, R_exit is a function of entry time t_in (which also determines R_enter = r_exp(t_in)), particle velocity v, and the angle θ between particle velocity and the local bubble surface normal at t_in.

A.2.1. Integration over phase space

Once the energy/momentum gain or loss of a single particle entering the bubble at a given t and R_enter = r_b(t) is known, the next task is to determine the net transfer rate of a population of particles characterised by a phase space density f(v). For this we first need to find the number of particles crossing the bubble surface element in unit time. The inward particle flux in the frame comoving with the surface element is (defined as positive definite)

$$\begin{array}{c}\hfill d{F}_{\mathrm{in}}=-(v^{†}·\frac{v_{\exp }}{v_{\exp }})f^{†}(v^{†})d^3{v}^{†},\end{array}$$(A.15)

where v^† = v − v_exp is the particle velocity in the surface element frame, v is its velocity in the bubble frame (i.e. its centre), and v_exp is the local velocity of the surface element in the same frame. The scalar product in Eq. (A.15) can be written as −v^† · v_exp/v_exp = −|v| cos θ + v_exp > 0, where cos θ is the angle between the particle velocity and the local radial direction in the bubble frame. We note that this quantity has to be positive, which places a constraint on the angles at which particles of given a velocity can enter the bubble. The particle flux thus becomes

$$\begin{array}{c}\hfill d{F}_{\mathrm{in}}=(v_{\exp }-|v|cos\theta )f(v)d^{3}v,\mathrm{where}cos\theta <\frac{v_{\exp }}{|v|},\end{array}$$(A.16)

and we have again used the invariance of f and d³v by dropping the symbol †. The normalisation of f gives the number density of stars, or physical density if multiplied with the mass of a star

$$\begin{array}{c}\hfill m{\displaystyle \int f}(v)\mathrm{d}v={\rho }^{\ast }\end{array}$$(A.17)

Next, let us assume that a reference frame exists in which the particle distribution is isotropic (referred to below as the lab frame), and that the bubble propagates with a velocity v_b relative to this frame. The total rate of momentum transfer between the bubble and the particle population can then be written as

$$\begin{array}{c}\hfill \frac{\mathit{dP}}{\mathit{dt}}={\displaystyle \int m\mathrm{\Delta }}{v}_{\Vert }cos{\theta }_{v}d{F}_{\mathrm{in}}d{S}_{\mathrm{b}}.\end{array}$$(A.18)

Here mΔv_∥ is the transferred momentum of a single particle parallel to its initial direction of propagation (the net transfer of transverse momentum vanishes by symmetry in a homogeneous particle field and a spherically symmetric perturber), θ_v is the angle between the velocity v of the incoming particle and the bubble velocity v_b measured in the bubble frame, and dF_in is the particle flux defined by Eq. (A.16). The bubble surface element is

$$\begin{array}{c}\hfill d{S}_{\mathrm{b}}=R_{\mathrm{enter}}^{2}dcos\theta d\varphi ,\end{array}$$(A.19)

where the angles refer to the polar and azimuthal angles of the local surface normal vector in a coordinate system anchored to the incoming particle velocity v (in the bubble frame). This choice is motivated by the fact that in this coordinate system the integrand of Eq. (A.18) is independent of ϕ and the corresponding integral thus immediately yields 2π.

To make further use of the symmetries in the problem, the integral over particle velocities is best taken in the frame where the distribution is isotropic (i.e. lab frame, denoted by prime) by writing $f(\mathit{v})d^{3}v=f\prime (\mathit{v}\prime )d^{3}v\prime =f\prime (v\prime )v^{\prime 2}dv\prime dcos{\theta }_v\prime d{\varphi }_v\prime$, where θ_v′ denotes the angle between v′ and v_b. By symmetry, the integral over the corresponding azimuth ϕ_v′ is again trivial and gives 2π.

To proceed further, we need the relations between particle velocities and angles in the lab and bubble frames. Defining ${\beta }_{\mathrm{b}}\prime =v_{\mathrm{b}}/v\prime$, one obtains

$$\begin{array}{c}\hfill cos{\theta }_v=\frac{cos{\theta }_v^{\prime }-{\beta }_{\mathrm{b}}^{\prime }}{\sqrt{1-2cos{\theta }_v^{\prime }{\beta }_{\mathrm{b}}^{\prime }+{\beta }_{\mathrm{b}}^{\prime 2}}},v=v^{\prime }\sqrt{1-2cos{\theta }_v^{\prime }{\beta }_{\mathrm{b}}^{\prime }+{\beta }_{\mathrm{b}}^{\prime 2}}.\end{array}$$(A.20)

Putting all of the above together, Eq. (A.18) takes the form

$$\begin{array}{cc}\hfill \frac{\mathit{dP}}{\mathit{dt}}& =4{\pi }^{2}m{R}_{\mathrm{enter}}^2{{\displaystyle \int -\infty }}_{\infty }^{}{f}^{\prime }(v^{\prime })v^{\prime 3}d{v}^{\prime }{\displaystyle {\int }_{-1}^1}(cos{\theta }_v^{\prime }-{\beta }_{\mathrm{b}}^{\prime })dcos{\theta }_v^{\prime }\hfill \\ \hfill & \times {\displaystyle {\int }_{D_{\theta }}}(\beta -cos\theta )\mathrm{\Delta }{v}_{\Vert }dcos\theta ,\hfill \end{array}$$(A.21)

where

$$\begin{array}{c}\hfill D_{\theta }=\{cos\theta :cos\theta \ge -1,cos\theta \le min(1,\beta )\}\end{array}$$(A.22)

and β ≡ v_exp/|v|.

Appendix B: SB collapse analysis

One of the shortcomings of the current analysis is the unclear behaviour of a SB after the SNe cease depositing energy and momentum into the SB (see Sect. 5.2.2). Here we are testing how sensitive our results are to the specific form of the evolution of the SB after this cease, at t_max (Sect. 2.2 explains the fiducial behaviour). We are analysing three ways how the functional form of SB can be modified: separation of inner and outer edge collapsing speeds, then the vanishing of the SB by not changing radii in the collapsing phase, but the density amplitude; and last the influence of shearing to collapse. These modifications somewhat mimic different processes that may happen, we analyse them in the given order.

By relaxing our treatment of SB collapsing speeds of inner and outer radii (R_i and R_o), it allows us to mimic the outer shell stalling that is predicted for atmospheric conditions (Zel’dovich & Raizer 1967). In the main part of the paper we had one constant to describe both (v_col), but in the current analysis we split it into two separate speeds: inner (still v_col) and outer (v_o, col) shell. To understand how these changes may influence we remind that R_i determines what is the mass that is displaced to form the acceleration field correction and R_o effectively determines what is the shape of the acceleration field. When we keep v_col as constant, we can see how sensitive our results are to the precise acceleration field shape. The analysis of the set-up is same as in Sect. 4.1 – a homogeneous environment⁶. Figure B.1 shows what is the introduced de-acceleration changes when varying v_o, col. The varying range is from collapsing along with inner radius to ones retaining expansion speed outwards, even after the SNe cease supporting high temperature. The de-acceleration range is −0.27…−0.24, or about ± ≈ 6%. We conclude that the behaviour of outer shell does not alter de-acceleration (s) much.

The second aspect of the first test is aimed to see the sensitivity of R_i collapse speed. We kept the R_o, col as constant and calculated the de-acceleration dependence on the v_col. The corresponding result is shown in Fig. B.2 – they are qualitatively same as in Fig. 8: very slow collapses can provide about an order of magnitude larger contributions than the fast ones. There is little differences in faster outer shell collapse speeds. Overall, we conclude that the specific form of collapsing phase does not alter the results much.

Fig. B.1. The dependence of de-acceleration on the outer shell Ro collapse speed (vo, col, the negative value indicates continuous expansion, the positive collapse of radius). While keeping the inner radius collapse speed (vcol = 10 km s−1) as constant. We can see that the de-acceleration remains relatively constant throughout various speeds, indicating that the outer shell collapse speed alters the overall transfer of gas angular momentum only little.

Fig. B.2. The dependence of the de-acceleration on the inner radius (Ri) collapsing speeds (vcol). In the calculation, the inner radius reaches 0 at different times, hence the lifetime of SBs are different and overall contribution per SB is elevated. In the analysis we kept the rate of SB formation (fb) as constant. The lines show different outer radius collapsing speeds, and some do not reach small values as the outer radius of SB would shrink below inner radius, which we consider physically inconsistent.

Another scenario for how the SB may collapse is by intrusion of matter into the SB region, for example turbulence eddies or instability imprints entering the SB. These are not symmetric entities, hence they cannot be described in a symmetric manner, and we take a statistical approach: we re-scale the acceleration field after t_max with a variable that is 1 at t_max and will drop linearly to 0 upon the end of SB life t_b. So between t_max…t_b we have acceleration of SB as

$$\begin{array}{c}\hfill \mathbf{a}={\mathbf{a}}^{\mathrm{fid}}(\mathrm{t}={\mathrm{t}}_{\max })\frac{t_{\mathrm{b}}-t}{t_{\mathrm{b}}-t_{\max }}.\end{array}$$(B.1)

The outcome is shown in Fig. B.3. Since the collapsing and the vanishing versions are not directly comparable we chose the lifetime of SB as depending variable. We can see that the vanishing description of the SB has a stronger impact on the de-acceleration, hence our results are more conservative, or underestimated. The physical reason why the vanishing SB influences de-acceleration more strongly than the collapsing model does is the amount of mass that is displaced. In the case of the collapsing model it shrinks with radius in cubic form when the radius decreases linearly, but in the case of the vanishing model, it is a linear decrease.

Fig. B.3. The dependence of de-acceleration on the type of collapsing. Vanishing describes the situation where the SB density is altered according to Eq. B.1 and collapsing to the one where vcol = 10 km s−1 and vo, col = 0 km s−1.

The last possibility that we analyse for the SB collapse is the shearing of the SB due to differential rotation of galaxies. The amount of shearing depends on the SB location in the galaxy: inner regions have higher angular speeds than outer regions, and therefore shear the SB more. Since our used method breaks down upon distortions from sphericity (spherical symmetry of density was needed to evaluate SB acceleration field and to keep the stellar orbit integration bounds constrained to finite region) we are limited for only testing. The sheared SB is constructed by sampling particles according to the total gas density profile (both uniform and SB-containing), and shearing their location linearly according to the range required for the test. Hence we relied on evaluating the acceleration field of the sheared SN and compared it with non-sheared one. The chosen amount of shearing was chosen based on observations: Watkins et al. (2023) provided the SB axis lengths for NGC 628, we used their catalogue of SBs to select only SBs larger than 100 pc and found that the minimum axis ratio was 0.56, and mean 0.75. We sheared the SB to the value close to the extreme of observed ones: 0.6. The resulting acceleration field is shown in Fig. B.4. We concluded from the figure that maximally sheared SB can cause 6 − 39% offsets from homogeneous case at r = 0.5 kpc. The same numbers for SB with axis ratio of 0.75 provide 2 − 23% offsets.

Fig. B.4. Acceleration field of sheared SB. The top panel shows the acceleration field components of a close-to-maximally sheared SB mass distribution compared to a non-shared distribution. The corridors describe the 65 and 95% range of the components with different directions from the centre. The dashed line shows at what radius the bottom panel is measured. In the bottom panel we depict how the acceleration field changes with direction from the centre of the SB. The tangential component and the angular dependence of acceleration field can be observed.

All Tables

All Figures

	Fig. 1. Illustration of SBs in galaxy NGC 6946 from the intensity map of the THINGS survey (Walter et al. 2008). The inset in the top right corner is a zoomed-in image of a fraction of the galaxy where the locations of three SBs are shown (pink circles). A more detailed map of all the bubbles in the galaxy is shown in Brinks et al. (2007).
In the text

	Fig. 2. Density contrast profile approximation of the SB at t = 30 Myr inspired by the simulated SB behaviour in Y17 and described in Sect. 2.1.
In the text

	Fig. 3. Acceleration field of SB. The solid and dashed red lines are the outer and inner radius of the fiducial model from Y17, and are described with Eqs. (8) and (9). The colour scale shows the radial acceleration values described in Sect. 3.
In the text

	Fig. 4. Cumulative distribution of open cluster IMF from Bhattacharya et al. (2022) used to determine the open clusters able to produce SNe in Sect. 2.3. The vertical dashed line shows MOC, lim defined with Eq. (10).
In the text

	Fig. 5. Velocity changes for stars passing through the SB. The top part of the panel shows stars that enter the SB with v0 = (15, 0) km s−1, while in the bottom half v0 = (5, 0) km s−1. The colour depicts the time the star enters the SB (see colour scale at right).
In the text

Fig. 6. Momentum exchange rate between the stellar population and an idealised thin-edged SB with constant expansion and contraction velocities vexp = 20 km s−1 and vcol = 10 km s−1, respectively, and expansion phase duration of 30 Myr. The ambient gas density is 0.015 M⊙ pc−3. The stellar density is 0.1 M⊙ pc−3, with a Maxwellian velocity distribution with σ = 25 km s−1. The relative velocity between the SB and the stellar distribution is vb = 10 km s−1. The solid line shows the semi-analytical solution given by Eq. (A.21); the symbols correspond to the full numerical solution.

In the text

	Fig. 7. Average number of SBs that are in a random point based on simple numerical simulation.
In the text

Fig. 8. De-acceleration s in a homogeneous environment. The long-term effects can be estimated using the relation km s−1 per Gyr as it approximates km2 s−2 kpc−1. The top panel shows the dependence of s on the kinematical parameters of the homogeneous environment: the velocity and velocity dispersion. The magenta point shows the difference between MW gas and stellar disc values. The bottom panel shows the dependence on the SB parameters of collapse speed and scaling, with kinematic parameters v and σ fixed at the magenta point in the top panel. The de-acceleration field can be approximated using Eq. (24).

In the text

	Fig. 9. Kinematic solution of the test-galaxy from solving the Jeans equations. The rotation curve of the disc is β = 0.95 times the circular velocity curve.
In the text

	Fig. 10. Dependence of gas disc expansion due to SBs in terms of relative angular momentum change (see Eq. (18)) and of mass transfer (see Eq. (23)) in the top and bottom panel, respectively.
In the text

	Fig. 11. DM rotation influence to the transfer of angular momentum between DM and the gas disc (ζ). Different k values correspond to different rotations, with value ⟨v⟩ = 0.37σ, the dispersion of the DM halo is shown in Fig. 9.
In the text

Fig. 12. Comparison between the surface brightness from the JWST image and the surface brightness from the Y17 best fit of the model. The observational flux of the SB (gold line) is plotted between three sigma (light gold dashed lines) and the standard deviation plus the systematic error of the bootstrap analysis for the binned data. The best fit of the model incorporating the five free parameters is superimposed on them (blue dashed line).

In the text

Fig. B.1. The dependence of de-acceleration on the outer shell Ro collapse speed (vo, col, the negative value indicates continuous expansion, the positive collapse of radius). While keeping the inner radius collapse speed (vcol = 10 km s−1) as constant. We can see that the de-acceleration remains relatively constant throughout various speeds, indicating that the outer shell collapse speed alters the overall transfer of gas angular momentum only little.

In the text

Fig. B.2. The dependence of the de-acceleration on the inner radius (Ri) collapsing speeds (vcol). In the calculation, the inner radius reaches 0 at different times, hence the lifetime of SBs are different and overall contribution per SB is elevated. In the analysis we kept the rate of SB formation (fb) as constant. The lines show different outer radius collapsing speeds, and some do not reach small values as the outer radius of SB would shrink below inner radius, which we consider physically inconsistent.

In the text

	Fig. B.3. The dependence of de-acceleration on the type of collapsing. Vanishing describes the situation where the SB density is altered according to Eq. B.1 and collapsing to the one where vcol = 10 km s−1 and vo, col = 0 km s−1.
In the text

Fig. B.4. Acceleration field of sheared SB. The top panel shows the acceleration field components of a close-to-maximally sheared SB mass distribution compared to a non-shared distribution. The corridors describe the 65 and 95% range of the components with different directions from the centre. The dashed line shows at what radius the bottom panel is measured. In the bottom panel we depict how the acceleration field changes with direction from the centre of the SB. The tangential component and the angular dependence of acceleration field can be observed.

In the text