Size distribution of superbubbles
MMon. Not. R. Astron. Soc. , 000–000 (0000) Printed 23 April 2020 (MN L A TEX style file v2.2)
Size distribution of superbubbles
Biman B. Nath , Pushpita Das , , M. S. Oey Raman Research Institute, Sadashiva Nagar, Bangalore 560080, India Indian Institute of Science Education & Research Kolkata, Mohanpur 741246, India Dept of Astronomy, University of Michigan, 323 West Hall, 1085 South University Ave., Ann Arbor, MI 48109-1107, USA
23 April 2020
ABSTRACT
We consider the size distribution of superbubbles in a star forming galaxy. Previousstudies have tried to explain the distribution by using adiabatic self-similar evolutionof wind driven bubbles, assuming that bubbles stall when pressure equilibrium isreached. We show, with the help of hydrodynamical numerical simulations, that thisassumption is not valid. We also include radiative cooling of shells. In order to take intoaccount non-thermal pressure in the ambient medium, we assume an equivalent highertemperature than implied by thermal pressure alone. Assuming that bubbles stall whenthe outer shock speed becomes comparable to the ambient sound speed (which includesnon-thermal components), we recover the size distribution with a slope of ∼ − . Key words:
ISM:bubbles – HII regions – ISM: structure – galaxies: ISM
The distribution of sizes of superbubbles created by stellarwinds and supernovae of OB associations can be a diagnosticof the star formation process in a galaxy. The seminal pa-per by Oey & Clarke (1997) [hereafter OC97] showed thatif the mechanical luminosity function of OB associations isdescribed by a power law , φ ( L ) ∝ L − β , then for a con-stant star formation rate, the differential size distribution ofsuperbubbles is given by N ( R ) ∝ R − β , for bubbles whoseevolution is dominated by ambient pressure and have stalled.For typical parameters, OC97 estimated this stalling radiusto be (cid:54) − β makes ita useful diagnostic of the star formation process in a galaxy.This predicted distribution has been confirmed by The HINearby Galaxy Survey , which studied the properties of HIholes in nearby galaxies (Bagetakos et al. 2011). They foundthat the size distribution has a power-law index of ∼ − . β ≈
2, which in turn is consistent with inde-pendent observations of OB association (McKee & Williams1997).The robustness of this predicted size distribution, how-ever, raises the question whether or not it depends on pa-rameters such as ambient pressure, density and so on, andif so, how. Since it was derived assuming that ambient pres- sure does not affect the adiabatic expansion of a stellar windbubble (Weaver et al. 1977) other than to impose a stallcriterion, it is not possible to answer these questions with-out going beyond these assumptions. Although the averagedistribution in
THINGS came out to have a power law of ∼ − .
9, individual galaxies had distributions whose power-law ranged between − −
4. Early-type spiral galax-ies showed a steeper slope than late-type spirals and dwarfgalaxies. Bagetakos et al. (2011) explained this by invokingthe fact that scale heights of disks in early-type spirals aresmaller than in late-type spirals, and are easy to break outof. This would limit the size of the largest holes, and con-sequently steepen the size distribution. However, one couldalso envisage a scenario in which a truncated or steepenedOB association luminosity function in early type galaxiesproduce a steep size distribution of HI holes, as is known toalso manifest in the HII region luminosity function (Oey &Clarke 1998). In order to disentangle the effects, one wouldhave to calculate the size distribution beyond the assump-tions of adiabaticity, which is inherent in the self-similarevolution of bubbles. And then the question remains as tohow the power-law distribution with index ∼ − c (cid:13) a r X i v : . [ a s t r o - ph . GA ] A p r B. B. Nath, P. Das, M. S. Oey the size distribution. We assume a constant star formationrate for simplicity. The paper is structured as follows. Wefirst review the results of OC97 in light of Monte-Carlo sim-ulations in §
2. Then we discuss the results when ambientpressure is included in the dynamics of superbubbles, in § §
4, we discuss the effects of cooling and present the resultsin §
5. The implications are discussed in § We first review the essence of OC97 calculations. They as-sumed that the mechanical luminosity function of OB as-sociation is given by φ ( L ) ∝ L − β . The luminosity of eachcluster is assumed to be constant in time, until the lowestmass SN progenitors expire at t e ∼
40 Myr. The wind bubbleis assumed to evolve in a self-similar manner, so that the ra-dius scales as, R ∝ L / t / (Weaver et al. 1977). Note thatthis assumes adiabatic expansion of the bubble. The pres-sure inside then evolves as p i ∝ L / t − / . If the ambientpressure is p , then these bubbles are assumed to stall when p = p i . This implies a stalling time, for the bubble to reacha final radius, of t f ∝ L / . This in turn leads to a scalingof the final stalling radius R f ∝ L / t / f ∝ L / .If the superbubbles are produced at a constant rate ψ ,then after time t , the number of bubbles with radii in therange R to R + dR will be dependent on φ ( L ) dL , where L and dL are the luminosity and range corresponding to thisinterval in radius R . Therefore the differential size distribu-tion will be given by, N ( R ) ∝ ψφ ( L ) (cid:16) ∂R f ∂L (cid:17) − ∝ L − β +1 / ∝ R − β , (1)where the last proportionality follows from R f ∝ L / forstalled final radii derived above. Therefore the above scal-ings imply a size distribution N ( R ) ∝ R − β , for stalledbubbles. In case of expanding bubbles in the momentumconserving phase, one again has R ∝ L / t / . In this case,OC97 considered stalling of bubbles when bubble speed be-comes comparable to ambient sound speed (or, equivalently,ambient pressure being comparable to ram pressure of theouter shock). This again leads to t f ∝ L / , and R f ∝ L / ,leading to N ( R ) ∝ R − β . They showed that stalled bub-bles dominate the size distributions. They also estimatedthat for typical ISM parameters, the largest size of stalledbubbles is ∼ (cid:12) yr − . The luminosity function ofOB association is taken as φ ( L ) = AL − β , with β = 2, asinferred from HII region luminosity function (Oey & Clarke1998; McKee & Williams 1997). The mechanical luminosityis related to the mass of the cluster in the following way.Given a Kroupa initial mass function (Kroupa 2002), thereis one OB star for 100 M (cid:12) cluster mass, and each OB starcan be assumed to give 10 erg. The average luminosity ofthe cluster (given a lifetime of ∼
37 Myr, corresponding to the main sequence life time of a 8 M (cid:12) star), is L ≈ × erg s − (cid:16) M cl M (cid:12) (cid:17) . (2)This is very close (within a factor of a few) to the result ofmechanical luminosity of star clusters including the effect ofstellar winds and supernovae, as calculated by Leitherer etal. (1999). Therefore the mass function of the clusters canbe written as dNdM cl = A (cid:16) M cl M (cid:12) (cid:17) − β . (3)Assuming the minimum and maximum mass of clusters tobe 100 and 10 M (cid:12) respectively, the average cluster mass isfound to be ∼ (cid:12) , corresponding to an average me-chanical luminosity of ∼ . × erg s − . The luminosityrange is ∼ × –9 × erg s − .The size distribution of bubbles is shown in Figure 1,for different times after the onset of star formation, for bub-bles growing in a self-similar manner as mentioned above.The ambient pressure is assumed to be P = 2 . × − dynes cm − as was considered by OC97. As expected fromthe arguments in OC97, the dominant slope of the distribu-tion is roughly 1 − β = −
3, for radii below the maximumsize of stalled bubbles at a given epoch. OC97 estimated thestalling timescale and the stalled sizes. Initially, low lumi-nosity bubbles stall, and at a given epoch, bubbles up to acertain size stall, beyond which the bubbles are in a grow-ing phase. We show this upper limit of stalled sizes withsquare boxes on the red and blue lines, corresponding to 5and 10 Myr (using equation 31 of OC97). (At later times,the biggest stalled bubbles cross the limit of sizes consideredby us.) Figure 1 shows that the fitted power-laws until thismaximum size has an index close to −
3. However, there isalso another class of small bubbles that are growing at anygiven epoch, arising from either massive but young clustersor low mass clusters.We assume the bubbles to start forming at a continuousrate at t = 0 = t l . The inset in Figure 1 plots the size distri-bution at 30 Myr with smaller bin size, to clearly show thecase of small and growing bubbles. These bubbles show adistribution that scales as R / as predicted by OC97 (theireqn 39). These bubbles are dominated by clusters which havenot yet reached stalling phase. In this case, R ∝ L / t / ,which gives ∂R∂L ∝ L − / t / . This distribution is then in-tegrated from t l until a time, t u , for which the minimumluminosity cluster is still growing. This time scale is givenby t u ∝ R / L − / (from self-similar evolution). Therefore,the size distribution of growing bubble is given by N ( R ) ∝ (cid:90) t u L − β ( ∂R∂L ) − dt ∝ R − β t − βu ∝ R / (4)Our Monte-Carlo runs for self-similar case show these grow-ing bubbles with a size distribution as predicted.The increasing ( R / ) and decreasing R − ) regimes ofthe size distribution give rise to a peak in the size distribu-tion, and the peak shifts towards large radii with time, untilthe time when the oldest population of minimum luminosityclusters has achieved stalling radius. In the case depicted inFigure 1 the minimum luminosity cluster achieves this statusat ∼ c (cid:13) , 000–000 uperbubble size distribution Figure 1.
Size distribution of self-similar superbubbles for a con-stant star formation rate at different epochs, for an ambient pres-sure p = 2 . × − dynes cm − , and ambient temperature T amb = 4 × K. The square boxes on the blue and red linesshow the values of the maximum stall radius at 5 and 10 Myr.The distribution at each epoch has been fitted with a power-lawbelow the maximum size shown by the square boxes, and the fit-ted slope are all roughly close to −
3, below the maximum stallradius. The inset shows the distribution at 30 Myr in detail (withsmaller bin size), which has the R / rising part for small andgrowing bubbles (the straight line shows a power-law of index2 / small or negligible from an observational point of view. Forexample, after 5 Myr, the peak is found to be at ∼
40 pc.We note here that the assumption of constant lumi-nosity with time for low mass clusters is justified becausethe mechanical luminosity not only arises from supernovae,which in the case of low mass clusters would be few andfar between, but also from stellar winds of massive stars.Since the typical mechanical power of stellar winds from OBstars is ∼ erg s − (Seo et al. 2018), which is also coin-cidentally comparable to (10 erg /
35 Myr), the mechanicalpower of OB associations before and after the first super-novae events differ at the most by a factor of 2. This is alsoborne out by the estimates of mechanical power using Star-burst99 (Leitherer et al. 1999).
When we consider the growth of bubbles beyond the self-similar evolution, we need to fully account for the effect ofambient pressure and radiative cooling. The effect of radia-tive cooling cannot be included without resorting to hydro-dynamical simulation. However, the effect of ambient pres-sure can be calculated numerically, and before going to thefull solution of hydrodynamical simulation, we will discussthis effect next.The wind bubble has four regions: a free wind region inwhich energy and mass are injected, a shocked wind region, acontact discontinuity and a shocked ISM region. The radiusof the contact discontinuity evolves in the presence of anambient pressure p , as, (equations 54 and 56 in Weaver et Figure 2.
Evolution of bubbles for ambient pressure p =3 × − dyne/cm and p = 5 × − dyne/cm , for averageluminosity L = 1 . × erg s − . Figure 3.
Evolution of pressure inside bubbles for the same casesas in Figure 2 al. (1977)) ddt (cid:18) π R ρ (cid:19) = 4 πR ( p − p ) ,dpdt = 23 L (4 π/ R − p (cid:18) R dRdt (cid:19) (5)Note that OC97 did not include p in this relation, and onlyused p to as a criterion for stalling bubble growth. We con-tinue with the assumption of OC97 that the bubbles stallwhen the interior pressure becomes equal to the ambientpressure. The results of eqn 5 show that the interior pres-sure decreases less rapidly than that assumed in OC97, andconsequently, the bubble size grows to larger values thanare admitted in self-similar evolution. We show two exam-ples in Figure 2, where the evolution of the bubble sizes inthe self-similar case is compared with that determined fromequations 5, for two values of ambient pressure. The curvesare for an average luminosity (given the above mentioned lu-minosity function), L av = 1 . × erg s − . We also showin Figure 3 the evolution of interior pressure in these twocases. The curves show that the interior pressure evolvesslower than in the case of self-similar case, as mentionedabove.The corresponding distribution will be steeper than theself-similar case. This is shown in Figure 4 for three differentepochs. The slopes of the distribution is close to ∼ −
4, inthis particular example, much steeper than − c (cid:13) , 000–000 B. B. Nath, P. Das, M. S. Oey
Figure 4.
The size distribution of bubbles for adiabatic evo-lution for an ambient pressure p = 2 . × − dyne cm − ,at 10 , ,
30 Myr. The distributions have a slope ∼ −
4, muchsteeper slope than − bles grow to larger sizes than before. This will make the sizedistribution steeper than the self-similar case.The main reason for the steepness is the fact that grow-ing bubbles dominate the population of superbubbles at anytime, since the bubbles tend to grow for a longer time in thiscase. In the case of domination by growing bubbles, suppos-ing a power-law scaling relation between radius and lumi-nosity as R ∝ L x , the slope of the size distribution is givenby N ( R ) ∝ φ ( L ) (cid:16) ∂R∂L (cid:17) − ∝ R − x − βx . (6)For example, for self-similar solution of Weaver et al. (1977), x = 1 / N ( R ) ∝ R − β , as shown by OC97. In ourcase, since the power-law index is close to − .
1, it implies x ≈ / .
1, although there is no simple analytical scalingrelation for radius and luminosity for this case. We haveconfirmed this from our numerical results.However, the assumption of bubble stalling when pres-sure equilibrium is reached is not valid, since the bubble shellcontinues to move outward because of its momentum. Thiswill continue to make the bubbles bigger in size until theouter shock speed becomes comparable to the sound speedof the ambient medium. This can be demonstrated with thehelp of numerical hydrodynamical simulations, which we de-scribe below.
Here we include the effect of radiative cooling in the evolu-tion of superbubbles. We have used
PLUTO for 1-D numer-ical hydrodynamical calculations (Mignone et al. 2007). Wehave solved the following equations: ∂ρ∂t + (cid:79) · ( ρ v ) = S m ∂ ( ρ v ) ∂t + (cid:79) · ( ρ v ⊗ v ) = − (cid:79) p∂e∂t + (cid:79) · [( e + p ) v ] = S e − q − (7) Figure 5.
Evolution of radius bubbles from hydrodynamical sim-ulations is compared to the self-similar case. The upper twopanels are for L = 2 × erg s − and lower panels, for L = 10 erg s − . The left panels show the case for ambientpressure p = 2 . × − dynes cm − , and the right panels,for p = 1 . × − dynes cm − . Red curves plot the adiabaticself-similar evolution (assuming stalling at pressure equilibrium).Blue curves show the results of using eqn 5 (assuming stalling atpressure equilibrium). Cyan curves show the results of numeri-cal simulation without cooling, and green curves show the resultsof simulation with cooling, assuming stalling when outer shockspeed equals ambient sound speed. Here, total energy e = ρ ( (cid:15) + v ), (cid:15) is specific energy, ρ ismass density, p is pressure and v is the fluid velocity. Theterms S m and S e correspond to the mass-loss rate ˙ M/V andinput mechanical energy (
L/V ) respectively, which are re-lated to each other by the wind speed, v w , as, L = 0 . Mv w ,and we assume v w = 2000 km s − (Chevalier & Clegg 1985).We have introduced the source terms as per the model ofChevalier & Clegg (1985). We have kept the ambient parti-cle density ( n amb ) and ambient temperature( T amb ) constantfor each run. Equations 7 have been solved in 1D sphericalcoordinate using the HLLC solver (CFL Number = 0.3)Since the appropriate medium for embedding the su-perbubbles is the warm neutral medium of ISM, we use n amb = 0 . − for all our runs (Wolfire et al. 2003).However, we consider the case of different equivalent ISMpressures by assuming different T amb . It is known that non-thermal pressures in ISM, especially in neutral medium thatwe are concerned with here, can be substantial (Elmegreen& Scalo 2004) and even more than the thermal pressure.For example, in the solar neighbourhood, Jenkins & Tripp(2011) estimated that turbulent and magnetic pressures arecomparable in magnitude, and are roughly three times largerthan the thermal pressure. In order to take into accountthe effect of non-thermal pressure, we use a correspond-ing equivalent ambient temperature, keeping n amb fixed. Weuse T amb = 4 × K, for corresponding ISM pressures p = 2 . × − dynes cm − for our fiducial runs, but have c (cid:13) , 000–000 uperbubble size distribution also varied the value of T amb in order to study its influenceon size distribution.Mass and energy are continuously injected within asmall radius r c . We use r c = 1 pc, in order to minimize non-physical cooling losses at the early stages of shock formation(see eqn 4 in Sharma et al. (2014)). According to this cri-terion, for the lowest luminosity considered here ( L = 10 erg s − ), r c should be less than (cid:54) . S m = ˙ M/V c and S e = L/V c are non-zero for r < r c ,where V c = π r c . In the last equation for energy conserva-tion, q − = n i n e Λ( T ), Λ( T ) being the cooling function. Wehave used a tabulated cooling function for solar metallicity.We turn off cooling when temperature comes down to 10 K,to mimic the effect of photoionization heating in the bubble.We show the results of simulations with and withoutcooling in Figure 5. The luminosities used here is a low valueof L = 2 × erg s − (top panels), and a high luminosityof L = 10 erg s − (bottom panels) . The left panels are for T amb = 4 × K, and the right panels, for T amb = 2 × K.The self-similar case is show in red, and the result of usingeqn 5 is shown for comparison in dark blue, although withoutthe condition of stalling at equal pressure. The results ofsimulations with and without cooling are shown in greenand cyan, respectively.We first notice that in the case of no radiative cooling(cyan), the evolution of the shell is roughly similar to thatof eqn 5, except for high mechanical luminosity, in whichcase the PLUTO runs give a slightly smaller radius (Figure5). This is because of the fact that eqns 5 do not take intoaccount mass injection, which is higher in our calculationsfor higher luminosities (since the wind speed is consideredequal in all cases). The injected mass increases the inertiaof the shell and decelerates it to some extent.Secondly, increasing the ambient pressure (without in-creasing the gas density) makes the bubbles smaller, as ex-pected physically.
Let us first discuss the size distribution in the adiabatic case.As mentioned above, the results of the simulations confirmour analytical results in §
3, where the pressure gradient termis included in the dynamics of bubbles. We had seen in thatcase that size distribution is steeper than −
3, because ofdomination of growing bubbles. This is confirmed by thesimulation results in which cooling is turned off, and we get asize distribution for the fiducial case with a power-law index ≈ −
4. The scatter plot of bubbles in this case is shown inFigure 6 as a function of mechanical luminosity (up to 10 erg s − ) and bubble age, for a snapshot at 30 Myr. The sizeof the bubbles are shown in different colours according tothe colour palette shown on the right, the red ones being thelargest and blue ones being the smallest. It can be seen thatthe envelopes for different colours (which can be thoughtof iso-size contours) delineate curves lines in the parameterspace. Consider the self-similar case for a moment, in whichthe combination Lt appear together in the relation for size R ∝ ( Lt ) / . If bubbles were to grow with this scaling,then these envelopes would be curves of t ∝ L − / . Butthese curves in Figure 6 are more complicated, signifyingnon-self-similar evolution of bubbles, even in the adiabatic Figure 6.
Scatter plot of bubbles in the parameter space of me-chanical luminosity L and bubble age, without cooling, after 30Myr of the onset of star formation, in which the colour of thedata points refer to bubble sizes, as shown in the colour paletteon the right. case, because of the presence of pressure gradient term thatis neglected in the Weaver et al. (1977) solution. Let us now discuss the effects of cooling. We notice fromFigure 5 that the inclusion of radiative cooling leads to alarge difference in the evolution of the shell radius.The evolution of the bubbles of different luminosities inthe presence of cooling and the pressure difference betweenambient and gas inside the bubble, can be roughly describedby a single parameter η , which can take into account radia-tion loss, where η is defined in terms of R = 0 . ηLt /ρ ) / (the pre-factor (cid:0) π (cid:1) ≈ .
76 applies to the adiabatic caseWeaver et al. (1977)). We show in Figure 7 the ratio of radiusto L / versus time for bubbles of four different luminosities.The curves show that they are nearly superposed on one an-other, which indicates that a single value of η can describetheir evolution, whose value in this case (for the choices ofambient conditions) is inferred to be ∼ .
25. Similar con-clusions have been reached by previous workers, e.g., MacLow & McCray (1988); Krause & Diehl (2014). More re-cently, Sharma et al. (2014) showed that bubbles can retaina fraction of the input energy, of order 0 . . η ∝ n − / ) and is of order ∼
20% for ambient density of 0 . − (their Figure 8), consistent with our estimate in thepresent work. Figure 7 shows that at early epochs, the valueof η can vary with luminosity, since the curves for differentluminosities do not quite overlap. However, they roughly doso by the time of bubble stalling, which is more relevant toour present work. In general, bubbles with lower luminositysuffer more loss of energy due to radiation (or, η is smallerfor lower L ), although this trend is reversed below 4 × c (cid:13)000
20% for ambient density of 0 . − (their Figure 8), consistent with our estimate in thepresent work. Figure 7 shows that at early epochs, the valueof η can vary with luminosity, since the curves for differentluminosities do not quite overlap. However, they roughly doso by the time of bubble stalling, which is more relevant toour present work. In general, bubbles with lower luminositysuffer more loss of energy due to radiation (or, η is smallerfor lower L ), although this trend is reversed below 4 × c (cid:13)000 , 000–000 B. B. Nath, P. Das, M. S. Oey
Figure 7.
Logarithm of the ratio of radius to mechanical lumi-nosity L in bubbles (in units of cm/(erg s − )) is plotted againsttime (in Myr), for four different values of L , shown in differentcolours. The ambient pressure is P amb = 2 . × − dyne cm − ,and equivalent temperature is 4 × K. The near superpositionof the cases of all luminosities show that a single value of η is ableto describe the evolution of bubbles in the presence of radiativecooling. The factor η is found to be ≈ .
25, making the bubblesa fraction ≈ .
76 times smaller than their adiabatic case. erg. This is because of the fact that the outer-shock speedis lower for low luminosity shells, and the resulting post-shock temperature is also reduced, leading to a higher rateof radiative loss, since the cooling function in the relevanttemperature range is a decreasing function of temperature.At the lowest luminosities considered here, the post-shocktemperature is close to 10 K even at the earliest epochs,where the cooling function drops, reversing the trend.We also notice that, as argued previously, bubbles donot stall when pressure inside becomes equal to the ambi-ent pressure. However, the continuation of expansion seen inthe simulation is also misleading, because our use of equiv-alent temperature for non-thermal pressure of the ambientmedium has limitations. The non-thermal pressure in theambient medium due to turbulence and other process willfragment or distort the shell, robbing it of its momentumwhich would have otherwise make it expand further. This isnot easy to simulate without initially introducing density in-homogeneities and turbulence in the ambient medium in thenumerical set-up, which would increase the number of freeparameters in the calculations. The fact that shells wouldnot be able to retain their identity when their nominal speedbecomes comparable to the ambient medium sound speed isself-evident. For example, for supernovae blast waves, thiswas the condition imposed by McKee & Ostriker (1977) intheir three phase model of ISM.Therefore, we impose a condition of stalling the shellwhen the outer shock speed equals the (isothermal) soundspeed of the ambient medium. This is shown by the hori-zontal section of the green curves in Figures 5.
Figure 8.
Scatter plot of bubbles in the parameter space of me-chanical luminosity L and bubble age, after 30 Myr of the onsetof star formation, in which the colour of the data points refer tobubble sizes, as shown in the colour palette on the right. We first show in Figure 8 a scatter plot of the bubbles as afunction of mechanical luminosity (up to 10 erg s − ) andbubble age, for a snapshot at 30 Myr after the onset of starformation. The bubble sizes being shown in different coloursaccording the colour palette shown on the right, the red onesbeing the largest and blue ones being the smallest. The starformation rate is considered to be uniform and equal to 1 M (cid:12) yr − , and the luminosity function of clusters is assumed toobey a power-law with index β = 2. If we take a snapshot at30 Myr for bubbles triggered by mechanical luminosity up to10 erg s − , then these are the bubbles that would show up.They will have different ages, as shown by their distributionalong the vertical axis. the vertical colour stripes indicatethat most bubbles have stalled, being at the same radiusat different times, except for the bubbles in the bottom ofthe figure. This is in contrast to the case without cooling,in which the population of superbubbles is dominated bygrowing bubbles, in Figure 6.The scatter plot shows that at any given time (here,at 30 Myr), the smallest bubbles are produced mostly bylow luminosity clusters (blue circles) and they are predomi-nantly young. This is shown by the fact that blue dots mostlyappear at the left bottom corner of this plot. The largestbubbles are created by clusters at the high luminosity endand can be both young and old. Most of the points in thescatter plot, however, arise due to stalled bubbles. Leavingaside the smallest bubbles (blue), there are vertical columnsof different colours (different sizes) in the figure. This im-plies that the sizes of the bubbles mostly correspond to themechanical luminosity of the bubbles, and almost indepen-dent of the age. This is because of the stalling condition wehave imposed at the time of outer shock speed becomingcomparable to the ambient medium sound speed.The stalled radii can also be written in terms of thestalling time t s as (equating the ambient isothermal sound c (cid:13) , 000–000 uperbubble size distribution Figure 9.
Final bubble size as a function of mechanical luminos-ity, in the case of self-similar evolution (where stalling is assumedto occur at pressure equilibrium stage, OC97), and from simula-tions with cooling (in which stalling is imposed when outer shockspeed equals the ambient isothermal sound speed). speed of (cid:112) P /ρ with the outer shock speed, R/t ), R cool,e = 53 P / ρ − / t s , (8)and the corresponding luminosity is given by, L e = (cid:16) (cid:17) (cid:16) π (cid:17) − P / ρ − / t s η − / . (9)These expressions allow us to estimate the largest bubblesize R e that is reached after t = t s = t e ≈
40 Myr, afterwhich the OB stars in a cluster drop off the main sequenceand the mechanical luminosity ceases to power any furthergrowth of the bubble. For the same fiducial ISM parametersas above, we get R cool,e = 1388 . L e ≈ × ergs − . Clusters with luminosity larger than this will continueto grow and not stall even at t e (40) Myr. We find that thesevalues are similar to the ones considered in OC97 ( R e =1300 pc, L e = 2 . × erg s − ). These estimates also giveus an idea of bubbles that are still evolving and have notreached steady-state at a certain epoch. For example, forthe same ISM parameters, at 10 Myr (after the onset of starformation), bubbles smaller than 347 pc (corresponding to L ∼ . × erg s − ) have reached steady-state and largerones are still evolving.Figure 9 shows the relation between final bubble sizefrom our simulations (red points) as a function of mechani-cal luminosity and shows that the size scales as L / . Thisfollows from the fact that radius of a bubble scales as R ∝ ( ηL ) / t / , where η ∼ .
25 takes into account the en-ergy loss by radiation and pressure gradient, as mentionedearlier. This implies an outer shock speed v ∝ ( ηL ) / t − / .Since the bubbles are assumed to stall when this speed be-comes comparable to the ambient sound speed, the stallingtime scales as t s ∝ L / , and consequently, R f ∝ L / .This stalling condition is similar to that used by OC97,since the ram pressure of the outer shock is related to thepressure of the shocked wind region in a bubble (Weaveret al. 1977), which dominates the pressure inside a bubble.Therefore, it is not a surprise that we get a similar relation Figure 10.
Same as in Figure 6, but for luminosities in the range10 –6 × erg s − . Figure 11.
Size distribution of bubbles in our simulation, forambient pressure p = 2 . × − dynes cm − , at 10, 20, 30Myr. The distributions at different epochs are fitted with powerlaw, and the power-law indices are found to be ∼ . between radius and luminosity as in OC97. However, themagnitudes of these two types of pressure are different. Theinner pressure is roughly ∼ . ∼ .
2. However when radiative cooling is taken intoaccount, then this criterion gives a smaller radius becauseof radiation loss.In the OC97 case, the stalled size is given by (using c (cid:13) , 000–000 B. B. Nath, P. Das, M. S. Oey their equations 31 & 32), R f,OC ≈ × / √ π P / ( µm H n ) / L / ≈
305 pc (cid:16) P . × − dyne cm − (cid:17) − / × (cid:16) n . − (cid:17) / (cid:16) L erg s − (cid:17) / . (10)This relation is shown by the green line in Figure 9. For thecase of cooling in our simulation, equating the shock speedwith the ambient (isothermal) sound speed leads to, R cool ≈ (cid:16) (cid:17) − / (cid:16) π (cid:17) / P − / ( µm H n ) / ( ηL ) / ≈
180 pc (cid:16) P . × − dyne cm − (cid:17) − / × (cid:16) n . − (cid:17) / (cid:16) L erg s − (cid:17) / (cid:16) η . (cid:17) / . (11)The determination of the outer shock speed and thestalling time, however, involves the smoothening of the outershock speed versus time, and the final determination of thestalled radius has an uncertainty of order (cid:46)
20% owing tothis. The stalled radii in the case of cooling are shown withred points in Figure 9.The second reason why we get a similar scaling betweenradius and luminosity is the fact that radiation and ambientpressure affects the size evolution through a single factor η , leaving the dependency of size on luminosity the intact.These two facts contrive to make the results in OC97 andin the present work look similar, even in the presence ofcooling.The corresponding size distribution is, therefore, againexpected to of the type N ( R ) ∝ R − β , and it is shown inFigure 11. The fitted power-law indices at different epochs(10, 20, 30 Myr) are roughly ∼ . .
9. This is close to thevalue of 1 − β (= − − η weakly depends on luminosity. Had η been totally independent of L , then, the size distributionwould have been exactly 1 − β (= − η is some-what smaller for lower L than its fiducial value, the shellsfor low luminosity clusters at stalled phase are somewhatsmaller (see equation 11), leading to a somewhat flatter sizedistribution than the fiducial 1 − β value. Since the radii ofbubbles have decreased in general by a factor of 1 .
7, and thepeak radius has decreased by a similar factor compared tothe adiabatic case in Figure 1, to ∼
10 pc. While determin-ing the size distribution we did not consider holes smallerthan 10 pc, and therefore we do not see any rising part inthe distribution with radiative cooling.One notes in Figure 11 that the slope slowly changeswith time, and becomes flatter. This is because of the factthat bubbles with large luminosity take a while to stall, andin any given snapshot, there would be newly formed (largeluminosity) bubbles which would be in a growing phase. Thestalling timescale is given by t s ≈ . (cid:16) P . × − dyne cm − (cid:17) − / × (cid:16) n . − (cid:17) / (cid:16) L × erg s − (cid:17) / (cid:16) η . (cid:17) / . (12)This implies that over time, the number of large luminosity bubbles, or consequently, large size bubbles would increase,whereas the small bubbles would have stalled quickly andtheir numbers would more or less freeze. This is shown inthe scatter plot for the bubbles, in Figure 10 for the rangeof luminosities 10 –6 × erg s − . The large luminositybubbles are seen to be growing in this plot, and the timescale for attaining stalled phase for L ∼ × erg s − is ≈ ≈ − .
7) in the presence of cooling, orif this depends on the ambient pressure. In order to answerthe question, we have run our simulations for two differentpressures.Figure 12 shows the scatter plot of bubbles (similar toFigure 6 and Figure 8 for our fiducial ISM pressure) for twodifferent values of ISM pressures, at 30 Myr. Essentially, wehave decreased the non-thermal contribution in the ambientpressure, by reducing the equivalent temperature from 4 × K to 2 × K (left panel) and then to 10 K (rightpanel). Comparing with Figures 6 (adiabatic case) and 8, wefind that the bubbles in these two cases are not dominatedby stalled bubbles, as was the case for Figure 8. For the leftpanel, we find that bubbles with L ∼ × erg s − stallafter a time scale of ∼ − β , as we confirm below.We show in Figure 13 the size distribution of the bub-bles for these values of pressure at two different times afterthe onset of star formation (10 and 30 Myr). In general, wefind that lowering pressure steepens the size distribution,by allowing low luminosity bubbles to grow to larger sizes.This trend is shown in Figure 14, for two different epochs,10 Myr and 30 Myr. At lower pressures, there is an evolutionof the slope between these time scales, since bubbles keepgrowing in low ambient pressure, compared to high pressureenvironments.The distributions at lower (non-thermal) pressures alsoshow a positive part at small sizes, in addition to the fallingnumbers at large sizes as we have seen earlier. The rising partcomes from growing young bubbles whose age is smaller thanthe stalling time of the bubble with the lowest luminosity.This generates a peak in the distribution. The peak dependsupon the radius of the stalled bubble with the lowest lumi-nosity at the time being observed. Figure 13 shows that athigher pressures the peak radius is independent of time butstarts varying with time as soon as we start decreasing thepressure. This occurs due to the fact that at higher pressures(such as, at 2 . × − dynes cm − ), the bubbles with lowluminosities stall much before 10 Myr hence the peak occurs c (cid:13) , 000–000 uperbubble size distribution Figure 12.
Scatter plot for bubbles for ambient pressure p = 1 . × − dynes cm − (left), and 6 . × − dynes cm − (right), at30 Myr. Figure 13.
Size distribution of bubbles for ambient pressure p = 2 . × − dynes cm − (blue), 1 . × − dynes cm − (red),6 . × − dynes cm − (green), at 10 and 30 Myr. at the same radius even when we observe the distributionat 30 Myr. However at lower pressures like p = 6 . × − dynes cm − the bubbles with low luminosities stall muchafter 30 Myr, and hence the peak shifts towards larger radiiwith time. In other words, the deviation of the slope from − L ∝ R of stalled bubbles,which would have produced a − It is known that observations support the 1 − β power lawfor the superbubble size distribution. For example, Bage-takos et al. (2011) found a power-law index of ≈ − . The HI Nearby Galaxy Survey of 20 nearby galaxies. This isconsistent with β ≈ . ≈
4) than late type galaxies. This indicates,a pre-ponderance of small size holes in early type galaxies. Theyassociated this phenomenon with the level of star formationactivity in different galaxy types. One obvious connection isthat in early type galaxies, the observed mass function ofstar clusters is steeper than late type galaxies, as evidencedby the HII region luminosity function (e.g., Kennicutt et al.(1989)). While the luminosity function in Sb-Sc galaxies hasan index of ∼ −
2, the index varies from ∼ − . ∼ . c (cid:13) , 000–000 B. B. Nath, P. Das, M. S. Oey
Figure 14.
The slope of size distribution as a function of ISMpressure, for two different epochs, 10 Myr (green) and 30 Myr(red). (Joung et al. 2009). According to their results, the turbu-lent pressure scales with surface density of star formationas P turb ∝ Σ / ∗ . Since our results indicate that decreasingnon-thermal pressure steepens the size distribution, this re-mains another possibility. Future observations will be ableto point towards the right explanation.We note that our calculations do not take into accountthe merging of superbubbles and its effect on the size distri-bution. This was discussed by OC97 in terms of a porosityparameter Q (Cox & Smith 1974), which is the ratio of su-perbubble volume to total ISM volume. For the Milky Way,the SFR is near the critical point where Q ∼
1. It can beshown that the total volume occupied by superbubbles canexceed the volume of Milky Way ISM, considering a cylin-der of 10 kpc radius and 500 pc thickness, within (cid:54) ∼ (cid:12) yr − .Following Clarke & Oey (2002), if one takes the averagemass of clusters as ∼ (cid:12) , then the number of super-bubbles produced per Myr is ∼ (cid:12) yr − . Then with a size distribution with aslope of −
3, the total volume exceeds the Milky Way ISMvolume in le ∼
10% (Bagetakos et al. 2011). This implies that mergingof superbubbles is important. It is also evident from the ob-servations of Simpson et al. (2012) that roughly ∼
30% ofHI shells show signs of merging. OC97 also explore how theMilky Way Q varies with the assumed value of β . Mergingamong superbubbles likely flattens the size distribution tosome extent, by increasing the number of large bubbles atthe expense of smaller bubbles. However, it is difficult toestimate the magnitude of this effect without a simulationthat includes non-uniformity in the spatial distribution ofstar clusters. We have studied the form and evolution of the size distri-bution of HI holes in the ISM of galaxies owing to super-bubbles triggered by OB associations, taking into accountradiative cooling, with the help of numerical hydrodynami-cal simulations. Previous works had assumed bubble growthstalls when the inner pressure of adiabatic bubbles equalsthe ambient pressure, which is not valid since the bubblesmaintain momentum-driven growth. Assuming that bubblesstall when the expansion speed becomes comparable to theambient sound speed, we show that the inclusion of radiativecooling and ambient pressure results in a power-law index ofthe size distribution with slope ∼ − . p = 2 . × − dyne cm − and density n = 0 . − .This is consistent with observations by THINGS . We havefurther shown that decreasing the ISM pressure can makethe population of growing bubbles dominating over stalledones, consequently making the size distribution steeper. Ourresults imply that the size distribution can help interpretthe evolution of bubbles, with the slope being ∼ − . We would like to thank K. S. Dwarakanath, Nirupam Roy,Aditi Vijayan, Siddhartha Gupta and Ranita Jana for use-ful discussions. The authors thank the anonymous refereefor helpful comments. PD would like to thank Raman Re-search Institute for local hospitality during her visits for theproject.
REFERENCES
Bagetakos, I., Brinks, E., Walter, F., De Block, W. J. G.,Usero, A., Leroy, A. K., Rich, J. W., Kennicutt, Jr, R. C.2011 AJ, 141, 23Banfi, M., Rampazzo, R., Chincarini, G., Henry, R. B. C.1993, A&A 280, 373Caldwell, N.,Kennicutt, R., Phillips, A. C., Schommer, R.A. 1991 ApJ, 370, 526Chevalier, R. A. & Clegg, A. W. 1985, Nature, 317, 44Clarke, C. & Oey, M. S. 2002, MNRAS, 337, 1299Cox, D. P. & Smith, B. W. 1974, ApJLett, 189, L105Elmegreen, B. G. & Scalo, J. 2004, ARAA, 42, 211Fierlinger, K. M., Burkert, A., Ntormousi, E., Fierlinger,S., Schartmann, M., Ballone, A., Krause, M. G. H., Diehl,R. 2016 MNRAS, 456, 710Gentry, E. S., Krumholz, M. R., Madau, P., Allessandro,L. 2019 MNRAS, 483,3647Gupta, S., Nath, B. B., Sharma, P. & Eichler, D. 2018MNRAS, 473,1537Jenkins, E. B. & Tripp, T. M. 2011, ApJ, 734, 65 c (cid:13) , 000–000 uperbubble size distribution Joung, M. R., Mac Low, M.-M., Bryan, G. L. 2009, ApJ,704, 137Kennicutt, R. C., Edgar, B. K, Hodge, P. W. 1989, ApJ,337, 761Krause, M. G. H., Diehl, R. 2014 ApJ, 794, L21Krause, M. G. H. et al.2015 ApJ, 794, L21Kroupa, P. 2002, Science, 295, 82Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999,ApJS, 123, 3Mac Low, M.-M., McCray, R. 1988 ApJ, 324, 776McKee, C. F., Ostriker, J. P. 1977 ApJ, 218, 148McKee, C. F., Williams, J. P. 1997 ApJ, 476, 144Mignone, A., Bodo, G., Massaglia, S., Matsakos, T. etal.2007, ApJS, 170, 228Oey, M. S., Clarke, C. J. 1997 MNRAS, 289, 570Oey, M. S., Clarke, C. J. 1998 AJ, 115, 1543Seo, J., Kang, H., Ryu D. 2018 JKAS, 51, 37Sharma, P., Roy, A., Nath, B. B., Shchekinov, Y. 2014,MNRAS, 443, 3463Simpson, R. J., Povich, M. S., Kendrew, S. et al.2012, MN-RAS, 424, 2442Weaver, R., McCray, R., Castor, J., Shapiro, P., Moore, R.,1977, ApJ, 218, 377Wolfire, M. G., McKee, C. F., Hollenbach, D., Tielens, A.G. G. M. 2003, ApJ, 587, 278Yadav, N., Mukherjee, D., Sharma, P., Nath, B. B. MN-RAS, 465, 1720
Here, we present the convergence tests of the numerical runsfor superbubble evolution for our fiducial case, for a differ-ent spatial resolution, with ∆ r = 0 . .
16 pcused earlier. Figure 15 shows the size distribution and fit-ted slopes at 10 , ,
30 Myr, for the fiducial case of P =2 . × − dyne cm − . The slopes ( − . , − . , − . − . , − . , − .
64) obtained fora coarser resolution (Figure 11). This confirms the conver-gence of our results with respect to numerical resolution.This implies that η has also reached convergence limit,and we show in Figure 16 the logarithm of ratio of radiusto L / vs. time for different luminosities, superimposed onthe expected evolution for η = 0 .
25. The curves show (as inFigure 7) that η depends on luminosity rather weakly.We should note that our runs do not simulate the effectsof turbulent mixing in ISM, which might render numericalconvergence ineffective (e.g., Gentry et al. (2019); Fierlingeret al. (2016)). However, those effects are beyond the scopeof the present work. Figure 15.
Similar to Figure 11 except that it is for ∆ r = 0 . Figure 16.
Similar to Figure 7 except that this is for ∆ r = 0 . (cid:13)000