© ESO, 2013

1. Introduction

Most massive stars are born in OB-associations in coeval starbursts on timescales of less than 1−2 Myr (Massey 1999). These associations can contain a few up to many thousand OB-stars, so-called super star clusters, but typically have 20−40 members (McCray & Kafatos 1987). Energy and mass are injected through strong stellar winds and subsequent supernova (SN) explosions of stars with masses above 8 M_⊙. The emerging shock fronts sweep-up the ambient interstellar medium (ISM) and, as the energy input in form of SN-explosions continues, superbubbles (SBs) are produced, which may reach dimensions of kiloparsec-size (e.g. Tenorio-Tagle et al. 2003). The swept-up ISM collapses early in the evolution of the SB into a cool, thin, and dense shell (Castor et al. 1975) present in HI and Hα observations. The bubble interior contains hot (>10⁶  K), rarefied material, usually associated with extended diffuse X-ray emission (Silich et al. 2005). Due to the stratification of the ISM in disk galaxies, the superbubbles can accelerate along the density and pressure gradient and blow out into the halo, appearing as elongated structures. Examples of huge bubbles and supergiant shells are the Cygnus SB with a diameter of 450 pc in the Milky Way (Cash et al. 1980) and the Aquila supershell extending at least 550 pc into the Galactic halo (Maciejewski et al. 1996). Our solar system itself is embedded in an HI cavity with a size of a few hundred parsecs called the Local Bubble (Lallement et al.2003), and is most likely generated by stellar explosions in a nearby moving group (Berghöfer & Breitschwerdt2002). Also in external systems like in the LMC (Chu & Mac Low 1990), in NGC 253 (Sakamoto et al. 2006) and M101 (Kamphuis et al. 1991) such bubbles, holes and shells are observed.

The acceleration of the shell promotes Rayleigh-Taylor instabilities and after it is fully fragmented, only the walls of the SB are observed. Through such a chimney the hot pressurized SN-ejecta can escape into the halo. The walls may be subject to the gravitational instability and interstellar clouds can form again, which triggers star formation (e.g. McCray & Kafatos 1987). This is observed, for example, on the border of the Orion-Eridanus SB (Lee & Chen 2009).

The knowledge of SB evolution is crucial for understanding the so-called disk-halo connection and it also gives us information about the chemical evolution of the galaxies, the enrichment of the intergalactic as well as the intracluster medium. The thick extraplanar layer of ionized hydrogen seen in many galaxies has probably been blown out of the disk into the halo by photoionization of OB-stars and correlated SNe (Tüllmann et al. 2006). With the high star formation rate of starburst galaxies, the energy released by massive bursts of star formation can even push the gas out of the galactic potential, forming a galactic wind. Outflow rates of 0.1 − 10   M_⊙/yr are common for starburst driven outflows (Bland-Hawthorn et al. 2007). If the hot and metal enriched material is brought to the surrounding intergalactic medium, it will mix after some time, increasing its metallicity. Galactic winds are observed in nearby galaxies (e.g. Dahlem et al. 1998), as well as up to redshifts of z ~ 5 (e.g. Swinbank 2007; Dawson et al. 2002), thus being an ubiquitous phenomenon in star forming galaxies. If the energy input is not high enough, the gas will fall back onto the disk due to gravity, after loss of pressure support, forming a galactic fountain (Shapiro & Field 1976; de Avillez 2000). In this case, the heavy elements released by SN-explosions are returned back to the ISM in the disk, presumably spread over a wider area, and future generations of stars will incorporate them.

Spitoni et al (2008), also using the Kompaneets approximation, have investigated the expansion of a superbubble, its subsequent fragmentation and also the ballistic motion of the fragments including a drag term, which describes the interaction between a cloud and the halo gas. In addition, these authors have analyzed the chemical enrichment of superbubble shells, and their subsequent fragmentation by Rayleigh-Taylor instabilities, in order to compare the [O/Fe]-ratios of these blobs to high velocity clouds (HVCs). They find that HVCs are not part of the galactic fountain, and even for intermediate velocity clouds (IVCs), which are in the observed range of velocities and heights from the galactic plane, the observed [O/Fe]-ratios can only be reproduced by unrealistically low initial disk abundances. The chemical enrichment of the intergalactic medium will be the subject of a forthcoming paper.

This paper is structured as follows: Sect. 2 shows how the evolution of superbubbles can be described. In Sect. 3 we present the results of this work, which is mainly the expansion of a bubble in space and time and also the onset of Rayleigh-Taylor instabilities in the shell. In Sect. 4 the models are used to analyze two Milky-Way superbubbles. A discussion follows in Sect. 5 and summary and conclusions are presented in Sect. 6.

2. Superbubble evolution

2.1. ISM stratification

For a homogeneous ISM, the propagation of a shock front originating from an instantaneous release of energy was described by Sedov (1946) and Taylor (1950), while the expansion of an interstellar bubble with continuous wind energy injection was studied by Castor et al. (1975) and Weaver et al. (1977). Yet, the description of SB evolution in such a uniform ambient medium is only valid for early stages of evolution. The vertical gas density distribution has a major effect on the larger superbubbles which have sizes exceeding the thickness of a galactic disk. On these scales, the ISM structure is far from homogeneous. In a Milky-Way-type galaxy, the warm neutral HI is plane-stratified and can be described by a symmetric exponential atmosphere with respect to the midplane of a galaxy (Lockman 1984). The scale height of this layer is ~500 pc. The Reynolds layer of warm ionized gas has a scale height of H = 1.5  kpc (Reynolds 1989) and the highly ionized gas of the hot (10⁶−10⁷  K) halo surrounds the galaxy extending to ~20  kpc (Bregman & Lloyd-Davies 2007) with a scale height of ~4.4  kpc (Savage et al. 1997). In the disk, the cold neutral and molecular gas are found to have H ~ 100  pc.

The expansion of superbubbles in a stratified medium was studied since many decades, both analytically (e.g. Maciejewski & Cox 1999) and numerically (Chevalier & Gardner 1974; Tomisaka & Ikeuchi 1986; Mac Low & McCray 1988, hereafter MM88). For an analytic description of superbubbles, Kompaneets’ approximation (1960) is a very good choice. Although it involves several simplifications (e.g. no magnetic field), it represents the physical processes involved very well (Pidopryhora et al. 2007).

2.2. Blow-out phenomenon and fragmentation of the shell

A focus of this paper is to analyze the blow-out phenomenon: in an exponential as well in an homogeneous ISM, a bubble decelerates first. But due to the negative density gradient of the exponentially stratified medium and the resulting encounter with very rarefied gas, the bubble starts to accelerate into the halo and even beyond, if the energy input is high enough, at a certain time in its evolution. MM88 call this process blow-out, which is also used by Schiano (1985) and Ferrara & Tolstoy (2000), but there, blow-out involves complete escape of the gas from the galaxy. From a more phenomenological point of view, Heiles (1990) distinguishes between breakthrough bubbles, which break out of the dense disk and are observed as shells and holes in the ISM, whereas blow-out bubbles break through all gas layers and inject mass and metals into the halo. In our definition (see Sect. 2.3), a SB will blow out of a specific gas layer at the time, when the outer shock accelerates, and if the shock stays strong all the time. In particular, we want to determine the energy input required for blow-out into the halo, and its dependence on ISM parameters (see Sect. 3.2.).

Shortly after the acceleration has started, Rayleigh-Taylor instabilities (RTIs) appear at the interface between bubble shell and hot shocked bubble interior. As the amplitudes of the perturbations grow, finger-like structures develop at the interface and vorticity of the flow increases due to shear stresses. Finally, in the fully nonlinear phase of the instability, turbulent mixing of the two layers starts, and the shell will eventually break-up and fragment. An azimuthal magnetic field in the shell will limit the growth rate of the instability due to magnetic tension forces, but not for all wave modes of the instability. In essence, only modes above a critical wave number will become unstable, giving a lower limit to the size of blobs (Breitschwerdt et al. 2000). This happens first at the top of the expanding bubble, where the acceleration is highest. The clumps generated this way and the hot gas inside the bubble – including the highly enriched material – are expelled into the halo of the galaxy or even into intergalactic or intracluster space, contributing to the chemical enrichment of galactic halo or intracluster medium.

After deriving the acceleration of the outer shock in Sect. 3.2, we analyze the timescales for the development of RTIs in the shell in Sect. 3.3.

2.3. Modeling superbubbles

Originally developed to describe the propagation of a blastwave due to a strong nuclear explosion in the Earth’s atmosphere, the approximation found by Kompaneets (1960) is also applicable to investigate the evolution of a superbubble in a disk galaxy analytically. We modify Kompaneets’ approximation (KA) in order to describe not only a bubble driven by the energy deposited in a single explosion or as a continuous wind (e.g. Schiano 1985), but to produce analytical models for the expansion using a time-dependent energy input rate due to sequential SN-explosions of massive stars according to a galactic initial mass function (IMF). The axially symmetric problem is described in cylindrical coordinates (r,   z). In order to use the KA for the investigation of SB evolution the following assumptions have to be made: (i) the pressure of the shocked gas is spatially uniform; (ii) almost all of the swept-up gas behind the shock front is located in a thin shell, and (iii) the outer shock has to be strong all over the evolution of the bubble. Using the third assumption we get our blow-out condition: if the outer shock has a Mach number M ≥ 3, i.e. the upstream velocity of the gas in the shock frame is at least 3 times the sound speed of the warm neutral ISM at the transition of deceleration to acceleration, then the bubble will blow-out into the halo and reach regions of higher galactic latitudes. Mac Low et al. (1989, hereafter MMN89) find via comparisons to their numerical simulations that the KA can be even used after the instabilities in the shell set in, because the pressure inside the bubble is not released very quickly. If the condition is not fulfilled (i.e. no strong shock), the system is obviously not energetic enough, hence cannot be described by the KA and will finally slow down and merge with the ISM like in the case of disk supernova remnants. Blow-out usually occurs when the shell reaches between one (Veilleux et al. 2005) and three scale heights (Ferrara & Tolstoy 2000). In the next section, this is confirmed and exact values for the height of the bubble are given, using different descriptions of the ambient ISM and different ways of energy input (see Table 2). Compared to other groups using a dimensionless dynamical parameter introduced by MM88 to decide if blow-out will happen, our criterion’s advantage is that the result is given in explicit numbers, and hence easier to use when compared to observations. Additionally, we investigate if the acceleration of the bubble at the time where fragmentation occurs will exceed the gravitational acceleration near the galactic plane. Using this simple comparison we can ensure that fragments of the bubble and the hot bubble interior will be expelled into the halo instead of falling back onto the disk.

A pure exponential atmosphere was used in the KA, which was already modified for a radially stratified medium (Korycansky 1992) and an inverse-square decreasing density (Kontorovich & Pimenov 1998). We adopted the original calculations to model the expansion of superbubbles in a more realistic fashion and investigate two cases of density distribution in our paper: the first one corresponds to a bubble evolving in an exponentially stratified medium symmetric to the midplane (1)where ρ₀ is the density in the midplane and H is the scale height of the ISM. The cluster is located at z = 0. Since it is an idealization that OB-associations are only found in the galactic midplane, but are rather offset in z-direction, we examine in our second case an off-plane explosion, where the density law is given by (2)with z₀ being the center of the explosion. The density at z = z₀ is either derived from the relation ρ₁ = ρ₀·exp [ − z₀/H] or can be a known value. Any real case encountered will be described by one of the two or lie in between, i.e. the offset is small enough for the midplane to be pierced by further explosions.

To get a spatial solution the Rankine-Hugoniot conditions for a strong shock for every point of the azimuthally symmetric shock surface have to be solved (e.g. Bisnovatyi-Kogan & Silich 1995). This gives the normal component of the expansion speed (3)where γ is the ratio of specific heats (γ = 5/3 for a monatomic ideal gas). We need to know the pressure P(t) in the bubble, (4)with Ω(t) being the volume confined by the shock and E_th the thermal energy in the hot shocked gas region. The volume is given by the integral (5)Introducing a transformed time variable (in units of a length) (6)makes it easier to solve the PDE which is obtained after rearranging and equating Eq. (3) and the time derivative of Eq. (6) and by assuming that the shock front is a time-dependent surface, defined as f(r,z,t) = 0 (7)By separation of variables one gets the solution r = r(y,   z), which describes the evolution of the half width-extension of the bubble parallel to the galactic plane. For the symmetric density law this is (8)and for an off-plane explosion we find (9)Due to the explicit relation between time t and y, the parameter y is used to represent the evolution of the bubble. At this point we introduce a dimensionless parameter . Figures 1 and 2 show the position of the outer shock at different values of for the symmetric and off-plane model, respectively.

Fig. 1 Position of the shock front in the symmetric model at certain values of the dimensionless time variable , 1.0, 1.4, 1.7, 1.9, 1.98, and 2.0 with the energy source in the galactic midplane.

Fig. 2 Same as Fig. 1, but for the off-plane model. The energy source is located at z0 = 0.7H above the plane (indicated by the dashed line).

A few points of the bubble’s surface are special, because explicit equations are available for them from the solution (Eqs. (8) and (9)). When , the top z_u and bottom z_d of the bubble can be evaluated. In the case of the symmetrically decreasing density, the top of the bubble as a function of is given by (10)and the bottom of the bubble is simply . In the other case, the coordinate z_u shows an offset of z₀ and we derive (11)and for the bottom of the bubble (12)Furthermore, the half-width maximum extension of the bubble, r_max, parallel to the galactic disk is found, where , (13)which is valid for both density laws.

Obviously, the top of the bubble reaches infinity by the time where approaches the value 2. This happens, because the newly swept-up mass asymptotically goes to zero, due to the exponential decrease in density. The solution then breaks down, because numerically the shock reaches infinity in a finite time. Also the remnant volume goes to infinity (see Figs. 1 and 2). In reality, a nonzero density at infinity leads to restricted values of the shock front velocity (Kontorovich & Pimenov 1998). Actually, blow-out occurs earlier in the evolution of a bubble with values of well below 2.

3. Results

We use different energy input schemes to calculate blow-out timescales, the energy input that is needed for blow-out, and the instability timescales. In the first and simplest case, a single, huge explosion forms the superbubble like in the original version of Kompaneets (1960). A time dependent energy input rate is implemented in the next model, where the number and sequence of SN explosions in a star cluster is given by an initial mass function. Finally, the wind model uses a constant energy input rate to drive the expansion of the bubble, which is just a special case of the IMF-model. In all calculations, a cluster is defined to have at least two member stars.

We apply these models to two cases of ISM density distribution. First, we study the evolution of a bubble in an ISM representing the Lockman layer of a Milky-Way-type galaxy (number density n₀ = 0.5  cm^-3, scale height H = 500  pc) with the center of the explosion placed in the midplane and exponentially, symmetrically decreasing density above and below the disk (Eq. (1)). In the second case, the star cluster is displaced by z₀ from the galactic midplane and the SB is expanding into a rather dense (n₀ = 10 cm^-3), low-scale height (H = 100  pc) medium described by a simple exponential law (cf. Eq. (2)). In this off-plane configuration the bubble shall not expand below the midplane before blow-out (i.e. roughly the time when RTIs start to appear in the shell), thus it is not affected by the increasing density below the plane until that time, which results in a one-sided blow-out. In other words, the value of z₀ has to correspond to the absolute value of the coordinate (Eq. (12)) in the case of an unshifted bubble (z₀ = 0) at the time when the acceleration sets in. Values of at range between 0.45H and 1.0H, thus we take an average value and put the explosion at z₀ = 0.7H¹ in order to be able to compare the models. This corresponds very well to MM88’s criterion for ’one-sided superbubbles’: bubbles for which the association is found above 0.6H blow out on one side of the disk only, and the bottom of the bubble should be decelerating more strongly than a spherical one would do.

3.1. Thermal energy

Basu et al. (1999) show that the thermal energy in the hot interior of a SB expanding into an exponentially stratified medium is very close to the value in the case of a homogeneous ISM until about four times the dimensionless timescale t/t_D, which is rather late in the evolution of the superbubble. We find that at the time of fragmentation ~2 − 3   t/t_D are reached (see Table 3). Thus, we can estimate the thermal energy in the region of hot shocked gas following the calculations of Weaver et al. (1977) for a wind-blown bubble in a uniform ISM. These take into account the equations of energy and momentum conservation, as well as the radius of a spherical bubble. The inner shock, where the energy conversion takes place is always roughly spherical, since it is close to the energy source. However, with increasing time, the dynamics of a bubble in an exponentially stratified medium will differ from that in a homogeneous medium, where no blow-out will happen.

3.1.1. SN-model

All SNe explode at the same time. We find that the fraction of the total SN-energy converted into thermal energy at the inner shock is E_th,   SN(t) = 2/3·E_SN·N_SN (see Appendix A). N_SN is the total number of SNe, with each explosion releasing a constant energy of E_SN = 10⁵¹  erg.

3.1.2. IMF- and wind-model

A time-dependent energy input rate L_SB(t) = L_IMF·t^δ is used, where the energy input rate coefficient L_IMF and the exponent δ are characterized by the slope of the IMF and the main sequence lifetime of massive stars. The calculations of the time-dependent energy input rate follow Berghöfer & Breitschwerdt (2002) and Fuchs et al. (2006) and shall be presented here briefly. The IMF describes the differential number of stars in a mass interval (m,m + dm) by a power law (14)where Γ is the slope of the IMF and m is always given in solar mass units. Integration from a lower mass limit m_l to an upper mass limit m_u gives the cumulative number of stars or – having m_l ≥ 8M_⊙ – the number of OB stars to explode as SNe in this stellar mass range (15)In our general treatment of the IMF either N₀ or m_u will be given and in order to get N_OB from the equation above, a correlation between N₀ and m_u is needed. We use integer mass bins and simply fix the number of stars in the last mass bin N(m_u − 1,m_u) = N_OB = 1, i.e. there is exactly one star in the mass bin of the most massive star. Hence the normalization constant is . When dealing with a real association, the number of stars N_OB in a certain mass range (m_l, m_u) can be deduced from observations and thus the normalization constant will be estimated. Moreover, using this information is statistically more relevant than using THE most massive star of the cluster because the distribution of stars in a real clusters may not follow integer mass bins.

Since the distribution of the stars by their mass is given by the IMF, and L_SB follows the time-sequence of massive stars exploding as SNe, we get the energy input rate L_SB(t) = d/dt[(N(m)·E_SN)]. We just need to express the function N(m) as a time-sequence, thus we treat the number of stars between (m,m + dm) as a function of mass and use the main sequence lifetime (t,t − dt) of massive stars (16)The main sequence lifetime and the mass of a star are connected through t(m) = κ·m^− α or m(t) = (t/κ)^− 1/α. The values of Fuchs et al. (2006) are used throughout this paper, where κ = 1.6 × 10⁸  yr and α = 0.932. The time-derivative of m(t) together with dN/dm from Eq. (14) are inserted into Eq. (16) (17)Since L_SB(t) ∝ t^δ, the exponent must be δ =  −⁽Γ/α + 1⁾. After summarizing the constants one obtains the energy input rate coefficient (18)The full equation of the thermal energy as a function of time for the IMF-model is (see Appendix A for details) (19)Using an IMF yields a more realistic framework for galactic SB expansion compared to a simplified point explosion. Therefore, we investigate the effect of changing the slope of the IMF in our calculations. We compare three different IMF-slopes for massive stars: Γ₁ =  −1.15 (Baldry & Glazebrook 2003), Γ₂ =  −1.35 (Salpeter 1955), and Γ₃ =  −1.7 (Brown et al. 1994; Scalo 1986) resulting in δ₁ = 0.23, δ₂ = 0.45, and δ₃ = 0.82.

With an exponent δ = 0, Eq. (19) respresents the thermal energy in case of a constant energy input rate, where the number of all SNe is averaged over the whole lifetime of the association, thus this model simply has a slope of Γ₀ =  −0.932 (wind-model). Now the larger time intervals between SN-explosions of stars with higher main sequence lifetime compensate the growing number of SNe per mass interval going to lower mass stars as it would be the case in the IMF-model. With a constant energy input rate L_w, our result (Eq. (19)) checks with Weaver et al’s relation of E_th,w(t) = 5/11·L_w·t.

3.2. Evolution of a superbubble until blow-out

The transformed time variable (Eq. (6)) contains not only the energy input and the ambient density, but also the volume of a bubble at a given time. Since this equation gives us an explicit correlation between time t and the time variable it is possible to describe the volume of the bubble as a function of instead of being dependent on time. To proceed with our analytical description, we need to find simple expressions of the volume for the two cases of density distribution.

Symmetric model

The bubble contour on the + z-side of the midplane can be approximated by being part of an ellipse with a semimajor axis of (20)and corresponding to the bottom of an unshifted bubble in a pure exponential atmosphere. The semiminor axis is equal to (see Eq. (13)). The center of the ellipse is located on the z-axis (21)such that the ellipse is given by (22)Rotating this curve around the z-axis and multiplying by two results in the total volume of the superbubble (23)The approximation works very well and we find a deviation of Eq. (23) from the numerical integration of the volume according to Eq. (5) of only ~2.2% at very late stages of evolution ().

Fig. 3 Blow-out timescales at the coordinate zu for the SN-model. Left: double exponential layer with H = 500  pc and n0 = 0.5  cm-3 (symmetric model); right: exponentially stratified medium with H = 100  pc and n0 = 10  cm-3 (off-plane model). The dashed lines indicate the number of SNe required for blow-out and the corresponding blow-out time.

Off-plane model

The volume will be approximated by treating the 3D shape of the bubble as a prolate ellipsoid (see also Maciejewski & Cox 1999). It is not dependent on the offset z₀, thus the bubble volume as a function of reads (24)At , the analytical result differs from the numerical value according to Eq. (5) by only ~1.8   %.

Obviously, the parameter is important for estimating the evolutionary status of a bubble, since it relates the size of the bubble to time. We can now calculate the age t of a bubble at any value of . Replacing and using the nondimensional volume instead of Ω   (t) helps to obtain a useful relation of the time derivative of Eq. (6) (25)When performing the integration, the thermal energy for each model, i.e. for different kinds of energy input, has to be replaced by the corresponding formula. Substituting and integrating yields the general expression (26)which will be solved separately for all three models in the remainder of this subsection. In order to simplify the integration on the right hand side, and especially to avoid a double integral in the case of in the calculations to follow, we use a series expansion of the bubble’s volume: In both cases of density distribution we expand until 43rd order to make sure the simplified integral has a deviation of ≪1   % from the numerical value at a time of .

For further investigation of the blow-out phenomenon we need to know the velocity and acceleration of the shock front. It is possible to derive these properties analytically at certain points of the bubble’s surface, where explicit equations exist. The velocity can be calculated for top, bottom and maximum radial extension of the bubble, and also at for the symmetric model as well as at for the off-plane model. We will concentrate on calculating the velocity at the top of the SB (which is equal to the absolute value at the bottom for the symmetric model) in this paper, since this is crucial in determining if blow-out happens or not. For all models, the velocity at the top of the bubble is given by (29)with the derivative of z_u with respect to (30)Whereas the velocity given by Eq. (3) depends on time t and the coordinate z, the velocity in Eq. (30) is only dependent on . This makes it easier to find general results in terms of the y-parameter at the time of blow-out for each model specification. The second derivative of z_u gives us the acceleration at this coordinate (31)The calculation of (Eq. (25)) and thus need to be done separately for each model, which is shown in Appendix B.

Once velocity and acceleration are obtained, we get the value of , where the velocity of the top of the bubble has its minimum, i.e. the transition from deceleration to acceleration along the density gradient. The value of the dimensionless variable is the same for all SBs of each model, but corresponds to a different time in the evolution of a bubble depending on the ISM parameters and the energy input. The time interval until is called blow-out timescale (see Figs. 3−5 for SN-, IMF-, and wind-model, respectively).

Fig. 4 Same as Fig. 3, but for the IMF-model (left: symmetric model; right: off-plane model). Thin lines indicate the number of OB-stars needed for blow-out. Γ is the IMF exponent (for details, see text).

Fig. 5 Same as Fig. 3, but as a function of a constant energy input rate (left: symmetric model; right: off-plane model). The number of SNe is obtained by converting the wind luminosity to an energy input averaged over a time interval of 20 Myr, from the explosion of the most massive star until the last SN.

By fixing the velocity of the outer shock at the time , the energy input required for blow-out of the disk can be derived. We assume that for M ≥ 3 at z_u the shock is sufficiently strong and the blow-out condition is fulfilled. Thus, the velocity of the shock has to be at least 3·c_s with respect to an ambient ISM at rest. Using a temperature of the surrounding medium of T ~ 6000  K typical for the Lockman layer of the Galaxy (Crawford et al. 2002), the minimum velocity corresponds to cm/s, where k_B is Boltzmann’s constant; the mean atomic mass is in a gas with mass density where the mean molecular weight of the neutral ISM of μ = 1.3 and the hydrogen mass of m_H = 1.7 × 10^-24  g are used. This criterion is valid for the low-scale height, high-density ISM as well, because one-sided SBs blow out of the dense disk and start to accelerate into the halo at ~2H = 200  pc (Table 2). At these distances from the midplane the presence of the warm gas layer already influences the evolution of the bubble. The expressions for the minimum energy input are derived in the following subsections.

3.2.1. SN-model

Since the evolution of the bubble depends on the scale height H, on the density ρ_0,1 of the ambient medium near the energy source, and on the energy E_th,SN, a timescale (in units of seconds) can be constructed from these quantities (32)Solving Eq. (26) yields the time t elapsed since the explosion for any chosen value of . For a constant value of energy released by a number of supernovae, the integration on the left hand side is simple. Expressing the time t as a function of , rearranging the equation and making use of the timescale (Eq. (32)) yields the following equation (33)For the calculation of the velocity, we first express Eq. (25) in terms of the timescale t_SN(34)Now, using the expressions from above, the velocity at the top of the bubble can be written as (35)The calculation of the acceleration at the top of the bubble is found in Appendix B. The critical value , where the acceleration of the outer shock starts, the corresponding dimensionless timescale and the coordinate of the top of the bubble at this time are summarized in Table 2 for both cases of density distribution. Fig. 3 shows the age of the bubble at the transition from deceleration to acceleration as a function of the number of SN-explosions in the range of 2 − 500 SNe. The dashed lines indicate the number of SNe that are required for blow-out for each kind of density distribution and the corresponding timescale for blow-out. In the case of the symmetric density law with H = 500  pc and n₀ = 0.5  cm^-3 about 56 SNe have to explode at once and the acceleration starts ~40 Myr after the initial explosion. Only three SNe are sufficient for an off-plane explosion in a pure exponential atmosphere at z₀ = 0.7H with H = 100  pc and n₀ = 10  cm^-3. In this case, blow-out happens after ~7.2 Myr.

In order to obtain the general dependence of the minimum number of SNe on the properties of the surrounding medium at the coordinate z_u, the thermal energy E_th,   SN appearing in t_SN (Eq. (32)) is replaced by 2/3·E_SN·N_SN as derived in Appendix A and the velocity at z_u (Eq. (35)) is solved for N_SN: (36)The velocity is fixed to be cm/s at the time of blowout and the constants k = 1 for the symmetric model and k = exp( − z₀/H) for the off-plane model, respectively, are used.

3.2.2. IMF- and wind-model

Next, we want to calculate the age, velocity and acceleration of a bubble driven by a time-dependent energy input rate, where different slopes of the IMF will be used. As a special case of this model, we can describe the evolution of a bubble powered by a constant energy input rate. Using a procedure similar as above, the characteristic timescale for the IMF-model in terms of the energy input rate coefficient L_IMF, the density and the scale height of the ISM is found to be (37)After inserting the thermal energy derived for the IMF-model (Eq. (19)) into Eq. (26) and solving the integral, one gets the time as a function of (38)Moreover, inserting the thermal energy for this model given by Eq. (19) into Eq. (25) and using the time t and the timescale t_IMF gives (39)where ϵ = 1/(δ + 3). Using further substitutions (40)and (41)yields a simplified expression of Eq. (39) (42)As it was done in the previous model, multiplying this equation by (Eq. (30)) results in the velocity at the top of the expanding superbubble (43)For the calculation of the acceleration , see Appendix B. Now that general expressions for velocity and acceleration of the bubble are known, we can estimate when a bubble starts to accelerate into the halo. The results for the IMF-model using different slopes and for the wind-model are listed in Table 2.

In order to present the blow-out timescales for superbubbles driven by a time-dependent energy input rate as a function of the number of OB-stars, N_OB, instead of the normalization constant N₀, which is contained in t_IMF (Eq. (38)), we have to make use of the fit N₀(N_OB) from Table 1. The coefficients of the power-law function found with a nonlinear least-squares fit provide an excellent fit to the relation with errors less than 1% for 2 ≤ N_OB ≤ 500 for all IMF-models. The timescales in Fig. 4 are shown for 2−500 association members, where the dashed lines represent the minimum number of OB-stars and the corresponding time until blow-out. Using IMFs with Γ₁ =  −1.15, Γ₂ =  −1.35, and Γ₃ =  −1.7, respectively, 27, 25, and 24 OB-stars are needed with an upper mass limit of 19, 18 and 17 M_⊙. The resulting blow-out timescales are 22.8, 20.1, and 15.2 Myr for a symmetric density distribution with H = 500  pc and n₀ = 0.5  cm^-3 (Table 2 and Fig. 4, left). In the case of an explosion at 70 pc above the plane in an ISM with H = 100  pc and n₀ = 10  cm^-3 (Table 2 and Fig. 4, right), the association needs to have at least ~8, 13, and 29 massive stars or an upper mass limit of 13, 15 and 18 M_⊙ (same order of IMF-slopes). Approximately 4.0, 3.5, and 2.6 Myr pass from the first SN-explosion until blow-out. We are interested again in obtaining an analytical expression for the minimum number of OB-stars to get blow-out. Equation (43) needs to be solved for the normalization constant, after L_IMF was replaced by Eq. (18), which yields (44)for any set of the parameters scale height and ISM-density. This has to be converted to a number of stars N_{OB,   blow} by using the fit presented in Table 1, which is obtained with the same fitting procedure as before. We find that our fits are very good for all IMF-slopes with an average deviation of ~1.4% in the range of 2 − 500 OB-stars. In Fig. 11 we compare the number of stars required for blow-out as a function of the scale height using certain values of the density for the SN- and the IMF-models. The discussion of the results is found in the next section.

Using a constant energy input rate (δ = 0), we obtain a dimensionless wind coefficient (45)instead of the normalization constant. Following Eq. (17), we can derive the minimum wind luminosity that fulfills the blow-out-criterion (46)in units of erg/s. We assume that the wind luminosity is the energy input of all SNe averaged over the period of Δt ≈ 20  Myr, the lifetime of the association, which is in general the time between the first and the last SN-explosion². Hence, the minimum wind luminosity can be converted into a total number of stars N_w,blow = L_w,blow·Δt/E_SN. Figure 5 shows the time elapsed since until the point of acceleration, for wind luminosities between 3 × 10³⁶ and 8 × 10³⁸ erg/s, corresponding to ~2−500 SNe calculated over a timespan of 20 Myr. The values for the minimum wind luminosity and corresponding ages of the SB are given by the dashed lines. For a symmetric expansion into the Lockman layer, an energy input rate of ~7.4 × 10³⁷  erg/s is necessary for blow-out which corresponds to ~47 SNe. Such a bubble needs 26.1 Myr until acceleration begins. Expansion of a SB into the low-scale height, high-density medium produced by an off-plane explosion at z₀ = 0.7H requires ~1.9 × 10³⁷  erg/s or about 12 SNe and takes ~4.7 Myr.

3.3. Rayleigh-Taylor instabilities in the shell

Infinitesimal perturbations at the interface between a denser fluid supported by a lighter fluid in a gravitational field generate waves with amplitudes growing exponentially with time in the initial phase. For an incompressible, inviscid, nonmagnetic fluid, the timescale τ_rti, characterizing the growth of the instability, results from a linear stability analysis combined with the conservation equations (47)where λ is the perturbation wavelength, g is the gravitational acceleration and ρ₁ and ρ₂ are the densities of the light and heavy fluid, respectively. In the case of a superbubble, where the dense shell is accelerated by the hot, tenuous gas, the most important perturbation wavelengths out of the Fourier spectrum are the ones, which are comparable to the thickness of the shell (which is itself a function of time), because these distort the shell so strongly that break-up can occur. Furthermore, identifying the gravitational acceleration with the acceleration at the top of the bubble results in an instability timescale at this coordinate (48)The density in the shell at the coordinate z_u is always given by for an adiabatic strong shock³. The density of the bubble interior as a function of time is with being the ejecta mass and being the volume of the bubble confined by the inner boundary of the shell. The new semimajor and semiminor axes simply reduced by the thickness of the shell are and . To calculate the thickness of the shell we use the fact, that the mass of the gas in a volume of undisturbed ISM is the same as that in the shell of swept-up ISM (if we neglect effects of evaporation, heat conduction or mass loading). In a symmetrically decreasing density distribution the bubble contour has an hourglass-shape, thus, due to symmetry, we need to include half of the bubble only in the calculation (49)Half of the shell mass is , which we obtain by integration over the density gradient. The bubble’s volume, confined by the surface A_ell of an ellipsoid, is included between the coordinates z = 0 and z_u. Using the new semimajor and semiminor axes – each itself a function of – in the equation of the ellipse means solving a double integral in the calculation of the volume (Eq. (23)). However, we can approximate by using an average thickness of the shell over the time in this case. The new semimajor axis is as large as ~91% of the regular value and the new semiminor axis ~90% (averaged over the time span ).

In the case of an off-plane explosion the thickness of the shell is given by (50)with the mass inside the shell with the surface A_ell of the complete ellipsoid.

Fig. 6 Fragmentation timescales for the SN-model; the dashed line marks the minimum number of SNe for a blow-out superbubble and the corresponding time until full fragmentation of the shell at the top of the bubble (for details see text and Tables 2 and 3). Left: double exponential layer with H = 500  pc and n0 = 0.5  cm-3 (symmetric model); right: exponentially stratified medium with H = 100  pc and n0 = 10  cm-3 (off-plane model).

Furthermore, we have to estimate the mass inside the hot bubble interior, which is the mass ejected by the SN-explosions. The total mass of each star belonging to the association goes into ejecta, except ~1.4   M_⊙ for a neutron star remnant, which is lower than the Oppenheimer-Volkoff limit, because most SN progenitors are lower mass stars among the massive stars. Thus, we have to consider the mass of per star with in units of solar masses and the total number of stars . In order to obtain the ejecta mass for the IMF-model as a function of time (i.e. time variable ), we fix the upper mass limit – which is related to the first SN-explosion – and introduce a variable lower mass limit . Again, connecting the mass of a star to its main sequence lifetime gives a time dependent ejecta mass. In order to account for the mass included in the mass interval (m_u − 1,m_u) for integral mass bins at time , we have to correct for the main sequence lifetime of stars with mass m_u − 1 (51)Integrating over the mass range of the association gives the ejecta mass of SNe exploded until some time (52)So we are able to derive the mass ejected until the time , where instabilities start to appear in the shell. Unfortunately, this formula does not hold for small associations, since all SNe may have exploded in a rather small time interval Δτ = τ_ms(8) − τ_ms(m_u − 1), possibly long before . In that case, the ejecta mass is replaced by the total mass of the ejecta after the last SN-explosion, otherwise the mass and thus the RTI-timescale would be overestimated. Including the ejecta mass (Eq. (52)) in the density ρ_in yields an instability timescale (Eq. (48)) as a function of the most massive star m_u in the IMF-model, but we prefer expressing it as a function of the total number of stars in an association. The relation between these two parameters is fitted with an approximation of the form (see Table 1) and has average errors of 1.2% over the range of 2−500 OB-association members.

In the case of the SN-model, all the mass of exploding stars (except of 1.4 M_⊙ per star) is released in the initial explosion. We find a useful power-law approximation for the total ejected mass of a star cluster as a function of the number of SNe (using an IMF with Γ =  −1.35) (53)Errors for this approximation are about 3% for small associations (2 OB-stars) and ~1% for large clusters (500 stars). It can be used to obtain the ejecta mass for the wind-model as well, but instead of the number of supernovae, N_SN, we have to include the energy input rate into the bubble until the point of acceleration in units of the standard SN energy . Calculating fragmentation this way only works for a wind, which is produced by averaging SNe and thus includes the mass of the exploded stars, but not for a true stellar wind.

Instabilities will dominate when the RTI-timescale at the top of the bubble becomes smaller than the dynamical timescale τ_dyn,I = a/ż_u (symmetric bubble) or τ_dyn,II = z_u/ż_u (off-plane model) of the system. In terms of the dimensionless time variable this happens at and we have to find the value, where for all models, which is shortly after the acceleration sets in. The exponentially growing instability is usually fully developed and the shell will break-up at when (see Table 2). This means that the exponentially growing amplitude of the perturbation has reached a size of e³ ≃ 20 times of the initial one, sufficiently large to assume full fragmentation.

Fig. 7 Same as Fig. 6, but for the IMF-model and indicating the minimum number of OB-stars in an association for different IMF slopes by the corresponding thin line (left: symmetric model; right: off-plane model).

Fig. 8 Same as Fig. 6, but for the wind model and a minimum constant energy input rate equal to a number of SNe distributed over 20 Myr (left: symmetric model; right: off-plane model).

Figure 6 shows the fragmentation timescale as a function of the total number of stars for the SN-model, Fig. 7 shows the same for the IMF-model for different slopes and finally, the timescale as a function of wind luminosity is presented in Fig. 8.

Due to the large fragmentation timescales for very small associations, the abscissa is chosen to range from 20−500 OB-stars for the symmetric model to highlight the behavior of this function. The off-plane model is shown again for 2−500 stars. For the wind-model these numbers correspond to 3.2 × 10³⁷−7.9 × 10³⁸  erg/s (symmetric model) and 3.2 × 10³⁶−7.9 × 10³⁸  erg/s (off-plane model).

In both cases of ISM density distribution, fragmentation is easily achieved for clusters with N_blow (i.e. those producing blow-out superbubbles) within 50 Myr, a reasonable timespan for fragmentation to take place before galactic rotation or turbulences have had a major influence on the bubble structure (see Figs. 6− 8, right panel) and no additional fragmentation criterion is needed. However, we want to check at this point whether the fragmenting SB can escape from the disk and accelerate into the halo. Thus, we compare the acceleration of the top of the bubble⁴ at with the gravitational acceleration near the galactic plane (we do not require the bubble to escape completely from the gravitational potential of the galaxy). We use g_z(R,z) =  −∂φ/∂z with φ, the disk potential in cylindrical coordinates (Miyamoto & Nagai 1975), given by (54)where R is the galactocentric radius, G is the gravitational constant, a = 7.258  kpc, b = 0.520  kpc, and the disk mass M = 2.547 × 10¹¹M_⊙ (Breitschwerdt et al. 1991) for a Milky-Way-type galaxy. This yields, e.g. for the galactocentric distance of the Sun at R_⊙ = 8.5  kpc a value of − g_z(R = R_⊙,z = 1   pc) = 3.47 × 10^-11  cm s^-2. We find for all models that SBs driven by the minimum blow-out energy, N_blow·E_SN, have an acceleration larger than − g_z(R,z = 1   pc) for all galactrocentric radii R. The complete set of results can be found in Table 3.

4. Application to the W4 superbubble

Winds of young massive stars of the cluster OCl 352 are supposed to be the energy source of the W4 superbubble. The cluster is located at a height of ~35  pc (Dennison et al. 1997) above the disk at a distance of ~2.35   kpc (West et al. 2007). In the following, we want to apply our off-plane wind-model to this superbubble.

According to West et al. (2007) only the structure above OCl 352 should be termed superbubble or chimney (G134.4+3.85), while the lower part is the W4 loop. The superbubble is in the process of evolving into a chimney, because the ionized shell of G134.4+3.85 already started to fragment at the top of the bubble, where the shell is expected to break-up first due to instabilities. From HI observations they derive a scale height of the ambient ISM of 140 ± 40  pc and they get bubble coordinates of pc (i.e. pc) and pc. With this ratio of z_u,I/r_max = 2.57 the bubble has reached an evolutionary parameter of in the Kompaneets model. At this time, the bubble’s extension from the star cluster to the top of the bubble is 6.66H and from the cluster to the bottom it is about − 1.35H. Comparing this to the physical dimensions yields a scale height of only ~32  pc in this region of the Milky Way. The scale height of 140  pc cannot be confirmed by the KA, since this would result in a value for z_d of almost 190 pc, i.e. the bottom of the shell reaching 150 pc below the midplane. Our model predicts that z_d should have reached only ~8  pc below the Galactic plane. The cluster itself is found at z₀ = 1.09H. We obtain an age of the SB of 1.8 Myr or 2.3 Myr, taking a density of n₀ = 5  cm^-3 and n₀ = 10  cm^-3, respectively, and using the energy input rate of 3 × 10³⁷  erg/s (Normandeau et al 1996).

Fig. 9 Left: color image of W4 (West et al. 2007) combining HI, Hα and infrared data. For comparison with our model, due to the faint structures and low contrast, the overlay is shown in a separate figure (right panel). The height of the star cluster above the plane is z0 = 35  pc, marked by the black dot. Right: Kompaneets bubble at overlaid on the same image. The upper part of the contour fits quite well (solid line), while the part of the shell below the cluster is not seen in the observations (dotted line).

According to West et al. (2007), Basu et al.’s (1999) way of treating the W4 SB and W4 loop as one entity (and fitting a Kompaneets model to that) seems inappropriate. However, we wanted to see if our model could improve their findings. Since Basu et al. (1999) do not use an off-plane model, we first checked their results with our model for an explosion in the midplane. We infer an age of ~2.3 Myr () when inserting all the values they use, which corresponds quite well to their derived age of 2.5 Myr. In this case, the distance of the top of the bubble to the star cluster of z₁ = 246  pc was used, which corresponds to the coordinate in the off-plane model including the offset z₀. In fact, the coordinate z_u,I = 211  pc should be used as the distance of the top of the bubble to the cluster when determining the aspect ratio, which is z_u,I/r_max = 211/74 = 2.85. With a slightly larger elongation than in the previous calculation, the bubble has reached the time parameter . With a height of 7.81H from z₀ to z_u,II, this yields a scale height of H ≅ 27  pc close to Basu et al.’s value of 25 pc. But, the fact that the bubble is shifted to a lower density environment changes the age of the bubble significantly. With an offset of z₀ = 1.30H, the age of the bubble is only ~1.7  Myr using n₀ = 10  cm^-3. This age is found within the previous estimates for the cluster of 1.3 − 2.5  Myr (cf. Dennison et al. 1997 and references therein) and is therefore consistent with the assumption that the bubble is blown by the wind of the O-stars in the cluster.

As the W4 loop and the W4 bubble are not dynamically connected, we concentrate in our approach on fitting only the SB shell above the cluster with our off-plane model to Fig. 10 of West et al. (2007). If we use and the coordinates given in West et al. (2007) we find that the contour of the model does not match the shell in the observations very well. The model would fit almost perfectly, if it was shifted upwards by about one scale height, but then the position of the star cluster would be located outside the contour like in Fig. 11 of West et al. (2007). The somewhat more elongated bubble at respresents nicely the shell above OCl 352 (Fig. 9, right) with the cluster – although not exactly in the center – matching the offset z₀ in the observations. We suggest that the part of the shell below OCl 352 (Fig. 9, right, dashed line) was decelerated due to the presence of the W4 loop and appears now flattened or has even merged with the upwards expanding part of the W4 loop.

Applying our criterion from Sect. 3, we find that blow-out of a bubble with the association at around one scale height in an ISM with H = 27  pc and n₀ = 10  cm^-3 is guaranteed for a wind luminosity as low as 8 × 10³⁵  erg/s. We thus support Basu et al.’s statement that the bubble is already on its way of blowing out into the Galactic halo. Even if the bubble was not shifted above the plane, around 4 × 10³⁶  erg/s would be sufficient for blow-out, which is well below the 1.67 × 10³⁸  erg/s found by MMN89 (as cited in West et al. 2007). The acceleration of the bubble, i.e. blow-out, has started already at , more than 1 Myr ago.

Fig. 10 Symmetric Kompaneets bubble at overlaid on the absorption map of neutral interstellar NaI of Lallement et al. (2003). Isodensity contours with an equivalent width of 20 m (inner contour) and 50 m (outer contour) represent the rarefied cavity.

4.1. Local Bubble

The Local Bubble (LB) was most likely produced by 14 − 20 SNe, which exploded 10 − 15 Myr ago as the Pleiades subgroup B1 was moving through our Galactic neighborhood (Berghöfer & Breitschwerdt 2002; Fuchs et al. 2006). Numerical simulations suggest that 19 SNe are responsible and the bubble is ~13.6  Myr old with the last SN having exploded 0.5 Myr ago (Breitschwerdt & de Avillez 2006; de Avillez & Breitschwerdt 2012). Observations (Welsh et al. 1999; Sfeir et al. 1999) show that the bubble is not confined at higher galactic latitudes and thus, should be termed “Local Chimney”, but an elongated structure (the chimney walls) still exists, extending to ~250 − 400  pc above and below the midplane (Lallement et al. 2003). Although the Local Bubble is tilted about 20° to the midplane and expands perpendicularly to the plane of Gould’s Belt, we apply here our symmetric model with a time-dependent energy input rate.

The walls of the Local Bubble (Fig. 5 of Lallement et al. 2003) can be fit very well with a Kompaneets model using an evolutionary parameter of . With the coordinate corresponding to ~370  pc, this yields a scale height of H = 80  pc. The radius of the bubble in the plane is almost 150 pc and r_max has an extension of 180 pc at in our model, corresponding to the observations (see Fig. 10). Since it was suggested that 19 SNe exploded in a time interval of 13.1 Myr and with a given lower mass limit of m_l = 8.2   M_⊙ (Fuchs et al. 2006), we can infer – using main sequence lifetimes of the stars – that the upper mass boundary should be m_u = 20.9   M_⊙ (independent of the IMF). So far, we assumed that there is exactly one star in the mass bin (m_u − 1,m_u) for general modeling. But this way of mass binning has to be modified and the mass interval of the most massive star needs to be adopted, because we know both m_u and N_OB in the case of the LB. Using an IMF with Γ₂ =  −1.35, the mass bin containing the most massive star is N(m_u − 1.9,m_u) = 1 and a normalization constant of N₀ = 613 is found. This simply means, that the remaining 18 stars are distributed within the mass interval from 8.2 − 19   M_⊙. A density of the undisturbed ISM of n₀ ≈ 7  cm^-3 is obtained to infer the presumable age of the LB. Acceleration of the shell started already at in this configuration, which was 3.3 Myr after the first SN exploded. The velocity of the bubble at this time was about 20 km s^-1; thus the bubble fulfills the Kompaneets criterion of blow-out into the halo. Rayleigh-Taylor instabilities started to appear at the top of the bubble at , about 3.4 Myr after the initial SN-explosion and full fragmentation took place at a time of 5.2 Myr. With an acceleration of ~2.5 × 10^-9  cm s^-2 at exceeding the gravitational acceleration near the galactic plane by two orders of magnitude, the fragmenting shell will not fall back onto the disk. Actually, it is even higher than the local vertical component of the gravitational acceleration − g_z(R = 8.5  kpc, z = 114   pc) = 2.37 × 10^-9 cm s^-2 such that blow-out of the Galaxy’s gravitational potential could be achieved. But since the LB is already disrupted on its poles before having reached one scale height of the Lockman layer, the bubble will not be able to expand to high-z regions, but the ejected material will probably mix with the lower halo gas, leaving behind a “Local Chimney”.

5. Discussion

5.1. Comparisons of the models: blow-out

Fig. 11 Numbers of SNe/OB-stars needed for blow-out are compared for SN-, IMF- and wind-model (energy input rate over 20 Myr gives the number of SNe for the wind-model). a) Symmetric case, midplane number density: 0.5 cm-3, scale height range: 220−500 pc; b) symmetric case, midplane number density: 5 cm-3, scale height range: 80−500 pc; c) off-plane explosion at z0 = 0.7H, midplane number density: 1 cm-3, scale height range: 270−500 pc; d) off-plane explosion at z0 = 0.7H, midplane number density: 10 cm-3, scale height range: 90 − 500 pc.

In order to compare the efficiencies of the three different models, we have to take a look at the number of SNe needed for blow-out and at the corresponding timescales.

For an ISM with properties of the Lockman layer of the Galaxy, the SN-model requires the largest number of SNe for blow-out. If taking into account an OB-associaton’s lifetime of 20 Myr, the wind model’s efficiency is between the IMF- and the SN-model. The IMF-model is the most efficient one and needs only half the number of SNe to explode over the whole lifetime of the SN-model. In fact, the IMF with a slope of Γ₃ =  −1.7 needs the lowest number of OB-stars to produce a blow-out SB, followed by the one with Γ₂ =  −1.35 and Γ₁ =  −1.15, but the numbers are approximately the same (see Table 2). Additionally, the blow-out-timescales for the IMF-models are shortest with a steeper slope yielding a shorter timescale.

For the low scale height-high density ISM used in combination with an off-plane explosion, the numbers of N_blow are lower in general. The IMF-model with Γ₃ =  −1.7 is an exception, since the number of OB-stars N_blow is now higher than in the symmetric case. The SN-model is the most efficient one with only three SNe needed for blow-out, and also the behavior of the IMF-slopes is the opposite as above. The minimum energy input for the wind-model is somewhere in between, similar to that of the IMF-model with Γ₂. But the timescales for blow-out are unaffected by this change and are still longest for the SN-model and shortest for the Γ₃-model. This is simply because the acceleration starts later in terms of the time variable for the SN-model, thus a larger volume has to be carved out by the SN-explosions.

We find the mathematical explanation for the blow-out efficiency of the models when looking at Eqs. (36), (44), and (46), which are used to describe N_{SN,   min}, N_{OB,   min} and L_w,   min, respectively. The dependence on the density is linear for the SN- and wind-model and slightly increasing with an exponent b₂ from Table 1 for the IMF-models. The SN-model is strongly influenced by the scale height, since N_{SN,   min} ∝ H³. For the wind-model, L_w,   min ∝ H² is found. The different IMF-slopes give us the following relations: N_{OB,   min} ∝ H^1.9 for Γ₁, H^1.6 for Γ₂, and finally H^1.2 for Γ₃.

In Fig. 11 we illustrate these facts by comparing the symmetric and the off-plane model for a low-density and high-density ISM each. The values were chosen such that the Lockman-layer density of n₀ = 0.5  cm^-3, (a), and a high-density layer (d) with n₀ = 10  cm^-3, were used as reference (see Sect. 3). According to that, we wanted to investigate what happens to N_min with a density 10 times higher, i.e. n₀ = 5  cm^-3, in the symmetric case (b) or a 10 times lower density, i.e. n₀ = 1  cm^-3, in the off-plane scenario (c). The plots start with a scale height on the x-axis, where at least N_min ≥ 2 are needed, which is the smallest association we use in this paper, and end at a scale height of 500 pc. We find that the models behave similarly throughout this scale height range irrespective of the density distribution, but the order of the efficiency of all models changes at H ~ 300−450  pc, which marks the transition of a low-scale height to a high-scale height medium. In the low-scale height regime (until ~350  pc) the wind-model is least efficient, followed by the Γ₃-model and the SN-model. The Γ₂-IMF is slightly better and an IMF with Γ₁ needs the lowest number of SNe. For a scale height above 400 pc, the sequence of the models is the same as the one calculated for the Lockman layer (described above), except for the off-plane model with n₀ = 1  cm^-3, where the IMF-models do not intersect before ~500  pc according to Fig. 11c.

5.2. Comparison of the models: fragmentation

Next, we want to compare the timescales until fragmentation and explore the efficiencies of the different models. The energy input needed for a bubble to achieve fragmentation is given by the minimum blow-out energy, i.e. the energy required for reaching a certain velocity v_acc at (see Table 2). The trend until the beginning of the fragmentation process (i.e. at ) is the same as that for blow-out for both ISM cases, in terms of and in absolute timescale: the SN-model is the least efficient one and an IMF with a steep slope yields the best results (shown in Table 3). Also, it takes the SN-model the longest total time ) until full fragmentation occurs for bubbles driven by N_blow. But the IMF-model with Γ₂ =  −1.35 now shows the lowest total timescale followed by Γ₁ =  −1.15 with no obvious correlation among the models. We find that the SN-model is the fastest model in completing the fragmentation process – which means from until – and a steep IMF (Γ₃) yields the highest values of this timespan. Since this is the opposite to the efficiency sequence until , no general trend for the fragmentation timescales can be seen.

Concerning the interpretation of these facts one has to look at the RT-instability timescale itself. It is decreasing more rapidly for the SN-model and a flat IMF-model, because the acceleration is increasing in a slightly shallower way than for steep IMFs. Thus, a value of 1/3   τ_rti is reached easier for a flat IMF than for a steep one.

6. Summary and conclusions

We have developed analytical models based on the Kompaneets approximation (KA) in order to derive in a fairly simple and straightforward manner the physical parameters of observed superbubbles and their ambient medium and to gain physical insight into the blow-out phenomenon associated with star forming regions.

In this paper we have deliberately refrained from building more complex models, which e.g. include stellar wind and Wolf-Rayet wind phases, because we have put our focus on the important dynamical phenomena of blow-out and fragmentation of the outer shell, and their dependence on the energy input source over time. A more detailed description of superbubble evolution models including further stellar evolutionary phases will be the subject of a forthcoming paper.

In our work the key aspect was to work out analytically the dynamics of (unfragmented) superbubbles for different energy input modes. We modified the KA to implement a more realistic way of energy input, i.e. modeling the time sequence of exploding stars in an OB association by including the main sequence lifetime of the massive stars and describing the numbers per mass interval by an initial mass function. We tested three different IMF slopes and also compared the IMF-model to a simple SN-model with an instantaneous release of energy and to a wind model with a constant energy input rate. Two different density distributions of the ISM were applied, a symmetric medium with parameters of the Lockman layer and a high-density, low-scale height pure exponential atmosphere with the star cluster dislocated from the galactic plane. Velocity and acceleration of the shock front can be calculated analytically and the question how many SNe are needed for blow-out into a galactic halo can be answered. The exact position of the outer shock in scale height units when the acceleration starts can be given. Furthermore, the timescale for the development of Rayleigh-Taylor instabilities in a SB shell is calculated, and thus, a fragmentation timescale can be derived. The overall pattern shows that at larger scale heights (H > 400  pc), independent of the ISM density, the SN-model needs the highest energy input, followed by the wind-model, whereas an IMF with a steep slope is the most efficient one. The same ranking applies to the blow-out timescales of these models (Figs. 3−5 and Table 2). At low scale heights (H ~ 100  pc) and moderate or high densities (n₀ ≥ 5  cm^-3), the picture changes completely, i.e. the SN-model requires the lowest energy input and the IMF-model with Γ₃ =  −1.7 is the least efficient one. The explanation is that a single release of energy is more powerful in sweeping up a thin layer of ISM, whereas it is easier for the IMF model with an increasing energy input with time (L_SB ∝ t^δ, 0 < δ < 1, Eq. (17)) to sustain the supply of energy over a larger distance. When comparing fragmentation timescales (Figs. 6−8 and Table 3), the IMF-models exhibit the lowest values, favouring flatter IMFs with Γ₁ and Γ₂. In terms of and absolute time, fragmentation happens first for the Γ₂-model for SBs driven by the minimum number of SNe for blow-out.

Still, the KA is a rather simple model. It does not account for magnetic fields, ambient pressure, inertia and evaporation of the shell. Also a galactic gravitational field and cooling of the shocked gas inside the cavity are neglected.

In our model, we include galactic gravity in a rudimentary way: if a bubble fulfills the blow-out criterion and additionally the shell’s acceleration at the top of the bubble at the time of fragmentation exceeds the vertical component of the gravitational acceleration in the disk, the SN-ejecta and fragmenting shell are expelled into the halo. We find that this is true for all bubbles created by an energy input of N_blow·E_SN.

Further analysis of a fragmenting superbubble would also involve to calculate the motion of the shell fragments, and ultimately the dynamics of a galactic fountain like e.g. in Spitoni et al. (2008). A simple ballistic treatment of the motion of blob fragments in a gravitational potential would, however, be too simplistic, as they would for some time experience a drag force due to the outflowing hot bubble interior. The drag would be proportional to the ram pressure of the hot gas and the cross section of a blob, which would, to lowest order, be proportional to the thickness of the fragmented shell; and finally it would depend on the geometry of the blob, which might be taken as spherical. However, due to compressibility, bow-shocks and head-tail structures might subsequently be formed, so that for a realistic treatment numerical simulations would be the best choice.

According to MM88 cooling of the bubble interior can be neglected for Milky Way type parameters of the ISM, but should be taken into account for smaller OB-associations or a dense and cool ISM. In general, cooling is not important as long as the timescale for radiative cooling is large compared to the characteristic dynamical timescale of the superbubble. Including a magnetic field would be quite important, but this goes beyond the scope of our analytical model and is left to numerical simulations. Ferriere et al. (1991) find that the presence of a magnetic field could slow down the expansion of a superbubble. Stil et al. (2009) argue that the scale height and age of a bubble are underestimated by ~50% when using a Kompaneets model without magnetic fields. However, they cannot produce such narrow superbubbles like W4 with their MHD simulations. Moreover, 3D high resolution numerical simulations (de Avillez & Breitschwerdt 2005) show that magnetic tension forces are much less efficient in 3D than in 2D in holding back the expanding bubble.

A slightly slower growth of a bubble would be also achieved by taking into account the inertia of the cold massive shell in the calculations. MM88 find a difference of ~10% in radius after comparison with the models of Schiano (1985), which neglect inertia. Also due to the ambient pressure of the ISM superbubbles should expand more slowly as it was suggested by Oey & García-Segura (2004). In our analytical calculations we find that these effects have to be compensated when we make comparisons to observed bubbles by including a rather high ISM density to prevent SBs from expanding too fast. Furthermore, a clumpy ambient medium cannot be considered by the KA.

We applied our models to the W4 superbubble and the Local Bubble, both in the Milky Way. It is certainly not easy to compare a simple model with an observed SB, which is not isolated. In the region of the W3/W4/W5 bubbles, several complexes and clouds are found and multiple epochs of star formation make it difficult to distinguish, which cluster has formed which bubble at what time. However, it is most likely that the cluster OCl 352 is responsible for driving the evolution of the bubble (West et al. 2007). Oey et al. (2005) suggest that winds or SNe of previous stellar generations are responsible for earlier clearing of this region, which could explain the low scale height of around 30 pc. Our calculations suggest that the bubble is younger than found by other authors, which is due to the offset of the association more than one scale height above the Galactic plane. Shifting to lower densities makes it easier to produce a blow-out superbubble in a shorter timescale. This is an important result and should be included in the models.

The Local Bubble is one of the rare cases for a double-sided bubble, which can be tested with our symmetric superbubble model. From geometrical properties, we estimate an ISM scale height of ~80  pc. In order to reproduce size and age of the bubble correctly, we deduce from our models that it was an intermediate density region (n₀ ~ 7  cm^-3) before the first SN-explosion around 14 Myrs ago. Also this place in the Milky Way is very complex (neighboring Loop I superbubble), which can’t be included in our modeling of superbubbles.

We conclude that blow-out energies derived in this paper are lower thresholds and might be higher if e.g. magnetic fields play a role. Accordingly, fitting the models to observed bubbles gives an upper limit for densities of the ambient ISM prior to the first SN-explosion. Observers are encouraged to use the model presented here for deriving important physical parameters of e.g. the energy input sources (number of OB stars, richness of cluster etc.), scale heights, dynamical timescales among other quantities. The solutions of the equations derived in detail here are easy to obtain by simple mathematical programs. Theorists may find it useful to compare our analytic results to high resolution numerical simulations in order to separate more complex effects, such as turbulence, mass loading, magnetic fields etc. from basic physical effects, incorporated in our model.

Acknowledgments

V.B. acknowledges support from the Austrian Academy of Sciences, the University of Vienna, and the Austrian FWF, and the Zentrum für Astronomie und Astrophysik (TU Berlin) for financial help during several short term visits.

References

Appendix A: Thermal energy in the bubble

The fraction E_th of the total energy, which is converted into thermal energy at the inner shock, is derived. The equation for the thermal energy inside the hot bubble is (A.1)where is the volume of a spherical remnant and R(t) is the radius of the outer shock. Since Basu et al. (1999) have shown that the thermal energy in the hot bubble is comparable to that of a homogeneous one until late evolutionary stages, the simplification of a spherical volume is used in the calculations. The pressure P(t) is obtained by taking into account momentum conservation of the bubble shell (see, for example, Castor et al 1975; Weaver et al 1977) (A.2)with homogeneous ambient density ρ₀. This equation is rearranged which gives (A.3)Combining this equation with Eq. (A.1) for the thermal energy results in (A.4)

Appendix A.1: SN-model

In the case of a bubble created by a single explosion, the radius of the outer shock is (e.g. Clarke & Carswell 2007) (A.5)Inserting R_SN(t) and its derivatives Ṙ and to Eq. (A.4) results in E_th,   SN = 2/3·E_SN·N_SN.

Appendix A.2: IMF-model

In order to describe the evolution of a superbubble driven by a time-dependent energy input rate – in addition to the equations above – energy conservation of the hot wind gas has to be considered for this problem (A.6)with the total thermal energy in this region given by Eq. (A.1). L_SB = L_IMF·t^δ is the energy input rate delivered by sequential SN-explosions according to an Initial Mass Function. Additionally it is assumed that the radius of the bubble in a self-similar flow scales with time like R_IMF(t) = A_s·t^μ. The constant is (A.7)and the exponent (e.g. Berghöfer & Breitschwerdt 2002). Now, the fraction E_th(t) of the total energy, which is converted into thermal energy at the inner shock, can be calculated. Inserting R_IMF(t) and its derivatives into Eq. (A.4) yields

Appendix B: Acceleration of the top of the bubble

The derivative of the bubble volume with respect to is needed in the following calculations, which is the same for all models presented below.

We use the derivative of the series expansion of the volume V_I in the symmetric case. For the offcenter model we find (B.1)

Appendix B.1: SN-model

In order to obtain the acceleration at the top of the bubble, we need to calculate the first term on the RHS of Eq. (31) which is the derivative of Eq. (35) with respect to (B.2)Multiplying the above equation by (Eq. (34)) yields the acceleration (in units of cm/s²) (B.3)

Appendix B.2: IMF- and wind-model

Next, we want to determine the acceleration at the top of the bubble for the IMF-model. Using δ = 0 gives the acceleration for the wind-model. Already, is known (Eq. (42)), therefore only has to be calculated. Determining the derivative of Eq. (43) with respect to yields (B.4)For the complete expression of the acceleration in the IMF-model, the equation above just has to be multiplied by (Eq. (42)).

Appendix C: List of variables and parameters

γ	...	ratio of specific heats (γ = 5/3)
Γ	...	slope of the IMF (Γ0 =  −0.932, Γ1 =  −1.15,
		Γ2 =  −1.35, Γ3 =  −1.7)
H	...	scale height
L IMF	...	IMF energy input rate coefficient
n 0	...	number density at galactic midplane
N 0	...	IMF normalization constant
N blow	...	minimum number of stars for blow-out
τ dyn	...	dynamical timescale
τ rti	...	Rayleigh-Taylor instability timescale
	...	bubble volume (symmetric model)
	...	bubble volume (off-plane model)
v acc	...	velocity at time of blow-out
y	...	transformed time variable
	...	transformed time variable in scale height units
	...	time of blow-out
	...	time of fragmentation of the shell
	...	time of onset of Rayleigh-Taylor instabilities
z 0	...	center of the explosion
z d	...	bottom of the bubble
z u	...	top of the bubble
	...	velocity of the bubble at zu
	...	acceleration of the bubble at zu

All Tables

All Figures

	Fig. 1 Position of the shock front in the symmetric model at certain values of the dimensionless time variable , 1.0, 1.4, 1.7, 1.9, 1.98, and 2.0 with the energy source in the galactic midplane.
In the text

	Fig. 2 Same as Fig. 1, but for the off-plane model. The energy source is located at z0 = 0.7H above the plane (indicated by the dashed line).
In the text

	Fig. 3 Blow-out timescales at the coordinate zu for the SN-model. Left: double exponential layer with H = 500  pc and n0 = 0.5  cm-3 (symmetric model); right: exponentially stratified medium with H = 100  pc and n0 = 10  cm-3 (off-plane model). The dashed lines indicate the number of SNe required for blow-out and the corresponding blow-out time.
In the text

	Fig. 4 Same as Fig. 3, but for the IMF-model (left: symmetric model; right: off-plane model). Thin lines indicate the number of OB-stars needed for blow-out. Γ is the IMF exponent (for details, see text).
In the text

	Fig. 5 Same as Fig. 3, but as a function of a constant energy input rate (left: symmetric model; right: off-plane model). The number of SNe is obtained by converting the wind luminosity to an energy input averaged over a time interval of 20 Myr, from the explosion of the most massive star until the last SN.
In the text

Fig. 6 Fragmentation timescales for the SN-model; the dashed line marks the minimum number of SNe for a blow-out superbubble and the corresponding time until full fragmentation of the shell at the top of the bubble (for details see text and Tables 2 and 3). Left: double exponential layer with H = 500  pc and n0 = 0.5  cm-3 (symmetric model); right: exponentially stratified medium with H = 100  pc and n0 = 10  cm-3 (off-plane model).

In the text

	Fig. 7 Same as Fig. 6, but for the IMF-model and indicating the minimum number of OB-stars in an association for different IMF slopes by the corresponding thin line (left: symmetric model; right: off-plane model).
In the text

	Fig. 8 Same as Fig. 6, but for the wind model and a minimum constant energy input rate equal to a number of SNe distributed over 20 Myr (left: symmetric model; right: off-plane model).
In the text

Fig. 9 Left: color image of W4 (West et al. 2007) combining HI, Hα and infrared data. For comparison with our model, due to the faint structures and low contrast, the overlay is shown in a separate figure (right panel). The height of the star cluster above the plane is z0 = 35  pc, marked by the black dot. Right: Kompaneets bubble at overlaid on the same image. The upper part of the contour fits quite well (solid line), while the part of the shell below the cluster is not seen in the observations (dotted line).

In the text

	Fig. 10 Symmetric Kompaneets bubble at overlaid on the absorption map of neutral interstellar NaI of Lallement et al. (2003). Isodensity contours with an equivalent width of 20 m (inner contour) and 50 m (outer contour) represent the rarefied cavity.
In the text

Fig. 11 Numbers of SNe/OB-stars needed for blow-out are compared for SN-, IMF- and wind-model (energy input rate over 20 Myr gives the number of SNe for the wind-model). a) Symmetric case, midplane number density: 0.5 cm-3, scale height range: 220−500 pc; b) symmetric case, midplane number density: 5 cm-3, scale height range: 80−500 pc; c) off-plane explosion at z0 = 0.7H, midplane number density: 1 cm-3, scale height range: 270−500 pc; d) off-plane explosion at z0 = 0.7H, midplane number density: 10 cm-3, scale height range: 90 − 500 pc.

In the text