Open Access
Issue
A&A
Volume 693, January 2025
Article Number A206
Number of page(s) 16
Section Extragalactic astronomy
DOI https://doi.org/10.1051/0004-6361/202450564
Published online 20 January 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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

thumbnail 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, Ri and Ro, 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 Rm = (Ri + Ro)/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−3. The profile itself is shown in Fig. A.1 of Y17. Our approximation assumes that the SNe have pushed out all the gas inside Ri. 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 Ri and Ro, where we describe the density distribution as a quadratic polynomial. Outside of Ro, the SB has had no influence and the density-correction is zero. Combining all the regions together, we have the formula

ρ ρ ρ env = { ρ env , R < R i B A ( R R m ) 2 , R i R R o 0 , R > R o . $$ \begin{aligned} \rho ^* \equiv \rho - \rho _{\rm env} = \left\{ \begin{array}{ll} - \rho _{\rm env},&R < R_{\rm i}\\ B - A(R - R_{\rm m})^2,&R_{\rm i} \le R \le R_{\rm o}\\ 0,&R > R_{\rm o} \end{array}{\!.}\right. \end{aligned} $$(1)

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

thumbnail 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 tmax.

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 Ri and Ro 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 Ri and Ro with a constant speed vcol. We denote the lifetime of the SB as tend, defined as the time when the Ri shrinks back to 0 kpc. The analytical form for Ri is

R i fid = { R u ( t / t i , e ) α i , e , t t max R max , i v col ( t t max , i ) t max < t t end 0 , t > t end $$ \begin{aligned} R_{\rm i}^\mathrm{fid} = \left\{ \begin{array}{lll} \!\!R_{\rm u}(t / t_{\rm i,e})^{\alpha _{\rm i, e}},&\,&t \le t_{\rm max}\\ \!\!R_{\rm max,i} - v_{\rm col}(t-t_{\rm max,i})&t_{\rm max} < \!\!&t \le t_{\rm end}\\ \!\! 0,&\,&t > t_{\rm end} \end{array}\right. \end{aligned} $$(2)

and for Ro

R o fid = { R o 0 + R u ( t / t o , s ) α o , t t max R max , o v col ( t t max , i ) t > t max . $$ \begin{aligned} R_{\rm o}^\mathrm{fid} = \left\{ \begin{array}{lll} \!\!R_{\rm o0} + R_{\rm u}(t / t_{\rm o,s})^{\alpha _{\rm o}},&t \le t_{\rm max}\\ \!\!R_{\rm max,o} - v_{\rm col}(t-t_{\rm max,i})&t > t_{\rm max} \end{array}\right.. \end{aligned} $$(3)

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

thumbnail 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) Ri and Ro 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: Ri and Ro are found from fiducial values by multiplying them with the term sr, or equivalently

s r = R i R i fid = R o R o fid · $$ \begin{aligned} s_{\rm r} = \frac{R_{\rm i}}{R_{\rm i}^\mathrm{fid}} = \frac{R_{\rm o}}{R_{\rm o}^\mathrm{fid}}\cdot \end{aligned} $$(4)

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

s r = ( N SN 100 ) 1 / 5 . $$ \begin{aligned} s_{\rm r} = \left(\frac{N_{\rm SN}}{100}\right)^{1/5}. \end{aligned} $$(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 (MOC) by introducing a normalising constant nIMF so that

0 m n IMF Ψ ( m ) d m = M OC . $$ \begin{aligned} \int \limits _0^{\infty }mn_{\rm IMF}\Psi (m)\mathrm{d}m = M_{\rm OC}. \end{aligned} $$(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 (mSN, minmSN, max):

N SN = m SN , min m SN , max n IMF Ψ ( m ) d m = M OC m SN , min m SN , max Ψ ( m ) d m 0 m Ψ ( m ) d m · $$ \begin{aligned} N_{\rm SN} = \int \limits _{m_{\rm SN,min}}^{m_{\rm SN,max}}n_{\rm IMF}\Psi (m)\,\mathrm{d}m = M_{\rm OC}\frac{\int \limits _{m_{\rm SN,min}}^{m_{\rm SN,max}}\Psi (m)\,\mathrm{d}m}{\int \limits _0^{\infty }m\Psi (m)\,\mathrm{d}m}\cdot \end{aligned} $$(7)

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

For the numerical evaluation of ϑ we need to specify Ψ, mSN, min, and mSN, max. For the Ψ we selected the Kroupa IMF (Kroupa 2001), the mSN, min = 8 M, but the mSN, max is not well determined. There have been dimming periods observed of massive stars, and hints of faint or absent SNe with masses above mSN, max = 17 M (Reynolds et al. 2015; Byrne & Fraser 2022), which shows that there is a lack of SNe with masses above mSN, 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 mSN, 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 mSN, max = ∞ M for integrations. In these cases, the values are ϑ17 = 0.0054 M−1 and ϑ = 0.0084 M−1. In the current paper we adopt the more classical ϑ = ϑ.

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

R i = R i fid ( ϑ M OC / 100 ) 1 / 5 , $$ \begin{aligned} R_{\rm i}&= R_{\rm i}^\mathrm{fid}(\vartheta M_{\rm OC}/100)^{1/5},\end{aligned} $$(8)

R o = R o fid ( ϑ M OC / 100 ) 1 / 5 . $$ \begin{aligned} R_{\rm o}&= R_{\rm o}^\mathrm{fid}(\vartheta M_{\rm OC}/100)^{1/5}. \end{aligned} $$(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 NSN, lim ≈ 100. The mass of the open cluster corresponding to this is

M OC , lim N SN , lim / ϑ 11 800 M . $$ \begin{aligned} M_{\rm OC, lim} \equiv N_{\rm SN, lim} / \vartheta \approx 11\,800\,\mathrm{M_\odot }. \end{aligned} $$(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 sr 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.

thumbnail 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 MOC, 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

s r = 1.19 . $$ \begin{aligned} s_{\rm r} = 1.19. \end{aligned} $$(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−6 M−1 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

a = G M ( < r ) r 3 r , $$ \begin{aligned} \boldsymbol{a} = -\frac{GM({ < }r)}{r^3} \, \boldsymbol{r}, \end{aligned} $$(12)

where

M ( < r ) = 0 r 4 π r 2 ρ d r = 4 π r 3 ρ env 3 · $$ \begin{aligned} M({ < }r) = \int \limits _0^r 4\pi r^{\prime 2} \rho ^* \mathrm{d}r\prime = - \frac{4\pi r^3 \rho _{\rm env}}{3}\cdot \end{aligned} $$(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 Rexit > Renter, 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 Rexit > Renter, 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 Rexit/Renter 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 Rexit/Renter, and a greater energy loss.

thumbnail 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 momentum1, since the SB will have had time to complete the expansion phase and contract to a size Rexit < Renter 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

v out 2 v in 2 = 1 + [ v ch 2 ( R exit ) v ch 2 ( R enter ) ] v in 2 , $$ \begin{aligned} \frac{v_{\rm out}^2}{v_{\rm in}^2} = 1 + \frac{\left[v_{\rm ch}^2(R_{\rm exit}) - v_{\rm ch}^2(R_{\rm enter})\right]}{v_{\rm in}^2}, \end{aligned} $$(14)

where vin and vout refer to the stellar velocities at the time of entry and exit, respectively, and v ch ( r ) = 8 π G ρ env r 2 / 3 $ v_{\mathrm{ch}}(r) = \sqrt{8\pi G \rho_{\mathrm{env}} r^2/3} $. Thus, stars with vin ≫ max[vch(Renter),vch(Rexit)] only obtain a minor correction to their initial energy and vice versa (unless, by chance, Renter ≈ Rexit). The same applies to their trajectories, and hence vch 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

Δ P = Δ p d F in d S b d t in , $$ \begin{aligned} \Delta \boldsymbol{P} = \int \Delta \boldsymbol{p} \, \mathrm{d}F_{\rm in} \, \mathrm{d}{S\!}_{\rm b} \, \mathrm{d}t_{\rm in}, \end{aligned} $$(15)

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

We note that in general the stellar flux Fin 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 tenter. 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 Rexit/Renter > 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 Renter, Rexit ≪ Rmax and their energy–momentum transfer is correspondingly small (see Eqs. A.13 and A.14).

thumbnail 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 (ΔLz′):

Δ L z = Δ L z d F in d S b d t in . $$ \begin{aligned} \Delta L_z = \int \Delta L_z\prime \, \mathrm{d}F_{\rm in} \, \mathrm{d}{S\!}_{\rm b} \, \mathrm{d}t_{\rm in}. \end{aligned} $$(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 ΔΘ = {ΔPx, ΔLz}.

3.2.1. Transfer of linear momentum

First, we consider the exchange of linear momentum, taking ΔΘ = ΔPx. 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 fb. 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 fbA, each SB contributing an exchange of momentum ΔPx. This yields a momentum gain rate of gas within the region −fbAΔPx and a corresponding rate for the stellar population (with opposite sign), or equivalently, an effective force per unit area, dFx/dA = −fbΔPx.

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

s d v gas d t = 1 ρ gas h d F x d A = f b Δ P x ρ gas h · $$ \begin{aligned} s \equiv \frac{\mathrm{d}v_{\rm gas}}{\mathrm{d}t} = -\frac{1}{\rho _{\rm gas}h} \frac{\mathrm{d}F_x}{\mathrm{d}A} = -\frac{f_{\rm b} \Delta \!P_x}{\rho _\mathrm{gas} h}\cdot \end{aligned} $$(17)

3.2.2. Transfer of angular momentum

The angular momentum gain or loss of the gas disc can be evaluated by taking ΔΘ = ΔLz; 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,

ζ L ˙ z L z = Δ L z f b Σ gas R v c , $$ \begin{aligned} \zeta \equiv \frac{\dot{L}_z}{L_z} = -\frac{\Delta L_z f_{\rm b}}{\Sigma _{\rm gas} R v_{\rm c}}, \end{aligned} $$(18)

where R is the galactocentric distance, and we have assumed that the rotational velocity is very close to the circular speed vc2 = −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,

t Σ gas + 1 R R ( R Σ gas v R ) = 0 $$ \begin{aligned} \frac{\partial }{\partial t} \Sigma _{\rm gas} + \frac{1}{R}\frac{\partial }{\partial R} \left(R \Sigma _{\rm gas} v_{\rm R}\right) = 0 \end{aligned} $$(19)

and

t ( Σ gas R 2 Ω ) + 1 R R ( R Σ gas v R R 2 Ω ) = T A , $$ \begin{aligned} \frac{\partial }{\partial t} \left(\Sigma _{\rm gas} R^2 \Omega \right) + \frac{1}{R}\frac{\partial }{\partial R} \left(R\Sigma _{\rm gas}v_{\rm R} R^2 \Omega \right) = T_{\rm A}, \end{aligned} $$(20)

respectively, where vR 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

v R R ( R 2 Ω ) = T S Σ gas T , $$ \begin{aligned} v_{\rm R} \frac{\partial }{\partial R} \left(R^2\Omega \right) = \frac{T_{\rm S}}{\Sigma _{\rm gas}} \equiv T, \end{aligned} $$(21)

where T is the torque per unit mass, and we assume ∂(R2Ω)/∂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 = ΔLzfbgas, whereby Eq. (21) yields

v R = Δ L z f b Σ gas ( R 2 Ω ) / R = ζ R [ 1 + d ln v c d ln R ] 1 , $$ \begin{aligned} v_{\rm R} = \frac{\Delta L_z f_{\rm b}}{\Sigma _{\rm gas} \,\partial (R^2\Omega )/\partial R} = \zeta R \left[1 + \frac{\mathrm{d}\ln v_{\rm c}}{\mathrm{d} \ln R}\right]^{-1}, \end{aligned} $$(22)

where we have used vc = 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 vR to mass flux by multiplying it with surface density and the area shifted over the time interval:

d M gas = Σ gas 2 π R v R d t . $$ \begin{aligned} \mathrm{d}M_{\rm gas} = \Sigma _{\rm gas}2\pi R v_{\rm R}\mathrm{d}t. \end{aligned} $$(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 kpc2. To convert this into the SB rate, we use the lifetime of SB tend = 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 fb = 127 kpc−2 Gyr−1.

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−2 (McKee et al. 2015) we obtain the rate of SBs fb = 42 kpc−2 Gyr−1.

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 fb measures the number of random points inside of SBs (defined as the distance from a random point to the SB centre is less than Rm) on average. At a SB formation rate of fb = 25 kpc−2 Gyr−1 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.

thumbnail 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: fb = 25 kpc−2 Gyr−1. Changing the value of fb 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.

Table 1.

Parameters for determining the level of the influence to the gas slab by the stars via SBs.

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−1, then during that time the gas can migrate from due to SB ≈ 10s/230, or 20% if s = 5 km2 s−2 kpc−1. 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−1, 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 sr) 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).

thumbnail 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 km2 s−2 kpc−1, 0.257 km2 s−2 kpc−1. The best model has the form

s approx = v bub + σ + a 1 a 3 ( v bub a 2 ) [ a 5 + exp ( σ + a 4 v bub ) ] , $$ \begin{aligned} s_{\rm approx} = \frac{v_{\rm bub} + \sigma + a_1}{a_3(v_{\rm bub} - a_2) \left[a_5 + \exp \left(\frac{\sigma + a_4}{v_{\rm bub}}\right) \right]}, \end{aligned} $$(24)

where a1 = −112 km s−1, a2 = −1.64 km s−1, a3 = 0.650 km2 s−2 kpc−1, a4 = 5.94 km s−1, and a5 = −0.217. This equation assumes that the velocities are in units of km s−1, distances in kpc, and the resulting de-acceleration in km2 s−2 kpc−1.

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 Θ = Lz, 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

Φ MN = G M MN R 2 + [ a MN + z 2 + b MN 2 ] 2 · $$ \begin{aligned} \Phi _{\rm MN} = \frac{-GM_{\rm MN}}{\sqrt{R^2 + \left[a_{\rm MN}+\sqrt{z^2 + b_{\rm MN}^2}\right]^2}}\cdot \end{aligned} $$(25)

The chosen values for the MN profile were aMN = 3.0 kpc, bMN = 0.4 kpc, and MMN = 7.11 × 1010 M. The NFW profile has density in the form of

ρ = ρ NFW r R NFW ( 1 + r R NFW ) 2 $$ \begin{aligned} \rho = \frac{\rho _{\rm NFW}}{\frac{r}{R_{\rm NFW}} \left(1 + \frac{r}{R_{\rm NFW}}\right)^2} \end{aligned} $$(26)

with the parameter values ρNFW = 0.001 × 1010 M kpc−3 and RNFW = 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

( ρ σ R 2 ) R = ( 1 β 2 ) ρ Φ R 1 α 2 R ( ρ σ R 2 ) , $$ \begin{aligned} \frac{\partial (\rho \sigma ^2_{\rm R})}{\partial R} = -(1-\beta ^2)\rho \frac{\partial \Phi }{\partial R} - \frac{1-\alpha ^2}{R}(\rho \sigma ^2_{\rm R}), \end{aligned} $$(27)

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

σ θ = α σ R , $$ \begin{aligned} \sigma _\theta&= \alpha \sigma _{\rm R},\end{aligned} $$(28)

v rot = β v circ . $$ \begin{aligned} v_{\rm rot}&= \beta v_{\rm circ}. \end{aligned} $$(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−1 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:

( ρ σ z 2 ) z = ρ Φ z · $$ \begin{aligned} \frac{\partial (\rho \sigma ^2_z)}{\partial z} = -\rho \frac{\partial \Phi }{\partial z}\cdot \end{aligned} $$(30)

We solved the equations using the Runge-Kutta-Fehlberg method by integrating ρσr2 or ρσz2 along the coordinate axis2 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 vrot distribution (prot,  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 switched3. The factor prot, DM describes the modification of the distribution function

f rot = 1 1 + exp ( k x ) 1 2 , $$ \begin{aligned} f_{\rm rot} = \frac{1}{1 + \exp (-kx)} - \frac{1}{2}, \end{aligned} $$(31)

with x ≡ |vrot|/σDM. The value of prot, DM is multiplied by 1 ± frot, depending on whether vrot 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.

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

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

thumbnail 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 = ± r 2 R 2 $ 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

I ( R ) = C l lim + l lim ρ [ r ( l ) ] d l + B . $$ \begin{aligned} I(R) = C \int _{-l_{\rm lim}}^{+l_{\rm lim}} \rho [r(l)] \, \mathrm{d}l + B. \end{aligned} $$(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 sr (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 105 iterations in an adaptive Monte Carlo scheme provided by Julia, as implemented by Vihola (2020):

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

Here Iobs 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 Imodel, and ΔIi represent the associated errors.

The constraints of the time of the expansion of the SB to tobs = 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 sr = 1.18 ± 0.04, translating to ≈230 progenitor SN explosions. For the rest of the parameters, we obtained the line-of-sight distance llim = 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.

thumbnail 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 v0, the SB is at galactocentric radius R0, the radius of SB is RSB, and the time of shearing is Δt. In that case, the far edge of the SB can shift by v0RSBΔt/R0 ≈ 225 pc from the centre when assuming values Δt = 30 Myr, R0 = 8 kpc, RSB = 300 pc, and v0 = 200 km s−14. 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 (vcol) 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 (fb) increases5. 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 104 M and radius 5 pc, and evaluated the relative acceleration originating from the SB: |aSB|/(|aSB|+|aOC|) = 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−1, 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−1, compared to SBs ≲ 0.02 M yr−1 (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−1, 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


1

Here we are referring to stars with velocities directed towards the SB interior, rather than those with outward-directed velocities that are caught up by the expanding SB (which constitute a minority compared with the former in terms of how many enter the SB).

2

The characteristic aligns with coordinates.

3

Technically, this produces a valid solution for the solved point, but loses the validity of isotropy in other radii. For the current application it is a suitable approximation.

4

The calculations assume there are no compensating forces (e.g. pressure gradients) to restore the round shape.

5

Here we are assuming that the limiting factor of SBs is the spatial saturation (see Sect. 3.2.3).

6

Including a constant value to the rate of the SBs fb

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

  1. Barnes, A. T., Chandar, R., Kreckel, K., et al. 2022, A&A, 662, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Barnes, A. T., Watkins, E. J., Meidt, S. E., et al. 2023, ApJ, 944, L22 [NASA ADS] [CrossRef] [Google Scholar]
  3. Benito, M., Criado, J. C., Hütsi, G., Raidal, M., & Veermäe, H. 2020, Phys. Rev. D, 101, 103023 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2017, SIAM Rev., 59, 65 [Google Scholar]
  5. Bhattacharya, S., Rao, K. K., Agarwal, M., Balan, S., & Vaidya, K. 2022, MNRAS, 517, 3525 [CrossRef] [Google Scholar]
  6. Binney, J. 2013, in Secular Evolution of Galaxies, eds. J. Falcón-Barroso, & J. H. Knapen, 259 [CrossRef] [Google Scholar]
  7. Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton: Princeton University Press) [Google Scholar]
  8. Brinks, E., & Portas, A. 2017, in Formation and Evolution of Galaxy Outskirts, eds. A. Gil de Paz, J. H. Knapen, & J. C. Lee, 321, 235 [NASA ADS] [Google Scholar]
  9. Brinks, E., Bagetakos, I., Walter, F., & de Blok, E. 2007, in Triggered Star Formation in a Turbulent ISM, eds. B. G. Elmegreen, & J. Palous, 237, 76 [NASA ADS] [Google Scholar]
  10. Byrne, R. A., & Fraser, M. 2022, MNRAS, 514, 1188 [CrossRef] [Google Scholar]
  11. Chandrasekhar, S. 1943, ApJ, 97, 255 [Google Scholar]
  12. Cowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135 [NASA ADS] [CrossRef] [Google Scholar]
  13. Cranmer, M. 2023, ArXiv e-prints [arXiv:2305.01582] [Google Scholar]
  14. de los Reyes, M. A. C., & Kennicutt, R. C. 2019, ApJ, 872, 16 [NASA ADS] [CrossRef] [Google Scholar]
  15. Della Croce, A., Dalessandro, E., Livernois, A., & Vesperini, E. 2024, A&A, 683, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Efron, B. 1979, Ann. Stat., 7, 1 [Google Scholar]
  17. Ehlerová, S., & Palouš, J. 2013, A&A, 550, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Finn, M. K., Johnson, K. E., Indebetouw, R., et al. 2024, ArXiv e-prints [arXiv:2401.01451] [Google Scholar]
  19. García-Conde, B., Antoja, T., Roca-Fàbrega, S., et al. 2024, A&A, 683, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Girichidis, P., Walch, S., Naab, T., et al. 2016, MNRAS, 456, 3432 [NASA ADS] [CrossRef] [Google Scholar]
  21. Goldbaum, N. J., Krumholz, M. R., & Forbes, J. C. 2015, ApJ, 814, 131 [NASA ADS] [CrossRef] [Google Scholar]
  22. Goldbaum, N. J., Krumholz, M. R., & Forbes, J. C. 2016, ApJ, 827, 28 [NASA ADS] [CrossRef] [Google Scholar]
  23. Herpich, J., Stinson, G. S., Dutton, A. A., et al. 2015, MNRAS, 448, L99 [CrossRef] [Google Scholar]
  24. Huško, F., Lacey, C. G., & Baugh, C. M. 2023, MNRAS, 518, 5323 [Google Scholar]
  25. Jeans, J. H. 1915, MNRAS, 76, 70 [Google Scholar]
  26. Joshi, R., & Widrow, L. M. 2023, ArXiv e-prints [arXiv:2312.01206] [Google Scholar]
  27. Kataria, S. K., & Shen, J. 2022, ApJ, 940, 175 [NASA ADS] [CrossRef] [Google Scholar]
  28. Keller, B. W. 2017, Superbubble Feedback in Galaxy Evolution, https://api.semanticscholar.org/CorpusID:125997457 [Google Scholar]
  29. Keller, B. W., Wadsley, J., Benincasa, S. M., & Couchman, H. M. P. 2014, MNRAS, 442, 3013 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kim, C.-G., & Ostriker, E. C. 2017, ApJ, 846, 133 [NASA ADS] [CrossRef] [Google Scholar]
  31. Kim, C.-G., Ostriker, E. C., & Raileanu, R. 2017, ApJ, 834, 25 [Google Scholar]
  32. Kinugawa, T., Horiuchi, S., Takiwaki, T., & Kotake, K. 2024, MNRAS, 532, 3926 [NASA ADS] [CrossRef] [Google Scholar]
  33. Kipper, R., Tenjes, P., Hütsi, G., Tuvikene, T., & Tempel, E. 2019, MNRAS, 486, 5924 [NASA ADS] [CrossRef] [Google Scholar]
  34. Kipper, R., Benito, M., Tenjes, P., Tempel, E., & de Propris, R. 2020, MNRAS, 498, 1080 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kipper, R., Tamm, A., Tempel, E., de Propris, R., & Ganeshaiah Veena, P. 2021, A&A, 647, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Kipper, R., Tenjes, P., Benito, M., et al. 2023, A&A, 680, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Kormendy, J. 2013, in Secular Evolution of Galaxies, eds. J. Falcón-Barroso, & J. H. Knapen, 1 [Google Scholar]
  38. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  39. Kuzmin, G. G. 1961, Publications of the Tartu Astrofizica Observatory, 33, 351 [NASA ADS] [Google Scholar]
  40. Kuzmin, G. G. 1963, Publications of the Tartu Astrofizica Observatory, 34, 9 [NASA ADS] [Google Scholar]
  41. Laplace, P. S. 1774, Memoir on the Probability of Causes of Events, 621 [Google Scholar]
  42. Lara-López, M. A., Pilyugin, L. S., Zaragoza-Cardiel, J., et al. 2023, A&A, 669, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. Li, X., Shlosman, I., Heller, C., & Pfenniger, D. 2023, MNRAS, 526, 1972 [NASA ADS] [CrossRef] [Google Scholar]
  44. Lian, J., Zasowski, G., Hasselquist, S., et al. 2022, MNRAS, 511, 5639 [NASA ADS] [CrossRef] [Google Scholar]
  45. Mac Low, M.-M., & McCray, R. 1988, ApJ, 324, 776 [NASA ADS] [CrossRef] [Google Scholar]
  46. Maxwell, A. J., Wadsley, J., & Couchman, H. M. P. 2015, ApJ, 806, 229 [NASA ADS] [CrossRef] [Google Scholar]
  47. Mayya, Y. D., Alzate, J. A., Lomelí-Núez, L., et al. 2023, MNRAS, 521, 5492 [CrossRef] [Google Scholar]
  48. McKee, C. F., Parravano, A., & Hollenbach, D. J. 2015, ApJ, 814, 13 [Google Scholar]
  49. Minchev, I., & Famaey, B. 2010, ApJ, 722, 112 [Google Scholar]
  50. Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533 [NASA ADS] [Google Scholar]
  51. Nath, B. B., Das, P., & Oey, M. S. 2020, MNRAS, 493, 1034 [NASA ADS] [CrossRef] [Google Scholar]
  52. Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563 [Google Scholar]
  53. Nuñez-Castiñeyra, A., Nezri, E., Mollitor, P., Devriendt, J., & Teyssier, R. 2023, JCAP, 2023, 012 [CrossRef] [Google Scholar]
  54. Patil, A. A., Bovy, J., Jaimungal, S., Frankel, N., & Leung, H. W. 2023, MNRAS, 526, 1997 [NASA ADS] [CrossRef] [Google Scholar]
  55. Reynolds, T. M., Fraser, M., & Gilmore, G. 2015, MNRAS, 453, 2885 [Google Scholar]
  56. Ryden, B., & Pogge, R. 2021, Interstellar and Intergalactic Medium (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
  57. Sarbadhicary, S. K., Wagner, J., Koch, E. W., et al. 2023, ArXiv e-prints [arXiv:2310.17694] [Google Scholar]
  58. Schönrich, R., & Binney, J. 2009, MNRAS, 396, 203 [Google Scholar]
  59. Sellwood, J. A., & Binney, J. J. 2002, MNRAS, 336, 785 [Google Scholar]
  60. Struck, C., & Elmegreen, B. G. 2017, MNRAS, 469, 1157 [NASA ADS] [CrossRef] [Google Scholar]
  61. van der Wel, A., Chang, Y.-Y., Bell, E. F., et al. 2014, ApJ, 792, L6 [NASA ADS] [CrossRef] [Google Scholar]
  62. Vera-Ciro, C., D’Onghia, E., Navarro, J., & Abadi, M. 2014, ApJ, 794, 173 [NASA ADS] [CrossRef] [Google Scholar]
  63. Vihola, M. 2020, Wiley StatsRef: Statistics Reference Online, 1 [Google Scholar]
  64. Walter, F., Brinks, E., de Blok, W. J. G., et al. 2008, AJ, 136, 2563 [Google Scholar]
  65. Watkins, E. J., Barnes, A. T., Henny, K., et al. 2023, ApJ, 944, L24 [NASA ADS] [CrossRef] [Google Scholar]
  66. Wu, J., Struck, C., D’Onghia, E., & Elmegreen, B. G. 2020, MNRAS, 499, 2672 [NASA ADS] [CrossRef] [Google Scholar]
  67. Yadav, N., Mukherjee, D., Sharma, P., & Nath, B. B. 2017, MNRAS, 465, 1720 [NASA ADS] [CrossRef] [Google Scholar]
  68. Zel’dovich, Y. B., & Raizer, Y. P. 1967, Physics of Shock Waves and High-temperature Hydrodynamic Phenomena (New York: Academic Press) [Google Scholar]

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πr3ρ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 v2 ≫ ISMR2, 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

v ˙ G M ( < r ) r 2 cos ϕ = 4 π G ρ ISM b 3 , $$ \begin{aligned} \dot{v}_{\perp } \approx \frac{G M( < r)}{r^2} \cos {\phi } = \frac{4\pi G \rho _{\rm ISM} b}{3}, \end{aligned} $$(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ϕ/cos2ϕ and observing that r < R requires R 2 b 2 / b < tan ϕ < R 2 b 2 / b $ -\sqrt{R^2 - b^2}/b < \tan \phi < \sqrt{R^2 - b^2}/b $, one obtains

Δ v = ϕ min ϕ max 4 π G ρ ISM b 3 b d ϕ v cos 2 ϕ = 4 π G ρ ISM b 2 3 v tan ϕ | R 2 b 2 / b R 2 b 2 / b = 8 π G ρ ISM 3 v b R 2 b 2 . $$ \begin{aligned} \Delta {v}_{\perp }&= \int _{\phi _{\rm min}}^{\phi _{\rm max}} \frac{4\pi G \rho _{\rm ISM} b}{3} \frac{b d\phi }{v\cos ^2\phi } = \left. \frac{4\pi G \rho _{\rm ISM} b^2}{3 v} \tan \phi \right|_{-\sqrt{R^2 - b^2}/b}^{\sqrt{R^2 - b^2}/b}\nonumber \\&= \frac{8\pi G \rho _{\rm ISM}}{3 v} b \sqrt{R^2 - b^2}. \end{aligned} $$(A.2)

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

Δ v = v 2 v 2 v v 2 ( Δ v v ) 2 = 32 π 2 G 2 ρ ISM 2 9 v 3 b 2 ( R 2 b 2 ) . $$ \begin{aligned} \Delta v_{\parallel } = \sqrt{v^2 - v_{\perp }^2} - v \approx -\frac{v}{2} \left(\frac{\Delta v_{\perp }}{v}\right)^2 = -\frac{32\pi ^2 G^2\rho _{\rm ISM}^2}{9 v^3} b^2 \left(R^2 - b^2\right). \end{aligned} $$(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:

F , coll = Δ p n v 2 π b d b = 64 π 3 G 2 ρ ISM 2 9 v 2 ρ 0 R b 3 ( R 2 b 2 ) d b = 16 π 3 G 2 ρ ISM 2 ρ R 6 27 v 2 . $$ \begin{aligned} F_{\parallel , \mathrm {coll}}&= \int \Delta p_{\parallel } \,\, nv \, 2\pi b db = -\frac{64\pi ^3 G^2\rho _{\rm ISM}^2}{9 v^2} \rho \int _0^R b^3 \left(R^2 - b^2\right) \,db \nonumber \\&= -\frac{16\pi ^3 G^2 \rho _{\rm ISM}^2 \rho R^6}{27 v^2}. \end{aligned} $$(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

ρ m f ( v ) d 3 v , $$ \begin{aligned} \rho \rightarrow m \int f({v}) d^3 v, \end{aligned} $$(A.5)

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

F = 16 π 3 G 2 ρ ISM 2 R 6 27 m f ( v ) v 2 ω d 3 v , $$ \begin{aligned} {F} = -\frac{16\pi ^3 G^2 \rho _{\rm ISM}^2 R^6}{27} m \int \frac{f({v})}{v^2} {\omega } \, d^3 v, \end{aligned} $$(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 v0 relative to the perturber (cavity). By symmetry, the dynamical friction force must be parallel to v0, hence we consider only the Fv0 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 d3v, one obtains

F = 16 π 3 G 2 ρ ISM 2 R 6 27 m f ( v ) v 2 cos θ d 3 v , $$ \begin{aligned} F = -\frac{16\pi ^3 G^2 \rho _{\rm ISM}^2 R^6}{27} m \int \frac{f^{\prime }(v^{\prime })}{v^2} \cos \theta \, d^3 v^{\prime }, \end{aligned} $$(A.7)

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

F = 16 π 3 G 2 ρ ISM 2 R 6 27 m f ( v ) v 2 d v × 1 + q cos θ v 0 2 ( 1 + q 2 + 2 q cos θ ) 3 / 2 d ϕ d cos θ . $$ \begin{aligned} F&= -\frac{16\pi ^3 G^2 \rho _{\rm ISM}^2 R^6}{27} m \int f^{\prime }(v^{\prime }) \, v^{\prime \,2} dv^{\prime }\nonumber \\&\times \int \frac{1 + q \cos \theta ^{\prime }}{v_0^2 (1 + q^2 + 2q \cos \theta ^{\prime })^{3/2}} \, d\phi ^{\prime } d\cos \theta ^{\prime }. \end{aligned} $$(A.8)

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

F = 64 π 3 G 2 ρ ISM 2 R 6 m 27 v 0 2 f ( v ) v 2 d v . $$ \begin{aligned} F&= -\frac{64\pi ^3 G^2 \rho _{\rm ISM}^2 R^6 \,m}{27 v_0^2} \int f^{\prime }(v^{\prime }) \, v^{\prime \,2} dv^{\prime }. \end{aligned} $$(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

d v r dt = d dt ( v · r r ) = F r + L 2 r 3 , $$ \begin{aligned} \frac{dv_r}{dt} = \frac{d}{dt} \left({v}\cdot \frac{{r}}{r}\right) = F_r + \frac{L^2}{r^3}, \end{aligned} $$(A.10)

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

d dr ( v r 2 2 + Φ + L 2 2 r 2 ) = 0 , $$ \begin{aligned} \frac{d}{dr}\left(\frac{v_r^2}{2} + \Phi + \frac{L^2}{2r^2}\right) = 0, \end{aligned} $$(A.11)

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

For the cavity we can write

Φ ( r ) = 2 π G ρ ISM 3 r 2 = 1 4 v ch 2 r 2 R enter 2 , $$ \begin{aligned} \Phi (r) = -\frac{2\pi G \rho _{\rm ISM}}{3} r^2 = -\frac{1}{4} v_{\rm ch}^2 \frac{r^2}{R_{\rm enter}^2}, \end{aligned} $$(A.12)

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

v ch = 8 π G ρ ISM R enter 2 3 . $$ \begin{aligned} v_{\rm ch} = \sqrt{\frac{8\pi G \rho _{\rm ISM} R_{\rm enter}^2}{3}}. \end{aligned} $$(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,

v out 2 = v in 2 + v ch 2 2 [ ( R exit R enter ) 2 1 ] . $$ \begin{aligned} v_{\rm out}^2 = v_{\rm in}^2 + \frac{v_{\rm ch}^2}{2} \left[\left(\frac{R_{\rm exit}}{R_{\rm enter}}\right)^2 - 1\right]. \end{aligned} $$(A.14)

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

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 Renter = rb(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)

d F in = ( v · v exp v exp ) f ( v ) d 3 v , $$ \begin{aligned} dF_{\rm in} = -\left({v}^{\dagger }\cdot \frac{{v}_{\rm exp}}{v_{\rm exp}} \right) f^{\dagger }({v}^{\dagger }) d^3 v^{\dagger }, \end{aligned} $$(A.15)

where v = vvexp is the particle velocity in the surface element frame, v is its velocity in the bubble frame (i.e. its centre), and vexp is the local velocity of the surface element in the same frame. The scalar product in Eq. (A.15) can be written as −v · vexp/vexp = −|v| cos θ + vexp > 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

d F in = ( v exp | v | cos θ ) f ( v ) d 3 v , where cos θ < v exp | v | , $$ \begin{aligned} dF_{\rm in} = (v_{\rm exp} -|v|\cos \theta ) \, f({v}) d^3 v, \quad \mathrm{where} \quad \cos \theta < \frac{v_{\rm exp}}{|v|}, \end{aligned} $$(A.16)

and we have again used the invariance of f and d3v 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

m f ( v ) d v = ρ $$ \begin{aligned} m\int f({v})\mathrm{d}v = \rho ^* \end{aligned} $$(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 vb relative to this frame. The total rate of momentum transfer between the bubble and the particle population can then be written as

dP dt = m Δ v cos θ v d F in d S b . $$ \begin{aligned} \frac{dP}{dt} = \int m\Delta v_{\parallel } \cos \theta _v \, dF_{\rm in} \, dS_{\rm b}. \end{aligned} $$(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 vb measured in the bubble frame, and dFin is the particle flux defined by Eq. (A.16). The bubble surface element is

d S b = R enter 2 d cos θ d ϕ , $$ \begin{aligned} dS_{\rm b} = R_{\rm enter}^2 d\cos \theta \, d\phi , \end{aligned} $$(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 ( v ) d 3 v = f ( v ) d 3 v = f ( v ) v 2 d v d cos θ v d ϕ v $ f(\boldsymbol{v}) d^3 v = f\prime(\boldsymbol{v}\prime) d^3 v\prime = f\prime(v\prime) \, v^{\prime\,2} dv\prime \, d\cos\theta_v\prime \, d\phi_v\prime $, where θv′ denotes the angle between v′ and vb. 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 β b = v b / v $ \beta_{\mathrm{b}}\prime = v_{\mathrm{b}}/v\prime $, one obtains

cos θ v = cos θ v β b 1 2 cos θ v β b + β b 2 , v = v 1 2 cos θ v β b + β b 2 . $$ \begin{aligned} \cos \theta _v = \frac{\cos \theta _v^{\prime } - \beta _{\rm b}^{\prime }}{\sqrt{1 - 2\cos \theta _v^{\prime }\beta _{\rm b}^{\prime } + \beta _{\rm b}^{\prime \,2}}}, \quad v = v^{\prime } \sqrt{1 - 2\cos \theta _v^{\prime }\beta _{\rm b}^{\prime } + \beta _{\rm b}^{\prime \,2}}. \end{aligned} $$(A.20)

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

dP dt = 4 π 2 m R enter 2 f ( v ) v 3 d v 1 1 ( cos θ v β b ) d cos θ v × D θ ( β cos θ ) Δ v d cos θ , $$ \begin{aligned} \frac{dP}{dt}&= 4\pi ^2 m R_{\rm enter}^2 \int _{-\infty }^{\infty } f^{\prime }(v^{\prime }) \, v^{\prime \,3} dv^{\prime } \int _{-1}^1 (\cos \theta _v^{\prime } - \beta _{\rm b}^{\prime }) \, d\cos \theta _v^{\prime } \nonumber \\&\times \int _{D_{\theta }} (\beta - \cos \theta ) \, \Delta v_{\parallel } \, d\cos \theta , \end{aligned} $$(A.21)

where

D θ = { cos θ : cos θ 1 , cos θ min ( 1 , β ) } $$ \begin{aligned} D_{\theta } = \{\cos \theta : \cos \theta \ge -1 , \cos \theta \le \min (1,\beta )\} \end{aligned} $$(A.22)

and β ≡ vexp/|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 tmax (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 (Ri and Ro), 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 (vcol), but in the current analysis we split it into two separate speeds: inner (still vcol) and outer (vo, col) shell. To understand how these changes may influence we remind that Ri determines what is the mass that is displaced to form the acceleration field correction and Ro effectively determines what is the shape of the acceleration field. When we keep vcol 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 environment6. Figure B.1 shows what is the introduced de-acceleration changes when varying vo, 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 Ri collapse speed. We kept the Ro, col as constant and calculated the de-acceleration dependence on the vcol. 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.

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

thumbnail 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 tmax with a variable that is 1 at tmax and will drop linearly to 0 upon the end of SB life tb. So between tmaxtb we have acceleration of SB as

a = a fid ( t = t max ) t b t t b t max . $$ \begin{aligned} \mathbf a = \mathbf{a }^\mathrm{fid}(\mathrm{t = t_{\rm max}}) \frac{t_{\rm b} - t}{t_{\rm b} - t_{\rm max}}. \end{aligned} $$(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.

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

thumbnail 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

Table 1.

Parameters for determining the level of the influence to the gas slab by the stars via SBs.

All Figures

thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail Fig. 7.

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

In the text
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.