ABSTRACT
Interstellar bubbles appear to be smaller in observations than expected from calculations. Instabilities at the shell boundaries create three-dimensional (3D) effects and are probably responsible for part of this discrepancy. We investigate instabilities and dynamics in superbubbles using 3D hydrodynamics simulations with time-resolved energy input from massive stars, including supernova explosions. We find that the superbubble shells are accelerated by supernova explosions, coincident with substantial brightening in soft X-ray emission. In between the explosions, the superbubbles lose energy efficiently, approaching the momentum-conserving snowplow limit. This and enhanced radiative losses due to instabilities reduce the expansion compared to the corresponding radiative bubbles in pressure-driven snowplow models with constant energy input. We note generally good agreement with observations of superbubbles and some open issues. In particular, there are hints that the shell velocities in the X-ray-bright phases are underpredicted.
Export citation and abstract BibTeX RIS
1. INTRODUCTION
The interstellar medium is commonly seen to show cavities delimited by arcs and shells in the Milky Way (e.g., Churchwell et al. 2006; Paladini et al. 2012) as well as in other star-forming galaxies (e.g., Dunne et al. 2001; Bagetakos et al. 2011). Bubble centers seem to hold individual stars and stellar explosion sites (e.g., Gruendl et al. 2000; Green 2004) as well as stellar groups and associations (superbubbles, e.g., Cash et al. 1980; Maciejewski et al. 1996; Breitschwerdt & de Avillez 2006; Jaskot et al. 2011). Interstellar bubbles are evidence of local energy input due to high-energy photons (e.g., Oort & Spitzer 1955; Dale & Bonnell 2011; Krasnobaev et al. 2014) and/or mechanical energy from winds and explosions (e.g., Freyer et al. 2003; van Veelen et al. 2009; Ntormousi et al. 2011; van Marle et al. 2012; Georgy et al. 2013; Krause et al. 2013; Rogers & Pittard 2013).
The expansion of interstellar bubbles was initially explained through self-similar, analytic descriptions (e.g., Sedov 1959, point explosion; Castor et al. 1975; Weaver et al. 1977, constant power wind with and without cooling). Since then, much effort has been devoted to incorporate additional complexity in such descriptions, such as non-uniform ambient density distributions (e.g., García-Segura & Mac Low 1995) or mass loading by entrained clouds (e.g., Pittard et al. 2001); see Ostriker & McKee (1988) for a review and more applications. The structure of bubble interiors in terms of density and composition (e.g., van Veelen et al. 2009; Georgy et al. 2013), as well as the shell structure (e.g., Vishniac 1983; Mac Low & Norman 1993; Ntormousi et al. 2011; van Marle & Keppens 2012; Krause et al. 2013; Pittard 2013), have been inferred from further analytical studies and multi-dimensional simulations.
Observations of interstellar bubbles and constraints on their stellar content showed significant deviations from theory: the bubbles appear to be smaller and thus appear to lose more energy during their evolution than predicted in the models (e.g., Drissen et al. 1995; García-Segura & Mac Low 1995; Oey 1996; Oey & García-Segura 2004; Butt & Bykov 2008; Harper-Clark & Murray 2009; Bruhweiler et al. 2010). X-ray-bright superbubbles are observed to have higher shell velocities than expected for their size and evolution time (e.g., Oey 1996).
Possible solutions to the first discrepancy include a correction for the underestimate of the ambient pressure (Oey & García-Segura 2004) and blow out in the case of superbubbles (e.g., Mac Low & McCray 1988; Breitschwerdt & de Avillez 2006). The latter accounts for an observational bias because shells are best seen in the direction where the density is the highest and, consequently, where the expansion is slowest. Harper-Clark & Murray (2009) have suggested that energy may be lost from bubbles via the leakage of hot gas through holes in the shell. In simulations, such holes have not been found, however, even in strong density inhomogeneities (Pittard 2013).
Time variability of the energy input has been suggested as an explanation for the high velocity of X-ray-bright superbubbles. Oey & García-Segura (2004) show in a one-dimensional (1D) hydrodynamical simulation including radiative cooling and heating due to photoionization that the supershells may be accelerated to the observed velocities after a supernova has exploded in a bubble created previously by the massive-star winds.
Time variability, however, also has another effect, which is not captured in the 1D models. It is well known that shells are unstable to the Rayleigh–Taylor instability whenever the shell accelerates, and to the Vishniac instability (Vishniac 1983) whenever it decelerates. The Vishniac instability leads to shell clumping, and thus quite possibly influences star formation. The Rayleigh–Taylor instability causes filaments of shell gas to enter the superbubble, where it subsequently mixes with the bubble gas, thereby reducing its temperature and increasing its density. We recently demonstrated both effects in three-dimensional (3D) hydrodynamics simulations of superbubbles. Even in the more complex case of three closely spaced massive stars, the filamentary network caused by the Vishniac instability is established on the shell early on, and is present up to the end of the simulations (Krause et al. 2013, hereafter Paper I).
The instabilities lead to an intermediate density region inside the shell (Krause et al. 2014, hereafter Paper II). The density is high enough to produce an X-ray luminosity comparable to the one seen in X-ray-bright superbubbles. We calculated the spectrum in detail, assuming equilibrium ionization, and showed that it is more complex than a thermal spectrum for a single temperature. This is in good agreement with observations, as is the edge-brightened X-ray appearance. Here, we show that time variability of the energy input also leads to a fundamental change in bubble dynamics: most of the time, superbubbles are in a momentum-conserving rather than a pressure-driven snowplow phase. Additional energy loss occurs in the mixing region because densities and temperatures favorable for atomic line cooling are reached in certain locations. In Section 2, we show that when also taking this into account, the momentum-conserving snowplow loses energy most efficiently among the power-law solutions, i.e., also faster than the pressure-driven snowplow. In Section 3, we demonstrate that our simulated superbubbles are almost always close to a momentum-conserving snowplow phase. They therefore lose more energy and have smaller bubble radii than the classical pressure-driven snowplows. In Section 4, we show that the X-ray luminosity correlates strongly with shell acceleration, and that a supershell may be accelerated beyond the velocity expected from even an adiabatic Weaver model for short time intervals.
2. BUBBLE DYNAMICS IN THE THIN-SHELL APPROXIMATION
Cooling in the shocked ambient medium of interstellar bubbles may lead to the formation of a thin and cool shell (snowplow phase). Even when radiative energy losses are negligible, the strong decline of the density inward of the outer shock in the classical solutions (e.g., Sedov 1959) justifies the thin-shell approximation for dynamical purposes (Kompaneets 1960; Bisnovatyi-Kogan & Silich 1995). The thin-shell model describes interstellar bubbles as systems of two homogeneous regions: the spherical interior carries the thermal energy, Et(t), which is subject to stellar energy input and adiabatic losses; its pressure, p(t) = (γ − 1)Et/V(r), where V(r) is the bubble volume, regulates the dynamics of the surrounding shell via the momentum equation
Here, Mv denotes the shell's momentum and A(r) is the spherical surface area. The shell carries the kinetic energy:
This situation is also called a pressure-driven snowplow (Ostriker & McKee 1988, their Section VI.A), in contrast to the momentum-driven snowplow, where p(t) = 0. For constant energy injection (Weaver model), the pressure-driven snowplow dissipates a constant fraction of the input power at the leading radiative shock, and thus energy accumulates with time. The momentum-conserving snowplow loses energy at a rate ∝t−3/4 (Ostriker & McKee 1988, their Section VI.A).
We modified this description allowing for arbitrary energy loss in the leading shock, as well as the bubble interior by requiring only consistent dynamics regarding bubble pressure and shell acceleration, rather than summing up kinetic energy and adiabatic expansion energy loss (Krause 2003). Therefore, only Equations (1) and (2) are solved.
In Krause (2003), we have shown that these equations may be integrated for arbitrary E(t) and spherically symmetric density ρ(r):
where M(r) = 4πρ(r) r3/3. Inserting a constant energy accumulation rate, E(t) = Lt, and a constant ambient density, ρ(r) = ρ0, recovers the classical solution by Weaver et al. (1977): r = α(L/ρ0)1/5t3/5, but with α = 0.83. This compares to α = 0.76 and α = 0.88 3 in Weaver et al. (1977) for radiative and non-radiative shells, respectively.
For further analysis of the simulation results, we now discuss suitable special solutions of Equation (3). For constant ambient density and an energy injection power law of E(t) = ctd, the shell radius evaluates to
which is valid for d > −1. The shell velocity v(t) follows by differentiation with respect to t. From this, we obtain the kinetic energy of the shell, Ek(t) = 0.5M(r(t))v(t)2, and the bubble's kinetic energy fraction:
Figure 1 shows the kinetic energy fraction as a function of the power-law index d for energy. Standard cases are marked: isolated supernova, 
k(d = 0) = 0.4, and constant luminosity wind, k(d = 1) = 0.3. The total energy in a bubble may decay, e.g., due to radiation, and therefore d may be negative. The thin-shell approximation, however, restricts it to d > −3/4 (also marked in Figure 1), because otherwise the total energy would be above 100%. This corresponds to the momentum-conserving snowplow. The lower limit to the kinetic energy fraction is 0.2, achieved for a strong increase of the bubble energy.
Figure 1. Kinetic energy fraction in the thin-shell approximation, k, as a function of the total energy evolution index d (compare Equation (5)). In the thin-shell approximation, k never falls below 0.2 (horizontal axis). For a wind with constant input power (d = 1), it is k = 0.3; for an energy-conserving explosion (d = 0), it is k = 0.4. Both cases are marked in the plot (dotted lines and filled circles). The lower limit for d is −3/4 (dashed line) because the kinetic energy fraction may not exceed unity.
Download figure:
Standard image High-resolution imageThis model is particularly useful for comparison to simulations, because it allows for a direct comparison of the evolution of the bubble radius with E(t) given by the stellar input to the corresponding one with E(t) measured from the simulation, which takes into account the complex heating and cooling history due to the 3D hydrodynamical evolution.
3. COMPARISON TO 3D HYDRODYNAMICS SIMULATIONS
We now compare these analytic solutions to the hydrodynamics simulations presented in Papers I and II. We show here a simulation with three massive stars (25, 32, and 60 M☉) at tens of parsecs distance from each other, where the individual bubbles merge early (3S1-hr in Paper I).
The kinetic energy fraction varies in different phases of the superbubble evolution, but is generally well described by the expectation from the thin-shell model in the respective phases (Figure 2). In particular, in between supernovae the kinetic energy fraction strongly dominates over the thermal one. We show the energy evolution for 1 Myr after each supernova for all three supernovae on a log–log plot in Figure 3. A power-law behavior with a slope of −3/4 is also indicated in Figure 3. Obviously, the energy loss of the bubble is limited by, and close to, the expectation from the thin-shell approximation. We note that the energy tracks of Thornton et al. (1998, their Figure 3, top left) for the radiative phase of an isolated supernova show the same effects as discussed here.
Figure 2. Energy evolution for a 3D hydrodynamical simulation of the interstellar medium around three massive stars. Top: time evolution of the kinetic (red dashed) and thermal (black solid) energy. Bottom: corresponding fractions of the total energy. The blue dotted lines mark interesting values for the kinetic energy fraction derived in the thin-shell approximation, namely, the lower limit (20%), the case of a constant luminosity wind (30%), and the case of the isolated, adiabatic supernova (40%).
Download figure:
Standard image High-resolution image
Figure 3. Evolution of the total energy after each supernova, i.e., for 1 Myr after the one at 4.6 Myr (solid black), at 7.0 Myr (dashed red), and at 8.6 Myr (dotted blue), respectively. The dash-dotted thick line is the maximum energy loss rate which we derive for the thin-shell approximation, and it is proportional to t−3/4.
Download figure:
Standard image High-resolution imageThe thin-shell approximation also adequately reflects the link between total energy and shell radius (Figure 4). The total energy differs from the input energy due to radiative losses. While the individual supernovae (4.6 Myr, 7 Myr, 8.6 Myr) are just discernible as mild steepenings of the curve, the overall evolution with time t is very well described by a t3/5 power law, as applicable for a bubble around an energy source with constant power. The same conclusion was reached previously by Mac Low & McCray (1988). We discuss this below in Section 5. This supports the use of the 3/5 law in analysis of observations in the literature (compare Section 1).
Figure 4. Evolution of the bubble radius. The top part compares the thin-shell approximation to the mean bubble radius from the simulation ("+" symbols) in a double logarithmic scale. We show a thin-shell model using the full input energy (solid line (A)), and one using the energy measured from the simulation (dotted line (B)). The ratio of the two different thin-shell approximations, A/B, is shown in the bottom part. The radiative energy loss reduces the bubble radius by 38% on average.
Download figure:
Standard image High-resolution imageHowever, the shell radius is now found to be substantially smaller, only 62% ± 1.6% of what it would be for an adiabatic bubble, i.e., without radiative losses (Figure 4, bottom part). This compares to a reduction by only 10% in the case of constant wind luminosity for the model of Weaver et al. (1977). The smaller radius in our superbubble models from 3D hydrodynamic simulations reflects the enhanced radiative energy losses compared to the Weaver model caused by the strong time dependence of the energy input in realistic cases.
4. CORRELATION BETWEEN SHELL ACCELERATION AND X-RAY EMISSION
We calculate the superbubble X-ray emission spectrum using a model for hot plasma emission ("Mekal"; Mewe et al. 1985, 1986; Liedahl et al. 1995). General properties of the X-ray emission are discussed in Paper II. Entrainment and mixing lead to strong X-ray emission in our simulations, close to the observed luminosity. The simulation results are converged in X-ray-bright phases. Figure 5 shows that the simulated X-ray luminosity strongly correlates with shell acceleration. The reason is that the supernova shock, which heats the intermediate density gas entrained due to 3D hydrodynamic instabilities, also accelerates the shell. As the shell accelerates, more instabilities entrain colder gas, the mixing of which combines with adiabatic expansion to reduce the X-ray luminosity. At the same time, the shell velocity starts to decline, quickly approaching the behavior in the momentum-conserving phase (compare Section 3 above). We find that the first supernova accelerates the shell beyond the velocity expected from an adiabatic Weaver model with constant, averaged energy input.
Figure 5. Evolution of shell velocity. The symbol size is chosen to be proportional to the logarithm of the X-ray luminosity in the full ROSAT band (0.1–2.4 keV); the color also encodes X-ray luminosity. The dotted line gives the shell velocity for an adiabatic comparison model, where the energy input rate is constant in time and set to the average energy-loss rate of the massive stars within 10 Myr. Vertical lines indicate the times of the supernova explosions.
Download figure:
Standard image High-resolution image
5. DISCUSSION
5.1. Caveats
In the earliest stage of the creation of an interstellar bubble, which corresponds to the main-sequence phase of the most massive star, the shell always decelerates. Thus, the contact surface between the shocked wind and ambient gas is stable and should therefore display a sharp change in density and temperature. The finite resolution of simulations such as in our code unavoidably produces intermediate temperature cells that artificially enhance radiative cooling beyond what one would expect from a thermal conduction model. Consequently, we discard this phase from further analysis, as in Paper II, i.e., up to about 4 Myr of simulation time.
From the first Wolf–Rayet phase onward, the shell accelerates significantly (compare Figure 5). Consequently, Rayleigh–Taylor instabilities disturb the contact surface and lead to entrainment of shell gas into the bubble interior. The complex interplay between the Vishniac instability, the Rayleigh–Taylor instability, and mixing then determines the evolution. The outcomes of these processes are naturally resolution dependent, as, for example, the smallest Rayleigh–Taylor mode grows fastest. Global properties may still converge with increasing resolution, and thus be robust, if dominated by resolved modes. Radiative dissipation (Paper I) and X-ray properties (Paper II) have essentially converged in this phase. Therefore, the results presented here should be robust.
Two reasons lead to enhanced energy loss in the simulated bubbles compared to the standard solution by Weaver et al. (1977). First, radiation losses are larger in a mixing layer of a few parsecs width. Densities are typically an order of magnitude higher than expected from thermal conduction, ρ(x) = ρc(1 − x)−2/5 (Mac Low & McCray 1988), where ρc is the central density and x is the radial coordinate in units of the shell radius. Second, the superbubbles are found in our simulations to almost always be close to the momentum-conserving snowplow phase, which, as shown above, dissipates energy more efficiently than the standard pressure-driven snowplow. Both processes contribute similarly to increase the overall energy loss.
5.2. Link to Observations
The momentum-conserving snowplow phases are a consequence of the time-dependent energy input, and are therefore included in, e.g., the 1D models of Oey (1996) or Oey & García-Segura (2004). The latter authors also accounted for photoionization by massive stars. While the exact run of the density differs, the effect of photoionization is quite comparable to mixing in creating an intermediate density layer inside of the shell. Jaskot et al. (2011) analyzed the X-ray properties of two superbubbles, including DEM L50, using the code of Oey & García-Segura (2004). They found too low X-ray emission, probably related to a very narrow photoionized layer, and suggested that this was due to unaccounted for thermal conduction. We would expect that the mixing layer proposed here can explain the X-ray emission at least equally well. We find a peak X-ray luminosity of 1.3 × 1036 erg s−1, whereas Jaskot et al. (2011) give 2–4.5 × 1036 erg s−1. Similar X-ray-bright superbubbles in Oey (1996) have 1.8 and 5.4 × 1035 erg s−1 (DEM L25 and DEM L301, respectively). Within uncertainties, there is thus reasonable agreement for the accelerating phase of our simulations, up to about 0.2 Myr after a supernova. Independent confirmation of recent supernova activity is difficult. 26Al, radioactive with a half-life of 0.7 Myr and ejected by supernovae, may provide an additional constraint in nearby (≲kpc) superbubbles. It is indeed detected (Diehl 2002) in the high-velocity, X-ray-bright (Oey 1996) Orion-Eridanus superbubble.
For the superbubbles DEM L25 and DEM L301, in both of which the most massive star should have had around 60 M☉, which should have exploded a fraction of a Myr ago, Oey & García-Segura (2004) find velocities of 40 and 20 km s−1, respectively, in the photoionized layer of their models. This compares to measurements of 60 and 40 km s−1 from Hα, respectively, in broad agreement, given unaccounted for details of ionization and geometry. Our simulated superbubbles have similar parameters to DEM L25 and DEM L301. At 4.5 Myr, their size is roughly 50 pc which is comparable to the observed one. Oey & García-Segura (2004) use an ambient density of 17 cm−3 to model these superbubbles, which is very comparable to the 10 cm−3 used in our simulations. Still, our simulated superbubbles reach only 8 km s−1 for the bulk shell velocity. This discrepancy might be caused by the neglect of photoionization in our models, as the shells in Oey & García-Segura (2004) are only partially ionized so that some low-velocity parts would not contribute to the Hα emission.
Mac Low & McCray (1988) showed that the shock waves of supernovae in superbubbles always become subsonic well before they reach the supershell because the thermal energy of the superbubble is larger than the energy increase due to the supernova. In our 3D simulations, superbubbles cool much more quickly and have only a fraction of a supernova energy immediately before each explosion. Consequently, the shock wave reaches the shell even if the supernova is central. Off-center supernovae, which are necessary for explaining X-ray-bright superbubbles in the self-similar framework (Chu & Mac Low 1990), are no longer required in our 3D models, but produce a similar X-ray luminosity as the central ones (Paper II).
6. CONCLUSIONS
We have analyzed the dynamical evolution of emerging superbubbles in 3D simulations, and compared it to a thin-shell model which allows for additional energy loss in the bubble interior. We found that the latter describes the global aspects of the 3D simulation results very well. The total bubble energy in this thin-shell model cannot decay faster with time t than t−3/4, as in the classical momentum-driven snowplow. In the simulations, the cooling rate after energy injection by a supernova approaches this dynamical limit. Superbubbles are therefore much closer to a momentum-conserving snowplow than to the pressure-driven solution of the Weaver model. While the bubble radius in our simulations still follows a 3/5 law, there is an enhanced energy loss due to the predominance of momentum-conserving phases and enhanced cooling in a mixing layer. The small bubble sizes that we find are in good agreement with the observational results mentioned in Section 1. Shells in our simulations accelerate beyond the velocity of a Weaver solution with equivalent, but time-averaged, power input from the stars and supernovae. This appears to be coincident with X-ray brightening, as expected from observations. Our model does not explain the high Hα velocities observed in X-ray-bright superbubbles, possibly due to a lack of radiative transfer in the code.
We thank the referee for very useful suggestions that significantly improved the manuscript. This research was supported by the cluster of excellence "Origin and Structure of the Universe" (www.universe-cluster.de).
Footnotes
- 3
We actually calculate 0.82 from the expression they give.