The Impact of Type Ia Supernovae in Quiescent Galaxies. I. Formation of the Multiphase Interstellar Medium

We use the run n0.02-T3e6 as an example to illustrate the 3 × 10⁶ K case. The cooling time for the medium at the beginning of the simulation is ${t}_{{\rm{c}},0}$ = 90 Myr. We define the time when any gas parcel cools below 2 × 10⁴ K as ${t}_{\mathrm{multi}}$, which is 75 Myr for this run.

The upper four panels of Figure 2 show slices of density, temperature, pressure, and color fraction at 70 Myr, right before the formation of cool gas. The SN color fraction, ${f}_{\mathrm{color}}$, is defined as the ratio between the color density and gas density (similar to metallicity), normalized to that of the SN ejecta. The medium is visually inhomogeneous, with signs of SN bubbles. The spatial variation of density and temperature is about a factor of 10. The pressure varies much less than density or temperature, although spherical blast waves/sound waves driven by SNe are visible. It is most obvious in the ${f}_{\mathrm{color}}$ panel, where part of the gas has been mixed with SN ejecta: areas with darker shades have been polluted by SN "colors," whereas regions with light shades have not been impacted by SNe directly. Comparing ${f}_{\mathrm{color}}$ with the other panels, one can see that low ${f}_{\mathrm{color}}$ corresponds to higher densities and lower temperatures, and vice versa.

Zoom In Zoom Out Reset image size

The lower four panels of Figure 2 show slices at t = 90 Myr, after the cool phase forms. The cool gas is seen as tiny clumps that have the highest densities and lowest temperatures (centered on white circles). Comparing to the plots at t = 70 Myr, one can see that the cool phase forms in regions with low ${f}_{\mathrm{color}}$. Note that the cool clumps do not form at the shock front of individual SNRs.

To better see the correlation of physical properties of the inhomogeneous gas, we plotted in Figure 3 two phase diagrams at t = 70 Myr, before the cool phase formation (same time as the first snapshot of Figure 2). The top panel shows the isochoric cooling time ${t}_{{\rm{c}},0}$ versus ${f}_{\mathrm{color}}$, color-coded by the total mass. At the beginning of the simulation, all gas resides in a single point on the diagram: ${t}_{{\rm{c}},0}$ = 90 Myr and ${f}_{\mathrm{color}}$ = 0. The dashed line shows the initial ${t}_{{\rm{c}},0}$. At t = 70 Myr, ${t}_{{\rm{c}},0}$ spans about 3 orders of magnitude, and ${f}_{\mathrm{color}}$ spans 5 orders of magnitude. There is a fairly tight correlation between ${t}_{{\rm{c}},0}$ and ${f}_{\mathrm{color}}$: gas with long ${t}_{{\rm{c}},0}$ tends to have higher ${f}_{\mathrm{color}}$. This is easy to understand: the heating and rarefaction effect of the blast waves make hot, tenuous bubbles, which arrive at a rough pressure balance with the ambient medium at ${R}_{\mathrm{fade}}$. For T ≳ 10⁵ K, at fixed pressure, ${t}_{{\rm{c}},0}$ is shorter for lower temperatures. New SN bubbles have the highest ${f}_{\mathrm{color}}$ and longest ${t}_{{\rm{c}},0}$. The tightest correlation of ${t}_{{\rm{c}},0}$ and ${f}_{\mathrm{color}}$, seen in the upper right corner of the plot, arises from these new SN bubbles. The majority of the gas has ${t}_{{\rm{c}},0}$ > 90 Myr, but some gas has ${t}_{{\rm{c}},0}$ ≪ 90 Myr. The short-${t}_{{\rm{c}},0}$ tail will form the cool phase within 10 Myr.

Zoom In Zoom Out Reset image size

The bottom panel of Figure 3 shows density versus temperature, color-coded by the isochoric ${t}_{{\rm{c}},0}$. The star symbol indicates the initial condition. The distributions of temperature and density lie roughly along a constant pressure. Gas now has temperatures and densities both above and below the original value. Gas with the shorter cooling time (darker color shade) has higher densities and lower temperatures. These features are consistent with the expectations described in Section 2. The medium becomes inhomogeneous under the impact of SNe Ia. The gas that is going to become part of the cool phase (i.e., has the shortest ${t}_{{\rm{c}},0}$) has the lowest ${f}_{\mathrm{color}}$, largest density, and lowest temperature. This gas has been outside of any SN bubbles. Also note that the overall gas pressure is somewhat higher than the initial condition, so the cooling process is not entirely isobaric. This can cause the cooling time to slightly deviate from ${t}_{{\rm{c}},1}$.

Now we turn to more quantitative analysis, by examining when cool gas starts to form and how much cool gas forms with time. Figure 4 shows the fraction of mass in the box that is in the cool phase (T < 2 × 10⁴ K), as a function of time. The cool phase occurs at about 0.8${t}_{{\rm{c}},0}$. After that, the mass of the cool phase appears to accumulate rapidly: by ${t}_{{\rm{c}},0}$, 13% of the mass in the box is cool; by 2${t}_{{\rm{c}},0}$, about 55% of the mass is cool. We warn, however, that the accumulation rate of the cool phase may not be accurate, since (1) the cool phase, once formed, is not well resolved; (2) some relevant physics—such as thermal conduction and sound-wave heating (sound waves can steepen into shocks after entering the cool gas)—are not included or captured; and (3) we have neglected the larger context and gravitational field of the elliptical. The dashed line indicates the theoretical estimate from Equation (8). Plugging in relevant parameters, Equation (8) shows that the onset of cool gas formation is predicted at 0.30${t}_{{\rm{c}},0}$, and 67 of the mass would be in the cool phase after that. The simulation results, on the other hand, show a later start of multiphase gas formation and a slower accumulation of cool mass over time. This suggests that the simple expectation where all gas unheated by SNe cools at ${t}_{{\rm{c}},1}$ has something missing, and possibly that hydrodynamical effects are playing a role. We will discuss this further in Section 5.

Zoom In Zoom Out Reset image size

Now we discuss a run with a higher initial temperature, 10⁷ K. We find that the multiphase gas is, once again, produced, but not until (2–3)${t}_{{\rm{c}},0}$ has elapsed, in contrast to the lower-temperature case, where the cool phase occurs before ${t}_{{\rm{c}},0}$. We use the run n0.08-T1e7 as an example. Such gas has an isochoric ${t}_{{\rm{c}},0}$ = 133 Myr. Figure 5 shows the fractional mass in the cool phase as a function of time (similar to Figure 4). As one can see, the cool phase does not form until t ≈ 2.7${t}_{{\rm{c}},0}$, which is about 375 Myr. By the end of the simulation, only 27% of mass is in the cool phase. It is expected, from Equation (8), that the cool phase formation time would be later, and the mass fraction of cool gas is less, compared to the lower-temperature run of n0.02-T3e6. For n0.08-T1e7, ${t}_{{\rm{c}},1}$ ≈ 0.85${t}_{{\rm{c}},0}$. Therefore, according to Equation (8), the cool phase would form at 0.85${t}_{{\rm{c}},0}$ with a mass fraction of 0.42, compared to 0.39${t}_{{\rm{c}},0}$ and 0.67% mass fraction for n0.02-T3e6. Quantitatively, however, the actually formation time of cool gas in the simulation is much more delayed.

Zoom In Zoom Out Reset image size

To better understand what is causing the delay, we plotted in Figure 6 the phase diagram of temperature versus density, color-coded by the isochoric ${t}_{{\rm{c}},0}$ at two times in this simulation. The top panel is for t = 120 Myr, slightly before ${t}_{{\rm{c}},0}$ (133 Myr), and the bottom panel is for t = 340 Myr, somewhat before the cool phase forms (375 Myr). The star symbols indicate the initial condition of the gas. At 120 Myr, the distributions of gas density and temperature do broaden over time, similar to what was seen for the n0.02-T3e6 run. That said, the shortest ${t}_{{\rm{c}},0}$, which is around 90 Myr, does not fall much below the initial value of 133 Myr. This means that not much gas has lost sufficient thermal energy to form a cool phase soon. The majority of the gas has a higher temperature and lower density, which corresponds to a longer ${t}_{{\rm{c}},0}$. Only at a much later time (bottom panel) does some gas reach a short enough ${t}_{{\rm{c}},0}$. At 340 Myr, the distributions of density and temperature lie above the initial condition. This implies that the overall gas has been heated considerably, and thus the thermal pressure is above the initial value.

Zoom In Zoom Out Reset image size

Besides the illustrated cases with T = 3 × 10⁶ K and 10⁷ K above, we extend the parameter space to a wide range of n and, in particular, have two additional simulations with T = 6 × 10⁶ K and 3 × 10⁷ K with the same n = 0.02 cm⁻³. All these simulations have the fiducial H/C = 1.02. We find that all these runs eventually develop multiphase gas, demonstrating the universality of such an outcome, even under the uneven heating of SNe Ia.

The time at which the cool gas forms, though, varies. Figure 7 shows ${t}_{\mathrm{multi}}$/${t}_{{\rm{c}},0}$ for all the simulations. The ratio ${t}_{\mathrm{multi}}$/${t}_{{\rm{c}},0}$ increases monotonically with increasing T. For the case with the highest temperature $T=3\times {10}^{7}$ K, cool gas forms at about 4.4${t}_{{\rm{c}},0}$. For the same T, the dependence of ${t}_{\mathrm{multi}}$/${t}_{{\rm{c}},0}$ on n is quite small. For example, all runs with T = 3 × 10⁶ K develop a cool phase within ${t}_{{\rm{c}},0}$, even when n varies by a factor of 40. This is also shown by the four runs with T = 10⁷ K: when n increases by a factor of 16, ${t}_{\mathrm{multi}}$/${t}_{{\rm{c}},0}$ only decreases very mildly from 3.2 to 2.6. We will discuss the possible reason for the delay of cooling in the next section.

Zoom In Zoom Out Reset image size

The fiducial run assumes a rough balance of heating and cooling, i.e., H/C = 1.02. How robust is the formation of the cool phase to an increase in H/C? To investigate that, we vary the SN rate for n0.02-T3e6 and n0.32-T1e7. We vary H/C from 0.8 to 1.8. Figure 8 shows the time when the multiphase gas forms, ${t}_{\mathrm{multi}}$ as a function of H/C. For n0.02-T3e6, the multiphase is persistent when H/C is $\leqslant 1.4$. The time when the cool phase first occurs is delayed as H/C increases. When H/C = 1.4, ${t}_{\mathrm{multi}}$ is almost twice the value for H/C = 1.02. When H/C = 1.8, the cool phase does not occur within the duration of simulation, which is 4${t}_{{\rm{c}},0}$. For n0.32-T1e7, a similar trend with H/C is seen. When H/C = 1.4, ${t}_{\mathrm{multi}}$ > 4${t}_{{\rm{c}},0}$. Note that even when H/C = 0.8, i.e., the cooling rate is higher than the heating, the cool phase still does not develop until after 2${t}_{{\rm{c}},0}$.

Zoom In Zoom Out Reset image size

The correlation between ${t}_{{\rm{c}},0}$/${t}_{\mathrm{multi}}$/ and H/C suggests that forming the cool phase is more difficult with more intense heating, which makes intuitive sense. However, from Equation (8), the onset of the cool phase should not depend on H/C. We will discuss the delay of cooling in more detail in Section 4.4.

We found in the previous section that the formation time of the cool phase is delayed from the idealized calculation of Equation (8). The delay is longer when the medium has a higher temperature. In this section, we will discuss turbulent mixing as an important cause of this delay.

Feedback from SNe not only heats the gas but also induces motions. In a static, uniform medium, SN explosions will drive spherical, outward shocks. Due to the pressure of the ambient medium, the blast waves will decay into sound waves. As the blast and/or sound waves from explosions at different locations interact with each other, the medium becomes turbulent. Turbulence tends to mix material and energy, therefore transporting heat from SN bubbles to the ambient medium. Gas that is not directly heated by SN-driven blast waves can be heated by this turbulent heat transport. With this extra heating, it takes longer for the gas to cool down. Indeed, turbulent mixing has been suggested as a way to suppress thermal instability (e.g., Banerjee & Sharma 2014). In our cases, the cool phase does form eventually, but we argue that mixing postpones the onset of it.

To quantify turbulent mixing, we use "color" to trace where SN ejecta goes, which is a proxy for the region that is covered by SN bubbles.⁷ If there is no mixing at all, and the domain of SN heating is approximated by a sphere with a radius ${R}_{\mathrm{fade}}$, then according to Equation (6), the volume fraction that is not heated by any SN bubble, ${f}_{\mathrm{unh}}$, decays exponentially with time.

With turbulent motions, the actual ${f}_{\mathrm{unh}}$ declines faster than what is predicted by Equation (6). To find gas that is not heated by any SN bubbles, we have chosen a sufficiently small value of SN color density as a cutoff,

$$\begin{eqnarray}&&{\rho }_{\mathrm{color}}=\kappa {m}_{\mathrm{color}}/{V}_{\mathrm{SN}}=\kappa {m}_{\mathrm{color}}P/{E}_{\mathrm{SN}},\end{eqnarray} \tag{ 9 }$$

where κ = 4.4 × 10⁻⁴. This chosen value of κ is sufficiently small that the results are insensitive to the precise value used. In Figure 9, the dashed line shows ${f}_{\mathrm{unh}}$ as a function of time according to Equation (6), while the solid lines show the actual ${f}_{\mathrm{unh}}$, for the two example runs, n0.02-T3e6 and n0.02-T3e7. The solid lines are well below the dashed line. Furthermore, the way that ${f}_{\mathrm{unh}}$ declines deviates from a simple power law and declines much faster at later time. We use the timescale t_d to empirically quantify the decline of ${f}_{\mathrm{unh}}$, which we take to be the time when ${f}_{\mathrm{unh}}$ equals 1/e.

Zoom In Zoom Out Reset image size

For n0.02-T3e6, t_d = 20 Myr, and ${t}_{{\rm{c}},0}$/t_d ≈ 4.5; for n0.02-T1e7, t_d = 61 Myr, and ${t}_{{\rm{c}},0}$/t_d ≈ 9. The ratio ${t}_{{\rm{c}},0}$/t_d indicates the effectiveness of mixing as a way to spread SN heating, relative to radiative cooling. Larger ratios indicate more efficient mixing, leading to a longer delay in the development of the cool phase.

Figure 10 shows ${t}_{{\rm{c}},0}$/t_d for all simulations with H/C = 1.02. The value of ${t}_{{\rm{c}},0}$/t_d ranges from 4 to 12, indicating that t_d is much smaller than ${t}_{{\rm{c}},0}$ and that turbulent mixing is generally important. The ratio increases with increasing temperature and decreasing density, but it is more sensitive to the former. For example, when increasing the temperature by a factor of 10 for the runs with n0.08, from 3 × 10⁶ K to 3 × 10⁷ K, the ratio changes from 3.7 to 10.6. When increasing the density by a factor of 16 for the run with T = 10⁷ K, the ratio only decreases mildly from 8.8 to 5.5. This trend is very similar to what is found in Figure 7. The similarity of the patterns in Figures 7 and 10 indicates that turbulent mixing is effective in delaying formation of the cool phase for a higher-temperature medium. The increase of ${t}_{{\rm{c}},0}$/t_d is also seen when H/C is larger, confirming that mixing contributes to the delayed cooling shown in Figure 8.

Zoom In Zoom Out Reset image size

In this paper, we have established the relation between the delayed cooling and the turbulent mixing by measuring t_d empirically. In Paper II, where we investigate the turbulence structure of hot ISM in detail, we will present a theoretical model for t_d based on simple physical arguments.

As discussed in Section 4.4, turbulent mixing can transport heat to gas that is not directly heated by SN bubbles, therefore smoothing the uneven heating by random SNe. In our idealized simulation setup, the only driver of turbulence is SNe themselves. In real galaxies, other processes can also contribute to the turbulence, such as galaxy mergers (Paul 2009; Iapichino et al. 2017) and AGN feedback (Gaspari et al. 2012a; Valentini & Brighenti 2015; Wang et al. 2019). These processes may contribute significantly to mixing depending on the magnitude of the turbulence they drive. Unfortunately, turbulent velocities and the relevant driving scales are poorly constrained by current observations of the hot gas in elliptical galaxies (Ogorzalek et al. 2017). In Paper II, we will discuss the SN-driven turbulence in these systems in more detail and compare them to the theoretical estimates of the turbulence driven by other processes.

Thermal conduction is another mechanism of heat transport (Spitzer 1956). Conduction may thus delay the formation of the cool phase. After the formation of the cool phase, depending on the size of cool clumps, thermal conduction can result in evaporation of the clouds or condensation of hot gas onto the cloud (Cowie & McKee 1977). However, our resolution is far coarser than the Field length; therefore, thermal conduction cannot be self-consistently modeled. Additionally, the cool clumps are as small as a few cells, thus not well resolved in these simulations. So the amount of cool gas once multiphase gas appears may not be accurate. Dedicated simulations resolving the conductive layer and cool clumps are needed to further study this problem (e.g., Borkowski et al. 1990; Ferrara & Shchekinov 1993; Brüggen & Scannapieco 2016; Liang & Remming 2020).

Our simulations only model a patch of the ISM in elliptical galaxies. Gravitational fields and gas stratification are not included. In a realistic environment with gravitational stratification, SN bubbles have a lower density than the ambient medium and tend to rise because of the buoyancy force; similarly, the unheated gas has a relatively higher density and tends to sink (Mathews 1990). Figure 11 shows the phase diagram of gas pressure and specific entropy for n0.02-T3e6. The snapshot is taken at t = 50 Myr, which is about 0.56${t}_{{\rm{c}},0}$, when the cool phase has not formed yet. The initial condition is indicated by the star. The pressure has a very narrow range within a factor of a few. But the specific entropy already spans a factor of >30. The range is similar to the observed range for radii over two logarithmic scales: R = 0.1–10 kpc (see Figure 1 of Voit et al. 2015). When stratification is present, a gas parcel that has a different specific entropy from the environment will move in the gravitational potential, until its entropy matches that of the ambient medium (if not considering mixing and cooling/heating).

Zoom In Zoom Out Reset image size

The effect of buoyancy damping on thermal instability in a stratified medium has been studied extensively in recent years (see Voit et al. 2017, and references therein). Due to the lack of gravity in our simulations, the conditions that we find for thermal instability to develop can be considered necessary but not sufficient in real galaxies with gravity. That is, if buoyancy damping is effective (i.e., when the cooling timescale is significantly longer than the dynamical timescale), the instability that can grow in our idealized box will likely be damped via buoyancy oscillation in a stratified medium.

One implication of the stratification is the metal transport. These high-entropy SN bubbles are metal enriched, and if mixing is not efficient, the rising of SN bubbles may preferentially distribute metals outside the central regions of the galaxy (Tang et al. 2009; Voit et al. 2015), and even to the circumgalactic medium (Zahedy et al. 2019). In addition, by the time we stopped the simulations, the multiphase ISM is still evolving and does not reach a steady state. Therefore, the quantitative properties of the simulated ISM may not be representative of real systems. This also poses an obstacle for comparing the simulations directly to the observations. These effects cannot be captured in our simulations with periodic boundary conditions but would be very interesting to investigate in future studies that include gravity and gas stratification.

A cool phase can form owing to the random heating of SNe Ia. The cool clumps, once formed, will fall toward the center of the galaxy/galaxy cluster. Cool gas produced by this mechanism may be an important source for the observed multiphase systems. The cool gas could then form stars and/or be accreted to the central supermassive black hole and fuel the AGNs (Ciotti & Ostriker 2007; Choi et al. 2012; Li et al. 2015b; Yuan et al. 2018). These processes consume the cool phase and in turn inject energy back into the system.

Since SNe Ia heat and metal-enrich the medium at the same time, the chemical abundance of the cool phase would be smaller than that of the hot medium, since the cool phase develops as a result of being missed by any SNe. The mean abundance difference between gas covered by SN bubbles and the cool phase is ${\rm{\Delta }}Z={M}_{\mathrm{met},\mathrm{Ia}}/M({R}_{\mathrm{fade}})$, the ratio between the metal mass ejected by an SN Ia and the enclosed mass within ${R}_{\mathrm{fade}}$. Assuming that the cool phase has a metallicity (metal mass fraction) of Z₀ and using Equation (3) for ${R}_{\mathrm{fade}}$, the percentage difference of metallicity is thus

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{\Delta }}Z}{{Z}_{0}} & = & \displaystyle \frac{(\gamma -1){E}_{\mathrm{SN}}{M}_{\mathrm{met},\mathrm{Ia}}}{{c}_{s}^{2}{Z}_{0}}\\ & = & 0.27\left(\displaystyle \frac{T}{{10}^{7}\,{\rm{K}}}\right)\left(\displaystyle \frac{0.01}{{Z}_{0}}\right)\left(\displaystyle \frac{{M}_{\mathrm{met},\mathrm{Ia}}}{1{M}_{\odot }}\right),\end{array}\end{eqnarray} \tag{ 10 }$$

where E_SN = 10⁵¹ erg is assumed for the second equality. This equation can apply to either the metallicity including all elements or a certain element (such as iron) when Z and ${M}_{\mathrm{met},\mathrm{Ia}}$ are replaced by the corresponding quantities for that element. Equation (10) emphasizes the abundance difference between the hot and the cool and indicates that the difference is more prominent when the ISM is hotter and/or has a smaller Z₀. In addition, the cool phase forms at different epochs or different locations of galaxies and would have different chemical abundances. This at least partly explains the observed diversity of metallicity of cool gas in early-type galaxies (e.g., Sarzi et al. 2006).

In current cosmological simulations, due to the coarse resolution, SN feedback can only be included as a subgrid model. As a result, the above mechanism for forming a cool phase due to uneven heating is not captured. The modeling of feedback from SNe Ia in these simulations is quite simple, usually as a heating term. The amount of heating for each volume (or mass) element in one time step is just the total energy from SNe that would have exploded within that time step, which is related to the stellar density and age. If the added SN energy is larger than the energy loss due to radiative cooling (n²ΛΔt), then the thermal energy of this element will increase and will not form any cool gas. Therefore, cosmological simulations miss an important way to produce the cool phase. Furthermore, due to the large inhomogeneity of the hot gas (Figures 3 and 6), the actual cooling rate of the medium deviates from ${n}^{2}{\rm{\Lambda }}(T)$, where n and T are the averaged quantities (averaged over the large length and math scales resolved in such simulations). This effect is investigated in more detail in Paper II.

To capture the real impact of SN feedback, the resolution requirement is that ${R}_{\mathrm{fade}}$ be resolved by at least a few computational cells. Note that for a typical warm/cool ISM in spiral galaxies, one needs to resolve the cooling radius R_cool to capture SN feedback (Simpson et al. 2015; Kim & Ostriker 2015; Li et al. 2015a). But in a hot ISM in elliptical galaxies, ${R}_{\mathrm{fade}}$ ≪ R_cool (Appendix C); thus, a general condition for the resolution would be resolving min(${R}_{\mathrm{fade}}$, R_cool) by 3–10 cells. For typical conditions in elliptical galaxies, ${R}_{\mathrm{fade}}$ ranges from 20 to 150 pc (Table 1), meaning that the spatial resolution has to be a fraction of that, which is beyond current computing capabilities. That said, doing so in smaller systems like galactic bulges and dwarf elliptical systems is possible (e.g., Tang et al. 2009); and thanks to the rapid development of computing facilities, simulating the whole massive elliptical system with resolved SN Ia feedback would be possible within a few years.