Multiphase Gas in Elliptical Galaxies: The Role of Type Ia Supernovae

Massive elliptical galaxies harbor large amounts of hot gas (T ≳ 106 K) in their interstellar medium (ISM) but are typically quiescent in star formation. The jets of active galactic nuclei (AGNs) and Type Ia supernovae (SNe Ia) inject energy into the ISM, which offsets its radiative losses and keeps it hot. SNe Ia deposit their energy locally within the galaxy compared to the larger few ×10 kiloparsec-scale AGN jets. In this study, we perform high-resolution (5123) hydrodynamic simulations of a local (1 kpc3) density-stratified patch of the ISM of massive galaxies. We include radiative cooling and shell-averaged volume heating, as well as randomly exploding SN Ia. We study the effect of different fractions of supernova (SN) heating (with respect to the net cooling rate), different initial ISM density/entropy (which controls the growth time t ti of the thermal instability), and different degrees of stratification (which affect the freefall time t ff). We find that SNe Ia drive predominantly compressive turbulence in the ISM with a velocity dispersion of σ v up to 40 km s−1 and logarithmic density dispersion of σ s ∼ 0.2–0.4. These fluctuations trigger multiphase condensation in regions of the ISM, where min(tti)/tff≲0.6exp(6σs) , in agreement with theoretical expectations that large density fluctuations efficiently trigger multiphase gas formation. Since the SN Ia rate is not self-adjusting, when the net cooling drops below the net heating rate, SNe Ia drive a hot wind which sweeps out most of the mass in our local model. Global simulations are required to assess the ultimate fate of this gas.


INTRODUCTION
The interstellar medium (ISM) in massive galaxies is hot, with temperatures T ∼ 10 6 -10 7 K and the dense central regions (number density n ∼ 0.01-1 cm −3 ) can cool in less than 1 Gyr (Werner et al. 2012).However, the ISM in most of these early-type galaxies is not undergoing a cooling flow (Peterson et al. 2003).They show little to no current star formation and their stellar populations are typically≳ 10 Gyrs old, barring a few exceptions (e.g., see Calzadilla et al. 2022).Hence the ISM of these massive ellipticals needs to be heated constantly to keep them in rough thermal balance, see Donahue & Voit (2022) for a recent review.
Heating of the ISM due to active galactic nuclei (AGNs) jets powered by the accretion of matter onto the supermassive black hole (SMBH) is one significant energy source (e.g.Allen et al. 2006;Fabian 2012;Main et al. 2017;Gaspari et al. 2017;Guo et al. 2023).AGN jets are typically found to deposit their energy further out (few×10 kpc) from the central regions of the galaxy (Bîrzan et al. 2004).However, on ∼ kpc scales, the jets are highly collimated and the cavities/bubbles associated with them are not volume-filling.While the mechanical energy input from AGN jets matches the net cooling of the halo gas in several systems (McNamara & Nulsen 2007;Rafferty et al. 2008;Olivares et al. 2023), the distribution of this heating throughout the ISM is still an open question.
Heating due to supernovae is another attractive model.In early-type galaxies with old stellar populations, type Ia supernovae (SNIa) dominate the stellar energy output channel and are a natural outcome of stellar evolution.The spatial distribution of the SNIa is expected to follow the stellar distribution in the galaxy.Hence, these supernovae deposit their energy more uniformly and within the ISM of the galaxy compared to the AGN jets.Studies such as Voit et al. (2015) have shown that the energy injected by the SNIa exceeds the energy lost due to radiative cooling in high-entropy outer regions of the elliptical galaxies.They propose that the mass deposited into the ISM by AGB stars can be swept out by an SNIa-driven wind, helping to keep the ISM density low and its cooling time (t cool ) long.Li et al. (2018) find that SN heating can suppress star formation in low mass ellipticals, while AGN feedback is more efficient for systems with stellar mass≳ 10 11 M ⊙ (also see Ciotti et al. 1991;Conroy et al. 2015;Voit et al. 2020;Molero et al. 2023).
Since the net energy output of SNIa is lower than that of the AGN outbursts, many individual galaxy/galaxy cluster-scale simulations have ignored their effect.Cosmological simulations have generally included the mass and energy injection due to supernovae using sub-grid models (for eg.Crain et al. 2015;Pillepich et al. 2018).Alternatively, in simulations that do include discrete SNe (eg.Su et al. 2019) it is not clear that they have the resolution to resolve individual supernova remnants of size a few×10 pc in the low-density ISM of massive ellipticals.While the energy injection due to supernovae is accounted for by such models, they cannot account for the small-scale turbulence driven by the remnants (on scales of 50-100 pc) and the density fluctuations they generate.These density fluctuations can trigger multiphase condensation in the halo gas if the ratio of the thermal instability growth time to the free-fall time (t ti /t ff , Sharma et al. 2012;Choudhury et al. 2019) or the ratio of t ti to the mixing time (t ti /t mix , Gaspari et al. 2018;Mohapatra & Sharma 2019) falls below a critical density dispersion amplitude (σ s )-dependent ratio (Voit 2021;Mohapatra et al. 2023).Observational studies such as Olivares et al. (2019) report multiphase filamentary structures around the brightest cluster galaxies (BCGs) in low-entropy and short cooling time regions of clusters.
Numerical studies on small scales (≲ 10 kpc regions) can resolve individual remnants, either by implementing mesh-refinement techniques around the supernovae injection regions or modeling even smaller ≲ kpc-scale regions of the ISM with uniform resolution.For example, Tang et al. (2009a,b) have looked at the role played by SNIa in the Milky Way's bulge.They find that the SNIa power a galactic bulge wind with a filamentary structure.This wind sweeps out mass from the inner regions and keeps the circumgalactic medium (CGM) hot.
In the context of massive elliptical galaxies, Li et al. (2020a,b) (hereafter MLi20a and MLi20b, respectively) have modeled small patches (of size ∼ few×100 pc) of the ISM at different distances from the galactic cen-ter.They show that the non-uniform heating by supernovae produces multiphase gas even when the net heating rate moderately exceeds the net cooling rate.They also report that the heating by SNIa is different from a uniformly injected volume-heating-the localized heating caused by SNIa drives compressive turbulence and produces large density fluctuations in the ISM.However, MLi20a; MLi20b have not included an external gravitational field in their setup.In a stably stratified atmosphere, buoyancy effects can suppress the formation of cold gas if t ti /t ff is large enough.The supernova-heated bubbles would rise against gravity and generate convective instabilities (such as the Rayleigh-Taylor instability) which can generate solenoidal turbulence and increase the mixing rate of the supernovae ejecta with the ambient ISM.Further, an overheated atmosphere would drive an out-flowing wind (as reported by Tang et al. 2009b, in the context of the Milky Way's bulge).Such a wind can transfer energy, mass and metals to outer regions of the galaxy.Thus an external gravitational field is expected to play a critical role in understanding the impact of Type Ia SNe on the ISM of massive galaxies.
A parallel set of calculations has studied the onset of thermal instability in stratified plasmas and the impact of ambient turbulence (and/or seed perturbations) on this process (e.g.Parrish et al. 2010;Choudhury & Sharma 2016;Voit 2018;Choudhury et al. 2019;Mohapatra et al. 2023).This is important for understanding the origin of multiphase gas in massive galaxies, groups and clusters, and its role in fueling star formation and black hole growth.Such calculations typically either seed initial density fluctuations of a fixed amplitude or drive turbulence on large scales.Li et al. (2020b) have shown that the properties of turbulence driven by small-scale SNIa remnants are quite different from the large-scale forcing typically employed in turbulence simulations.The remnants typically drive turbulence on sub-100 pc scales with a large compressive to solenoidal ratio and seed much larger density perturbations than expected from σ s -M (rms Mach number) scaling relations in forced turbulence (Federrath et al. 2008(Federrath et al. , 2010;;Konstandin et al. 2012).These differences in driving can affect the mixing rate between the hot remnants and the ambient ISM, as well as the occurrence of multiphase condensation for a given set of turbulence parameters.
In this work, we aim to study the properties of SNIa energy injection in a density-stratified, kpc-scale patch model of the ISM.We conduct hydrodynamic simulations including randomly injected supernovae, radiative cooling and shell-by-shell volume heating.In our fiducial set, we study the effect of increasing supernova injection rate on multiphase condensation, properties of turbulence in the ISM and the development of an SNIapowered wind.We also check the dependence of our results on the ambient stratification (which control t ff ), different initial gas density (which affect t ti ) as well as the effect of including mass ejected from AGB stars.We note from the outset that our study in this paper is restricted to local simulations that model a patch of the ISM/ICM.These calculations are analogous to the extensive literature on simulations of the impact of corecollapse SNe on the ISM of disk galaxies (e.g.Kim & Ostriker 2015;Martizzi et al. 2015;Fielding et al. 2017).One important difference is that the local approximation is much less well motivated for the Type I/elliptical galaxy application because of the absence of a thin disk.Given, however, the extensive literature on local simulations in studying the onset of thermal instability it is valuable to assess the role of Type Ia supernovae in that context before proceeding to global simulations.
In the next section, we introduce some theoretical background and estimate the different relevant energy sources and sinks.We state our model and numerical methods in Section 3. In Section 4, we present the results from our fiducial set of simulations, where we progressively increase the SNIa heating rate keeping all other parameters constant.We study the effect of other simulation parameters and summarize the astrophysical implications of this work in Section 5. We discuss our caveats and prospects in Section 6 and finally conclude in Section 7.

Analytic estimates
In this section, we begin by providing some estimates of the different energy sources and sinks for the typical massive elliptical galaxies studied in Werner et al. (2012Werner et al. ( , 2014)).

Radial profiles of stars, gas density and entropy
We first obtain the radial profiles of gas number density n, entropy S from Fig. 6 of Werner et al. (2012).We normalize the gas number density n = 0.07 cm −3 and entropy S 0 = 6.8 keV cm 2 at r = 2 kpc, in agreement with the X-ray observations.The entropy profile flattens in the center and follows a power-law scaling ∝ r 1 at large r.In the upper panel of Figure 1, we show the radial profiles of gas number density n, stellar mass density ρ * and entropy S.
For calculating the stellar mass distribution, we assume a Hernquist profile (Hernquist 1990) with m * = 2×10 11 M ⊙ (see Table 1 in Merritt & Ferrarese 2001, for typical Bulge masses of these galaxies) and r * = 2 kpc (r * ≈ r eff /1.8, we obtain r eff using mass-size relations from Fig. 1 of Trujillo et al. 2011).

Cooling rate of the hot ISM
The ISM of an early-type galaxy is generally hot, with temperatures around 10 6 -10 7 K.It cools through free-free Bremsstrahlung emission and far UV/X-ray line-emission.In order to calculate the cooling rate, we set the gas metallicity Z = 0.5Z ⊙ and generate a temperature-dependent cooling table Λ(T ) using Grackle (Smith et al. 2017).The cooling rate density L is given by: We show L as a function of r in the lower panel of Figure 1.

Heating due to SNIa
For a stellar population of age ∼ 10 Gyr, the SNIa energy injection rate is given by ĖSN ≃ 10 51 ergs × 300(Myr) −1 (M * /10 10 M ⊙ ), (1b) assuming that each supernova injects E SN = 10 51 ergs of energy and that SNIa go off at a rate 300(Myr) −1 (M * /10 10 M ⊙ ), consistent with the delaytime distribution of SNIa in Maoz & Graur (2017) (also see Scannapieco & Bildsten 2005;Barkhudaryan et al. 2019).We show ėSN ( ĖSN per unit volume) in the lower panel of Figure 1.While the exact values of the cooling rate and the SN heating rate will vary across different systems, the comparison in Figure 1 motivates that the SNe heating rate is likely within an order of magnitude of the cooling rate density within the central 10 kpc of the ISM in massive galaxies.This emphasizes the potential importance of Type Ia SNe as a heating source that can compensate for cooling losses and affect the development of thermal instability.This is one of the main motivations for studying the effect of SNIa heating in local patches of the ISM in this paper.

Heating due to thermalisation of stellar winds
Stars on the AGB eject their outer envelopes into the ISM, at a rate assuming a stellar age ∼ 10 Gyrs and using models of Conroy et al. (2009).Since the stars typically have a velocity dispersion σ * ∼ 200-300 kms −1 with respect to the ISM, the wind material can become thermalized with the ISM and contribute to its heating. 1 The heating rate density due to the stellar winds is given by: We show the radial profile of ė * in the lower panel of Figure 1, assuming σ * = 250 kms −1 .The stellar heating rate density is roughly an order of magnitude smaller than the cooling rate density and the supernova heating rate density.However, it could be more important in compact galaxies with large σ * , as discussed in Conroy et al. (2015).

Condition for absence of shell formation
In the standard picture of a supernova going off in the ISM (see chapter 39 of Draine 2011), the evolution of the remnant can be divided into four stages-(1) free-expansion, (2) Sedov-Taylor phase (energy conservation), (3) Snowplow phase (where a dense cooling shell forms from the swept-up material) and (4) Fadeaway (where the shock wave fades into a sound wave in pressure equilibrium with the ambient medium).In the early-type galaxies considered in this study, the ISM is hotter (temperature T ≳ 10 6 K) and has a lower density (n ≲ 0.1 cm −3 ) compared to typical conditions in a galactic disk.Under such conditions, the shock wave may reach pressure equilibrium before the snowplow phase, i.e., before the onset of significant cooling losses (MLi20a).
The fade radius R fade , at which the shock-wave evolves into a sound wave is given by: where P and T are the pressure and temperature of the ambient ISM, respectively and α is a free parameter of order unity.The fade radius is less than the radius at end of the Sedov-Taylor phase (when strong cooling commences) when the ambient density n falls below a critical number density n crit , which is given by MLi20a: The strongest dependence in equation 2b is on the temperature of the ambient medium.Higher temperatures -when the temperature is above the peak of the atomic cooling curve -correspond to less efficient cooling and higher ambient pressure, both of which lead to SNR reaching pressure equilibrium prior to significant cooling losses.The condition n < n crit is satisfied for most regions of the ISM in the early-type galaxies that we study.

Model equations
We use Euler equations to model the ISM as an ideal gas with an adiabatic index γ = 5/3.We evolve the following equations: where ρ is the gas mass density, v is the velocity, P = ρk B T /(µm p ), where µ ≈ 0.6m p is the mean particle weight, m p is the proton mass and k B is the Boltzmann constant.In the continuity equation (3a), n SN denotes the rate of supernovae injection per unit volume and M SN denotes the mass injected by each supernova.In addition to the mass injected by supernovae, we also include the contribution of AGB winds, denoted by ρ * (only in a subset of our simulations).We denote the acceleration due to gravity as g in the momentum equation (3b).In the energy equation ( 3c), the total energy density is given by E and the gravitational potential is denoted by Φ, with g = −∇Φ.We include L, ėSN and ė * defined by eqs.(1a) to (1d), as well as a shell-by-shell heating rate density ėshell .

Numerical methods
We perform hydrodynamic simulations using a performance portable version of the Athena++ (Stone et al. 2020) code implemented using the Kokkos library (Trott et al. 2021).We use the second-order Runge-Kutta (RK2) time integrator, the Harten-Lax-van Leer-Contact (HLLC) Riemann solver and the piece-wise linear spatial reconstruction method (PLM).For cells with unphysically large velocities or temperatures, we use a first-order flux correction algorithm described in the appendix of Lemaster & Stone (2009).

Domain size and decomposition
We simulate a cuboidal patch of the ISM with dimensions L × L × 1.5L using a Cartesian grid, where L = 1 kpc and the longer dimension is oriented along the z-direction.By default, we divide the simulation domain into a grid of size 512 × 512 × 768, such that each resolution element is cubical in size (i.e.dx = dy = dz).Our simulation box is centered at the origin (0, 0, 0).

Boundary conditions
Similar to Mohapatra et al. (2023) (hereafter M23), we use periodic boundary conditions along the x and y directions.The proper boundary conditions in the vertical direction for this problem are non-trivial in that if a SNe driven outflow develops, one would like boundary conditions that allow such an outflow (and allow the density and pressure at the boundary to change).If an outflow does not develop, the boundary conditions should maintain hydrostatic equilibrium with densities and pressures similar to the initial condition.After experimentation, we settled on a compromise in which we use constant boundary conditions for ρ and P at the z boundaries.At the upper z boundary, we implement diode boundary conditions for the z-direction velocity v z .At the lower boundary, we set v z = 0 for gas with T ≳ T floor and diode otherwise.Setting v z = 0 at the lower boundary prevents the atmosphere from undergoing a global collapse and the second condition prevents cold infalling gas from gathering at the negative z-direction.
To minimize the impact of the boundary conditions on our results, we treat the regions beyond |z| > 0.5L as boundary regions and do not include them while analyzing the results of our simulations.

Initial density and pressure profiles
We set up a gravitationally stratified atmosphere with a constant g oriented along the − ẑ direction.At time t = 0, pressure and density follow exponentially decreasing profiles along the + ẑ direction and the gas is at hydrostatic equilibrium.The initial z-profiles are given by: H is the scale height of pressure/density and P 0 , ρ 0 (= P 0 /gH) are the initial values of pressure and density at z = 0, respectively.

Supernova injection
The supernovae heating rate density ėSN is given by: ėSN = ṅSN E SN (5a) In our local ISM patch simulations, we set ṅSN = ṅ0 exp(−z/H SN ), such that the net heating rate due to supernovae and f SN is a parameter that sets the normalization of the SN heating rate.In the above set of equations, H SN denotes the scale height of the stellar/supernova distribution, L 0 is the cooling rate at t = 0 and f SN is a parameter that we vary across simulations.We set H SN = 0.5H so that the heating rate by supernovae has the same variation with z as L (L ∝ ρ 2 ∝ exp(−z/0.5H)),motivated in part by the observations in Figure 1.As we shall show, even with an initial condition balancing heating and cooling both globally and as a function of z, the simulations with significant Type Ia heating tend to rapidly blow a significant amount of gas out of the box.We tried other values of H SN /H as well and found that in those cases the transition to rapid blowout happened even more quickly.We fix f SN at t = 0, which sets ṄSN .The expected number of supernovae in a time step dt is given by ṄSN dt.We draw the total number of supernovae in a given time-step from a Poisson distribution with this mean.We obtain the x and y coordinates of the center of the remnant from a uniform distribution and use an exponential distribution for the z coordinate.
The SNIa remnants in the ISM of an elliptical galaxy are not expected to undergo a snowplow phase (see sec.2.2) and remain in the energy-conserving phase until they fade away into sound waves.Once the center of an individual remnant is determined, we inject E SN = 10 51 ergs of energy as thermal energy and M SN = M ⊙ of mass inside a sphere of radius 0.6 × R fade .Since our simulations have a spatial resolution of ∼ 2 pc, each remnant is resolved by ≳ 10 cells.

Heating due to AGB winds
In a few of our simulations, we have included the effect of mass and energy added to the ISM due to ejecta from AGB stars.We assume that the AGB-wind material is thermalized with the ambient ISM.After setting ṄSN using eq.(5c), we obtain Ṁ * and ė * using eqs.(1b) to (1d) and σ * = 250 kms −1 .We inject mass and energy into the ISM with the same z-dependence as the supernova injection (∝ exp(−z/H SN )).Although we do not include the effect of AGB winds in our fiducial set, we discuss their impact in Appendix A.

Modifications to the cooling function
To control the code-time step set by the cooling function, we reduce L to zero for T < T cutoff .The modified cooling function is given by: where we set T floor = 3 × 10 4 K.

Thermal heating rate and shell-by-shell energy balance
To study the impact of different fractions of supernova heating with respect to the net cooling, as well as to prevent a runaway cooling flow, we implement an additional heating source ėshell which replenishes (1 − f SN ) fraction of the energy lost due to cooling of hot gas (T ≥ T hot = 10 5.5 K) in every z-shell at each timestep.This method is a useful toy model that is agnostic to the origin of other sources of heating (such as AGN, radiation, cosmic rays, conduction, etc.) and has been used in many previous studies, such as Sharma et al. (2012); Choudhury et al. (2019) and M23.Despite its simple prescription, it leads to similar results as feedback heating (Gaspari et al. 2012;Prasad et al. 2015).Mathematically, this is given by: ėshell (5f)

Important time-scales
The different relevant time-scales of our setup are- The Brunt-Väisälä oscillation time: the turbulent mixing time: the gravitational free-fall time: the cooling time: the thermal instability growth time: In the above set of equations H S is the scale height of entropy, k is the wave number, E sol (k) is the power spectrum of the divergence-free (or solenoidal) component of the velocity field, v int,sol is the solenoidal component of velocity on the integral scale ℓ int,sol , α heat is defined such that ėshell ∝ ρ α heat , and n bub and v bub are the number density and the velocity of the bubble inflated by a supernova remnant, respectively.

Initial conditions
We conduct a variety of simulations to model small patches of the ISM of an elliptical galaxy at different locations away from the galactic center.For all our simulations, we initialize the gas at a constant initial temperature T 0 = 3×10 6 K, which is consistent with typical ISM temperatures in elliptical galaxies (e.g., see Fig. 1 in Voit et al. 2015).For our fiducial set of simulations, we set n 0 = 0.08 cm −3 (at roughly 2 kpc out in Fig 1), motivated by observations of gas a few kpc from the center of massive galaxies, see for eg.Fig. 1 in Voit et al. (2015).At t = 0, we seed density fluctuations δρ/ρ in the gas, with 20% amplitude and a flat power spectrum (i.e.no scale dependence).These are consistent with measurements of X-ray brightness fluctuations (e.g., see Zhuravleva et al. 2018, albeit on somewhat larger scales compared to our setup).Previous studies such as Choudhury et al. (2019); Voit (2021); M23 have shown that a stratified, thermally unstable medium is expected to become multiphase if the ratio t ti /t ff is becomes smaller than a δρ/ρ-dependent threshold.We choose the density/pressure scale height H = 12 kpc such that the initial t ti /t ff is close to the condensation threshold proposed by M23.Since we expect the supernova-driving to generate further density fluctuations in the ISM (as seen in MLi20b), our choice of initial conditions lets us directly study their impact on multiphase condensation.Our choice of H and T are also consistent with the observed profiles of density and temperature (at a few kpc from the center) in massive galaxies that host extended multiphase filaments (see blue scatter points in Fig. 1 of Voit et al. 2015, for reference).

List of simulations
For this study, we conducted a total of 16 simulations which are listed in Table 1.We mainly study the effect of varying the fraction of supernova heating rate with respect to the net cooling rate.For our fiducial set, we conduct four simulations with f SN = 0.01, 0.1, 0.5 and 0.99, as indicated in their respective simulation labels.In addition to understanding the effects of resolved SNIa heating, our parameter choices are motivated by the inferred ratio between the heating and cooling rates observations of massive elliptical galaxies (Werner et al. 2012(Werner et al. , 2014;;Voit et al. 2015).For example, we find that in massive ellipticals such as NGC5846, NGC5044 which have a flat gas entropy profile, beyond the inner 10 kpc the stellar density (and as a result the SNIa heating rate) drops off much faster with increasing radius compared to the net cooling rate.So heating due to the SNIa could dominate in the core, but other heating sources are expected to contribute in the outskirts.
To check the effect of a faster cooling rate, we repeat the fiducial set with double the initial density (n 0 = 0.16 cm −3 ) and smaller initial seed density fluctuations (δρ/ρ = 0.05).These runs are indicated by 'hdens' in the simulation label.
For all our other parameter choices, we conduct a pair of simulations-one with a small supernova heating fraction (f SN = 0.01) and another with a large supernova heating fraction (f SN = 0.99).The observed galaxies are have different densities, temperatures and strengths of stratification, which also vary across the inner and outer regions of the ISM.The 'ldens' runs are useful to study the effect of a weaker cooling rate due to their smaller initial gas density.We conduct a pair of simulations with smaller density/pressure scale heights-'H6hdens' to look for the effect of stronger stratification on the 'hdens' runs.Our 'H3' runs have initial density profiles that are directly comparable to that of a massive elliptical galaxy (see Fig. 1), noting that we use a slightly cooler T 0 .Finally, the ' ṁ' runs include mass and energy injection due to mass loss from AGB stars, as described in Section 3.3.3(all other runs do not account for AGB winds).
We evolve all the simulations till t end = 500 Myr.We define t mp as the time at which the fraction of gas mass at T < 10 4.2 K exceeds 1% for the first time.For all simulations that have more than 1% of their total mass in the cold phase, we list t mp in column 5 of Table 1, otherwise we denote t mp > 500 Myr.In columns 6-12 we present different statistical properties of the simulations, averaged over the last 15 Myr of their evolution before t mp (or t end if they never have M cold /M tot > 0.01 during the simulation).For simulations where t mp ≤ 30 Myr, we present the statistical properties for 15 Myr ≤ t < t mp , to allow enough time for the system to evolve before we analyze them.
We rerun all 16 simulations listed in Table 1 at half the resolution (using 256 2 ×384 resolution elements) and find all of the simulation properties listed in columns 5-12 to be within a few % of their high-resolution counterparts.

RESULTS -FIDUCIAL RUNS
We start by presenting results from our fiducial runs, where we progressively increase the fraction of the SNIa heating rate with respect to the net cooling rate.We study the effect of this parameter (f SN ) on the distributions of thermodynamic quantities such as density and temperature, time evolution of mass and energy, formation of cold gas (T ≲ 10 4.2 K) and vertical shell-averaged profiles of energy fluxes.As we shall see, an important feature of our results is that when Type Ia SNe heating is important (larger f SN ), there is no local equilibrium between SN heating and cooling possible.SNe eventually drive gas out of the box, decreasing the cooling rate, leading to stronger SNe driven wind, etc.The Notes: Column 1 shows the simulation label.The number following f SN denotes the ratio of the supernova energy injection rate to the net cooling rate at t = 0.The default initial number density n 0 is set to 0.08 cm −3 unless specified as 'hdens' with n 0 = 0.16 cm −3 or 'ldens' with n 0 = 0.04 cm −3 , respectively and are listed in column 2. The runs with ' ṁ * ' in the label include mass and energy injection into the ISM by AGB winds.In column 3, we show the density/pressure scale height H which is set to 12 (= 12 kpc) for most of our runs (H is denoted in the simulation label otherwise).We list ṄSN , the total number of supernovae injected in the simulation per Myr in column 4. In column 5, we show the multiphase gas formation time tmp (> 500 if the gas in the simulation remains in a single phase).We show mass-weighted volume averages of the standard deviation of velocity v in column 6 and dimensionless turbulence properties-the Froude number Fr, the Mach number M and the compressive component of the Mach number Mcomp in columns 7, 8 and 9, respectively.
In columns 10 and 11, we list the ratio of thermal instability growth time-scale t ti to the free-fall time-scale t ff and the turbulent mixing time-scale t mix , respectively.Finally, we show the amplitude of logarithmic density fluctuations in the hot phase in column 12.All averaged quantities/statistics (columns 6-12 ) are calculated for gas with 10 6 K ≤ T ≤ 3 × 10 7 K within the inner cube (|x|, |y|, |z| ≤ 0.5), and averaged over a 15 Myr duration before tmp (or t end if tmp > t end = 500 Myr).
longer-timescale outcome of this inevitable instability cannot be studied in our local approximation and will ultimately require global simulations.

Projection maps
In Figure 2 we show snapshots of the x-projections of the logarithms of gas number density log 10 (n x ) (column 1), mass-weighted temperature log 10 (T x ) (column 2), normalized density log 10 (ρ x ) (column 3) and normalized pressure log 10 ( Px ) (column 4), respectively2 .For the f SN 0.01 and f SN 0.1 runs that do not form multiphase gas (i.e.M cold /M tot < 0.01 at all times), we present the snapshots at t = t end , whereas we show the snapshots at t = t mp for the f SN 0.5 and f SN 0.99 runs.
Among the single phase runs (f SN 0.01 and f SN 0.1), we find that the density of the gas decreases with increasing f SN .The gas is also hotter.Since the SNIa rate is fixed throughout the course of the simulation unlike the shell-by-shell heat prescription (which self-adjusts to (1 − f SN ) fraction of the total cooling rate at each time-step), whenever the net cooling rate falls below the initial cooling rate the system overheats.The amount of overheating increases with increasing f SN .
Higher SNIa rates are also associated with larger density and pressure fluctuations.This is because the supernovae remnants are small (a few ×10 pc-scale regions) and get overheated compared to the remainder of  the ISM.The overheated regions expand and rise buoyantly (as seen in fig. 3) and drive sound waves in the ISM.As these bubbles interact with the ISM, they form mushroom-shaped clouds characteristic of the Rayleigh-Taylor instability and drive turbulence in the ISM.The amplitudes of density and pressure fluctuations increase with increasing f SN , as seen in the right two panels of Figure 2.
In snapshots of the multiphase runs (f SN 0.5 and f SN 0.99), we observe dense small-scale clouds which correspond to gas at or below the cooling cutoff temperature (T ≲ 10 4.2 K).These clouds are not in hydrostatic equilibrium with the ambient ISM and rain down to the bottom of the box, forming a trail along the direction of gravity.The structure and distribution of cold gas in supernova-driven turbulence are quite different from that of the cold clouds in turbulence driven on large box scales studied in M23 (see their Fig. 1).In our simulations, the clouds are not volume-filling (in contrast to their natural driving runs) and do not form large, boxscale filaments (in contrast to their compressive driving runs).

Time-evolution
Here we present the time evolution of some key properties of our fiducial set of simulations.In rows 1-4 of Figure 4, we show the time-evolution of net mass and energy (normalized to their values at t = 0), the amplitude of velocity dispersion of hot gas (σ v,hot ) and the mass fraction of cold gas (M cold /M tot ), respectively.We  (second row) normalized to their initial values, velocity dispersion (σ v,hot , column 3) of hot gas (T ≳ 10 5.5 K) and mass fraction of cold gas (T ≲ 10 4.2 K, column 4) for our fiducial set of simulations.A higher SNIa rate (larger fSN) leads to a larger σ v,hot .Stronger driving triggers multiphase condensation in the fSN0.5 and fSN0.99 runs.Imbalanced heating due to the SNIa also gives rise to a wind.In steady state, the mass and energy within the box decrease with increasing fSN due to the increasing impact of these two processes.
first discuss the evolution of these quantities at initial times (t < 100 Myr).We find that the net mass decreases (almost monotonically) as a function of time.The mass loss rate increases with increasing f SN .The net energy also shows a similar behavior, although the energy loss rate is much smaller and the decrease is not monotonic.We observe a clear effect of different SNIa rates on σ v,hot , where the runs with larger SNIa rates drive stronger motions in the ISM.This is consistent with the velocity measured in ob-served Hα filaments on 100 pc scales in Li et al. (2020c).However, the cold and hot phases are not co-moving in our simulations.The condensed cold gas exits our simulation domain before it becomes entrained in the hot phase, so it does not trace the velocity structure of the hot phase.
Only the f SN 0.5 and f SN 0.99 runs trigger multiphase condensation3 .The cold gas mass fraction shows an initial increase, followed by a steep decrease.The fraction drops below 10 −4 in roughly t ff ≈ 80 Myr after its peak value as the cold gas rains down through the bottom of the box.
The f SN 0.99 run converts a larger fraction of its mass into the cold phase at an earlier time compared to the f SN 0.5 run.This is because the higher SNIa rate generates stronger density fluctuations where the dense regions cool faster, so the threshold t ti /t ff for the simulation to remain single phase is higher (Choudhury et al. 2019;Voit 2021,M23).
By t = 500 Myr, the net mass and energy reach a rough steady state and show a slow decrease with time.Their values are smaller for larger f SN runs, since they lose more gas due to multiphase condensation as well as to outflows (which we will show in the next section).The motions in the hot phase have roughly the same σ v,hot as initial times.The system remains in a single phase after the initial round of multiphase condensation.

Mass injection by SNIa and boundary fluxes
Type Ia supernovae inject both mass and metals into the ISM.In addition to forming metals during the SNIa explosions, the wind driven by the SNIa can transport these metals to the outer regions of the galaxy/into the CGM/ICM (Tang et al. 2009b).Accounting for the contribution from SNIa is important to explain the observed radial trends of metallicity in galaxy clusters, especially for the observed amounts of Fe and Ni (Gatuzz et al. 2023a,b). 4 Here we discuss the evolution of the mass fluxes and sources.The relevant sources/sinks of mass are the mass outflow rate Ṁout , and the mass injection due to the supernovae.They are defined as:  The mass fluxes at the domain boundaries are an order of magnitude larger than the mass injected by the SNIa, even for the fSN0.99 run.Thus gas inflow/outflow at the boundaries dominates the evolution of net mass.
In Figure 5 we show the evolution of Ṁout at the upper z-boundary Ṁout,upp (row 1) and the lower z-boundary Ṁout,low (row 2), respectively.In the third row, we present the mass injection rate due to the SNIa remnants ṀSN (row 3).All mass fluxes are divided by M 0 , the total mass at t = 0.The mass injected due to the SNIa remnants is dependent on f SN as it sets the SNIa rate.However, its value is an order of magnitude (or more) lower than the mass fluxes at the boundaries and M 0 / ṀSN ≪ 500 Myr.So ṀSN is too small to have a significant impact on the gas density during the course of a simulation.On the other hand, M 0 / Ṁout,upp and M 0 / Ṁout,low range between 80-500 Myr for f SN ≥ 0.5 and thus affect the overall mass distribution significantly during the course of a simulation.
Figure 6.The z-profile of number density of gas of gas with 10 6 K < T < 10 7.5 K at t = 0 Myr, t = min(tmp, 100 Myr) and t = 500 Myr.Increasing the SNIa rate leads to multiphase condensation and a strong outflow at late times.These factors contribute to a drop in n0 as well as a change in its functional form by 500 Myr.
Before the formation of multiphase gas at t = t mp , Ṁout,upp is positive and Ṁout,low is negative for the f SN 0.5 and f SN 0.99 runs, implying that the system is losing mass due to outflows even before multiphase condensation.At t = t mp , we observe a sharp, correlated dip in Ṁout,upp and Ṁout,low for both runs, which indicates cold gas raining down and hot gas flowing in from the upper boundary to replace it.This is also associated with a sharp decrease in the net mass and energy of the box.The amplitude of this drop is larger for the larger f SN run which forms more cold gas.
The f SN 0.01 run does not show strong outflows or inflows at any given time.The evolution of the f SN 0.1 runs is similar to the late-time evolution of the f SN 0.5 and f SN 0.99 runs, where most of the material is flowing in the positive z-direction with a similar outflow rate.

Vertical profile of density
In Figure 6, we show the vertical profiles of the shell-averaged number density ⟨n⟩ for the f SN 0.01 and f SN 0.99 runs from the fiducial set at t = 0, min(t mp , 100 Myr) and t end .At t = 0, the squiggles in the profiles are due to the seed white noise of amplitude σ s = 0.2 that we add.By t = min(t mp , 100 Myr), we find that ⟨n⟩ drops while maintaining a similar dependence with z, as a result of the weak outflows driven at early times.We also note that most of the small-scale perturbations have disappeared, likely due to viscous dissipation.
By t = t end , stronger outflows in f SN 0.99 run lead to a further drop in ⟨n⟩.The number density profile is also steeper.This implies that the SNIa heating drives an outflow, which can affect both the total amount of gas in the ISM as well as its radial distribution.So the gas density profiles found in observations (e.g.Werner et al. 2012Werner et al. , 2014) ) can evolve with time due to heating by the SNIa.

Energy sinks, sources and boundary flux
The heat injected by the SNIa is highly anisotropic on small ∼ 10 pc scales, where overheated remnants are expected to rise buoyantly and drive a convective flow.Further, the SNIa rate does not depend on the net cooling rate unlike the AGN feedback loop, where the cooling and heating are expected to be coupled through a delayed cycle (Prasad et al. 2015;Tremblay et al. 2018).Although the energy injection due to SNIa and the radiative cooling rate are of similar magnitude (see Fig. 1), no local equilibrium is possible: SNe over/underheat the ISM when the net cooling rate increases or decreases.During episodes of weaker cooling, the energy injected by the SNIa powers an outflow, at least on the scale of our local box.
We define the energy outflow rate Ėout and the net energy injection/loss rate densities below 5 : Ėshell = ėshell dV.
In Figure 7 we present the evolution of the net outflowing energy flux ( Ėout,upp − Ėout,low ) (first row), ĖSN (second row), Ėcool (third row) and Ėshell (fourth row).All quantities are normalized by E 0 , the net energy at t = 0.All the fluxes, source and sink terms are significant during the course of a simulation (i.e.E 0 / Ė ≲ 500 Myr).By construction, the energy injection rate due to the SNIa increases with increasing f SN and does not show much variation with time.Before t = t mp , we observe a net energy outflow at the boundaries.Around t = t mp , Ėout is briefly negative for the f SN 0.99 run, possibly because cold gas exits the box from the bottom and hot gas enters from the top at this time.At late times, Ėout is positive and increases with increasing f SN , with an amplitude comparable to the SNIa heating rate.
Since the cooling rate is proportional to the square of the gas density, its evolution follows a similar trajectory as the evolution of net mass in Figure 4, barring the peak in the cooling rate at t = t mp .The shell-by-shell heating rate is (1 − f SN ) fraction of the radiative cooling losses at all times except for t mp ≲ t ≲ t mp + t ff in the f SN 0.5 and f SN 0.99 runs.Because it only replenishes cooling losses for gas with T ≳ 10 5.5 K (see eq. 5f).Hence the absence of a peak in the net Ėshell corresponds to the formation of gas at intermediate and low temperatures (T ≲ 10 5.5 K).
For the f SN 0.01 and f SN 0.1 runs, ĖSN and Ėout are much smaller compared to the cooling and shell-by-shell heating rates at all times.For the two multiphase runs (f SN 0.5 and f SN 0.99), we find that at t ≲ t mp , the cooling rate has a similar amplitude as the SNIa heating rate and is larger than the outflow rate.However, at late times the cooling rate drops significantly due to mass loss and most of the energy injected by the SNIa is accounted for by a larger Ėout .

Vertical profiles of energy fluxes
Here we present the contribution of different flux components to the net energy flux at different z-shells.This is important to understand the mechanism through which supernovae transfer their energy to the ISM as compared to the heat transfer mechanism by AGN (e.The total energy flux Ėout defined in eq.7c can be decomposed into different physical components including wind (advective), convective, and wave as follows (see Parrish et al. 2009, for a similar analysis): Ėwave = dxdy ⟨δP δv z ⟩ , where ⟨⟩ represents the average over a z-shell and the fluctuations such as δn are defined as n − ⟨n⟩ for each shell.We break down this discussion into two  subsections-(1) at early times (t ≲ min(t mp , 100 Myr)), shown in the upper panel of Figure 8 and (2) at t = t end , shown in lower panel.Since the evolution of the zprofiles of the energy fluxes does not vary much across the fiducial set, we present our analysis for the f SN 0.5 run only.For t = t mp (or at t = 100 Myr), we find that the net energy flux is dominated by the contribution of the thermal wind flux.The contribution of the kinetic wind flux is low, which is expected since Ėwind,kin / Ėwind,therm = (γ − 1)M 2 /2 and M ≲ 0.1 for all our simulations (see column 8 in Table 1).The flux due to convection is also generally positive since the supernova remnants are overheated compared to their surroundings and rise due to buoyancy (δT × δv z > 0).Tang & Churazov (2017) show that for an outburst in a uniform medium under spherical symmetry, the fraction of energy carried away by sound waves is ≲ 12%.Bambic & Reynolds (2019) showed that this fraction could be larger if a cocoon of shocked plasma generates small-scale waves, boosting their fraction to ∼ 25%.We do not find much energy in the form of sound waves, except small bumps in their z-profiles which likely corre-spond to recently injected SNIa.On average, the energy flux fraction in sound waves is ≲ 1% of ĖSN .
At t = t end (the lower panel of Figure 8), we find the total energy flux Ėtot is positive at all z-shells for f SN ≥ 0.1, implying an outflow in the positive zdirection.It is completely dominated by the contribution of the thermal wind flux.Since the system loses mass due to outflows or condensation, the gas density decreases as a function of time, hence Ėconv becomes weaker at late times.

Distribution functions
In Figure 9, we present the mass-weighted probability distribution functions (PDFs) of gas number density (column 1), temperature (column 2) and Mach number (column 3) .These help us understand the general characteristics of turbulence driven by supernovae in a low-density ISM.All the PDFs are averaged over 10 Myr before min(t mp , t end ).
The two single-phase runs show a single peak, corresponding to the hot phase.The two multiphase runs show two distinct peaks, corresponding to the hot and cold phases, respectively.Unlike large-scale driven turbulence, the gas density distribution is not lognormal.For the f SN 0.01 and f SN 0.1runs show a low-density tail and a corresponding high-temperature tail in the temperature PDF, which are likely to be associated with gas inside the hot bubbles directly heated by the remnants.The gas motions in the hot phase are subsonic, with the peak M ≲ 0.1 for all runs.The cold phase in the f SN 0.5 and f SN 0.99 runs is mildly supersonic, with M ∼ 2-3.
With increasing f SN , the density and temperature PDFs become broader and the peak M increasesdenoting stronger turbulence driven by a higher SNIa rate.We also find more gas at intermediate temperatures (10 4.2 K ≲ T ≲ 10 5.5 K).The hot, intermediate and cold phases are observed to coexist in filaments of the M87 galaxy, with temperatures ranging from 100 K to 10 7 K (Werner et al. 2013;Anderson & Sunyaev 2018).Since the intermediate temperature gas cools fast and has a short expected lifetime, turbulent mixing with the hot phase is one of the proposed mechanisms to continuously generate them.Heating of the cold regions (at T ≲ 10 4 K) by the SNIa and condensation from the dense regions of the hot phase can also generate gas at these intermediate temperatures.Although the filaments in M87 are likely associated with gas uplifted by the AGN jet, the SNIa can play a role in producing gas at these temperatures through the above processes.Here we discuss the effects of SNIa driving on the density-temperature phase diagram.

Density-temperature phase diagram
The different trend-lines shown in the plots are useful for tracking the nature of the thermodynamic perturbations in the gas.The phase diagrams let us distinguish between different physical mechanisms that are responsible for energy transfer in the ISM -for example, adiabatic fluctuations can be due to compressive turbulence and sound waves, whereas isobaric fluctuations are caused by weak subsonic motions in a stratified medium (Mohapatra et al. 2020).Further, the growth rate of different modes of thermal instability t ti (for example isobaric, isochoric) is expected to be different (Das et al. 2021).We can determine the relevant mode for our simulations from the phase diagram and calculate the correct t ti .
In Figure 10, we show the phase diagram for the fiducial set.For the f SN 0.5 and f SN 0.99 runs, we observe the hot phase at T ≳ 10 6 K and the cold phase below the cutoff temperature at 10 4.2 K.At T ≳ 10 7 K, we find adiabatic fluctuations that are likely to be associated with recently injected supernovae.The perturbations in the intermediate temperature gas are mostly isobaric (10 4.2 K ≲ T ≲ 10 5.5 K).Unlike previous studies such as Mohapatra et al. (2023), we do not observe a sharp isochoric drop at T ∼ 10 5.5 K, which corresponds to the peak of the cooling curve.We note a slight steep- Figure 11.The time-evolution of the minimum value of the ratio tti/t ff plotted against σs for the fiducial set.The dark gray points show the value of this ratio at t = 0. Simulations that do not form multiphase gas are shown as unfilled colored data points at t = t end and simulations that form multiphase gas are shown as filled colored data points at t = tmp.The colored dashed lines show the evolution of this ratio from t = 0 to t = min(tmp, t end ).The dash-dotted trend line is a scaled version of the condensation criterion proposed in M23.Although all four simulations start with similar min(tti)/t ff , their evolution is strongly dependent on the SNIa rate.
ening of the PDF at T ∼ 10 4.5 K, which could be due to the under-resolution of the coldest clumps, for which we do not resolve the cooling length (ℓ cool = c s t cool , see Fielding et al. 2020).
With increasing f SN , we observe more contribution of adiabatic perturbations-which can be associated with compressive turbulence and sound waves directly driven by recently injected supernovae.Arévalo et al. (2016) have analyzed the nature of X-ray brightness fluctuations for M87 using the hardness ratio of X-ray spectra in Chandra observations (also see Zhuravleva et al. 2018).They find the filamentary regions to be isobaric, weakly shocked regions to be adiabatic and the remaining regions to be a combination of the two.Our findings suggest that the subsonic turbulence driven by the SNIa can contribute to either of the two observed modes.

Condition for multiphase gas formation-fiducial runs
For the gas around galaxies and galaxy clusters, the ratio t ti /t ff is an important criterion to determine whether a thermally unstable system forms multiphase gas (Sharma et al. 2012;McCourt et al. 2012).A small t ti /t ff (≲ 10) has been associated with the existence of multiphase gas in both simulations (Prasad et al. 2018) and observations (Olivares et al. 2019).More recent theoretical studies have indicated that the threshold value of t ti /t ff for multiphase gas formation depends on the amplitude of density/entropy fluctuations (Choudhury et al. 2019;Voit 2021;M23).
In Figure 11, we plot the minimum value of the z shell-averaged value of t ti /t ff versus the amplitude of logarithmic density fluctuations for gas with T > 106 K.The dash-dotted black line shows the condensation curve proposed by M23, scaled down 6 by a factor of 0.6.
We observe a clear impact of the role played by the SNIa-driving on the occurrence of multiphase condensation.In the f SN 0.01 and f SN 0.1 runs, the seed density fluctuations at t = 0 are quickly damped by viscous forces and the weak SNIa driving does not generate large density fluctuations.As a result, the simulations no longer satisfy the condensation criterion and remain in a single-phase till t = t end .The stronger SNIa driving in the f SN 0.5 and f SN 0.99 runs increases the amplitude of density fluctuations.However, the larger number of injected SNIa also drive an outflow at initial times (see Figs 5 and 7) and overheat small regions of the gas where the remnants deposit their energy.These processes decrease the gas density and as a result increase min(t ti )/t ff .After accounting for the combined effect of the two, the system still satisfies the condensation criterion and forms multiphase gas.

EFFECT OF VARYING OTHER PARAMETERS AND SUMMARY
In the previous section, we discussed the effects of varying the level of SNIa heating while keeping all other parameters constant.In this section, we vary different simulation parameters such as the mean gas density, the strength of gravity, inclusion of heating due to AGB winds and assess the impact of switching SNIa heating on/off in these systems.These help us understand the local effects of SNIa heating on the ISM in different regions of early-type galaxies.

Scaling of fluctuations with the rms Mach number
Observational studies of the hot gas around giant elliptical galaxies or galaxy clusters often rely on indirect techniques such as using the amplitude of surface brightness fluctuations (which are dependent on the gas density, e.g.Zhuravleva et al. 2014Zhuravleva et al. , 2015Zhuravleva et al. , 2018) ) or fluctuations in the thermal Sunyaev-Zeldovich effect to infer the amplitude of the turbulent gas velocities (which depend on the thermal pressure of the gas, e.g.Khatri & Gaspari 2016;Romero et al. 2023).On the other hand, direct measurements of the turbulent velocities using the high spectral-resolution telescope XRISM 7 can be used to infer the state of perturbations in the hot ISM.In this subsection, we present the scaling relation between the gas density/pressure fluctuations and the turbulent Mach number for SNIa-driven turbulence.
In Figure 12, we show the joint evolution of σ 2 s,hot (row 1) and σ 2 ln( P ),hot (row 2) versus M (column 1) and its compressive component M comp (column 2).The dashed lines in each panel show empirical fits to the data.The dotted line in the top right panel shows the σ s -M comp scaling relation for homogeneous turbulence forced on large scales from Konstandin et al. (2012).
In general, we find that both σ s,hot and σ ln( P ),hot increase with increasing M or M comp .The amplitude of pressure fluctuations is much smaller than that of density fluctuations.The fluctuations are larger for runs with higher f SN , in line with our expectations.The multiphase runs (represented by the filled markers) are 7 https://heasarc.gsfc.nasa.gov/docs/xrism/about/typically associated with larger f SN and thus are found in the top right of each panel, whereas the single-phase regions are typically found in the bottom left, with a few exceptions.Among the runs with the same f SN but different densities, the 'ldens' ('hdens') runs have a lower (higher) cooling rate, hence a weaker (stronger) SNIa driving and smaller (larger) σ s,hot , σ ln( P ),hot compared to the fiducial set.Variations in the other parameters such as 'H' or the inclusion of stellar heating do not have any significant effect on these results.
Both σ s,hot and σ ln( P ),hot follow power-law scaling with M and M comp with some scatter around the empirical fits.We find a tight relation between σ ln( P ),hot and M comp .Similar to MLi20b, we find both fluctuations to be much larger than predicted by scaling relations in literature, such as Konstandin et al. (2012); Mohapatra et al. (2021Mohapatra et al. ( , 2022)).Unlike the proposed scaling relations, which were based on idealized turbulence driven by exciting particular modes in k-space, the supernovae overheat small regions of the gas which expand, rise buoyantly and drive turbulence in the ISM.Heating by the SNIa and radiative cooling of the ISM are also associated with isobaric density perturbations, as seen in fig.10.We expect the differences between such idealized driving and driving by the SNIa to be due to the heating and cooling of the gas, as well as the differences in driving scale, modes, etc.
Compared to MLi20b, who study SNIa driven turbulence in an unstratified box (i.e.without external gravity), we observe some key differences.The amplitude of perturbations that we observe are smaller by a factor of 4 or more compared to that of MLi20b for the same M.This is likely due to the smaller ratio between the power in compressive and solenoidal modes in our simulations, even though the contribution from compressive modes still dominate the total power.In our simulations, we find that the bubbles inflated due to the energy deposited by the SNIa rise buoyantly and form mushroom-shaped clouds, which are a typical characteristic of the Rayleigh-Taylor (RT) instability (see fig. 3).The motion of the hot bubble gas with respect to the ambient medium also generates the Kelvin-Helmholtz (KH) instability.Both the RT and KH instabilities can drive solenoidal turbulence in the ISM and increase its contribution to the total power.
Although the amplitude of density and pressure fluctuations generated by SNIa-driven turbulence is large, they may be difficult difficult to detect since they occur on ≲ 100 pc scales.First of all, small-scale fluctuations would be canceled out due to averaging along the line of sight.Secondly, most current X-ray and microwave telescopes lack the resolution to resolve small scales even for nearby massive elliptical galaxies (R fade ∼ 46 pc ∼ 0.5 arc-second for the M87 galaxy, which is close to the resolution limit of the Chandra telescope).Future X-ray telescopes such as AXIS 8 would be useful to measure these ISM properties with their higher angular resolutions and better sensitivities.The synchrotron emission from the electrons accelerated in the SNIa-shocks could be detected in radio wavelengths.However, it would be challenging to distinguish them from the radio emission from other sources, such as shocks due to AGN activity and mergers.

Condition for multiphase gas formation-all runs
In this subsection, we revisit the criteria for multiphase gas condensation for all our simulations.In addition to min(t ti )/t ff we discussed earlier, we also check the importance of turbulent mixing in suppressing mul-8 https://axis.astro.umd.edu/  13.First The time taken by a simulation form multiphase gas normalized by the initial thermal instability time-scale tmp/avg(tti)t=0 plotted against the amplitude of logarithmic density fluctuations in the hot phase σ s,hot .When the simulation forms multiphase gas, we show the data points as filled markers.When they remain single phase till t = t end , we show the lower-limit of this ratio and mark them as unfilled data points and an upward arrow.Second row: The minimum value of the ratio tti/t ff plotted against σs.The dark grey points show the value of this ratio at t = 0 and the dashed lines show the evolution of this ratio till t = min(tmp, t end ).The dash-dotted curve shows the tti/t ff − σs condensation criterion, which separates the plot into single and multi-phase regions.Third row: Similar to the second row, but showing the ratio tti/tmix instead.Fourth row: Here we show the larger of the two ratios shown in rows 2 and 3.The dot-dashed line shows the condensation criterion proposed by M23.The solid line is a re-scaled version of this criterion which fits our results better.
will let us address the fate of the SNe heated gas.This will also resolve the challenge of boundary condition sensitivity highlighted in the previous paragraph.

HEATING MODEL
We include an additional heating term that replenishes a (1 − f SN ) fraction of the radiative losses in each z-shell.This modification also improved the stability of the atmosphere against a large-scale cooling flow.This overly idealized heating model is motivated by the observational evidence for global thermal stability in massive galaxies but in detail it cannot be correct and the formation of multiphase gas and the long-term evolution of the hot ISM may be sensitive to the details of the correct heating function.
Many massive elliptical galaxies have AGN and the AGN jets can drive large-scale motions in the ISM.The in-fall of satellite galaxies can also add mass and drive turbulence in the ISM.We have not considered their effects in this study.

RESOLUTION
We have performed all simulations listed in Table 1 at two resolutions-512 2 × 768 and 256 2 × 384.We find the results of our single-phase runs are convergent with increasing the resolution.Among the multiphase runs, we find that the results diverge for t ≥ t mp .The lower resolution runs form more cold gas compared to their high-resolution counterparts.This may be due to the effects of excessive averaging at the boundary layers between the hot and cold regions which forms more gas at intermediate temperatures (10 5 K < T < 10 6 K).The intermediate temperature gas cools fast and forms more cold gas in the lower resolution simulations.

PHYSICS
We have ignored the effect of important physics such as magnetic fields (Hopkins et al. 2020;Wang et al. 2021;Buie et al. 2022), cosmic rays (Kempski & Quataert 2020;Butsky et al. 2020;Beckmann et al. 2022) and conduction (Parrish et al. 2009;El-Badry et al. 2019), see Faucher-Giguère & Oh (2023) for a review.Magnetic fields and cosmic rays are expected to be energetically important and can affect the properties of the thermal energy-driven outflow.Conduction can affect the energy exchange between the SNIa inflated bubbles and the ISM.All three can affect the formation of cold gas and its kinematics.In future studies, we plan to include some of these physical properties and study their impact.

CHEMICAL EVOLUTION
We have ignored the evolution of the chemical composition of the ISM due to metals injected into the ISM by the SNIa.SNIa are one of the main sources of elements such as Fe, Co, Ni, etc.The gas cooling rate at T ∼ 10 6 K is sensitive to the chemical composition of the ISM.We plan to include heavy element injection and transport in future works and study their properties such as their radial extent, spatial variations, etc.

CONCLUSION
In this paper we model a 1.5 kpc 3 local stratified patch of the hot ISM of a massive elliptical galaxy.We study the effect of different strengths of heating due to type Ia supernovae (SNIa), motivated by the observational fact that the SNIa heating rate is of order the radiative cooling rate in many massive galaxies (see, e.g., Fig. 1).We fix the heating rate due to the randomly injected SNIa to a fraction f SN of the net cooling rate of the ISM and compare the hot ISM properties against a uniform shell heating model typically used in idealized simulations.We have conducted a total of 16 simulations, where we vary ISM properties such as gas density and the strength of stratification.Here we summarize some of our key findings: • The SNIa deposit their energy in small ∼ 20 pc size regions of the ISM.These regions expand and rise buoyantly, driving turbulence in the ISM.The turbulence is associated with large density fluctuations (Figs. 2 and 3).
• The high-density regions have a short cooling time.
Since the ISM at these temperatures is thermally unstable, if the ratio of the thermal instability growth time to the free fall time (t ti /t ff ) is small enough, the dense regions cool down to the cooling cutoff temperature at 10 4.2 K, they are out of hydrostatic equilibrium and rain down through the bottom of the box.Since the SNIa seed these density fluctuations (also seen in Li et al. 2020a), a larger SNIa injection rate is more likely to trigger multiphase condensation (Figs. 11).Much of the literature has focused on the formation of multiphase gas by thermal instability (for eg.Sharma et al. 2012;McCourt et al. 2012;Choudhury et al. 2019) or AGN driven-uplift of gas (eg.Pulido et al. 2018;Huško & Lacey 2023).Our results show that SNIa are also efficient sources of multiphase gas production in the hot ISM/ICM of massive galaxies, groups, and clusters.A spatially and temporally resolved treatment of SNIa is likely critical for understanding the formation of multiphase gas in massive galaxies and its role in fueling star formation and black hole growth.
• As the SNIa rate is fixed, once the gas density drops due to multiphase condensation, the SNIa overheat the ISM and drive an outflow.The net mass remaining in the simulation decreases with increasing SNIa rate (figs. 4, 5).Even in the simulations that do not form multiphase gas, whenever the net heating exceeds the net cooling, the SNIa drive an outflow which further decreases the net mass and the cooling rate.This evolution reflects the fact that heating by SNIa is inevitably unstable: unlike core-collapse SNe or (plausibly) AGN there is no connection between the SNIa rate and the radiative cooling of the hot gas.Global simulations will be required to understand the ultimate outcome of this instability.
• In the initial phases of our simulations (prior to multiphase gas condensation), the total energy flux is set by the sum of convective and thermal wind fluxes, whereas sound waves and the wind kinetic energy carry negligible amounts of energy.At late times, after multiphase gas condensation, the convective energy flux also drops and the wind becomes the dominant energy transport mechanism (Fig. 8).Multiphase condensation and outflows also strongly alter the z-profile of gas density (Fig. 6).
• The SNIa drive subsonic turbulence in the ISM which causes isobaric perturbations in it (figs.9, 10).The amplitudes of the density and pressure fluctuations are proportional to the compressive component of the rms Mach number.However, the fluctuations generated by SNIa are much larger than predicted by the density fluctuationsrms Mach number scaling relation in homogeneous turbulence (Konstandin et al. 2012), in agreement with Li et al. (2020b).We expect that these differences are due to the additional thermodynamics in our simulations, i.e., perturbations associated with the direct heating of the ISM by the SNIa and the radiative cooling of the ISM.
• For all 16 simulations, we find that multiphase condensation occurs if the criterion min(t ti )/t ff ≤ c 1 exp(c 2 σ s ) is satisfied.We obtain c 1 = 0.6 and c 2 = 6 using an empirical fit (Fig. 13).In contrast with the simulations of idealized driven turbulence in Mohapatra et al. (2023), we find turbulent mixing does not suppress multiphase condensation even when t ti /t mix ≫ 10.This is because the SNIa drive compressive turbulence which is poor at mixing perturbations in the gas.Thus one needs to be careful in interpreting the importance of t ti /t mix in observations, where it is difficult to determine the source/driving mode of turbulence.
8. ACKNOWLEDGEMENTS RM would like to thank Minghao Guo, Patrick Mullen and Jim Stone for help with setting up the simulations using Athenak and post-processing the results.We thank Chang-Goo Kim, Prateek Sharma, Romain Teyssier and Eve Ostriker for useful discussions.We thank the anonymous referee for their comments, which helped improve this paper.This work was supported in part by a Simons Investigator award from the Simons Foundation (EQ) and by NSF grant AST-2107872.The analysis presented in this article was performed in part on computational resources managed and supported by Princeton Research Computing, a consortium of groups including the Princeton Institute for Computational Science and Engineering (PICSciE) and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University.We also used the Delta GPU machine at National Center for Supercomputing Applications, Illinois, United States through allocations PHY230106 and PHY230045 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services and Support (ACCESS) program, which is supported by National Science Foundation grants 2138259, 2138286, 2138307, 2137603, and 2138296  .Time-evolution of net mass (row 1, col 1), energy (row 1, col 2), velocity dispersion of hot gas (row 2, col 1) and mass fraction of cold gas (row 2, col 2) for our fiducial fSN0.99 run and the fSN0.99 ṁ * run implementing mass and energy injection from AGB winds.We find no major differences between the two, except a larger fraction of cold gas formation for the fSN0.99 ṁ * run.
• 'H6hdens' runs: https://youtu.be/Ti4wZ0F x8 In Section 2.1.4,we discussed that the energy injected into the ISM due to the thermalisation of the material ejected by AGB stars can contribute to its heating.The mass injected by these stars can also change the density of the ISM and its cooling rate.We have described our implementation of mass and energy input due to AGB wind ejecta in Section 3.3.3.Here, we study their effect on the time-evolution of important parameters in our simulations.
In Figure A1, we show the evolution of the net mass, the net energy, the hot gas dispersion velocity and the mass fraction of cold gas for our fiducial f SN 0.99 run and the f SN 0.99 ṁ * run.Comparing the two, we do not find any major

Figure 1 .
Figure 1.Upper panel: Radial profiles of gas number density n, stellar mass density ρ * and entropy S for a typical massive elliptical galaxy, normalized to their values at 2 kpc.Lower panel: Radial profiles of the cooling rate density L, supernova heating rate density ėSN and stellar wind heat rate density ė * , assuming the stellar velocity dispersion σ * = 250kms −1 .

Figure 2 .
Figure 2. Snapshots of x-projected quantities for our fiducial set of simulations with different fractions of heating due to supernovae at t = tmp (t = t end for fSN0.01 and fSN0.1 runs that do not form multiphase gas)-Column 1: Logarithm of gas number density n; Column 2: Logarithm of mass-weighted temperature T , Column 3: Logarithm of x-projected density normalized by the z-profile of density; Column 4: Logarithm of x-projected pressure normalized by the z-profile of pressure.The runs with larger fSN show larger density and pressure fluctuations and are more likely to undergo multiphase condensation.The cold gas forms in small dense clumps.

Figure 3 .
Figure 3. Snapshots of x-slices of the gas number density showing the evolution (from left to right) of a supernova inflated bubble in one of our simulations.The sound wave driven by the supernova can be seen in the first column.The under-dense bubble rises buoyantly and gets disrupted due to the Rayleigh-Taylor instability.These motions drive turbulence in the ISM.

Figure 4 .
Figure4.Time-evolution of net mass (first row), energy (second row) normalized to their initial values, velocity dispersion (σ v,hot , column 3) of hot gas (T ≳ 10 5.5 K) and mass fraction of cold gas (T ≲ 10 4.2 K, column 4) for our fiducial set of simulations.A higher SNIa rate (larger fSN) leads to a larger σ v,hot .Stronger driving triggers multiphase condensation in the fSN0.5 and fSN0.99 runs.Imbalanced heating due to the SNIa also gives rise to a wind.In steady state, the mass and energy within the box decrease with increasing fSN due to the increasing impact of these two processes.

Figure 5 .
Figure5.Time-evolution of the mass fluxes at z = 0.5 (upp, first row) and z = −0.5 (low, second row) and the supernova mass injection rate ( ṀSN, third row) for our fiducial set of simulations.The mass fluxes at the domain boundaries are an order of magnitude larger than the mass injected by the SNIa, even for the fSN0.99 run.Thus gas inflow/outflow at the boundaries dominates the evolution of net mass.

Figure 8 .
Figure 8. Vertical profiles of energy fluxes for the fSN0.5 run at t = tmp (upper panel) and t = t end (lower panel), normalized by the mean SNIa rate over the entire volume.Convection and thermal wind flux are the two main contributors to the net energy flux at early times.At late times, the net energy flux is mostly positive and dominated by the thermal wind flux.

Figure 9 .
Figure9.Mass-weighted probability distribution functions of gas number density, temperature and Mach number for our fiducial set at t = min(tmp, t end ).The two peaks in the distribution for the fSN0.5 and fSN0.99 runs correspond to the hot and cold phases.With increasing fSN, the width of the PDFs of the hot phase show an increase, consistent with increasing strength of turbulence.

Figure 10 .
Figure 10.Phase diagram showing the joint mass-weighted distribution of the gas temperature and density (in log-scale) for our fiducial set of simulations.The perturbations in the hot phase are mostly isobaric in nature.Recently injected SNIa appear as adiabatic fluctuations.The peak near T ∼ 10 4 K for 'fSN0.5'and 'fSN0.99'runs correspond to cold gas that has reached T floor (= 10 4.2 K).
Figure 12.First row : The square the amplitude of logarithmic density fluctuations (σ 2 s ) plotted against the Mach number (M) (left column) and its compressive component Mcomp.The larger symbols represent the data-points at t = min(tmp, t end ) and the smaller symbols show their values are at earlier times.The filled symbols denote simulations that form multiphase gas.(right column).The dotted line on the right panel shows the predicted scaling relation for subsonic turbulence with compressive forcing from Konstandin et al. 2012.Second row : Similar to the first row, but for the square of logarithmic pressure fluctuations σ 2 ln(P/⟨P ⟩) instead.The dashed lines in all four panels show our empirical fits to the data.The pressure fluctuations show a tight scaling with Mcomp.
Figure A1.Time-evolution of net mass (row 1, col 1), energy (row 1, col 2), velocity dispersion of hot gas (row 2, col 1) and mass fraction of cold gas (row 2, col 2) for our fiducial fSN0.99 run and the fSN0.99 ṁ * run implementing mass and energy injection from AGB winds.We find no major differences between the two, except a larger fraction of cold gas formation for the fSN0.99 ṁ * run.