Evolution of supernovae-driven superbubbles with conduction and cooling

TL;DR

This study uses hydrodynamic simulations to explore how thermal conduction and cooling influence the evolution and internal structure of supernova-driven superbubbles, revealing significant deviations from classical models.

Contribution

It introduces a new analytic model incorporating cooling effects, improving predictions of superbubble evolution over classical similarity solutions.

Findings

Evaporation reduces interior temperature and increases density by over an order of magnitude.

Most conductive heat flux is lost to cooling, lowering evaporation rates compared to classical models.

Supernova blast waves remain supersonic in the interior for the first ~30 SNe.

Abstract

We use spherically symmetric hydrodynamic simulations to study the dynamical evolution and internal structure of superbubbles (SBs) driven by clustered supernovae (SNe), focusing on the effects of thermal conduction and cooling in the interface between the hot bubble interior and cooled shell. Our simulations employ an effective diffusivity to account for turbulent mixing from nonlinear instabilities that are not captured in 1D. The conductive heat flux into the shell is balanced by a combination of cooling in the interface and evaporation of shell gas into the bubble interior. This evaporation increases the density, and decreases the temperature, of the SB interior by more than an order of magnitude relative to simulations without conduction. However, most of the energy conducted into the interface is immediately lost to cooling, reducing the evaporative mass flux required to balance conduction. As a result, the evaporation rate is typically a factor of $\sim$3-30 lower than predicted by the classical similarity solution of Weaver et al. (1977), which neglects cooling. Blast waves from the first $\sim$30 SNe remain supersonic in the SB interior because reduced evaporation from the interface lowers the mass they sweep up in the hot interior. Updating the Weaver solution to include cooling, we construct a new analytic model to predict the cooling rate, evaporation rate, and temporal evolution of SBs. The cooling rate, and hence the hot gas mass, momentum, and energy delivered by SBs, is set by the ambient ISM density and the efficiency of nonlinear mixing at the bubble/shell interface.

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{4\pi }{3}{R}^3{\rho }_0\frac{\mathrm{d}R}{\mathrm{d}t}\right)=4\pi R^{2}P.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{4\pi }{3}{R}^3\frac{P}{\gamma -1}\right)={\dot{E}}_{\mathrm{in}}-4\pi R^2\frac{\mathrm{d}R}{\mathrm{d}t}P,$$

$$R\left(t\right)$$

$$P\left(t\right)$$

$$E_{\mathrm{SB}}$$

$$=\frac{1}{2}{M}_{\mathrm{shell}}{v}_{\mathrm{shell}}^2+\frac{4\pi }{3}{R}^3\frac{P}{\gamma -1}.$$

$$E_{\mathrm{kin}}=\frac{15}{77}{\dot{E}}_{\mathrm{in}}t;E_{\mathrm{th}}=\frac{5}{11}{\dot{E}}_{\mathrm{in}}t,$$

$$E_{\mathrm{SB}}\left(t\right)=\frac{50}{77}{\dot{E}}_{\mathrm{in}}t.$$

$${\dot{E}}_{\mathrm{cooling},\mathrm{shocked}\mathrm{ism}}=\frac{27}{77}{\dot{E}}_{\mathrm{in}}.$$

$$E_{\mathrm{kin}}=\frac{3}{10}{E}_{\mathrm{SB}};E_{\mathrm{th}}=\frac{7}{10}{E}_{\mathrm{SB}},$$

$$4\pi r^2(v_{\mathrm{shell}}-v)(\mathcal{E}+P)\approx \frac{5}{2}\dot{m}{c}_{s,\mathrm{iso}}^2=-4\pi r^2\kappa \frac{\mathrm{d}T}{\mathrm{d}r},$$

$$\dot{m}=\frac{-8\pi \mu m_{\mathrm{p}}{r}^2\kappa \left(\mathrm{d}T/\mathrm{d}r\right)}{5k_{B}T}.$$

$${\dot{m}}_{\mathrm{nc}}=\frac{16\pi \mu m_{\mathrm{p}}RC{T}_{\mathrm{int}}^{5/2}}{25k_B},$$

$$T_{\mathrm{hot},\mathrm{nc}}(t)$$

$$M_{\mathrm{hot},\mathrm{nc}}(t)$$

$$\propto {\rho }_0^{-2/35}{\dot{E}}_{\mathrm{in}}^{27/35}{C}^{2/7}{t}^{41/35}.$$

$$\frac{\partial \mathcal{E}}{\partial t}+\nabla \cdot \left(\left(\mathcal{E}+P\right)\mathbf{v}\right)=\nabla \cdot \left(\kappa \nabla T\right)-n\left[n\Lambda \left(T\right)-\Gamma \left(T\right)\right]$$

$${\kappa }_S=\frac{1.70\times 10^{11}{T}_7^{5/2}}{1+0.029\ln \left(T_7{n}_{\mathrm{e},-2}^{-1/2}\right)}\mathrm{erg}{\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1},$$

$${\kappa }_P=2.5\times 10^5{T}_4^{1/2}\mathrm{erg}{\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1},$$

$$q=-\kappa \nabla T.$$

$$\frac{1}{|q|}=\frac{1}{\left|-\kappa \nabla T\right|}+\frac{1}{\left(3/2\right)\rho c_{s,\mathrm{iso}}^3}.$$

$${\kappa }_{\mathrm{eff}}^{-1}={\kappa }^{-1}+\left|\nabla T\right|/((3/2)\rho c_{s,\mathrm{iso}}^3).$$

$${\kappa }_{\mathrm{mix}}=\frac{(\lambda \delta v)\times \rho k_{\mathrm{B}}}{\mu m_{\mathrm{p}}}.$$

$$\lambda \delta v$$

$$=0.1\mathrm{pc}\mathrm{km}{\mathrm{s}}^{-1}{\left(\frac{v_{\mathrm{rel}}}{10\mathrm{km}{\mathrm{s}}^{-1}}\right)}^2{\left(\frac{{\rho }_{\mathrm{high}}/{\rho }_{\mathrm{low}}}{100}\right)}^{-1}\left(\frac{\Delta t_{\mathrm{SNe}}}{0.1\mathrm{Myr}}\right).$$

$$T_{\mathrm{hot}}\equiv \frac{2}{3}\frac{\mu m_p}{k_b}\frac{E_{\mathrm{hot},\mathrm{th}}}{M_{\mathrm{hot}}}.$$

$${\hat{E}}_{\mathrm{SB}}\left(t\right)\equiv \frac{E_{\mathrm{SB}}\left(t\right)}{N_{\mathrm{SNe}}\left(t\right)}=\frac{E_{\mathrm{SB}}\left(t\right)}{\lfloor t/\Delta t_{\mathrm{SNe}}\rfloor +1},$$

$${\dot{E}}_{\mathrm{diff}}$$

$${\dot{E}}_{\mathrm{advec}}$$

$${\mathcal{L}}_{\mathrm{int}}\equiv {\int }_0^{r(T=T_{\min ,\mathrm{shell}})}n\left[n\Lambda \left(T\right)-\Gamma \left(T\right)\right]4\pi r^2\mathrm{d}r.$$

$$\dot{m}=\frac{-4\pi r^2\kappa \left(\mathrm{d}T/\mathrm{d}r\right)-{\mathcal{L}}_{\mathrm{int}}}{\left(5/2\right)c_{s,\mathrm{iso}}^2}.$$

$$\begin{gathered} R(t)=83\mathrm{pc}\times {\left(1-\theta \right)}^{1/5}{\left(\frac{E_{\mathrm{SN}}}{10^{51}\mathrm{erg}}\right)}^{1/5}{\left(\frac{\Delta t_{\mathrm{SNe}}}{0.1\mathrm{Myr}}\right)}^{-1/5}\times \\ {\left(\frac{{\rho }_0}{1.4m_p\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1/5}{\left(\frac{t}{1\mathrm{Myr}}\right)}^{3/5} \end{gathered}$$

$$\begin{gathered} {\hat{p}}_{\mathrm{rad}}(t)=4\times 10^5{M}_{\odot }\mathrm{km}{\mathrm{s}}^{-1}{\left(1-\theta \right)}^{4/5}{\left(\frac{E_{\mathrm{SN}}}{10^{51}\mathrm{erg}}\right)}^{4/5}\times \\ {\left(\frac{\Delta t_{\mathrm{SNe}}}{0.1\mathrm{Myr}}\right)}^{1/5}{\left(\frac{{\rho }_0}{1.4m_p\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{1/5}{\left(\frac{t}{1\mathrm{Myr}}\right)}^{2/5}. \end{gathered}$$

Full text

††pubyear: 2018††pagerange: Evolution of supernovae-driven superbubbles with conduction and cooling–E

1 Introduction

Correlated supernovae (SNe) sweep out large cavities in the interstellar medium (ISM) around young star clusters. These “superbubbles” (SBs) are filled with hot gas and often survive for tens of Myr, with repeated SN explosions providing energy to keep their interiors pressurized. Hot gas within SBs is observed to emit X-rays, and recombination lines are observed to originate at the interface between hot SB interiors and cooled shells of swept-up ISM (e.g. Cash et al., 1980; Dunne et al., 2001; Cecil et al., 2002; Townsley et al., 2006; Jaskot et al., 2011). Resonant lines such as O VI, usually observed in absorption, trace warm-hot gas in the interface (e.g. Savage & Lehner, 2006; Bowen et al., 2008; Li et al., 2017). The effects of SBs can also be seen in the distributions of dust and gas in the Milky Way and external galaxies, where SBs that have broken out of galactic disks leave holes in the neutral gas distribution (e.g. Heiles, 1979; Brinks & Bajaja, 1986; Puche et al., 1992; Weisz et al., 2009; Bagetakos et al., 2011) and dust filaments tracing chimneys that vent hot gas from the disk into the halo (e.g. Norman & Ikeuchi, 1989; Strickland et al., 2004). At much earlier stages, expanding SBs may be important to the destruction of the molecular clouds that host nascent super star clusters (e.g. Leroy et al., 2018).

Most stars form in clusters, so a large fraction of SNe explode inside SBs. SBs are therefore the dominant source of hot gas in galaxies and are thought to play a critical role in the launching of galactic winds (e.g. Chevalier & Clegg, 1985; Nath & Shchekinov, 2013; Sharma et al., 2014; Vasiliev et al., 2015; Kim & Ostriker, 2017, 2018; Fielding et al., 2017, 2018). The momentum injected into the ISM by SBs is a key driver of turbulence, which regulates star formation (e.g. Norman & Ferrara, 1996; Elmegreen & Scalo, 2004; McKee & Ostriker, 2007; Ostriker & Shetty, 2011; Kim et al., 2013).

When SNe explode in isolation, the energy is often deposited in high-density gas that can rapidly cool. Radiative losses thus remove $\sim$90% of the initial energy of isolated SNe (e.g. Cioffi et al., 1988; Blondin et al., 1998; Thornton et al., 1998), while at the same time producing an expanding shell with terminal momentum of order $10^{5}M_{\odot}{\rm kms^{-1}}$ (e.g. Kim & Ostriker, 2015; Iffrig & Hennebelle, 2015; Walch & Naab, 2015; Martizzi et al., 2015). In contrast, when clustered SNe explode within SBs, a larger fraction of the energy may be retained as internal hot gas, while simultaneously accelerating the surrounding warm/cold ISM. Once a bubble has been excavated by early SNe and/or stellar winds, the SNe that explode at later times deposit energy in a hot, low-density bubble interior. Continued SN explosions maintain a strong termination shock and keep SBs overpressurized long after the leading shock becomes radiative. Cooling is inefficient in the low-density hot gas in SB interiors, with most radiative losses occurring in the thin interfaces between hot and cool components of the SB.

At the outer edge of a SB, there is an abrupt transition between low-density, hot gas ($T\sim 10^{6-7}$ K; Townsley et al. 2006) and the shell of dense gas that has been swept up by the expanding shock and subsequently cooled. In this region, which is typically quite thin ($\lesssim$ 5 pc; e.g. Kregenow et al. 2006; Sankrit & Dixon 2007), thermal conduction carries energy from the hot interior into the cool shell, and cool gas evaporates into the bubble. This interface, as well as the interfaces around cold clouds embedded in the bubble, is where most cooling losses occur (e.g. Li et al., 2017), since it is there that the temperature passes through the peak of the cooling curve. Cooling is aided by hydrodynamic instabilities in the interfaces that drive mixing of hot and cool gas (Miles, 2009; Sano et al., 2012; Krause et al., 2013; Inoue et al., 2013; Fierlinger et al., 2016), and by the propagation into the shell of shocks that re-heat cooled gas.

Many theoretical studies have investigated the structure and dynamics of SBs. Castor et al. (1975) and Weaver et al. (1977) developed an analytic similarity solution to predict the expansion history and internal structure of bubbles driven by an energy source with constant mechanical luminosity. Their model was designed to describe a bubble cleared out by stellar winds rather than SNe, but subsequent studies (e.g. McCray & Kafatos, 1987; Mac Low & McCray, 1988) showed that, because the SNe rate predicted for a simple stellar population with a typical IMF is roughly constant over $\sim$40 Myr, their solution (which we hereafter refer to as the “Weaver solution”) can be straightforwardly adapted to describe SBs driven by discrete, repeated SNe. Of course, treatment of discrete SNe explosions as a continuous energy source becomes less accurate when the time interval between SNe is long. An assumption made in adapting the Weaver solution for SBs driven by SNe is that the blasts produced by individual SNe become subsonic and share their energy fully with the hot bubble interior, rather than propagating as shocks to deposit energy in the shell.

The Weaver solution has since been widely used to interpret observations of SBs in the Milky Way and nearby galaxies. Many subsequent studies have carried out one-, two-, and three-dimensional simulations of SBs driven by discrete SNe (e.g. Tomisaka et al., 1981; Tomisaka & Ikeuchi, 1986; Silich et al., 1996; Sharma et al., 2014; Gentry et al., 2017; Kim et al., 2017; Yadav et al., 2017; Vasiliev et al., 2017; Fielding et al., 2018; Gentry et al., 2019). These studies have broadly validated some of the predictions of the Weaver solution when its assumptions apply, particularly the predicted scaling relations of the SB expansion rate, energy, and momentum with the SN rate, ambient density, and SB age. They have also highlighted some effects that are not captured in the spherical similarity solution, such as the collision of shocks from discrete SNe with the cooled shell, the development of dynamical instabilities at the bubble/shell interface, the departure from spherical symmetry arising from interaction with the inhomogeneous ISM, the ablation over time of cold clouds left behind within the bubble, and the vertical extension and breakout of SBs from galactic disks. Most notably, the Weaver solution does not account for cooling losses in the shell/bubble interface, so simulations with significant radiative cooling generally produce less energetic SBs than predicted by the analytic model.

Most previous numerical studies have not explicitly modeled thermal conduction. In the Weaver solution, conduction has no effect on the bubble energetics and expansion rate – it only changes the structure of the SB interior, increasing its density and decreasing its temperature. Some numerical studies have used the evaporative mass flux into the bubble interior that is predicted by the Weaver solution as a “subgrid” model for effects of thermal conduction (e.g. Mac Low & McCray, 1988; Silich et al., 1996; Keller et al., 2014). A few recent studies have tested the effects of explicitly including conduction as a source term in the energy equation in some of their simulations (Sharma et al., 2014; Gupta et al., 2016; Yadav et al., 2017), but the effects of conduction on the shell/bubble interface and on the interior structure of SBs have yet to be examined in detail. Moreover, the prediction of Weaver et al. (1977) for the evaporative mass flux from the shell into the bubble – which is based on the assumption that the energy fluxes from conduction and evaporation exactly cancel – has not been tested.

In this paper, we present results from idealized, one-dimensional simulations of SBs in a uniform medium for a wide range of SB parameters. We extend previous studies by explicitly solving for both thermal conduction and cooling, treating turbulent mixing at the shell/bubble interface via a parameterized term in the energy equation. We study in detail how the evaporative mass flux due to conduction varies with model parameters such as thermal conductivity, ambient density, SN rate, and nonlinear mixing efficiency. We also identify transitions in the SB character as the bubble ages, demonstrating that the continuous energy injection approximation is only valid at late times (which we quantify). We develop a modified analytic solution for the SB expansion rate and interior properties that takes into account both conduction and cooling in the shell, with the latter induced by small-scale mixing at the hot/cool interface. This analytic solution including cooling can be considered a successor to the no-cooling Weaver solution. As we shall show, cooling not only slows SB expansion and reduces momentum injection, but also alters the structure of the conductive bubble/shell interface and reduces the evaporation rate. Throughout the paper, we compare our simulations and analytic models to the classical Weaver solution.

The remainder of this paper is organized as follows. In Section 2, we review the classical theory of SB dynamics and the effects of conduction on the bubble interior in the absence of cooling. We summarize the numerical methods and model parameters we use for our SB simulations in Section 3. Section 4 presents and analyzes our fiducial simulations, focusing on a single choice of SN rate and ambient density. Here we present radial profiles (Section 4.1), time-evolution of integrated SB properties (Section 4.2), and the energy balance of fluxes through the shell/bubble interface (Section 4.3). Section 5 develops an analytic model to predict the internal structure and evolution of integrated properties of SBs with cooling. Here we compare the model predictions to a suite of simulations with a range of ambient densities, SNe rates, and nonlinear mixing efficiencies. Section 6.1 considers the evolution of SBs at early times, when the effects of discrete SNe cannot be accurately modeled under the approximation of continuous energy injection. Here we also compare to the early-time dynamical evolution that would be driven by stellar winds. We summarize our results, compare to previous work, and discuss astronomical implications in Section 7. The Appendices provide further details about the numerical models used in the simulations.

2 Classical Theory

2.1 Dynamical evolution

In this section, we summarize the classical theory of spherical SB evolution driven by continuous energy injection.

After a short period of adiabatic expansion that is analogous to the Sedov solution for a single SN, high-density shocked gas at the leading edge of a SB cools radiatively, forming a thin, dense shell that is separated from the hot bubble interior by a contact discontinuity. The timescale for shell formation is set by the cooling time in the shocked gas, which is generally short ($\ll 0.1$ Myr) compared to the lifetime of the SB (e.g. Kim et al., 2017).

After shell formation, the structural properties of a hot, pressurized bubble expanding into a homogeneous ISM are set by energy and force balance. Weaver et al. (1977) derived a similarity solution for the evolution of such a bubble under the assumption that cooling losses within the bubble and inner interface are negligible (although the leading shock and outer layer of the shell are assumed to be radiative). Their solution assumes a continuous energy injection luminosity from stellar winds, but it can be adapted to describe a bubble driven by a series of SN explosions with fixed time-averaged power (e.g. McCray & Kafatos, 1987; Mac Low & McCray, 1988; Koo & McKee, 1992; Kim et al., 2017).

Neglecting gravity,111For both the classical SB solution and the models developed in this work, the global force of gravity in countering shell expansion is negligible compared to pressure over the full range of SB parameters we consider. the expansion of the pressure-driven spherical SB is described by a momentum equation:

[TABLE]

Here $R$ and $P$ represent the radius and internal pressure of the SB, and $\rho_{0}$ is the density of the ambient medium. It is assumed that (a) the SB radial momentum is dominated by the cooled shell, $p_{{\rm rad}}\approx M_{{\rm shell}}v_{{\rm shell}}$, where $v_{{\rm shell}}={\rm d}R/{\rm d}t$ is the shell velocity, and the shell contains most of the mass that has been swept up, so $M_{{\rm shell}}\approx\left(4\pi/3\right)R^{3}\rho_{0}$; (b) the external pressure is negligible; and (c) the SB interior is isobaric due to the short sound crossing time within the bubble.

Assuming that the internal energy of hot gas in the SB interior is reduced by adiabatic expansion with negligible cooling losses, the internal energy of the SB follows

[TABLE]

where $\dot{E}_{\rm in}=E_{\rm SN}/\Delta t_{\rm SNe}$ is the time-averaged energy input rate from SNe, $E_{{\rm SN}}\approx 10^{51}\,{\rm erg}$ is the mechanical energy injected by one SN, $\Delta t_{\rm SNe}$ is the time interval between SNe explosions, and $\gamma=5/3$ is the ratio of specific heats.

Equations (1) and (2) can be directly solved for $R(t)$ and $P(t)$, yielding

[TABLE]

This solution is only valid up to the time when the pressure in Equation (4) becomes comparable to that of the ambient ISM.

The total energy of the SB is the sum of the thermal energy of the hot interior and the kinetic energy of the dense shell:

[TABLE]

Substituting Equations (3) and (4) into Equation (5) yields

[TABLE]

so that

[TABLE]

That is, $50/77\approx 65\%$ of the energy input from SNe is retained in the SB, and $27/77\approx 35\%$ is lost. The energy loss represents the portion of the $P\,{\rm d}V$ work done by the expanding bubble that goes into heating the ISM gas at the leading shock where it is swept up by the expanding shell. While the kinetic energy of the shocked ISM gas is retained, it is assumed that the post-shock thermal energy is quickly radiated away in the exterior layer of the shell:

[TABLE]

The relative proportions of thermal to kinetic energy in the classical SB solution are

[TABLE]

which is very similar to the energy proportions for the Sedov solution for a single non-radiative supernova blast ($E_{{\rm kin}}=0.28E_{{\rm SN}};\,E_{{\rm th}}=0.72E_{{\rm SN}}$).

To summarize, the classical Weaver solution assumes that there is negligible cooling in the bubble interior and shell/bubble interface. As a result, the dynamics and structural properties of the SB (i.e., radius, pressure, momentum, total energy, and the ratio of thermal to kinetic energy) are simply set by force balance and energy conservation, after accounting for radiative losses in the leading shock. Thermal conduction has no effect on these properties; at a given pressure, conduction controls the density and temperature of the bubble interior (Section 2.2).

In Sections 5.1 and 5.2 we provide a re-analysis of the predicted dynamical evolution allowing for cooling.

2.2 Bubble interior

Weaver et al. (1977) also proposed a similarity solution to describe the interior structure of a wind-driven bubble. In their solution, radial profiles of temperature, density, and velocity are determined by simultaneously solving the energy conservation and continuity equations within the bubble, under the assumptions that cooling losses in the bubble and interface are negligible and the bubble expansion is described by Equations (3)-(4). In this model, the conductive energy flux from the bubble into the shell is exactly balanced by the evaporation of cool gas from the shell into the bubble, such that the shell does not gain or lose thermal energy.

In fact, even without making assumptions regarding the detailed interior structure of the SB, the evaporation rate of the shell into the bubble can be estimated using the approach adopted by Cowie & McKee (1977) to compute the evaporation rate of a cold cloud embedded in a hot medium. The total conductive heat flux through a surface at radius $r$ is $-4\pi r^{2}\kappa\left({\rm d}T/{\rm d}r\right)$, where $\kappa$ and ${\rm d}T/{\rm d}r$ are the local thermal conductivity and temperature gradient, respectively. The inward mass flux in the frame of the expanding shell is $\dot{m}=4\pi r^{2}\left(v_{\rm shell}-v\right)\rho$, where $\rho$ and $v$ are the local density and flow velocity at the radius $r$ through which $\dot{m}$ is measured. In the absence of cooling, energy balance requires

[TABLE]

where ${\cal E}$ is the energy density, which can be approximated as $P/(\gamma-1)$ under the assumption that thermal energy dominates in the SB interior. Here we have defined the isothermal sound speed $c_{s,{\rm\,iso}}=\sqrt{k_{\rm B}T/\mu m_{\rm p}}$, where $\mu$ is the mean molecular weight, $k_{\rm B}$ is Boltzmann’s constant, and $m_{\rm p}$ is the proton mass. The mass flux required to balance the conductive heat flux is then

[TABLE]

The mass flux profile $\dot{m}(r,t)$ is thus fully specified when the temperature profile and thermal conductivity are known. Taking $\kappa=CT^{5/2}$ for the isotropic thermal conductivity (Spitzer 1962; “$C$” is a constant) and considering the “no cooling” limit, a local expansion of the temperature profile near the interface (e.g. Cowie & McKee, 1977) can be used to show that $r{\rm d}T/{\rm d}r\to\left(2/5\right)T_{\rm int}$, where $T_{{\rm int}}$ is the asymptotic interior temperature inward from the conductive interface.

The evaporative mass flux at the shell/bubble interface in the no-cooling limit thus becomes

[TABLE]

where $R$ is the SB radius. This expression is identical to that derived by Cowie & McKee (1977) for the evaporation rate of a cold cloud in a hot medium.

Applying Equation (12) to compute the evaporation rate into the SB interior requires knowledge of $T_{\rm int}$. A simple estimate can be obtained by assuming a uniform interior with density $\rho=P/(kT_{\rm int}/\mu m_{\rm p})=\left(\int\dot{m}\,{\rm d}t\right)/\left(4\pi R^{3}/3\right)$ and combining Equation (12) with Equations (3) and (4). We denote the resulting value as $T_{\rm int}=T_{\rm hot,nc}$ to make clear that it depends on the no-cooling assumption:

[TABLE]

Substituting into Equation (12) and integrating yields

[TABLE]

Of particular note is the weak dependence of $T_{\rm hot,\,nc}$ on time, the input energy rate, and the ambient density. With a constant energy input rate, $\dot{M}_{\rm hot,nc}$ has the same weak scalings as $T_{\rm hot,\,nc}^{-1}$.

Weaver et al. (1977) did not assume a constant-density interior, but instead calculated the interior structure by numerically integrating the equations of continuity and energy conservation in the SB interior under the assumption of self-similarity. This yielded predictions similar to Equations (13) and (14), except that our temperature and mass normalizations are respectively 22% and 46% lower.

Several previous studies (e.g. Mac Low & McCray, 1988; Silich et al., 1996; Keller et al., 2014) have used Equation (12) as a subgrid model to approximate the effects of thermal conduction in simulations of SBs. We note, however, that Equations (11) and (12) are not expected to hold if there is significant cooling in the shell/bubble interface. Cooling will remove some of the energy conducted into the shell, reducing the amount of evaporation needed to balance conduction. We provide an analysis of the effects of cooling on the evaporation rate and temperature in Section 5.3.

3 Numerical simulations

We simulate a 1D domain in spherical coordinates using the Athena++ code222https://princetonuniversity.github.io/athena/ (Stone et al., 2019), a new version of Athena (Stone et al., 2008; Stone & Gardiner, 2010). We use reflecting and outflow boundary conditions, respectively, at the inner and outer boundaries. The background medium is characterized by a mean density $\rho_{0}$ and pressure $P_{0}$; we adopt $P_{0}=3000\,k_{B}\,{\rm K\,cm^{-3}}\times(\rho_{0}/(\mu_{\rm H}\,m_{{\rm p}}{\rm cm^{-3}}))$ as a fiducial value. We assume a mean molecular weight $\mu=0.62$, appropriate for ionized gas,333Note that the mean molecular weight $\mu$ should in principle vary throughout the simulation domain, depending on the ionization state of the gas. For simplicity, we fix $\mu=0.62$. where $n=\rho/(\mu m_{\rm p})$, and $\mu_{\rm H}=1.4$, with $n_{\rm H}=\rho/(\mu_{\rm H}m_{\rm p})$.

We solve the equations of hydrodynamics on a 1D grid with uniform linear cell resolution $\Delta r$. Thermal conduction contributes a flux term to the energy equation:

[TABLE]

where $\kappa$ is thermal conductivity and $\cal{E}$ is the energy density. The gas temperature is $T=P/(nk_{B})$, and $\Lambda\equiv(\mu/\mu_{\rm H})^{2}\Lambda_{\rm H}$ and $\Gamma\equiv(\mu/\mu_{\rm H})\Gamma_{\rm H}$ are cooling and heating rate coefficients. Subscripts of “H” denote the coefficients with respect to the hydrogen number density.

We use an operator-split implementation of cooling and heating adapted from Fielding et al. (2018). The cooling curve $\Lambda_{\rm H}(T)$ is a 10-segment piecewise power law fit to the cooling curve from Oppenheimer & Schaye (2013) for $T>10^{4}\rm K$ and the one from Koyama & Inutsuka (2002) at lower temperatures; see Appendix E for details. We adopt an absolute temperature floor of 100 K. We modify the cooling curve from Fielding et al. (2018) slightly to drop off steeply over $100<T/{\rm K}<500$. This modification prevents the temperature in the shell from reaching the floor, which can lead to numerical problems, but it has little effect on the dynamics of the SB.

We use a constant, spatially uniform photoelectric heating rate of $\Gamma_{\rm H}=n_{{\rm H,0}}\Lambda\left(P_{0}\mu m_{p}/(k_{B}\rho_{0})\right)$ for $T<10^{4}\,\rm K$; here $n_{{\rm H,0}}$ is the hydrogen number density of the ambient medium. The scaling of $\Gamma_{\rm H}$ with $n_{\rm H,0}$ crudely accounts for the increased photoelectric heating rate when $n_{\rm H,0}$ and the star formation rate are higher. The practical consequence of adopting this heating rate is that heating and cooling exactly balance in the ambient medium. At $T>10^{4}$ K, we take $\Gamma_{\rm H}=0$, since we expect dust grains to be destroyed at higher temperatures.

3.1 Thermal conduction

Conduction in an ionized gas is primarily mediated by electrons. The thermal conductivity for a hydrogen plasma is (Spitzer, 1962):

[TABLE]

where $T_{7}\equiv T/(10^{7}\,\rm K)$ and $n_{\rm e,-2}\equiv n_{\rm e}/(10^{-2}\,\rm cm^{-3})$, with $n_{\rm e}$ the number density of free electrons. We take $n_{\rm e}=1.2n_{\rm H}$, appropriate for an ionized gas with solar metallicity. Equation (16) includes a correction for the self-consistent electric field produced by the diffusing electrons.

Because the electron density only enters Equation (16) in the logarithm, the density-dependence of Spitzer conductivity is very weak. To a good approximation, $\kappa$ can thus be modeled as being a function of temperature only, $\kappa_{S}\approx C\times(T/{\rm K})^{5/2}$, with $C=6\times 10^{-7}\,{\rm erg\,s^{-1}cm^{-1}K^{-7/2}}$. Our simulations, however, use the full expression from Equation (16).

Equation (16) does not hold at low temperatures, where the gas becomes neutral. At low temperatures, we take the conductivity for neutral atomic collisions from Parker (1953):

[TABLE]

where $T_{4}\equiv T/(10^{4}\,\rm K)$. The Spitzer and Parker conductivities are equal at $T=6.6\times 10^{4}$ K, so we transition between the two forms at this temperature. The effects of including Parker conductivity at low temperatures are very small, since values of $\kappa_{P}$ in the shell are orders of magnitude lower than $\kappa_{S}$ in the bubble and interface.

In the classical limit in which many collisions occur over the scale length of the temperature gradient, the conductive heat flux is

[TABLE]

However, Equation (18) overestimates the true heat flux when the electron mean free path is large compared to the temperature scale height, $T/|\nabla T|$. When the temperature gradient is very steep, the diffusive energy flux predicted by Equation (18) exceeds the energy flux that can actually be carried by electrons, which is $q_{{\rm max}}\sim(3/2)\rho c_{s,\,{\rm iso}}^{3}$ (Parker, 1963; Cowie & McKee, 1977).444The saturated heat flux is $q_{{\rm max}}\sim n_{e}k_{{\rm B}}T_{e}v_{e}$, where $T_{e}$ and $v_{e}$ represent the temperature and thermal speed of the electrons. The value of $v_{e}$ depends on the electric field set up by the electron current and distribution function of the electrons. Cowie & McKee (1977) calculated $q_{{\rm sat}}\approx 5\phi\rho c_{s,\,{\rm iso}}^{3}$, where $\phi$ is an uncertain dimensionless constant of order unity. Laser fusion experiments suggest $\phi\approx 0.3$ (Max et al., 1980; Balbus & McKee, 1982), leading to our adopted $q_{{\rm sat}}\approx\frac{3}{2}\rho c_{s,\,{\rm iso}}^{3}$. Measurements of the heat flux in the solar wind suggest slightly larger values (e.g. Bale et al., 2013). Following Balbus & McKee (1982), we use a conductive flux that smoothly interpolates between the classical and saturated values:

[TABLE]

This expression recovers Equation (18) in the limit where $\left|-\kappa\nabla T\right|\ll\left(3/2\right)\rho c_{s,\,{\rm iso}}^{3}$, and it asymptotes to $q=(3/2)\rho c_{s,\,{\rm iso}}^{3}$ in the opposite limit. Saturation is implemented via an effective conductivity, $\kappa_{\rm eff}$, that includes the effects of saturation:

[TABLE]

The actual conductive heat flux is then $q=-\kappa_{\rm eff}\nabla T$.

We examine the effects of saturation of conduction in detail in Appendix D. Accounting for saturation slightly decreases the average conductive energy flux from the bubble into the shell, and thus also decreases the evaporative mass flux into the bubble. At early times, the instantaneous heat flux immediately after the arrival of a shock from a new SNe at the shell is more than an order of magnitude larger when saturation is not included than when it is. However, in a time-averaged sense, the effects of saturation on integrated bubble properties are relatively minor, because at late times the magnitude of the classical conductive flux, $|\kappa\nabla T|$, is comparable to $(3/2)\rho c_{s,\,{\rm iso}}^{3}$. For example, including saturation decreases the hot gas mass produced per SN at late times by $\sim 15\%$ for our default SB parameters (Figure 18).

Our simulations become prohibitively expensive when $\kappa$ is very large. For computational tractability, we impose a ceiling on $\kappa$ (see Appendix C). We ensure that the ceiling is sufficiently high that it has negligible effects on our results.

3.2 Mixing due to nonlinear instabilities

In 3D simulations (and in real SBs and SN remnants), the thickness of the shell/bubble interface in which cooling occurs is enhanced by nonlinear development of Kelvin-Helmholtz (KH), Rayleigh-Taylor (RT), Richtmyer-Meshkov (RM) and related instabilities (e.g. Vishniac & Ryu, 1989; Blondin et al., 1998; Bucciantini et al., 2004; Folini & Walder, 2006; Ntormousi et al., 2011; Michaut et al., 2012; Sano et al., 2012; Pittard, 2013; Badjin et al., 2016). Such instabilities and the resulting turbulent mixing do not occur in 1D simulations, leading the interface region to become extremely thin and diminishing cooling losses (see Appendix A and Gentry et al. 2017), but their effects can be approximated through subgrid models (e.g. Duffell, 2016).

To model consequences of turbulence that are not directly captured in one-dimensional simulations, we include an effective diffusivity in the energy equation to represent nonlinear mixing. Because we wish to represent the mixing of material of different temperature phases due to instability, we implement the model through the conduction operator with a temperature-independent conductivity, $\kappa_{\rm mix}$. The actual conductivity, $\kappa$, in each simulation zone is set to the larger of $\kappa_{\rm mix}$ and thermal conductivity ($\kappa_{S}$ or $\kappa_{P}$, with appropriate modifications for saturation of conduction) described in Section 3.1.

The energy diffusion associated with mixing may be characterized by a length scale and effective nonlinear velocity amplitude. We adopt a parameterization

[TABLE]

We treat the product $(\lambda\delta v)$ as an adjustable, spatially uniform constant that parameterizes the efficiency of mixing. The value of $(\lambda\delta v)$ will be different for each unstable mode, but we expect mixing to be dominated by fluctuations of the largest-scale modes.

We estimate a possible range of values of $(\lambda\delta v)$ as follows. Instabilities driven by shear, such as the KH and RM instabilities, have a linear growth rate $\omega\sim kv_{{\rm rel}}(\rho_{{\rm high}}/\rho_{\rm low})^{-1/2}$, where $\rho_{\rm high}$ and $\rho_{\rm low}$ are the densities of the high- and low-density media (in this case, those of the shell and the bubble interior), $v_{\rm rel}$ is the tangential velocity difference across the corrugated shell-bubble interface, and $k=2\pi/\lambda$ is the wavenumber. The most effective mixing is likely to be caused by modes with growth rate comparable to the SNe rate, $1/\Delta t_{\rm SNe}$, since modes that grow more slowly will be suppressed by the arrival of blast waves from subsequent supernovae, and modes that grow more quickly are associated with smaller spatial scales and a thinner nonlinear mixing layer. Setting the growth rate equal to $1/\Delta t_{\rm SNe}$ yields a wavelength $\lambda\sim\Delta t_{\rm SNe}v_{\rm rel}(\rho_{{\rm high}}/\rho_{\rm low})^{-1/2}$. We can multiply by the typical nonlinear velocity difference for shear-driven instabilities, $\delta v\sim\omega\lambda\sim v_{{\rm rel}}(\rho_{{\rm high}}/\rho_{{\rm low}})^{-1/2}$, yielding

[TABLE]

Note that in the above, we focus on parametric scalings and omit factors of $2\pi$. For SBs expanding in the multiphase ISM, corrugation and tangential velocities may develop as a result of the shock overrunning cold, dense clouds. With hot gas flowing in to the interface at velocities of several hundred $\rm km\,s^{-1}$, depending on the level of corrugation in the front, the tangential velocity difference could be $10\lesssim v_{\rm rel}/(\rm km\,s^{-1})\lesssim 300$. Our simulations with conduction typically have $\rho_{\rm high}/\rho_{\rm low}\sim 100$ at the bubble/shell interface, so for $\Delta t_{\rm SNe}=0.1$ Myr, a plausible range of $\lambda\delta v$ is $0.1\lesssim\lambda\delta v/({\rm pc\,km\,s^{-1}})\lesssim 10$. We also note that the value of $\lambda\delta v$ could in principle change over time, since at late times (a) individual SNe blast waves no longer arrive at the interface as discrete shocks, (b) the dominant instability may change, and (c) a significant fraction of the mixing and cooling may occur at the edges of cool clouds within the SB interior, not at the shell-interior interface.

We emphasize that the above is intended solely to provide a sense of the plausible level and parametric scalings for the coefficient in mixing-driving energy exchange, and should not be taken as more than a conceptual guide. The “correct” value of $\lambda\delta v$, as well as any parameter dependencies and time dependencies, is presently quite uncertain. Although instabilities do aid mixing in 3D simulations of SBs, they are likely only marginally resolved in existing global simulations, and/or may be seeded by grid-scale numerical noise rather than density and velocity perturbations that are representative of the real turbulent, multiphase ISM. To constrain the mixing-driven diffusivity, focused high-resolution three-dimensional simulations of the interface that allow for a realistic inhomogeneous ambient medium to seed structure are essential. As we shall show, the net cooling rate is related to $\lambda\delta v$, so a possible path forward in future work is to use the cooling rates from 3D simulations as a calibration point and adjust $\lambda\delta v$ to match their energetics. We will explore such calibration further in future work.

Given the uncertainty in the mixing rate, we compare simulations with a range of values of $\lambda\delta v$, representing different efficiencies of nonlinear mixing, throughout this paper. We find that the cooling rate in the simulations is not fully converged at the resolution we can achieve for $\lambda\delta v<1\,\rm pc\,km\,s^{-1}$, so in most cases we vary $\lambda\delta v$ between 1 and 10 pc km s*-1*. For these values of $\lambda\delta v$, $\kappa_{S}\gg\kappa_{{\rm mix}}$ in the bubble interior, and $\kappa_{{\rm mix}}\gg\kappa_{S}$ in the cool shell. $\kappa_{{\rm mix}}$ begins to dominate over $\kappa_{S}$ in the region of the interface where the temperature falls below $\approx 2\times 10^{5}\,{\rm K}$ and $n_{{\rm H}}$ exceeds $\approx 0.2\,{\rm cm^{-3}.}$ We show results for a wider range of $\lambda\delta v$ in Appendix B and compare the numerical convergence of simulations with and without mixing in Appendix A.

3.3 Energy injection

In our simulations, SNe explode with a uniform time spacing $\Delta t_{\rm SNe}$. When a SN occurs, we inject $E_{\rm SN}=10^{51}$ erg of kinetic energy near the inner boundary of the simulation domain. This energy is added to an ejecta mass of $M_{\rm ej}=10\,M_{\odot}$, approximately representing the IMF-averaged total mass of SNe ejecta and swept-up circumstellar material. Within the injection region, we use a linear ejecta radial velocity profile, such that the radial velocity is 0 in the innermost cell and increases linearly to the outer edge of the injection region. For a given injection region size, the normalization of the velocity profile is set such that the total kinetic energy added per SN is $E_{\rm SN}$. Our fiducial simulations only include energy injection from SNe, which are expected to dominate the feedback energy budget over the lifetime of a young cluster. We consider the effects of stellar winds at early times in Section 6.2.

We use a default injection region of $2<r/{\rm pc}<4$, with a reflecting inner boundary at $r=2$ pc. We experimented with a range of injection region sizes, radial velocity profiles, and ejecta temperatures, and with injecting thermal rather than kinetic energy. In all cases, we find that although the properties of the bubble interior near the center of the SB are sensitive to the injection scheme, the structure of the SB near the shell/bubble interface and the integrated properties of the SB (i.e. radius, momentum, energy, and hot gas mass) are not.

3.4 Resolution

We discuss convergence with resolution in Appendix A but briefly summarize the results here. For 1D simulations without conduction or mixing, convergence is difficult to achieve because cooling losses decrease with increasing resolution. In this case, most of the cooling losses within the SB occur in the thin interface between the cool shell and hot interior, where the temperature falls from $T\sim 10^{6}\,{\rm K}$ to $T\sim 10^{3}\,{\rm K}$, passing through the peak of the cooling curve at $T\sim 10^{5}\rm K$. Without some diffusive process in the problem, at late times nothing in the problem sets a physical scale for the interface.555At early times, shocks from individual SNe propagate into the shell; this is discussed in Section 6.1. As a result, the region in which significant late-time cooling occurs becomes narrower and narrower with increasing resolution, such that cooling losses become negligible at very high resolution. In the real ISM, the transition from high to low temperature occurs over a macroscopic scale due to a combination of thermal conduction and physical mixing. To represent the latter in 1D simulations, we apply an effective diffusion model. As a result of including explicit terms for conduction and mixing, our simulations converge at relatively low resolution.

Because the simulations with conduction are computationally more expensive (see Appendix C), we compare simulations with and without conduction at different resolution levels. The fiducial simulations we show with conduction and mixing have $\Delta r=0.024$ pc and are well converged. Those without conduction have a higher resolution of $\Delta r=0.0016$ pc. Even at this resolution, the no-conduction case is not completely converged, but we show in Appendix A that the integrated bubble quantities in the no-conduction simulations asymptotically approach the predictions of the classical SB theory without cooling (Section 2). Except where explicitly noted, all the simulations are converged at the resolution we present.

3.5 Definitions

We use the following definitions for different components of the SB in analyzing the simulations. We define the “bubble” component to include all gas with $T>10^{5}$ K or $|v_{r}|>1\,\rm km\,s^{-1}$. This effectively selects all gas that has been perturbed; we define the ambient medium to include all gas that is not part of the bubble. We define the “hot” component of the bubble as all cells with $T>10^{5}$ K, and the “shell” component as all gas with $|v_{r}|>1\,\rm km\,s^{-1}$ and $T<2\times 10^{4}$ K.

For each component, we calculate the integrated mass, radial momentum, and kinetic, thermal, and total energy by summing over the relevant densities for all cells within that component, multiplied by each cell’s volume. We define the average temperature of gas in the hot component in terms of its total thermal energy:

[TABLE]

We define $R_{\rm bubble}$ as the outer edge of the outermost cell in the bubble component and $R_{\rm shell}$ as the mass-weighted mean radius of cells in the shell component. The expansion velocity of the shell is $v_{{\rm shell}}={\rm d}R_{{\rm shell}}/{\rm d}t$.

We often normalize the SB energy, radial momentum, and hot gas mass by the number of SNe that have exploded prior to a given time. This is denoted with a caret; e.g.,

[TABLE]

where $\left\lfloor x\right\rfloor$ denotes the floor function of $x$.

4 Fiducial Simulation Results

We present our simulation results in three parts.

First, in this section, we examine in detail the dynamics and internal structure of a single “fiducial” simulation: one with ambient density $n_{\rm H,0}=1\,\rm cm^{-3}$ (comparable to the mean ISM density in the solar neighborhood) and a SNe interval $\Delta t_{\rm SNe}=0.1\,\rm Myr$. Assuming a SNe rate of 1 SN per 100 M*⊙* of stellar mass formed, and a constant average SNe rate over $\sim$ 40 Myr, similar to the predictions of stellar population synthesis models for a Kroupa (2001) IMF (e.g. Leitherer et al., 1999), the fiducial SNe rate of $\Delta t_{\rm SNe}=0.1\,\rm Myr$ corresponds to a relatively large star cluster with $M_{\rm star}\sim 4\times 10^{4}\,M_{\odot}$. For our fiducial model we adopt $\lambda\delta v=1\,\rm pc\,km\,s^{-1}$ when mixing is turned on.

Then, in Sections 5 and 6.1, we present results for simulations with a range of SNe rate, ambient density, and nonlinear mixing efficiency.

4.1 Radial profiles

Figure 1 shows radial profiles of density, temperature, pressure, and radial velocity for simulations with our fiducial problem setup at $t=1$ Myr after SNe commence. The fiducial version (black) includes both conduction and nonlinear mixing. We also compare to a simulation without conduction or mixing (blue), and to one with only mixing but no conduction (gold).666We do not show a simulation with conduction and no mixing, because in this case the cooling rate is not converged at the resolution level our simulations can currently achieve, and the asymptotic behavior at high resolution is uncertain (see Appendix A).

In the simulation without conduction or mixing, there is a large degree of substructure in the radial profiles of all quantities. Such substructure is absent in simulations of a single SN explosion; it arises from the interference of the forward and reverse shocks generated by different SNe within the SB. Even at late times, contact discontinuities caused by discrete SNe can be seen in the density and temperature profiles. These structures would very likely not survive in the interior of real superbubbles, as they would be smoothed out by both instabilities and conduction.

We note that in the model without conduction but with mixing turned on, the sharpest density and temperature spikes are smoothed out. While some turbulent mixing is expected to occur in SB interiors, in the real case the parameters characterizing mixing would presumably differ from those at the outer reaches of SBs, unlike our simple prescription. Even though the detailed interior structure is unrealistic for both this case and the no-mixing, no-conduction case, the important physical point is that the mean interior density and temperature in these cases are similar.

For all quantities, conduction smooths out much of the structure in the bubble interior, erasing the memory of individual SNe. The density and temperature in the bubble interior are significantly higher and lower, respectively, in the simulation with conduction than in those without it, because conduction precipitates an evaporative mass flux from the shell into the bubble, diluting the hot gas in the interior.

The simulation that includes nonlinear mixing but not thermal conduction has a hot, low-density interior similar to the one without conduction or mixing. However, both it and the simulation with conduction and mixing produce a smaller SB at fixed time than the simulation with no conduction or mixing. As we shall show, the size and energy of the SB are set primarily by mixing, which determines the efficiency of cooling losses. Once the cooling rate is set, conduction sets the density and temperature of the SB interior. We note that the internal temperature in the interior is insensitive to cooling in the shell, however.

Figure 1 includes (in red) the profiles predicted by the Weaver similarity solution. The detailed internal structure of all simulations is different from the Weaver prediction, but much of this disagreement is a result of the discreteness of SNe explosions in the simulations; particularly at early times, the radial profiles in the bubble interior fluctuate significantly on short timescales whenever a new SN explodes.

As expected, the Weaver solution predicts a higher-density and lower-temperature interior than is produced by the simulations without conduction. The temperature in the fiducial simulation with conduction and mixing is very close to the Weaver solution, but the density is lower. As we will discuss in Section 4.3, this is primarily a consequence of cooling losses in the shell/bubble interface, which decrease the evaporative mass flux into the bubble. The Weaver solution does not account for cooling within the bubble interior or bubble/shell interface, so the bubble size in the simulation without conduction or mixing, which also has very small cooling losses, is similar to that predicted by the Weaver solution.

4.2 Integrated SB properties

We investigate the temporal evolution of the fiducial SB model in Figure 2, which compares bubble-integrated properties for the same set of simulations and analytic Weaver model shown in Figure 1.

The total bubble energy and radial momentum are a factor of $\sim$2 lower in the simulations with mixing than in the one without it. Comparing the two simulations that include mixing, the radius, total energy, and momentum are essentially identical; i.e., once mixing is included, thermal conduction has a negligible effect on the cooling efficiency and dynamics of the SB. However, the difference between the simulations with and without conduction becomes clear when comparing the mass and temperature of hot gas in the bubble interior. In the simulations without conduction, $\lesssim 10M_{\odot}$ of hot gas are added to the bubble per SN; this can be explained entirely by the ejecta mass that is added at each SN injection. With conduction, $M_{\rm hot}$ increases by $\sim 150\,M_{\odot}$ per SN at late times, primarily as a result of evaporation into the bubble from the shell.

The dynamics of the simulation without conduction or mixing (blue) are in good agreement with the predictions of the Weaver similarity solution. The internal energy, size evolution, and ratio of thermal to kinetic energy in the Weaver solution are calculated without any consideration of the effects of conduction on the bubble interior: they follow from the assumptions of force balance and adiabatic expansion (Section 2.1). At high resolution, cooling losses in the interface become very small in the simulation without conduction or mixing (Appendix A), leading to similar SB energetics to what is predicted by the Weaver solution. Indeed, in the limit where cooling within the bubble and interface vanishes (or in the limit of an infinitesimally thin interface; see Gentry et al. 2017), the energy, size, and momentum evolution should reproduce the Weaver solution almost exactly, modulo any modifications that result from the discreteness of SNe in the limit of large $\Delta t_{\rm SNe}$.

In terms of $M_{\rm hot}$ and $T_{\rm hot}$, the results from the simulation with conduction and mixing are in better agreement with the Weaver solution than the simulations without conduction, but still have slightly higher temperature and a factor $\sim$ 4 less hot gas mass than the Weaver solution predicts. As we will show in Section 4.3, both deviations are the result of cooling in the interface.

The ratio of kinetic to thermal energy at late times is similar in all simulations and is comparable to the value $E_{\rm kin}/E_{\rm SB}=3/10$ predicted by the Weaver solution. For a pressure-driven shell, this quantity depends only on the bubble’s size evolution $R(t)$; the classical power law of $R\left(t\right)\sim t^{3/5}$ that is predicted for a constant energy input rate leads to $E_{\rm kin}/E_{\rm SB}=3/10$. As we show in Section 5, including cooling in the interface does not change this prediction.

4.3 Conduction, evaporation, and cooling in the interface

We now examine the interface between the cool shell and hot bubble interior. Figure 3 shows a zoom-in of the shell at $t=3.5$ Myr, once shocks from discrete SNe have subsided and the interface has settled into quasi steady-state. We compare the fiducial simulation with conduction and nonlinear mixing (black) to the version with mixing but no conduction (gold). Within the shell and cooling region, the two simulations are very similar. In both cases, the interface between the bubble and shell is narrow, with the transition from $T\sim 10^{6}$ K and $n_{\rm H}\sim 10^{-2}\,\rm cm^{-3}$ in the hot interior to $T\sim 10^{2}$ K and $n_{\rm H}\sim 10^{2}\,\rm cm^{-3}$ in the shell occurring over a scale of $<1$ pc. Most of the cooling occurs in a narrow region in the interface on the inner edge of the shell, where the temperature passes through $T\sim 10^{4-5}$ K. Some additional cooling occurs in the outer edge of the shell, where shocked ISM being swept up by the expanding bubble is heated above the temperature of the ambient ISM.

In the simulation with conduction, the diffusive energy flux (the sum of conductive and mixing terms; 5th panel of Figure 3) is nearly flat interior to the shell’s inner edge. Conduction sets the shape of the temperature profile near the edge of the bubble such that the conductive heat flux, $q\sim-r^{2}T^{5/2}{\rm d}T/{\rm d}r$, is approximately constant with radius, at least over scales that are small compared to the size of the SB. Most of the conductive energy flux is deposited directly in the shell/bubble interface; i.e., the conductive flux drops rapidly in the region where cooling losses are significant and goes to 0 deep in the shell.

The dotted vertical line delimits the region to the left where (Spitzer) conduction dominates the diffusive energy flux, and to the right where mixing dominates the diffusive flux. For both simulations, $\dot{E}_{\rm diff}$ is similar in the high-density, low-temperature regime, where mixing dominates over Spitzer conduction and cooling becomes strong. In the simulation without conduction, the heat flux becomes negligible in the bubble interior, where $\kappa_{\rm mix}$ is small.

The advective energy flux from evaporation of the shell into the bubble (bottom panel) is significantly smaller than the diffusive flux at late times; as we discuss below, this is a consequence of radiative losses in the interface. In the simulation without conduction, the advective flux is small and often negative (i.e. mass and energy are advected into the shell from the interior). For the simulation snapshot shown at $t=3.5$ Myr, the advective energy flux in the simulation with conduction rises gradually toward the SB interior. However, the magnitude and shape of the advective flux profile varies with time as sound waves from individual SNe arrive at the interface.

We now compare the relative contributions of conduction, evaporation, and cooling in the interface region. We calculate the energy fluxes as follows. At each snapshot, we define the outer edge of the hot interior as the outermost zone with $T>10^{5}\,\rm K$. We then measure the fluxes through a surface 2 pc interior to this zone. Because the temperature and pressure gradients in the shell/bubble interface are steep and the fluxes vary rapidly there, measuring the conductive and advective fluxes through a surface slightly interior to the bubble leads to more stable results. We note that in this region, the conductive flux is much larger than the flux from mixing.

The total energy fluxes are:

[TABLE]

representing the energy diffusion (conduction plus mixing terms) into the shell and the advection out of the shell, respectively. Note that the sign of the velocity difference in $\dot{E}_{\rm advec}$ is such that an inward flow relative to the expanding shell (evaporation) is a positive advection flux. Correspondingly, a bulk flow that carries energy into the shell (condensation) is a negative advection flux.

Our results are not sensitive to the particular choice of measuring energy fluxes through a surface 2 pc interior to the shell boundary. As Figure 3 shows, the energy fluxes (particularly the advection term) are not strictly constant in radius at any particular time, but “ripples” from individual SNe explosions cancel out in a time-averaged sense, such that the time-averaged fluxes are consistent through any surface reasonably close to the shell boundary.

In calculating cooling losses, we wish to separate the cooling of shocked ISM exterior to the contact discontinuity (i.e., on the outer edge of the shell) from cooling in the interface between the hot interior and the cool shell. Only cooling in the interface on the inside of the shell, which removes energy that was carried in by conduction, is relevant for the energy balance between conduction and evaporation. Cooling in the outer edge of the shell primarily removes thermal energy from the $P{\rm d}V$ work done by the expanding bubble, so it should be discounted in calculations of the energy flux through the interface.

Gas can in principle mix through the shell, and some cooling does occur within the shell, so it is not possible to cleanly separate cooling of shocked ISM on the shell exterior from cooling in the interface. However, because the temperature in the shell at late times is always well below the peak of the cooling curve, the fraction of cooling that occurs within the shell is sub-dominant. We define the boundary between cooling in the interface and cooling of shocked ISM on the shell exterior at the first local minimum in the temperature profile within the shell component. We thus define the cooling in the bubble/shell interface as:

[TABLE]

$\mathcal{L}_{\rm int}$ in principle includes some cooling within the bubble interior but Figure 3 shows that this is negligible compared to the cooling in the interface.

Figure 4 shows the energy flux through the interface as a function of time for the fiducial simulation with conduction and mixing. We compare the conductive energy flux into the shell, the advective energy flux out of the shell due to evaporation, and the cooling losses in the shell/interior interface. All quantities are normalized by the total time-averaged energy input rate from SNe. The top and bottom panels show the same fluxes, with the smoothing timescale longer in the bottom panel.

The classical Weaver model, which neglects cooling in the interface, predicts that energy fluxes from conduction and advection should exactly cancel. If this were true, the gold and green lines in Figure 4 would be equal. In fact this prediction does not hold true for the simulation: at late times, the conductive heat flux into the shell is a factor $\sim$5 larger than the advective energy flux from shell evaporation. Cooling losses are quite significant, with $\sim$4 times more energy radiated away from the interface than is evaporated from the shell. The relative importance of cooling and advection is relatively constant; i.e., the cooling efficiency does not vary significantly with time. The dotted line in the bottom panel of Figure 4 shows the total interface energy losses from the combination of evaporation and cooling. This tracks the total energy conducted into the interface quite closely. This means that the basic picture of energy balance in the shell does hold, once cooling in the interface is accounted for.

We next consider the mass flux into the bubble interior from evaporation. In the classical model without cooling, this is predicted by Equation (11), but cooling in the shell/bubble interface significantly reduces the evaporation rate. In fact, if cooling is sufficiently strong, evaporation can be reversed and replaced with condensation of hot gas onto the shell (e.g. McKee & Cowie, 1977).

Figure 5 shows the mass flux from the shell into the bubble. The blue line shows the true mass flux measured in the simulation, $\dot{m}=4\pi r^{2}(v_{\rm shell}-v)\rho$, while the gold line shows the mass flux predicted if advection and conduction exactly cancel and cooling in the interface is negligible (Equation 11), as is assumed in classical models. The true time-averaged mass flux is a factor of $\sim$5 lower than this prediction at late times, because cooling rather than advection removes most of the energy from the interface.

When cooling in the shell-bubble interface is important, a term $\mathcal{L}_{\rm int}$ must be added to the left-hand side of Equation (10), representing the losses from cooling in the interface. In this case, less evaporation from the shell is required for energy balance, because some of the energy conducted into the shell is lost to cooling. The revised mass flux is then

[TABLE]

The green line in Figure 5 shows the mass loss predicted from energy balance in the shell after cooling is accounted for (Equation 29). This matches the true mass flux within 10% at late times; any residual disagreement can likely be attributed to ambiguity in the separation of cooling in the interface from cooling of shocked ISM on the outside of the shell, which can lead the total cooling in the interface to be slightly over- or under-estimated. The hot gas mass in the bubble can thus be understood purely as a consequence of the evaporation required to carry away “surplus” conduction energy that radiation is unable to remove.

The red line in Figure 5 shows the mass flux rate predicted by the Weaver solution (Equation 12). It is not identical to the gold line because the shape and normalization of the temperature profile inside the simulated SB differ somewhat from the predictions of the Weaver solution (see Figure 1 and 2). The higher temperature in the simulation is itself a result of cooling, which leads to less mass in the bubble interior.

Figures 4 and 5 also show that at early times ($t\lesssim 3$ Myr for the fiducial simulation parameters), all the energy and mass fluxes fluctuate significantly on a timescale comparable to $\Delta t_{\rm SNe}$. In particular, the mass and energy advection fluxes become negative for a brief period once every $\Delta t_{\rm SNe}$, when the blast wave from a new SNe reaches the shell. During these periods, hot gas is being pushed into the shell instead of flowing out of it. While some minor fluctuations persist at later times, the advective fluxes in the bubble remain positive after $t\approx 3$ Myr. As we will show in Section 6.1, the time after which blast waves from individual SNe cease pushing mass into the shell corresponds to the time when they become subsonic before reaching the shell.

5 Modified solution with cooling

5.1 Dynamics

Our simulations show that at late times, the fraction of the total energy input that is lost to cooling is approximately constant in time, varying primarily with the ambient ISM density. We now provide a simple analytic argument to predict the bubble dynamics in the presence of cooling.

We define the fraction of the energy input that is lost to cooling in the interface to be $\theta$, such that the cooling rate in the interface is $\mathcal{L}_{{\rm int}}\equiv\theta\dot{E}_{{\rm in}}$. The SB momentum equation (Equation 1) is unchanged from the case without cooling, but in the internal energy equation (Equation 2), $\dot{E}_{\rm in}$ is replaced with $(1-\theta)\dot{E}_{\rm in}$. All the dynamical predictions of classical SB theory (Section 2.1; Equations 3-9) still hold, after the substitution $\dot{E}_{\rm in}\to(1-\theta)\dot{E}_{\rm in}$.

We can then write the radius in terms of dimensional quantities as

[TABLE]

The radial evolution from Equation (30) can also be translated to a prediction for the radial momentum per supernovae, $\hat{p}_{{\rm rad}}=\left(4\pi R^{3}\rho_{0}/3\right)\left({\rm d}R/{\rm d}t\right)/\left(t/\Delta t_{{\rm SNe}}\right)$:

[TABLE]

In Figure 6, we compare the evolution of $R(t)$ and $\hat{p}_{\rm rad}(t)$ predicted by Equations (30) and (31) to a suite of simulations with a range of ambient densities (top panel) and a range of SNe rates (bottom panel). Computing the analytic prediction requires an estimate of $\theta$ consistent with the input parameters for each simulation; we calculate this as described in the next section. We shall show that (at late times) $\theta$ is expected to depend on ambient density $\rho_{0}$ but not on $\Delta t_{\rm SNe}$. At late times, the agreement between the simulations and analytic prediction is good. At early times, the simulated SBs are larger and grow more slowly with time than predicted by Equation (30).

The top panels of Figure 6 show that the time at which the simulations reach the behavior predicted by Equations (30) and (31) does not vary with ambient density (and thus, it does not vary with $\theta$): in all cases, the simulations converge to the predicted behavior after $1-2$ Myr, corresponding to $10-20$ SNe. The bottom panels show that the time to reach the predicted radial evolution does vary with $\Delta t_{\rm SNe}$: simulations with lower SNe rates reach the analytically predicted behavior at later times.

For simulations with a low ambient density and/or high SN rate, the SB radius can reach several hundred pc by the end of the simulation. We run all simulations for 10 Myr to facilitate straightforward comparison. However, we note that our idealized simulations with a uniform ambient medium become less realistic as the bubble grows large ($R\gg 100$ pc), since real galactic disks are vertically stratified and SBs become elongated and break out once the bubble radius becomes comparable to the disk scale height (Mac Low & McCray, 1988; Mac Low et al., 1989; Kim et al., 2017; Kim & Ostriker, 2018; Fielding et al., 2018).

For the set of simulations at fixed $\Delta t_{\rm SNe}$, the predicted and actual values of $\hat{p}_{\rm rad}$ are quite insensitive to density. This is true in both early stages dominated by discrete SN events, and in later stages that approach the cooling-modified pressure-driven bubble solution. While Equation (31) includes an explicit positive dependence of $\hat{p}_{\rm rad}$ on $\rho_{0}$, this is more than compensated by the decrease of $1-\theta$ with $\rho_{0}$. Given the dependence of $\theta$ on $\rho_{0}$ that we derive in Section 5.2, we expect $\hat{p}_{{\rm rad}}\propto\rho_{0}^{-1/5}$ at $n_{\rm H,0}\gg 1\,\rm cm^{-3}$. For the set of simulations at fixed $n_{\rm H,0}$, the increase of $\hat{p}_{\rm rad}$ with $\Delta t_{\rm SNe}$ is only due to the explicit dependence on $\Delta t_{\rm SNe}$ in Equation (31), as $\theta$ is expected to be independent of $\Delta t_{\rm SNe}$. We note that agreement between the prediction and simulation also occurs later at larger $\Delta t_{\rm SNe}$.

5.2 Predicting the cooling efficiency

The interface cooling efficiency, $\theta={\cal L}_{\rm int}/\dot{E}_{\rm in}$, is expected to depend on the ambient density and on the efficiency of nonlinear mixing as parameterized by $\lambda\delta v$. To estimate $\theta$, we assume that cooling occurs only in a nonlinear mixing layer, where the conductivity is given by Equation (21), and that the pressure is constant within the interface.777We have verified in our simulations that once equilibrium is reached at late times, pressure is always constant (within $\sim$ 5%; see e.g. Figure 3) through the interface and that most cooling losses ($\gtrsim 95$ %) indeed occur in the nonlinear mixing layer. The evaporative energy flux becomes negligible in the cooling layer, where the temperature and enthalpy of the gas are low, so the advection term in the energy equation may be neglected (see Figure 3). For simplicity, we neglect photoelectric heating so the volumetric net cooling rate is $n^{2}\Lambda(T)$.

We introduce a local coordinate $\xi=r/R$, where $R$ is the shell radius. Energy flux conservation then requires ${\rm d}\left(-\kappa_{\rm mix}{\rm d}T/{\rm d}r\right)/{\rm d}r=-n^{2}\Lambda$, or

[TABLE]

Note that in the above, we have only worked to lowest order in $\xi$, neglecting curvature of the front since the cooling/mixing layer is very thin.

Based on the sharply peaked cooling profile shown in Figure 3, we adopt the ansatz that most of the cooling occurs at a characteristic temperature $T_{\rm pk}$, so that we can approximate $\Lambda(T)/T^{2}\approx\Lambda(T_{\rm pk})/T_{\rm pk}^{2}$. From dimensional analysis, we assume that in the mixing region, the curvature of the temperature profile follows $d^{2}\ln T/d\xi^{2}\approx\alpha/(\Delta\xi)^{2}$ within a localized layer of thickness $\Delta\xi$, where $\alpha$ is a dimensionless constant. We can then write

[TABLE]

The total cooling rate in the interface is

[TABLE]

Finally, we substitute for $\Delta\xi$ from Equation (33) and use Equations (3), (4), and (6) for $R$, $P$, and $\dot{E}_{\rm th}$ with the modification $\dot{E}_{\rm in}\to(1-\theta)\dot{E}_{\rm in}$ . With the net cooling ${\cal L}_{\rm int}\equiv\theta\dot{E}_{\rm in}$, we obtain

[TABLE]

Because $T_{\rm pk}$ depends primarily on the adopted cooling function (and does not vary significantly with $\lambda\delta v$ or $\rho_{0}$; see Appendix E), we thus expect

[TABLE]

where $A_{\rm mix}$ is a dimensionless constant. If we adopt $\alpha=1$ and $T_{\rm pk}=2\times 10^{4}$ K (in rough agreement with our simulations; see Appendix E), we obtain $A_{\rm mix}\approx 1.7$. From our simulations, we can directly measure the ratio $\mathcal{L}_{{\rm int}}/\dot{E}_{{\rm th}}$, and we find $A_{\rm mix}\approx 3.5$ to provide a good fit to the numerical results (Figure 7).

The cooling efficiency, $\theta$, can be expressed in terms of the ratio of cooling luminosity to the rate of thermal energy increase as

[TABLE]

where the dependence of the quantity $\mathcal{L}_{{\rm int}}/\dot{E}_{{\rm th}}$ on physical parameters is predicted by Equation (37).

Finally, we note that the cooling efficiency is independent of time (i.e. ${\cal L}_{\rm int}\propto t$) because the time dependence of the width of the mixing region exactly compensates for that of $n^{2}\Lambda(T)R^{2}$.

5.3 Predicting the mass flux rate

With the cooling rate specified, we now compute the mass flux rate into the bubble interior in the presence of cooling. This quantity does not follow immediately from energy balance, because the conductive heat flux into the interface also changes when cooling becomes important.

Following Cowie & McKee (1977), we begin by setting energy fluxes out of the shell, which result from the combination of evaporation and cooling, equal to the heat flux into the shell from conduction:

[TABLE]

Following Weaver et al. (1977) (see also Section 2.2), we define a temperature scale

[TABLE]

for the SB interior that depends on the coefficient $C$ in the conductivity (see Section 3.1). This has the same parameter and time dependencies as the “no-cooling” expression in Equation (13), with an additional factor $(1-\theta)^{8/35}$ due to cooling. We can now write the temperature profile in terms of a dimensionless function $\tau$ as

[TABLE]

we use $\xi=r/R$ for the dimensionless radial coordinate as before. Dividing Equation (39) by $4\pi R^{3}P/t=2\dot{E}_{\rm th}$, we can write

[TABLE]

In Equation (42) we have introduced the parameter $w\equiv\dot{m}/\dot{m}_{\rm nc}$ to characterize the hot gas evaporation rate relative to the classical mass flux rate $\dot{m}_{\rm nc}$ predicted in the absence of cooling, where we substitute $T_{s}$ for $T_{\rm int}$ in Equation (12) for $\dot{m}_{\rm nc}$.

Equation (42) is an ordinary differential equation for $\tau$ as a function of $\xi$ with parameters $w$ and ${\cal L}_{\rm int}/\dot{E}_{\rm th}$; this determines the behavior of $\tau\left(\xi\right)$ near the shell/interior interface. Physically, we know that $\tau\to 0$ while $|d\tau/d\xi|\gg 1$ when $\xi\to 1$, so that the first term on the left-hand side of Equation (42) may be neglected. We thus have $\tau^{5/2}d\tau/d\xi=(2/7)d\tau^{7/2}/d\xi\approx-(1/2)\mathcal{L}_{\rm int}/\dot{E}_{\rm th}$. The definition in Equation (35) allows us to write

[TABLE]

Using the $\tau(1)=0$ boundary condition, we find that the limiting behavior of $\tau\left(\xi\right)$ near $\xi\to 1$ is

[TABLE]

This asymptotic behavior is different from that in the classical case with no cooling, where the 2nd term on the left-hand side of Equation (42) vanishes, yielding $\tau\left(\xi\right)\sim\left(1-\xi\right)^{2/5}$ (Weaver et al., 1977; Cowie & McKee, 1977).

In Figure 8, we show the temperature profiles for simulations with five different ambient densities at late times ($t\approx 10$ Myr). The solid lines show the actual profiles from the simulation, while the dashed lines show the prediction of Equation (44). In calculating the analytic prediction, we use Equations (37) and (38) to predict $\theta$, assuming $A_{\rm mix}=3.5$. We also multiply the prediction of Equation (44) by a constant factor of 1.24 in all cases. As we describe in Appendix D, saturation of conduction in the simulations, which was not accounted for in the calculations above, slightly reduces the mass flux in the simulations and increases the temperature relative to the analytic prediction. Once this constant multiple is included, the scaling with $\theta$ and spatial dependence on $\xi$ predicted by Equation (44) is in good agreement with the simulations.

The scalings of $\tau$ and $T_{s}$ with $\theta$ and other parameters of the problem can be combined to predict the temperature in the interior of the SB. From Equation (44), we can parameterize the temperature profile near $\xi\to 1$ as $T=T_{{\rm sb}}\left(1-\xi\right)^{2/7}$, implying a predicted interior temperature $T_{{\rm sb}}\sim 1.8[\theta/(1-\theta)]^{2/7}T_{s}$ (including the factor 1.24 to account for saturated conduction). As a scaling relation, we have

[TABLE]

where $T_{0}\approx 7.7\times 10^{6}\,{\rm K}$ is a calibrated constant. The Spitzer conductivity coefficient $C$ is a constant in our simulations, but we include the predicted scaling with $C$ as an indication of how bubble properties are expected to change if the conductivity were reduced by the presence of a magnetic field.

Finally, an estimate for the mass flux from the shell into the bubble can be calculated by combining Equation (45) with the cooling-modified version of Equation (6). We note that because the temperature is not constant throughout the interior, the outer value $T_{\rm sb}$ should only be considered an approximate mean value. With this, we have:

[TABLE]

Adopting the scalings in Equation  (45), we solve to obtain an expression for $\dot{m}$:

[TABLE]

Calibrating to the simulations, we find $\dot{m}_{0}\approx 2900\,M_{\odot}\,\rm Myr^{-1}$. Integrating Equation (47) yields the total mass of hot gas in the SB. The hot gas mass produced per SNe is:

[TABLE]

This can be combined with the radial evolution from Equation (30) to predict the mean number density of the interior, $n_{{\rm sb}}=3M_{{\rm hot}}/(4\pi R^{3}\mu m_{p})$. We find

[TABLE]

We note that our analytic predictions for $T_{\rm SB}$ and $\dot{m}$, and thus, for $\hat{M}_{\rm hot}$ and $n_{\rm SB}$, break down at small $\theta$ ($\theta\lesssim 0.05$, corresponding to $n_{\rm H,0}\lesssim 0.001\,\rm cm^{-3}$), where they predict $\dot{m}\to\infty$ and $T_{\rm sb}\to 0$, rather than recovering the scaling predicted by the Weaver solution. This occurs because in the limit of small $\theta$, the cooling term in Equation (42) no longer dominates the left-hand side, leading to a different asymptotic form of $\tau(\xi)$ and no dependence of $\tau$ on $\theta$.

Figure 7 compares the results of a suite of simulations with different $n_{\rm H,0}$ and $\lambda\delta v$ to the scaling relations we predict for $\mathcal{L}_{{\rm int}}/\dot{E}_{{\rm th}}$ (top panel; Equation 37), $\theta$ (second panel; Equation 38), $\dot{m}$ (3rd panel; Equation 47), and $T_{\rm sb}$ (4th panel; Equation 45). The bottom panel shows the ratio of the total energy flux out of the interface from advection and cooling to the heat flux into the interface; this is predicted to be $\approx 1$. All the simulations have $\Delta t_{\rm SNe}=0.1$ Myr. For context, we also show the predictions of the Weaver solution in the bottom three panels.

We measure all the quantities plotted at $t\approx 10$ Myr. At late times, shocks due to individual SNe have usually subsided near the outer regions of the bubble, causing the energy fluxes to reach quasi-equilibrium. However, as we discuss in Section 6.1, this does not occur in simulations with very efficient cooling (large $n_{\rm H,0}$ and $\lambda\delta v$), because the mass of hot gas in the bubble interior is reduced, causing SNe exploding within the bubble to sweep up less mass as they expand and decelerate more slowly than they would without cooling. In such cases, the quantities plotted in Figure 7 can fluctuate significantly between the arrival of blast waves. We minimize the effects of these fluctuations by averaging all quantities over a timescale of $\Delta t_{\rm SNe}$.

In calibrating our analytic predictions for the cooling efficiency (Equation 37), interior temperature (Equation 45), and evaporation rate (Equation 47), we choose coefficients such that the prediction exactly matches the simulation with $n_{\rm H,0}=1\,\rm cm^{-3}$ and $\lambda\delta v=3\,\rm pc\,km\,s^{-1}$. All other analytic curves shown in Figure 7 for differing $n_{\rm H,0}$ and $\lambda\delta v$ adopt the same value of the respective coefficient.

Comparison of the simulated and predicted values of $\dot{m}$ and $T_{\rm sb}$ with the Weaver solution shows that cooling significantly reduces the evaporation rate, and modestly changes the temperature. More cooling (larger $n_{\rm H,0}$ and/or larger $\lambda\delta v$) is associated with greater reductions in evaporation, and greater increases in the temperature of the hot gas. We further note that a decrease in the Spitzer conductivity coefficient $C$ would decrease $\dot{m}$ and increase $T_{\rm SB}$.

We compare the prediction of Equation (48) to the simulations in Figure 9, where we show simulations with a range of ambient density (top) and supernovae rate (bottom). There is significant disagreement with the simulation prediction at early times, both because $M_{\rm hot}$ is enhanced prior to shell formation, and because the evaporative mass flux is reduced at early times when shocks from discrete SNe propagate into the shell (Figure 5).

The bottom panel of Figure 9 shows that the simulation with the highest SN rate, $\Delta t_{\rm SNe}=0.01$ Myr produces a somewhat lower $\hat{M}_{\rm hot}$ than predicted by Equation (48), even though shocks from individual SNe become negligible quite early. The primary reason for this is that the effects of saturation of conduction become more significant when the SB interior is hotter, and $T_{\rm SB}$ increases with increasing SNe rate. Saturation suppresses evaporation into the SB and is not included in our analytic models, leading to a modest overprediction of $\hat{M}_{\rm hot}$ at smaller $\Delta t_{\rm SNe}$.

Overall, we regard the agreement between the simulations and analytic predictions in Figures 6 - 9 as quite satisfactory. We emphasize that while a normalization coefficient is calibrated using the simulations, the scaling relations with ambient density $n_{\rm H,0}$, SN interval $\Delta t_{\rm SNe}$, and mixing parameter $\lambda\delta v$ are not. Nevertheless, Figures 6 and 9 show that the predicted analytic scaling behaviors are only reached after the SB has evolved for some time. We return to the issue of early vs. late evolution in Section 6.1.

Some deviations from the analytic prediction for $\mathcal{L}_{{\rm int}}/\dot{E}_{{\rm th}}$ are expected due to issues such as (a) contributions of the photoelectric heating term to $\mathcal{L}_{\rm int}$, particularly at high $n_{\rm H,0}$, (b) cooling occurring at a range of temperatures rather than at a single temperature $T_{\rm pk}$, and (c) some cooling occurring outside the nonlinear mixing region.

Deviations from the analytic predictions for $\dot{m}$ and $T_{\rm sb}$ are expected to arise due to deviations from the predicted $\mathcal{L}_{{\rm int}}/\dot{E}_{{\rm th}}$ as well as from saturation of conduction and the discreteness of individual SNe. The largest deviation of $\dot{m}$ and $T_{\rm sb}$ in Figure 7 from the analytically predicted scaling relation occurs at large $\lambda\delta v$ and high ambient densities. This is the regime where cooling is most efficient and deviations from the classical Weaver solution are largest. The effects of the saturation of conduction are also largest when cooling is very efficient, because blast waves from individual SNe remain supersonic when they reach the shell until late times.

5.4 Evolution as a function of radius

The assumption of an unstratified homogeneous medium breaks down when the SB is sufficiently large, for $R_{\rm bubble}$ comparable to the scale height of the host galaxy’s ISM disk, $H\sim 100-300$ pc. The “feedback yields” delivered to the ISM are also limited by their values when $R_{\rm bubble}\sim H$. It is therefore useful to quantify the specific energy, momentum, and hot gas mass as a function of $R_{\rm bubble}$.

Figure 10 shows the time evolution of the specific radial momentum, hot gas mass, and thermal and kinetic energy as a function of SB radius for simulations with a range of ambient densities and SN rates. The SB properties at $R_{\rm bubble}\sim 100-200$ pc quantify approximately how much energy and momentum are available to drive turbulence prior to breakout. The value of $\hat{M}_{\rm hot}$ when the SB reaches a similar radius is indicative of the mass available to drive a hot galactic wind.

Solving Equation (30) for $t$ allows us to rewrite the analytically predicted scaling relations in terms of SB radius:

[TABLE]

Equations (50)–(52) have the same explicit dependence on $R$, $\rho_{0}$, and $\Delta t_{{\rm SNe}}$ as in the classical Weaver solution but include an additional dependence on $\theta$. In the Weaver solution, $\hat{M}_{\rm hot}(R)$, $\hat{E}_{\rm th}(R)$, and $\hat{E}_{\rm kin}(R)$ are independent of $\rho_{0}$. This is no longer the case when cooling is included, because $\theta$ increases with $\rho_{0}$. The comparison between the numerical results and analytic predictions show that the “late-time” solutions are generally reached by the time $R_{\rm bubble}\sim 100$ pc, so that the predicted scalings with density and SN rate are satisfied.

At $R_{\rm bubble}\sim 100-200$ pc, our models predict $10^{5}\lesssim\hat{p}_{{\rm rad}}/\left(M_{\odot}\,{\rm km\,s^{-1}}\right)\lesssim 6\times 10^{5}$, $40\lesssim\hat{M}_{{\rm hot}}/M_{\odot}\lesssim 300$, $6\times 10^{49}\lesssim\hat{E}_{{\rm th}}/{\rm erg}\lesssim 3\times 10^{50}$, and $3\times 10^{49}\lesssim\hat{E}_{{\rm kin}}/{\rm erg}\lesssim 1.5\times 10^{50}$. After breakout, a mass of hot gas not exceeding $\hat{M}_{{\rm hot}}$ per SN would be able to vent out of the disk. The venting of hot gas also effectively limits both the momentum injection to the ISM and to the galactic wind, as shown in simulations with clustered star formation in stratified disks (e.g. Kim & Ostriker, 2018; Fielding et al., 2018).

6 Early evolution

6.1 Effects of discrete shocks

The analysis in Section 5 primarily considered the limit in which repeated SNe can be approximated as a continuous energy source with luminosity $\dot{E}_{\rm in}=E_{\rm SN}/\Delta t_{\rm SNe}$. This approximation is not valid at early times. We now consider the effect of shocks from individual SNe.

Figure 11 shows the structure of the shell/bubble interface at an early time ($t\approx 0.3$ Myr), during and shortly after the arrival of a blast wave from a new SN at the shell. The left panels show the shell/interface region just as the shock arrives and begins to compress the interface. There is a strongly negative advective flux into the bubble interior (i.e., hot gas from the interior is being pushed into the interface), as well as a strong heat flux into the shell. The right panels show the same fluxes a short period later, after the shock has compressed the dense shell and the direction of the advective fluxes has reversed, with both mass and energy flowing from the interface into the bubble interior. The cooling rate in the interface is also significantly elevated over its value before the shock’s arrival. The energy fluxes into and out of the interface are more than an order of magnitude larger during and shortly after the arrival of the shock at the interface than their equilibrium values (see Figure 3). As we discuss below, this leads to rapid increases in $M_{\rm hot}$ followed by a period of enhanced cooling.

To further assess the effects of shocks from SNe within the SB, we trace the blast waves from individual SNe as they travel through the bubble interior. In each simulation snapshot following a new SNe explosion, we identify the shock from the most recent SN as the cell with maximum pressure. We then track the blast wave as it expands by identifying the local maximum in $P(r)$ in each snapshot that is nearest to the predicted location of the shock based on its position and velocity in the previous snapshot. We follow each blast wave until it reaches the shell or fades (i.e., no suitable local maximum can be identified in $P(r)$ that is associated with the location of the shock in the previous snapshot).

Figure 12 (top panel) shows the radial expansion of blast waves from individual SNe in the first 4 Myr of the fiducial simulation. We calculate the velocity of each blast wave (2nd panel) as $v_{{\rm shock}}={\rm d}r_{{\rm shock}}/{\rm d}t$. We compute the Mach number, $\mathcal{M}$, of each shock as a function of time (3rd panel) from the adiabatic sound speed, $c_{s}=\sqrt{\gamma P/\rho}$. Because gas at the pressure local maximum has already been heated by the shock, we calculate $c_{s}$ in front of the shock; i.e., in the medium into which it is propagating. In practice, we use the median sound speed in the region $(2-5)$ pc in front of the pressure local maximum to calculate $\mathcal{M}$.

Wiggles in the blast wave velocity curves at late times correspond to periods when the blast waves collide and interfere with reflected waves from previous SNe in the bubble interior. The sound speed in the bubble generally decreases with radius (see Figure 1), so the shock velocity drops faster than the Mach number. In some cases, the Mach number actually increases as the blast wave approaches the shell and the temperature falls; in such cases, the sound speed falls more steeply with radius than the shock velocity.

Shocks from individual SNe remain supersonic when they reach the shell at early times, only becoming subsonic after $\sim$3 Myr. This transition coincides with the time when $M_{\rm hot}$ ceases to increase stepwise, and SNe blast waves stop driving significant hot gas mass into the shell (Figure 5). It is important to distinguish when individual blast waves are supersonic because this determines (a) whether they sweep up significant mass from the interior, and (b) whether shocks propagate into and directly heat gas in the shell.

As long as the blast waves remain supersonic, they slam into the shell impulsively. This can result in (a) a brief period when a large amount of hot gas is carried in and/or kinetic energy is partly converted to thermal energy as the shock hits the shell, leading to enhanced conduction, and (b) a sudden increase in the kinetic energy and momentum of the shell, as well as $M_{\rm hot}$. It is also possible for swept-up hot gas deposited in the shell to condense and cool there, leading to a net loss of energy and hot gas mass from the bubble. Models with continuous energy injection are less applicable at early times, when blast waves within the bubble reach the shell supersonically. At later times, once blast waves become subsonic before reaching the shell, most of their kinetic energy has been thermalized, and the continuous luminosity approximation is appropriate.

We derive a scaling relation for the time after which blast waves from individual SNe become subsonic before reaching the shell as follows. We suppose that a strong shock propagating through the interior of the bubble can be represented as a front with velocity $v_{s}$ containing a fraction of the mass in the SB interior, such that $E_{\rm SN}=fM_{\rm hot}v_{s}^{2}$. Here the constant $f$ encapsulates uncertainty in the fraction of mass that is swept up and the fraction of the blast wave’s energy that is kinetic. We also have $M_{\rm hot}c_{s}^{2}\sim(50/99)(1-\theta)E_{\rm SN}(t/\Delta t_{\rm SNe})$ from integrating over the life of the bubble. We thus expect that shocks will become subsonic at a time $t_{\rm subsonic}\sim(2/f)\Delta t_{\rm SNe}/(1-\theta)$. In practice, we find

[TABLE]

where the constant “10” is found by calibrating to the simulations. Equation (53) can also be cast in terms of the total number of SNe that have exploded, as $N_{{\rm SNe,\,subsonic}}\approx 10/\left(1-\theta\right)$.

We compare this prediction to the simulations in Figure 13, where we plot the time after which SNe become subsonic for simulations with a range of $n_{\rm H,0}$, $\lambda\delta v$, and $\Delta t_{\rm SNe}$. We measure $t_{\rm subsonic}$ in the simulations as the first time at which the Mach number of an blast wave falls below $\mathcal{M}=1$ prior to reaching the shell. As in the previous sections, we compute $\theta$ for the analytic prediction using Equations (37) and (38). The overall agreement with Equation (53) is quite good. For a typical ISM density of $n_{\rm H,0}=1\,\rm cm^{-3}$, we thus predict that a relatively massive cluster with $M_{\rm star}\sim 4\times 10^{4}\,M_{\odot}$ (corresponding to $\Delta t_{\rm SNe}=0.1$ Myr) will transition from the discrete-SNe regime to the regime of continuous energy injection after $\approx 3$ Myr. Small clusters with $M_{\rm star}\lesssim 2\times 10^{3}M_{\odot}$ will never enter the continuous energy injection regime, because SNe will cease before $t_{\rm subsonic}$.

Given the radial evolution predicted by Equation (30), we can also formulate Equation (53) in terms of the radius of the SB:

[TABLE]

with $R_{0}\approx 83$ pc. Implicit in Equation (54) is the assumption that SNe do not cease before $t_{\rm subsonic}$. For the fiducial bubble parameters considered in Section 4, this yields $R_{\rm subsonic}\approx 120$ pc, which is comparable to the radius at which the bubble is expected to begin breaking out of the galactic disk (see Section 5.4).

Mac Low & McCray (1988) also considered whether shock waves from individual SNe within a SB should become subsonic. They assumed that the bubble interior was well-described by the Weaver et al. (1977) similarity solution (with no cooling) and found, contrary to the result of our simulations, that blast waves from individual SNe should essentially always become subsonic before reaching the shell. We can understand the differences in these conclusions as follows: The deceleration of blast waves within the SB is due to the sweeping up of mass in the bubble interior. Cooling reduces the mass in the interior and thus increases the velocity and Mach number of individual SNe blast waves. The scaling relations from Figure 13 predict the Mach number to scale as $\mathcal{M}\sim\left[\Delta t_{{\rm SNe}}/\left(\left(1-\theta\right)t\right)\right]^{1/2}$, so $\mathcal{M}$ is larger when $\theta$ is nonzero. This is the primary reason for the discrepancy between the prediction from Mac Low & McCray (1988) and our simulation results.888Equation (53) appears to predict that blast waves within the SB should remain supersonic for the first $\sim$ 10 SNe when $\theta=0$. As discussed in Section 5, the solution for the interior structure of the SB that we derive is not expected to hold in the limit of $\theta\to 0$. Using the SB properties predicted by the Weaver solution and reproducing the calculations from Mac Low & McCray (1988) leads to the same scaling relation with $\Delta t_{\rm SNe}$ that is predicted by Equation (53). However, in this case the analytically predicted coefficient to $\Delta t_{\rm SNe}$ is $\approx 1.3$, while we found it to be $\approx 10$ when we calibrated to the simulation results with cooling.

6.2 Effects of stellar winds

Our calculations thus far have only included energy injected by SNe, neglecting the effects of stellar winds and radiation from young stars. Integrating over the whole lifetime of a young star cluster, the energy injected into the ISM by SNe dominates over that injected by stellar winds by an order of magnitude (e.g. Leitherer et al., 1999), and the dynamical impact from expansion of H ii regions (driven by both radiation pressure and ionized gas pressure) is also an order of magnitude or more lower than that of SNe (e.g. Kim et al., 2018). However, stellar winds and radiation are dominant for the first few Myr of a cluster’s life, before SNe begin (e.g. Geen et al., 2015; Agertz et al., 2013). We therefore now consider the effects of stellar winds on the early-time evolution of SBs.

For a solar-metallicity cluster with a Kroupa (2001) IMF, the mechanical wind luminosity is roughly constant, $L_{{\rm winds}}\approx 9\times 10^{33}\,{\rm erg\,s^{-1}\times\left(M_{{\rm star}}/M_{\odot}\right)}$, over the first $\sim$3 Myr, after which it declines steeply (e.g. Leitherer et al., 1999). This is, coincidentally, quite similar to the mean mechanical SNe luminosity at later times, which is $L_{{\rm SNe}}\approx 8\times 10^{33}\,{\rm erg\,s^{-1}\times\left(M_{{\rm star}}/M_{\odot}\right)}$. This means that, for any SN rate, the effects of winds can be approximated by adding a $3\,\rm Myr$ period of continuous energy injection prior to the first SNe.

We have carried out a set of simulations with $n_{\rm H,0}=(0.1-10)\,\rm cm^{-3}$ and $\Delta t_{\rm SNe}=(0.03-0.3)\,\rm Myr$ that include a 3-Myr long period of continuous energy injection prior to the explosion of the first SN. The time-averaged wind luminosity is taken to be $\dot{E}_{\rm in}=10^{51}\,\rm erg/\Delta t_{\rm SNe}$, equal to the energy input rate once SNe commence, and winds cease as soon as SNe begin. We neglect any additional feedback due to ionization or radiation pressure in these calculations, since its character is quite different from that of mechanical feedback due to winds and SNe. Following Leitherer et al. (1999), the injected mass flux associated with the wind is a factor of 3.7 higher than the value for SNe; i.e., $37M_{\odot}/\Delta t_{\rm SNe}$ instead of $10M_{\odot}/\Delta t_{\rm SNe}$.

Figure 14 compares the evolution of $\hat{p}_{\rm rad}$ and $\hat{M}_{\rm hot}$ in the fiducial simulation to one in which energy injection is through winds for the first 3 Myr. For both simulations, the SB’s age is measured from the time when energy injection begins. For the simulation with winds, normalized quantities at $t<3$ Myr are obtained by dividing by $(t/\Delta t_{\rm SNe})$ as opposed to $\left\lfloor t/\Delta t_{{\rm SNe}}\right\rfloor+1$. That is, for both simulations, we normalize $\hat{p}_{\rm rad}$ and $\hat{M}_{\rm hot}$ relative to the amount of energy injected so far.

Overall, the global evolution the SB is quite similar in the two simulations, with $\hat{p}_{\rm rad}$ and $\hat{M}_{\rm hot}$ generally agreeing within a factor of two. At $t>3$ Myr, once SNe have commenced in the simulation with winds, the two simulations agree within a few percent. This is perhaps unsurprising, because the total energy injected at fixed time is the same after $t=3$ Myr, with or without winds. Once several SNe have exploded, the SB’s global properties are well described by the continuous energy injection approximation even in the case of discrete SNe explosions (e.g. Figure 6). We find good agreement in global SB properties between simulations with and without winds over the full range of SNe rates and ambient densities we tested. Differences between simulations with and without winds are generally larger for SBs with lower SNe rates or higher ambient densities.

Although the SB’s global properties are relatively insensitive to whether energy is injected continuously or through discrete SNe, inclusion of winds can diminish the importance of shocks at early times, because there is more time for hot gas to evaporate into the wind-driven bubble before SNe begin. If the SB’s age is measured from the time of the first SN, then including winds effectively shifts the SB’s evolution in time, as if SNe had begun 3 Myr earlier. This additional hot gas mass decelerates shocks when individual SNe begin, such that they can become subsonic at earlier times relative to the onset of SNe.

We find that for cases where $t_{\rm subsonic}$ predicted by Equation (53) is less than 3 Myr, inclusion of winds causes all SN blastwaves to become subsonic before reaching the shell. For cases where $t_{\rm subsonic}$ predicted by Equation (53) exceeds 3 Myr, we find that when winds are included (and $t$ is measured from the time when winds begin), the time at which SN blasts become subsonic always agrees with the prediction of Equation (53) within 20%. With the inclusion of winds, shocks from early SNe still reach the shell supersonically when the ambient density is high ($n_{{\rm H},0}\geq 1\,{\rm cm}^{-3}$ for $\Delta t_{\rm SNe}=0.1\rm\,Myr$) or the SN rate is low ($\Delta t_{{\rm SNe}}\geq 0.1\,{\rm Myr}$ for $n_{{\rm H},0}=1\,{\rm cm}^{-3}$), but not when the SNe rate is high or the ambient density is low.

The above results suggest that winds can diminish the importance of shocks from discrete SNe at early times. However, in practice ambient ISM densities are high enough and most cluster masses are low enough that the early SNe will still have supersonic blasts. Moreover, the stage when winds are important coincides with the time clusters are still embedded within dense ($n_{\rm H,0}>100\,{\rm cm}^{-3}$) molecular clouds. With the expected $\theta\gtrsim 0.9$ for gas at these densities (see Equations  37 and 38), Equation (30) suggests that SB evolution would in many cases transition from wind-dominated to supernova-dominated even before breaking out of the natal cloud. In such situations, SBs within the low-density diffuse ISM will be driven primarily by SNe.

Additionally, we note that observations of young clusters find feedback from stellar winds to be weaker than naively expected. That is, the X-ray luminosities observed from the regions surrounding young clusters ($\rm age\lesssim 3\,\rm Myr$, when no SNe are expected to have occurred yet) are much lower (typically by a factor of $20-50$) than predicted by models that take the wind luminosity from stellar population synthesis and assume the kinetic energy of winds is efficiently thermalized through shocks (e.g. Chu et al., 2003; Harper-Clark & Murray, 2009; Rosen et al., 2014; Leroy et al., 2018). These observations remain imperfectly understood theoretically: potential explanations include slower-than-predicted winds from massive stars, efficient turbulent mixing of shocked winds with the cool ISM, and escape of hot gas through holes in the shell at early times. Indeed, at the densities of molecular clouds, the high values of $\theta$ would be consistent (through Equation 52) with an order of magnitude or more reduction in the hot-gas energy that powers X-rays. In any case, the fact that the observable signatures of winds in real SBs are weaker than naively expected suggests a corresponding reduction of the dynamical effects of winds in real SBs, motivating the primary focus on the SN-driven stage of SB evolution in this paper.

7 Summary and Discussion

We have carried out idealized 1D simulations of superbubbles (SBs) driven by repeated supernovae (SNe), focusing in particular on the effects of thermal conduction and radiative cooling near the interface between the cool shell and the hot bubble interior. In order to account for the effects of interface instabilities that occur in real SBs but are not captured in 1D simulations, we also include a diffusion term in the energy equation to represent nonlinear mixing.

Unlike most previous studies, which either neglected conduction or assumed that conduction is balanced exactly by evaporation into the shell while neglecting cooling, we model conduction, cooling, and evaporation self-consistently. By modifying the model of Weaver et al. (1977) to include cooling losses, we have derived an analytic model to predict the thermal properties and dynamical evolution of SBs.

Our primary results are as follows:

Bubble expansion and momentum: With cooling, the bubble expansion rate is reduced compared to the classical Weaver solution, and the momentum injection per supernova is correspondingly reduced (Figure 2 and Equations 30-31). For our fiducial simulation with one SN every 0.1 Myr and an ambient density of $n_{\rm H,0}=\rm 1\,cm^{-3}$, the radial momentum per SN, $\hat{p}_{\rm rad}$, is reduced by a factor $\sim 2$ compared to the no-cooling case. For higher ambient densities, there is a monotonic increase in the cooling efficiency, and a decrease in the expansion rate and momentum. However, $\hat{p}_{\rm rad}$ of a few times $10^{5}M_{\odot}\,{\rm km\,s^{-1}}$ is still typical for the full parameter range we have considered (Figure 6 and 10). 2. 2.

The bubble interior: Thermal conduction substantially increases the mass of hot gas in SB interiors (Figures 1 and 2) relative to models without conduction, although this is still reduced below that of the classical Weaver solution due to cooling in the shell. For our fiducial simulation parameters, $\hat{M}_{\rm hot}\sim$150 $M_{\odot}$ of gas evaporates from the cool shell into the hot interior per SN. Values of $\hat{M}_{\rm hot}$ are somewhat lower for higher density and SN rate (see Figures 9,10). Without conduction, only $\sim$10 $M_{\odot}$ of hot gas is produced per SN; this is attributable entirely to the ejecta and circumstellar material. The temperature of the bubble interior is much cooler ($T\sim 10^{6-7}$ K) with conduction than without it ($T>10^{8}$ K). 3. 3.

Shocks from individual SNe: At early times ($t\lesssim 3$ Myr and $R\lesssim 120\,\rm pc$ for our fiducial simulation parameters), blast waves from individual SNe remain supersonic throughout their traversal of the bubble interior (Figure 12). This is not predicted for the classical SB solution (e.g McCray & Kafatos, 1987; Mac Low & McCray, 1988); it occurs because cooling losses reduce the gas mass swept up by expanding shocks in the interior, causing blast waves to decelerate more slowly than in the classical solution. The duration of the period dominated by discrete supersonic blasts depends on the time interval between SNe and the cooling efficiency (Figure 13 and Equations 53-54), such that SBs around lower-mass clusters and SBs in gas with higher ambient density will remain in the discrete-SNe limit longer.

When blast waves from individual SNe reach the shell before becoming subsonic, they remain strong shocks that carry energy and mass into the shell, accompanied by an enhancement in the conductive heat flux. This quickly heats the inner layer of the shell and precipitates an outflux of hot gas into the bubble (Figures 11 and 12). In all models, there is a brief period following the arrival of a shock when the energy flux into the shell exceeds what can be carried away by cooling, leading to an enhanced mass flux into the bubble. Nevertheless, the time-averaged energy advection is into rather than out of the shell at early times, and the time-averaged cooling is higher at early times compared to late times (Figure 4).

We have tested the effects of including a 3-Myr period where energy is injected continuously due to stellar winds before SNe begin (Section 6.2). This has little effect on the global evolution of the SB because the mechanical luminosity of winds is similar to that of SNe. However, because including winds allows more hot gas to accumulate in the SB interior before SNe begin to explode, winds tend to decrease the effects of shocks from individual SNe. Nevertheless,except for models with very low ambient density and high SN rate, the first SNe to occur will still produce shocks when they collide with the SB shell. 4. 4.

Cooling in the shell/bubble interface: Once the SB reaches a state in which transient blast waves can be neglected, energy is lost from the interior only by conduction to, and pressure work on, the shell. The conductive heat flux from the hot bubble into the shell is balanced by the combination of (a) mass evaporation from the shell into the bubble and (b) cooling at a rate ${\cal L}_{\rm int}$ in the shell/bubble interface (Figures 3, 4). For a typical ambient density of $n_{\rm H,0}=1\,\rm cm^{-3}$ and our fiducial mixing efficiency, cooling losses are dominant: $\mathcal{L}_{{\rm int}}/\dot{E}_{{\rm cond}}\approx 0.8$ and ${\cal L}_{\rm int}/\dot{E}_{\rm in}\approx 0.6$. In other words, roughly 80% of the energy conducted into the interface is cooled away, and $\sim$ 60% of the total bubble energy is lost to cooling in the interface. As a result, the mass flux from evaporation of the shell into the interior is only $\sim$20% of what it would be without cooling losses (Figure 5). Once cooling losses are subtracted, the mass flux can be fully accounted for by what is needed to balance the conductive heat flux. The classical similarity solution of Weaver et al. (1977) did not account for cooling in the interface, so we find significantly less mass in the bubble interior than was predicted by this earlier model. 5. 5.

New analytic models: We construct analytic models to predict the dynamical evolution of SBs in the presence of cooling (Section 5.1), the dependence of the cooling efficiency on SB parameters (Section 5.2), and the evaporation rate and interior temperature profile of SBs with cooling (Section 5.3). Our models build on the classical solution for bubbles driven by continuous energy injection (Weaver et al., 1977; McCray & Kafatos, 1987), but self-consistently account for cooling losses as well as effects of thermal conduction.

The effects of cooling in our analytic model are encapsulated in a dimensionless, time-independent cooling efficiency, $\theta\equiv{\cal L}_{\rm int}/\dot{E}_{\rm in}$, which depends only on the ambient density and the efficiency of turbulent mixing at the interface (Equations 37-38). For $n_{\rm H,0}=1\,\rm cm^{-3}$ and our fiducial mixing parameter $\lambda\delta v=1{\rm pc\ kms^{-1}}$, $\theta\approx 0.61$. For a given $\theta$, we provide scaling relations for the evolution of a SB’s radius, internal pressure, radial momentum, and total energy (see Equations 30-31 and related discussion in Section 5.1). In addition, we provide predictions for the interior temperature, hot gas mass, and density (Equations 45-49). We also predict the shape of the temperature profile at the edge of the bubble, which differs from the no-cooling case (Equation 44). Our simulations are in good agreement with the analytically predicted scaling relations (Figures 6 - 10 and 13). 6. 6.

Role of mixing and conduction in setting SB properties: Cooling losses are set primarily by turbulent mixing, not conduction (Figure 2). The dynamics and many global properties of the SB (size, momentum, pressure, energy) are uniquely specified by the cooling efficiency $\theta$. In turn, $\theta$ depends on the ambient ISM density and on the efficiency of nonlinear mixing (Figure 7). Thermal conduction sets the evaporation rate of the shell, and thus, the temperature and density of the bubble interior. Because the evaporation rate also depends on the amount of cooling losses in the interface, the interior properties of the SB depend on both the cooling efficiency and thermal conductivity (Section 5.3). However, for any nonzero level of conduction, the interior temperature is fairly insensitive to the values of parameters (Equation 45). 7. 7.

Numerical convergence: Without conduction or nonlinear mixing, it is very difficult to achieve convergence in 1D simulations. Most cooling losses occur in the thin interface between the hot bubble interior and cool shell, and the thickness of this interface becomes arbitrarily small at high resolution in simulations without conduction or mixing, leading to negligible cooling losses. Including thermal conduction alone broadens the interface and somewhat increases cooling losses but does not lead to full convergence at the resolution level we can achieve. Inclusion of diffusive mixing broadens the interface further and leads to more rapid convergence (Figure 16). 8. 8.

Saturation of conduction: Particularly at early times, conduction is often saturated (i.e., the classical conductive flux exceeds what can actually be transported by electrons). Failure to account for saturation of conduction leads to a modest but non-negligible ($\sim$ 20%) overestimate of the hot gas mass in the bubble (Figure 18).

7.1 Limitations of 1D models

Because our simulations are 1D and highly idealized, assuming that SNe explode at uniform time intervals in an unstratified homogeneous medium, several processes that affect real superbubbles are imperfectly modeled.

Most significantly, 1D simulations cannot explicitly model instabilities that lead to mixing of gas of different temperature phases. We represent the effects of turbulent mixing and heat diffusion using an effective conductivity (Equation 21). This model is motivated by shear-driven instabilities (Section 3.2) and has one free parameter $\lambda\delta v$; i.e, the wavelength of unstable modes times the turbulent velocity, which quantifies the efficiency of mixing. As shown in Section 5.2 and Appendix B, this parameter determines the cooling rate for a given ambient density and thus enters as a scaling factor in many of the SB properties predicted by our model.

Since cooling depends on $\lambda\delta v$, one possible approach is to calibrate $\lambda\delta v$ based on the cooling rates in 3D SB simulations. Presently available results are generally consistent with the range of $\lambda\delta v$ that we have considered (See Appendix B), but were not designed to enable this calibration. Some quantitative uncertainty remains regarding interface cooling. Previous works that included an inhomogeneous and/or turbulent ambient medium (Kim et al., 2017; Fielding et al., 2018) found converged cooling and momentum injection rates at moderate resolution in 3D. However, instabilities in such simulations may not be have been fully resolved. Other simulations in which the ambient medium is uniform and instabilities grow from the resolution scale have not found convergence (e.g. Gentry et al., 2019). Moreover, the amount numerical diffusion decreases rather slowly as the resolution improves (Agertz et al., 2007; Hopkins, 2015).

Cooling losses in SBs also affect the injection of the momentum that drives ISM turbulence (and regulates star formation), and the mass of hot gas that can escape as a wind. Thus, present limitations in understanding of interface mixing remains a persistent source of uncertainty in theoretical research in several areas. As such, obtaining more precise constraints on the “true” value of $\lambda\delta v$ represents an important open problem, which we will explore in future work. This will require very high resolution 3D simulations of SB expansion in models of the ISM with realistic ambient inhomogenous structure. Repeating the simulations presented here in 3D at the same resolution is currently prohibitively expensive computationally, but it may be possible to resolve the interface instabilities using adaptive mesh refinement. Given a calibration of mixing, it may be possible to implement it as a subgrid diffusion model in galactic and CGM simulations where resolution is inadequate to directly capture the relevant turbulent mixing.

Our simulations also do not include magnetic fields. Because electrons move along magnetic field lines, magnetic fields tend to suppress conduction perpendicular to the magnetic field (e.g. Braginskii, 1965). It is thus possible that our simulations somewhat overestimate the true conductive flux from the bubble to shell. However, the expulsion of gas from the shell into the bubble is expected to comb the magnetic field in the radial direction, such that any magnetic field frozen into the shell may have little effect on the conductive energy flux (e.g. Cowie & McKee, 1977). The RM instability also tends to produce radially aligned magnetic fields (Inoue et al., 2013). 3D simulations with conduction are required to quantify the effect of magnetic fields on the conductive heat flux. Previous work has also suggested that magnetic fields suppress the development of instabilities in the shell/interior interface (e.g. Ntormousi et al., 2017). In the context of the 1D models studied in this work, suppression of conduction due to magnetic fields can be approximated as decreasing the constant $C$ in the thermal conductivity. The scaling relations we provide (Equations 45-49) show that the effects of re-scaling thermal conductivity are relatively weak: for example, lowering $C$ by a factor of 10 only reduces $M_{\rm hot}$ by a factor of 1.9.

7.2 Comparison to previous work

Compared to the classical SB similarity solution (e.g. Weaver et al., 1977; McCray & Kafatos, 1987; Mac Low & McCray, 1988), our simulations predict SBs to be smaller, have lower energy and momentum, and contain less hot gas, usually with slightly higher temperature. The magnitude of the difference between our predictions and the classical solution depends on the ambient density and efficiency of mixing in the interface (Figure 7). For our fiducial problem setup, the momentum per SN $\hat{p}_{\rm rad}$ is a factor of $2-3$ lower than in the classical solution; in SBs with high ambient density and efficient mixing, cooling reduces $\hat{p}_{\rm rad}$ by up to a factor of $\sim$ 10.

Although they did not account for cooling in the interface in their similarity solution, Weaver et al. (1977) and McCray & Kafatos (1987) did note that at late times, cooling would eventually become important in the SB interior. They predicted that at this time, a SB’s expansion would cease to be pressure-driven, causing the momentum equation and predicted radial evolution to change. In particular, McCray & Kafatos (1987, their Equation 11) predicted that a SB’s radial evolution would depart from the $R\sim t^{3/5}$ behavior predicted for a pressure-driven SB after a cooling time of

[TABLE]

so that, for example, the interior of a SB with $\Delta t_{\rm SNe}=0.1$  and $n_{\rm H,0}=10\,\rm cm^{-3}$ would cool after $t_{c}\sim 5$ Myr.

None of our simulations reach such a regime within the first 10 Myr: Figure 6 shows that they all follow $R\sim t^{3/5}$. We also experimented with running the simulation with $n_{\rm H,0}=10\,\rm cm^{-3}$ (which is the most likely to cool according to Equation 55) for 40 Myr. Even in this case, we do not find any significant departure from the predicted $R\sim t^{3/5}$ behavior.

The reason for this discrepancy is that cooling becomes important in the interface long before it becomes important in the hot interior. As we show in Section 5.3, cooling in the interface (a) increases the temperature and (b) decreases the density in the interior. Both effects delay interior cooling. This delay is more severe at high ambient density and large $\theta$ (Equation 49), such that we do not expect cooling in the interior to become important for any of the range of SB parameters studied in this work.

Kim et al. (2017) recently carried out 3D simulations of SBs in an inhomogeneous medium consisting of cold clouds formed via thermal instability embedded in a warm intercloud medium. Their simulations included cooling and allowed for development of instabilities at the hot/cold interface, but did not include thermal conduction. The hot gas mass in the interior of their simulated SBs thus came primarily from the ablation of cold clouds within the SB, rather than from shell evaporation as it does in our simulations. Despite contributions to $M_{\rm hot}$ from ablation of cold clouds, they found relatively low values of $M_{\rm hot}$, with $\hat{M}_{\rm hot}=(10-100)M_{\odot}$ in most cases. This owed in part to the fact that clouds in the interior often cooled quickly even after being surrounded by the hot SB interior. Fielding et al. (2018) also found, from 3D simulations of SBs evolving in an inhomogeneous, turbulent medium, that $\hat{M}_{\rm hot}=(10-100)M_{\odot}$ prior to breakout from the ISM.

In our simulations and analytic model, the fraction of the total injected energy that is retained in the SB approaches a constant value (Figure 2). In contrast, Kim et al. (2017) found that the bubble energy per SN declined steadily at late times. This discrepancy is very likely due to the fact that the ambient medium in the simulations of Kim et al. (2017) was cloudy, while the ambient medium in our simulations is uniform. When comparing 3D simulations of SBs in a uniform medium to otherwise identical simulations in an inhomogeneous, turbulent medium, Fielding et al. (2018) found that those in a uniform medium approached a constant energy per SNe, while the energy per SNe decreased with time in the case of a turbulent medium. In the context of the models we develop, a decrease in $\hat{E}_{\rm SB}$ at late times is expected if the efficiency of mixing, parameterized by $\lambda\delta v$, increases with time.

Yadav et al. (2017) also found an approximately constant asymptotic energy per SNe in 3D simulations with a uniform medium. Which behavior is more realistic is not yet clear: the ISM is indeed inhomogeneous and turbulent, but numerical diffusion, which can lead to excess cooling losses, is also enhanced in an inhomogeneous medium. Yadav et al. (2017) also included conduction in one of their 3D simulations of SBs driven by young star clusters. In agreement with our results, they found that although conduction increases $M_{\rm hot}$, it does not significantly change the dynamics or total energy of the SB because the degree of cooling loss is primarily set by turbulent mixing. One question of interest for future work is whether this conclusion holds true in a turbulent medium, where conduction might be expected to accelerate the ablation of cold clouds within the SB interior.

While $\hat{M}_{\rm hot}$ in our models with conduction is significantly enhanced relative to simulations without conduction, the values of the radial momentum per SN at $t\lesssim 3$Myr are generally similar to recent results from 3D simulations. In particular, Kim et al. (2017) found $\hat{p}_{\rm rad}\sim(0.8-3)\times 10^{5}\,{\rm M_{\odot}\,km\,s^{-1}}$, relatively insensitive to ambient density but increasing with $\Delta t_{\rm SNe}$, comparable to the results shown in Figure 6 at similar times. This is expected, as the early times in our simulations are dominated by shocks from individual SN, similar to the results in 3D simulations without conduction. At late times, however, the pressure in the bubble interior becomes more important in driving expansion than individual SNe shocks in our simulations, and the values of $\hat{p}_{\rm rad}$ rise above those in 3D simulations without conduction. In the Kim et al. (2017) simulations, the pressure within SBs at late times was reduced by efficient mixing and cooling. An interesting question for future research is whether, during the period before SB breakout, $\hat{p}_{\rm rad}$ can be enhanced by conduction, when considering multiphase ISM disks with realistic scale heights and cluster masses (e.g. Figure 10).

7.3 Implications and prospects

One of the important conclusions of this work is that, compared to the levels produced by the combination of supernova ejecta and ablation by KH instabilities of embedded clouds by high-velocity hot flows, conduction can significantly enhance the mass of hot gas created in SBs. The enhancement of $\hat{M}_{\rm hot}$ by conduction has potentially important implications for phase balance in the ISM and circumgalactic medium, as well as for galactic winds.

In 3D simulations of SBs expanding in a two-phase ISM, Kim et al. (2017) found that both warm and hot gas were accelerated by SBs, but only the hot gas had high enough velocity that it might be expected to escape from an intermediate- or high-mass galaxy’s gravitational potential and drive a galactic wind. Fielding et al. (2018) similarly found little warm gas with $v>100{\rm km/s}$ in their simulations of SBs in a turbulent medium. Assuming that one SN explodes for every $\sim 100\,M_{\odot}$ of stars formed (and that only the hot gas escapes in a wind) the Kim et al. (2017) results correspond to a relatively low wind mass loading factor, $\eta=\dot{M}_{{\rm outflow}}/{\rm SFR}$, with $0.1\lesssim\eta_{\rm hot}\lesssim 1$. The results for mass loading from Fielding et al. (2018) were similar, even after breakout from the disk. Kim et al. (2017) thus concluded that SBs are only likely to launch winds with $\eta_{\rm hot}>1$ if the SNe rate is sufficiently high that SBs break out of the disk before shell formation, when cooling is not yet important and $\hat{M}_{\rm hot}$ is much higher. This is only expected to occur if the SNe rate is very high, as from a massive cluster (see also Fielding et al., 2018), and in particular if there is a burst of star formation above typical local values (this rate determines the disk scale height in the first place).

Due to the conduction-driven evaporation in our simulations, we predict values of $\hat{M}_{\rm hot}$ that are roughly an order of magnitude higher than those found by Kim et al. (2017) (see Equation 48 and Figure 9). This implies mass loading factors $\eta_{\rm hot}$ that are usually larger than 1. Such values of $\eta$ are more compatible with the values needed to explain the low cosmic star formation efficiency, though they are still lower than the values typically found in current cosmological simulations (e.g. Muratov et al., 2015; Anglés-Alcázar et al., 2017). We note, however, that there is not necessarily a simple correspondence between $\hat{M}_{\rm hot}$ and the mass that is actually removed from galaxies in winds: hot outflowing gas can both entrain additional gas in the galaxy or CGM after leaving a SB, increasing $\eta$, or be sufficiently decelerated by the material it sweeps up that it fails to escape, decreasing $\eta$ (e.g. Muratov et al., 2015; El-Badry et al., 2018; Fielding et al., 2018). Nevertheless, Kim & Ostriker (2018) found, using self-consistent simulations of star-forming three-phase ISM disks with uniformly high resolution in both midplane and outflow/fountain regions, that the value $\eta_{\rm hot}\sim 0.1$ in hot winds was comparable to that found in idealized SB simulations without conduction.999In an extended suite of simulations for different galactic conditions covering gas surface density up to $100\,M_{\odot}\,{\rm pc}^{-2}$ (Kim et al., in prep), the mass loading factors of hot winds also do not vary much, with typical $\eta_{\rm hot}\sim 0.1$. The interaction of hot, outflowing gas with warm fountain gas did not significantly accelerate the latter, nor reduce the energy flux of the former.

The possibility of a large increase of $\hat{M}_{\rm hot}$ in SBs due to thermal conduction has important implications for galactic winds from normal galaxies. In addition to increasing the hot-gas mass outflow loading to $\eta_{\rm hot}>1$, a greater momentum flux in the hot medium could boost the acceleration of warm gas out of galaxies, possibly increasing the wind loading factor in galaxies with deep potential wells.

Finally, we remark that our analytic solution for SBs with cooling can be easily implemented in any application that currently uses the Weaver et al. (1977) model. In addition to use in interpretation of observed supershells (e.g. Oey & Massey, 1995; Jaskot et al., 2011), other potential uses include setting the inputs to models of starburst-driven winds for varying environments (e.g. Thompson et al., 2016), comparing to outflows driven by nascent super-star clusters within dense molecular clouds (e.g. Leroy et al., 2018), and implementing subgrid models for galaxy formation simulations (e.g. Keller et al., 2014).

Acknowledgements

We thank the referee for a thoughtful report. We are grateful to Chris Mckee, Yuan Li, and Ken Shen for helpful discussions, and to Drummond Fielding for providing an implementation of radiative cooling in Athena++. KE was supported by an NSF graduate research fellowship. The work of ECO was supported in part by ATP grant NNX17AG26G from NASA and Investigator grant 510940 from the Simons Foundation. CGK was supported in part by grant 528307 from the Simons Foundation. EQ and KE are supported by a Simons Investigator Award from the Simons Foundation and by NSF grant AST-1715070. DRW is supported by a fellowship from the Alfred P. Sloan Foundation and acknowledges support from the Alexander von Humboldt Foundation. This work was initiated as a project for the Kavli Summer Program in Astrophysics held at the Center for Computational Astrophysics of the Flatiron Institute in 2018. The program was co-funded by the Kavli Foundation and the Simons Foundation. We thank them for their generous support.

Appendix A Numerical convergence

In Figure 15, we show the effects of varying the simulation resolution on integrated SB properties, for the fiducial model. Simulations with (without) conduction and diffusive mixing (adopting $\lambda\delta v=1{\rm pc\ km\ s^{-1}}$) are shown in the left (right) panels. Because the simulations with conduction are computationally more expensive than those without it, the four highest resolution levels ($\Delta r\leq 0.003$ pc) are only reached in the no-conduction, no-mixing simulations.

The right panels show that the simulations without conduction or diffusive mixing do not fully converge, even at the highest resolution level: increasing the resolution leads to higher SB energy and radial momentum, and a larger bubble radius. Radiative losses in the shell/interior interface are progressively reduced at higher resolution. Indeed, the highest-resolution simulation approaches the $\hat{E}_{\rm SB}/E_{\rm SN}\approx 0.65$ value predicted by the Weaver solution if there is no cooling within the bubble or interface. The simulations with conduction and mixing (left panels) converge much more rapidly: although there are some minor differences between resolution levels at early times, none of the integrated SB properties show systematic trends with resolution at late times. In addition, the simulations without conduction show that the hot gas mass is lower than with conduction, and does not converge with resolution.

The effects of conduction and mixing on the rate of cooling in the interface are investigated further in Figure 16, which shows the cooling rate within the interface (excluding cooling of shocked ISM on the shell exterior; see Section 4.3) as a function of resolution. This comparison is made at late time ($t\approx 10$ Myr), such that effects from individual SN shocks are no longer important. In the top panel, we compare simulations without conduction or diffusive mixing (solid gray), those with normal Spitzer conduction but no diffusive mixing (open black), those with mixing but no other conduction (open red), and those with both conduction and mixing (solid red).

Without conduction or mixing, the cooling rate decreases with increasing resolution, without any apparent limit. At late times, shocks do not propagate into the dense gas in the interface, and nothing else in the problem sets a physical scale for the transition between hot bubble and cool shell. The volume of the SB with $T\sim 10^{5}\,\rm K$ and efficient cooling thus drops with increasing numerical resolution, such that cooling losses becomes negligible at arbitrarily high resolution. At infinite resolution, $\mathcal{L}_{\rm int}$ is thus expected to go to zero (see also Gentry et al., 2017).

Including conduction makes the trend with resolution less steep, but does not lead to full convergence over the range of $\Delta r$ where our simulations with conduction are computationally tractable. This likely occurs because (temperature-dependent) thermal conductivity becomes very small at the cool outer edge of the interface bordering the shell. Conduction alone is expected to lead to convergence once $\Delta r$ is small compared to the minimum Field length in the problem (McKee & Cowie, 1977; Begelman & McKee, 1990), but the Field length varies substantially throughout the interface. At $n_{\rm H}=1\,\rm cm^{-3}$ and $T=2\times 10^{4}\,{\rm K}$ (representative of the value at the outer edge of the interface; see Figure 3), $\lambda_{\rm Field}\approx 10^{-3}$ pc, which is smaller than we can achieve in simulations with conduction. Our results are consistent with those of Fierlinger et al. (2016), who found that thermal conduction alone did not lead to convergence in 1D simulations of SNe-driven bubbles, except in an extreme case in which thermal conductivity was artificially enhanced by 14 orders of magnitude.

Irrespective of whether conduction is included, including diffusive mixing with $\lambda\delta v=1\,\rm pc\,km\,s^{-1}$ (red lines in the top panel of Figure 16) leads to fairly rapid convergence. This reflects the fact, already shown in Figure 2, that cooling losses are set by mixing, not conduction. Numerically, this is because the mixing coefficient is not temperature dependent, while the conductivity coefficient is. Thus, the thickness of the layer where there is mixing and therefore strong cooling does not continue to drop at high resolution, leading to convergence in the net cooling rate (see also discussion in Section 5.2).

The bottom panel of Figure 16 compares the convergence of simulations for different values of $\lambda\delta v$. We do not include thermal conduction in these simulations, as it increases the computational expense and does not significantly affect the cooling rate. As expected, simulations with smaller $\lambda\delta v$ (less efficient mixing) require higher resolution to be converged. At the resolution used for most of the simulations in this work ($\Delta r=(0.024-0.048)$ pc), only simulations with $\lambda\delta v\geq 1\,\rm pc\,km\,s^{-1}$ are fully converged. Simulations with lower, nonzero values of $\lambda\delta v$ do converge at higher resolution. As we showed in Section 5.2, the converged cooling rates increase with increasing $\lambda\delta v$.

Appendix B Efficiency of mixing

In Figure 17, we compare the evolution of the total energy, radial momentum, hot gas mass, and bubble radius in simulations with different values of $\lambda\delta v$ at fixed resolution. Consistent with the predictions in Section 5, increasing $\lambda\delta v$ leads to more efficient cooling and thus lower SB energy, momentum, hot gas mass, and size. The simulation with $\lambda\delta v=0.1\,\rm pc\,km\,s^{-1}$ appears almost identical to the one with no nonlinear mixing at all ($\lambda\delta v=0\,\rm pc\,km\,s^{-1}$). However, Figure 16 shows that the simulations with $\lambda\delta v<1\,\rm pc\,km\,s^{-1}$ are not fully converged at this resolution. For these two simulations, all quantities shown in Figure 17 would be somewhat larger at higher resolution.

Our results thus directly demonstrate and quantify how cooling depends on the efficiency of mixing across what would otherwise be a contact discontinuity (see discussion in Fierlinger et al. 2016; Yadav et al. 2017; Kim et al. 2017; Fielding et al. 2018; Gentry et al. 2019). In our 1D simulations, the mixing efficiency is realized via the effective diffusivity $\lambda\delta v$. In reality, this mixing is the result of nonlinear evolution of interface instabilities.

In principle, as discussed in Section 7.1, one way to constrain the effective value of $\lambda\delta v$ from instability-driven mixing would be to calibrate to 3D simulations. However, previous work provides only rather limited information on this. Some recent studies have been dominated by numerical diffusion at the resolution they can achieve. Other studies at higher resolution were not designed to study late-time evolution in the limit of small $\Delta t_{\rm SNe}$. We conclude that a focused study will be required to properly calibrate $\lambda\delta v$, as results from existing work are not sufficient to address this question.

Appendix C Conductivity Ceiling

High values of $\kappa$ require prohibitively small timesteps for numerical stability. In particular, for Spitzer conductivity the CFL condition (Courant et al., 1928) requires that the simulation timestep, $\Delta t$, satisfies

[TABLE]

At high temperatures and low densities, this condition can become extremely restrictive. For example, immediately after new SN explosions, $T\sim 10^{8}\,{\rm K}$ and $n_{\rm H}\sim 10^{-3}$ near the center of the SB, requiring a timestep of $\sim$15 seconds at our default resolution! Because a single global timestep is used for all cells in the simulation, high values of $\kappa$ in regions with low density will grind the simulation to a halt.

We prevent arbitrarily small timesteps by placing a density-dependent ceiling on thermal conductivity:

[TABLE]

This upper limit was chosen such that at a density of $n_{\rm H}=10^{-1}\,\rm cm^{-2}$, $\kappa_{\rm max}$ is the Spitzer value (Equation 16) for a temperature of $10^{7}\,\rm K$. The density dependence of $\kappa_{\rm max}$ effectively sets a minimum timestep at fixed resolution:

[TABLE]

The actual minimum timestep in the simulations is 0.3 times $\Delta t_{\rm min}$, because we use a CFL number of 0.3. We experimented with different normalizations of $\kappa_{\rm max}$ to determine the smallest value that can safely be used without significantly changing the simulation results. Conduction is important primarily in the interface between the cool shell and hot bubble interior, where the temperature is lower and the density is higher than near the center of the SB. The ceiling on $\kappa$ does change the conductive energy flux and temperature profile near the center of the bubble, but this has a negligible effect on the integrated properties of the SB. We have verified that the adopted ceiling on $\kappa$ changes the hot gas mass (as well as the SB radius, total energy, and radial momentum) by less than 1% from the values when no ceiling is used, while decreasing the typical computational expense by a factor of $\sim$ 100. Using a significantly smaller normalization in Equation (58) does begin to produce deviations from simulations without a ceiling; for example, decreasing the normalization by a factor of 5 leads to a $\sim$ 10% decrease in $M_{\rm hot}$.

Appendix D Saturation of conduction

Our fiducial simulations implement saturation as described in Section 3, adopting an effective conductivity such that the heat flux interpolates smoothly between the classical and saturated values.

We consider the effects of saturation in Figure 18. The left two panels compare the magnitude of the classical conductive heat flux (gold), the saturation limit (green), and actual flux implemented in the simulation (black), which interpolates from the classical value when $\kappa|\nabla T|\ll(3/2)\rho c_{s}^{3}$ to the saturation limit when $\kappa|\nabla T|\gg(3/2)\rho c_{s}^{3}$. We show two snapshots, one at early times, shortly after the arrival of the blast wave from a new SN, and one at later times, once the interface has reached quasi-equilibrium. Both snapshots are from the fiducial simulation, in which saturation is implemented. At early times (upper left panel), conduction is sometimes strongly saturated: after the arrival of a new SN blast wave at the interface, the classical conductive flux exceeds the saturation limit by more than an order of magnitude. At later times (lower left panel), the effects of saturation are modest, because the classical conductive flux is somewhat less than the saturation limit. The saturation prescription implemented in the simulation (Equation 20) does slightly reduce the conductive flux relative to the classical limit, but the change is less than a factor of two.

In the right panels of Figure 18, we compare the total SB energy and hot gas mass of the fiducial simulation (black) to one in which saturation has been turned off (i.e., the conductive heat flux always follows the classical prediction; gold). The differences in total energy between the two simulations are negligible at late times, implying that saturation does not significantly change the efficiency of cooling. On the other hand, there is a clear offset in $\hat{M}_{\rm hot}$ between the two simulations, in that the simulation without saturation produces somewhat more hot gas per SNe. We show in Section 4.3 that $\hat{M}_{\rm hot}$ is set by the difference between the conductive heat flux and the amount of cooling in the interface. Including saturation does not substantially change the amount of cooling, but it does reduce the conductive flux, leading to moderately less evaporation from the shell.

The Weaver solution does not include saturation; that is, it assumes that the conductive heat flux is always given by Equation (18). Because saturation reduces the evaporative mass flux into the bubble, it heightens one of the important differences between our simulations and the Weaver solution. However, the cumulative effects of saturation are modest overall, and the change in $\hat{M}_{\rm hot}$ due to saturation is much less than the change in the evaporation due to cooling in the interface.

Appendix E Cooling curve

In our derivation of the expected cooling efficiency as a function of $n_{\rm H,0}$ and $\lambda\delta v$ (Section 5.2), we made the approximation that most of the cooling in the interface occurs at a characteristic temperature $T_{\rm pk}$, and we assumed that $T_{\rm pk}$ does not vary much with $n_{\rm H,0}$ or $\lambda\delta v$.

Figure 19 examines these assumptions in more detail. The left panel shows the cooling curve used in our simulations (black line). The middle panel shows the temperature profile of a snapshot of the fiducial simulation at $t=10$ Myr, once the interface is in flux equilibrium and transient effects from individual SNe have died down. We color the temperature profile by the cell-integrated net cooling rate, $\dot{E}_{\rm cool}=n(n\Lambda-\Gamma)\times 4\pi r^{2}\Delta r$, in each cell. Most of the cooling losses occur in a narrow region where the temperature is between $T\approx 10^{4}$ K and $T\approx 2\times 10^{5}$ K. Some additional cooling occurs at $T<10^{3}$ K on both sides of the shell. In the center of the shell, the temperature reaches an equilibrium where cooling and heating exactly balance.

Most of the cooling occurs below the peak of the cooling curve. As discussed in Section 5.2, the pressure in the interface is approximately constant. Although $\Lambda(T)$ peaks at $T\sim 10^{5}\,\rm K$, the volumetric cooling rate, $\sim n^{2}\Lambda(T)\sim P^{2}\Lambda(T)/T^{2}$ peaks at $T\sim 2\times 10^{4}$ K for fixed pressure $P$, and this is the temperature where most of the cooling occurs. This is shown in the right panel, which plots the temperature where cooling is most efficient, for simulation with a range of $n_{\rm H,0}$ and $\lambda\delta v$. For each simulation, we calculate $T_{\rm pk}$ as the temperature where the cooling efficiency, ${\rm d}\dot{E}_{{\rm cool}}/{\rm d}\log T$, peaks. No significant trend in $T_{\rm pk}$ is evident with $n_{\rm H,0}$ or $\lambda\delta v$, validating the assumption in Section 5.2.