The small-scale dynamo in a multiphase supernova-driven medium

We use the Pencil Code (Brandenburg & Dobler, 2002; Brandenburg et al., 2020) to model SN-driven turbulence as described previously in Gent (2012), Gent et al. (2013a), and Gent et al. (2020). We solve the system of nonideal, compressible, nonadiabatic, MHD equations

with the ideal gas equation of state closing the system, assuming an adiabatic index (ratio of specific heats) $c_{\rm p}/c_{\rm v}=5/3$. Most variables take their usual meanings; a list of notations used is given in Table 2.

Terms containing $\nu_{6},~{}\chi_{6}$ and $\eta_{6}$ apply sixth-order hyperdiffusion to resolve grid-scale instabilities (see, e.g., Brandenburg & Sarson, 2002; Haugen & Brandenburg, 2004), with mesh Reynolds number set to be $\simeq 1$ at the scale of the zone size ${\delta x}$. Terms including $\zeta_{D},~{}\zeta_{\nu}$ and $\zeta_{\chi}$ resolve shock discontinuities with artificial diffusion of mass, momentum, and thermal energy, respectively. They depend quadratically on the local strength of the shock (see Gent et al., 2020, for details). Equations (2) and (4) include momentum and energy conserving corrections for the artificial mass diffusion $\zeta_{D}$ applying in Equation (1). Following Gent et al. (2021), resistive shock diffusion is omitted from Eq. (3).

SNe at rate $\dot{\sigma}$ inject $E_{\rm th}=10^{51}~{}{\rm erg}$ thermal energy. In dense regions up to 5% is instead injected as kinetic energy $E_{\rm kin}$ (see Kim & Ostriker, 2015; Gent et al., 2020). Nonadiabatic heating $\Gamma_{\rm UV}$ and cooling $\Lambda(T)$ (as detailed in Gent et al., 2013a) follow Wolfire et al. (1995) and Sarazin & White (1987).

We simulate a nonstratified, nonrotating domain with initial gas number density $n=1~{}{\rm cm^{-3}}$ covering (256 pc)³ and periodic across all boundaries. Resolution along each edge spans 0.5 to $4~{}{\rm pc}$, corresponding to grid sizes of $64^{3}$ to $512^{3}$ zones. Each model SN rate is given relative to the estimated rate in the Solar neighbourhood of the Milky Way, $\dot{\sigma}_{\rm Sn}=50~{}{\rm Myr}^{-1}~{}{\rm kpc}^{-3}$. SNe occur uniform randomly in space at times following a Poisson process. Models with common $\dot{\sigma}$ use the same SN schedule and locations. For ambient gas number density $n\simeq 1~{}{\rm cm^{-3}}$ and $\dot{\sigma}\simeq\dot{\sigma}_{\rm Sn}$ an estimated forcing scale of 60–100 pc (Joung & Mac Low, 2006; de Avillez & Breitschwerdt, 2007; Hollins et al., 2017) should support at least 2–4 turbulent cells. For higher SN rates the forcing scale reduces, increasing the number of turbulent cells (Joung et al., 2009), but the forcing scale remains unchanged for the lower SN rates applying here, where an individual SN at gas number density $1~{}{\rm cm^{-3}}$ merges with the local sound speed within a 70 pc radius (Gent et al., 2020)

In addition, Model LP5r8sm-S, is vertically stratified without rotation and has open vertical boundaries. Initially $n\simeq 1~{}{\rm cm^{-3}}$ at the midplane. With this we examine the effect that vertical expansion of SN remnants or superbubbles has on the SSD. Here $\Gamma_{\rm UV}$ is amplified to $3.5\times$ that of Wolfire et al. (1995) to support the thickness of the disk in the absence of ionization heating and cosmic ray pressure gradients (see Hill et al., 2018). The domain size of this model is $512~{}{\rm pc}$ along the disk and $\pm 1.534~{}{\rm kpc}$ perpendicular to it. SN are located uniform random horizontally and normal random with scale heights of $90~{}{\rm pc}$ and $325~{}{\rm pc}$ for Type I and Type II SNe, respectively (Ferrière, 2001). This model otherwise is as described in Gent et al. (2013a, b).

All models are listed in Table 1 where the model labelling convention is explained.

Angular brackets indicate the quantity is averaged over the volume, or also with a subscript $T$ when over the volume of individual temperature phases. An overbar indicates averaging over a domain and an interval of time, intervals which may vary, as explained in the text.

In the case of ${\overline{e_{\rm K}}}$ the interval is selected for each model to exclude initial transients and the period after the SSD approaches saturation. Even with saturation around 5% of equipartition energy, some models show damping in kinetic energy. Models with identical SN rates, have near identical kinetic energy density evolution, so direct comparisons are not affected by this choice. The magnetic energy density $\langle e_{\rm B}\rangle$ is then normalized by ${\overline{e_{\rm K}}}$ to ease model comparison.

The erratic growth of the volume averaged magnetic energy shown in Figure 1 and discussed in Sect. 3 shows that the kinematic stage of the SSD does not have a well-defined single growth rate (Gent et al., 2021). SSD growth typically follows an exponential of the form

In Figure 1 we fit such a function to specified time intervals for each model, but the SSD growth (or decay) at other times differs, with ${\overline{\gamma}}$ varying or no clear exponential behaviour occurring. Indeed within the domain, we can find very different growth patterns within and between thermal phases.

To interpret how the time and volume averaged growth rate ${\overline{\gamma}}$ is influenced by the phases and dynamical properties, we can use Eq. (3) to identify instantaneous changes to the magnetic energy. Hyperdiffusion, being purely numerical, is neglected. Taking the curl of Eq. (3) and contracting it with ${{\mbox{\boldmath$B$}}}\mu_{0}^{-1}$, we obtain an equation for the change of the magnetic energy density

Negative values represent decay of the magnetic field and positive its amplification. Dividing by $e_{\rm B}$, we obtain an equation for the relative growth rate exponent, of the form $e_{\rm B}\propto\exp(\Gamma t)$, at each instant in time and location in space

Here, $\Gamma({{\mbox{\boldmath$x$}}},t)$ is a function of position $x$, distinct from ${\overline{\gamma}}$ defined above, which is the volume and time averaged quantity. Statistical analysis of $\Gamma$ in relation to various physical properties will help determine how the varying growth depends on the multiphase structure of the ISM.

The growth rate $\Gamma$ does not indicate the absolute magnitude of the energy change. We also need to identify where the largest changes in magnetic energy density occur, since these need not correlate with the growth rates. We therefore also define the absolute growth rate

replacing $e_{\rm B}$ in the denominator with the time dependent volume averaged ${\langle e_{\rm B}\rangle}(t)$. Rescaling by time dependent $\langle e_{\rm B}\rangle$ assists comparison between all stages of the SSD. Both $\Gamma$ and $\tilde{\Gamma}$ have values that span orders of magnitude of both signs. To take advantage of logarithmic scales in our histograms we omit negative and negligible growth rates.

Similarly, a field of values for the magnetic Reynolds number Rm is calculated directly from the induction equation by taking the ratio of the advection to the diffusion terms,

where ${{\mbox{\boldmath$j$}}}=\mu_{0}^{-1}\nabla\times{{\mbox{\boldmath$B$}}}$. In a related approach, Evirgen et al. (2019) decompose the terms of Equation (2) to identify separately the spatial variation of each force.

There are various approaches to measuring the proportion of the ISM that contains gas in each phase, as discussed in detail within Gent et al. (2013a). Here, we measure the fractional volume of SSD activity for each phase $i$ defined as

where $V_{i}$ is the volume occupied by the gas in the temperature range defining phase $i$, as listed in Table 3, and $V$ is the total volume. Both the phase volume and the total volume omit locations with negative or negligible $\Gamma$ to focus on SSD growth. For each phase, $f_{V,i}$ is computed from the snapshot within the time intervals listed in Table 3.

To examine how properties of the ISM in each phase differ between times with rapid growth and slow growth or decay, for each model we compute the average growth rate $\gamma(t)$ as the rate of change in ${\langle e_{\rm B}\rangle}$ during the specified intervals from Table 3. The intervals are chosen at each resolution to exclude initial transient magnetic energy decay, prior to the onset of any SSD, and most of the subsequent saturated dynamo stage, so that we capture data of most relevance to SSD. The time series is then binned according to where $\gamma(t)$ is lowest, median and highest.¹¹1 Sampling only the lowest, median, and highest 5%, 15%, or 25% yielded results consistent with our method of simply splitting the growth rates into three bins.

Histograms for each phase from snapshots at times belonging to each time bin are then accumulated. Fractional volume of SSD activity $\overline{f}_{V,i}$ for each phase $i$ for a cumulative histogram is the mean of $f_{V,i}$ in snapshots contributing to that histogram. Summary statistics of mean $\log\Gamma$ or $\log\tilde{\Gamma}$ and relevant averaged physical quantities are calculated from cumulative histograms in each bin.

In Figure 3 we display the set of cumulative joint histograms of absolute growth $\tilde{\Gamma}$ and (twice the) enstrophy or the norm squared vorticity $|\omega|^{2}$ of cold, warm and hot gas from Model HP5r2sl. As is common across all models and variables, there is large variance and histograms overlap between the phases. However, a clear trend appears of increasing vorticity and growth rate with temperature. In the warm and hot phases there is discernible positive correlation between vorticity and growth rate within the phase. The phases and time intervals used for the histograms are listed in Table 3.

The most striking correlation we find is between relative growth rate $\Gamma$ and the norm squared vorticity $|\omega|^{2}$, which we show in Figure 4(a). We find that the growth rate $\log\Gamma$ is strongly proportional to $\log|\omega|^{2}$, and vorticity increases as we move from cold to hot gas at all times. At times belonging to the high growth rate $\gamma(t)$ bin (right panel) the hot gas increases in vorticity, while the vorticity in the cold and warm gas does not differ so much between bins. The fractional volume of the hot gas increases slightly for high resolution within the upper bin. So the efficiency of the SSD is linked to high vorticity and high local and average growth rates $\Gamma$ and $\gamma$ are correlated with increased vorticity in the hot gas, consistent with results from models using more idealised turbulence (Federrath et al., 2011; Achikanath Chirakkara et al., 2021)

To confirm whether this correlation is reflected in the absolute gains in magnetic energy we plot the same summary for $\tilde{\Gamma}$ in Figure 4(b). At high resolution the same trends are preserved as for $\Gamma$, with the changes in $\tilde{\Gamma}$ most associated with changes in the characteristics of the hot gas. However for low resolution $\tilde{\Gamma}$ is actually anti-correlated with vorticity. Increased mixing at low resolution dampens vorticity and inhibits the formation of hot gas, damping the SSD.

Energy growth at low resolution is better correlated with gas density, the cold and then warm phases exhibiting higher $\tilde{\Gamma}$, consistent with results from isothermal modelling, except at $\delta x=2~{}{\rm pc}$ for the upper bin of $\gamma(t)$ (right) when vorticity and fractional volume increases for the hot gas. Mean vorticity in the warm and cold phases does not change markedly between bins.

Results from models of isothermal SSDs indicate that highly compressible flows inhibit dynamo activity (Haugen et al., 2004; Federrath et al., 2011, 2014; Seta & Federrath, 2021), while high magnetic Reynolds number Rm makes SSDs more likely and is correlated with higher growth rates (Schekochihin et al., 2002, 2005; Iskakov et al., 2007; Brandenburg, 2011; Schober et al., 2012; Seta et al., 2020; Warnecke et al., 2022).

In Figure 5(a) we display the summary statistics for Mach number $\cal{M}$ against relative growth rate $\Gamma$ for each bin of volume-averaged magnetic energy growth rate. Overall the expected trend of $\Gamma$ reducing as $\cal{M}$ increases is visible, more so in the lowest growth rate regions (left panel). We observe that hot gas with high sound speed has low Ms and high relative growth rate $\Gamma$ and cold gas the inverse. What is, perhaps, unexpected is to see that within each phase there is a trend of increasing $\Gamma$ with increasing $\cal{M}$. This can be explained by the counter effect of increased velocities, which drive higher Mach numbers for a given sound speed, but also higher Rm. So increased $\Gamma$ with $\cal{M}$ should have complementary trends in Rm.

In this highly compressible system it is challenging to calculate Rm as would conventionally be defined for weakly compressible simulations. Instead, we use the ratio of advective to diffusive terms given by Eq. (9). With strongly fluctuating characteristic velocities and length scales applying between the continuously shifting phases, and the very sporadic and flexible scale of the SN forcing, Rm varies as a function of position and time. An illustration of the extent to which Rm varies is shown in the histograms of Figure 6. While the mean Rm is quite modest of order 16 and 45 in the cold and warm phases, there are regions in all phases where Rm $>10^{3}$, and in the hot phase Rm $>10^{4}$. For the summary statistics the log mean values are quite informative to identify the trends in SSD dependencies across simulations over time and between phases.

Results for Rm versus $\Gamma$ are shown in Figure 5(b). Velocities increase with temperature, increasing Rm and supporting faster SSD, an overall trend as would be expected for $\Gamma$. Within each phase and within each $\gamma(t)$ bin, however, $\Gamma$ is weakly anticorrelated with Rm, and across each bin the patterns are largely unaffected for cold and warm gas, suggesting higher velocities may be associated with higher Mach numbers. Results are similar for Re versus $\Gamma$ as shown in Figure 5(c) for samples from models that include explicit viscosity $\nu$. The higher resolution models ($\delta x=1$) show a good correlation between Re and growth rate. The lower normalization of that correlation seen in the lower resolution models may well be due to reduced kinetic energy from excess cooling in under-resolved density gradients. The comparison between the 1 pc and 0.5 pc models without viscosity suggests that the 1 pc model is well converged. For the hot gas Rm, Re, and $\Gamma$ all increase for the upper growth rate bin, suggesting increased velocities are associated with rotational rather than compressive flows.

If this were the case, we would expect to see this reflected in the structure of the flow, which we consider in Figure 7 (see analysis of hydrodynamic flow by Käpylä et al., 2018). We use a Helmholtz decomposition to separate the velocity field $u=u_{1}+u_{2}$ into purely solenoidal $u_{1}=\nabla\times{{\mbox{\boldmath$F$}}}$ and compressive flows $u_{2}=\nabla\Phi$, in which $F$ and $\Phi$ are, respectively, a vector potential and scalar potential pair from which orthogonal flows $u_{1}$ and $u_{2}$ can be derived.

Both flows appear strongly correlated with $\Gamma$, although not as clearly as vorticity in Figure 4(a). It is reasonable to expect that the magnitude of compressive and solenoidal flows would be correlated, and as such some common correlation to $\Gamma$ would be evident. The clearer correlations for (b) rotational rather than (a) compressive flows indicates that the latter is coincidental. Growth rate trends in the warm and cold gas do not alter significantly between bins of $\gamma(t)$ even as flow strength varies. In the hot gas $\Gamma$ clearly increases when $|\nabla\times{{\mbox{\boldmath$F$}}}|$ increases between low (left) and high (right) $\gamma(t)$ bins, especially for high resolution models.

In Figure 7(c) the summary statistics are plotted of the compressive ratio (Kritsuk et al., 2007) defined by

This shows that the ratio of solenoidal energy to total compressive and solenoidal energy $(1-R_{\rm cs})$ is typically between 60% (cold) and 90% (hot). At higher resolution the fraction of energy in solenoidal flow is higher. The compressive ratio reduces from cold to hot phases. The growth rate is directly correlated with the fraction of solenoidal energy.

The multiphase structure and intermittent forcing of the ISM causes the fluid Reynolds number Re to vary just as was shown for Rm in Figure 6. This results in a wide variation of magnetic Prandtl number Pm. In Gent et al. (2021, see Figures 5[b1] and [b2]) we examined the role of Pm on the SSD in ISM simulations, by varying both viscosity $\nu$ and resistivity $\eta$. For a fixed $\eta$, the SSD did not depend on Pm and the energy at which the dynamo saturated was consistently around 5% of equipartition. For a fixed $\nu$, on the other hand, SSD growth rate did increase with Pm for ${\rm Pm}\leq 5$ in the parameter space considered. The energy at saturation increased with Pm, remaining steady for ${\rm Pm}\gtrsim 5$. This suggested that it is the ratio of the Lorentz force to Rm rather than Pm that is the controlling parameter (see also Oishi & Mac Low, 2011, for a similar suggestion in magnetorotational dynamos).

Therefore, here we vary Pm by varying only $\eta$. Where the viscosity, $\nu$, and the resistivity, $\eta$, are explicitly defined, we can nominally express the magnetic Prandtl number, ${\rm Pm}=\nu/\eta$, which is what we shall use here. Alternatively, as a function of position Pm can also vary due to the flows and length scales characteristic at different regions and times altering the effective ratios of Rm to Re. In the Appendix we determine the effective magnetic Prandtl number ${\rm Pm}_{\rm eff}$ in various models to examine the relevance the explicit diffusivities to the dynamics (see Figure 12).

In Figure 8(a) we show the volume averaged magnetic energy density ${\langle e_{\rm B}\rangle}$ for models with ${\rm Pm}\in[1,5]$ with resolution $\delta x=1~{}{\rm pc}$. Around $20~{}{\rm Myr}$ there is a burst of SSD, common to all three runs, highlighted between the two leftmost vertical dotted lines. Growth rates ${\overline{\gamma}}$ applying for $18<t<20~{}{\rm Myr}$ depend on Pm in this range. Sensitivity to Pm reduces as Pm approaches 5, consistent with the idea that growth rates become asymptotic for higher Pm. We also fit ${\overline{\gamma}}$ for ${\rm Pm}=1$ at $41<t<44~{}{\rm Myr}$, where the growth rate exceeds that for ${\rm Pm}=5$ earlier. Under most circumstances we would expect ${\overline{\gamma}}$ to be lower for the larger $\eta$, but during this period the fractional volume of the hot gas $f_{V,i}$ is higher (Figure 8(b)) in all models than during the earlier subinterval. A surge aligned to the rightmost vertical dotted line is also evident in the high Pm models, although it is inhibited in those by the onset of saturation.

The critical Rm above which SSD can be excited and the typical scaling relations with Rm for the rate of magnetic energy growth are not well determined. For ${\rm Pm}=1$ SSD actually decays during two periods, before its final surge. This could be due to Rm varying or changing conditions resulting in a lower critical Rm.

We plot the hot gas fractional volumes $f_{V,h}$ for the models with ${\rm Pm}\in[1,5]$ in Figure 8(b). At first, the hot gas fractions $f_{V,h}$ for all models coincide. Despite the initial weakness of magnetic effects, the fractional volumes diverge slightly within 20 megayears. The resistive timestep differs between the models; small changes in timestep lead to the chaotic solutions adopting alternate trajectories. For example a delay of a few decades in scheduling a single SN, even at the identical location, alters the specific ambient conditions modifying the remnant evolution and subsequent dynamics, propagating into diverging trajectories between models.

Nevertheless the trends are consistent between models. In the subinterval between the two leftmost vertical dotted lines, $f_{V,h}$ is enhanced, corresponding to a burst in SSD activity. The dip in $f_{V,h}$ between the middle two vertical lines corresponds to decay in the ${\rm Pm}=1$ model and slower growth at higher Pm, followed by another SSD burst aligned with the third vertical line and a peak in $f_{V,h}$, least so for ${\rm Pm}=5$. The value of $f_{V,h}$ then rises again at 40 megayears, particularly for ${\rm Pm}=1$, with an accompanying boost in SSD. The described behavior supports the hypothesis that higher SSD activity is correlated with increasing fractional volume of the hot gas $f_{V,h}$.

To examine our hypothesis further, we split the evolution of the magnetic energy between phases ${\langle e_{\rm B}\rangle}_{T}$ and plot these in Figure 8(c) for the models with ${\rm Pm}=1$ and 5. These are normalised by ${\overline{e_{\rm K}}}$. Horizontal dotted lines of matching color for each phase indicate the time-averaged kinetic energy density by phase ${\overline{e_{\rm K}}}_{T}/{\overline{e_{\rm K}}}$. These lines, plotted for ${\rm Pm}=5$, are very similar for ${\rm Pm}=1$.

Usually ${\langle e_{\rm B}\rangle}_{T}$ is smallest in the hot gas, but during the SSD activity bursts, its growth is the most rapid and its energy density then exceeds that of the other two phases. The saturation energy of the hot gas is similar for both models, around $0.05{\overline{e_{\rm K}}}_{h}$ at 40 megayears for ${\rm Pm}=5$, and 55 megayears for ${\rm Pm}=1$. In contrast the saturation energy in the cold and warm phases is affected by ${\rm Pm}$. In the nonlinear phase the magnetic energy in the hot phase ${\langle e_{\rm B}\rangle}_{h}$ decays to around $0.01{\overline{e_{\rm K}}}_{h}$, while the magnetic energy in the cold phase ${\langle e_{\rm B}\rangle}_{c}$ grows, above ${\overline{e_{\rm K}}}_{c}$ in the case of ${\rm Pm}=5$.

Is the latter consistent with compression in the statistical steady MHD state? Under compression, magnetic field strength relates to density $|B|\propto\rho^{\alpha}$, where $\alpha\in[0,1]$ for compression along field lines up to disc or slab-like compression (Tritsis et al., 2015). In Figure 9, for the two models in Figure 8(c) we plot the mean gas number density in the cold gas (left axis) and the root mean square strength of the magnetic field normalised by cold gas density $\langle b_{\rm rms}\rangle_{C}=\langle B^{2}\rangle_{C}^{1/2}\langle\rho\rangle^{-\alpha}_{C}$ (right axis), taking $\alpha=1/2$. In the interval $50~{}{\rm Myr}<t<60~{}{\rm Myr}$ the dynamo has clearly saturated for the warm and hot gas (Figure 8). If the continued growth of the magnetic field in the cold gas were due to compression of the field with the gas, then we would expect the ratio of $b$ to remain steady or even decay. Instead, in both models this ratio continues to increase until $\sim 65~{}{\rm Myr}$, indicating that dynamo remains active in the cold phase. The increasing mean density in the cold gas toward the end of the simulation reflects the increasing fractional volume of hot gas and mean thermal energy density. The cooling in the hot gas is too inefficient to balance the supply of thermal energy from subsequent SNe in the closed system, and thermal runaway will eventually result. Applying $2/3<\alpha<1$ in the analysis does not qualitatively alter the development of $b$.

The strength of the SSD in SN-driven turbulence has been found to be sensitive to the SN rate, $\dot{\sigma}$ (Balsara et al., 2004; Gent et al., 2021). Balsara et al. (2004) used an SN rate of $8\leq\dot{\sigma}\leq 40\dot{\sigma}_{\rm Sn}$. Because of their confinement within the periodic boundaries, these models were dominated by thermal runaway, as described by Kim & Ostriker (2015), in which the fractional volume of hot gas reaches $f_{V,h}>0.9$. Therefore, SSD has a steady exponential growth, because the ISM has effectively become a homogeneous, single phase medium. Indeed, in similar experiments by Gent et al. (2021), thermal runaway was easily excited at high resolution ($\delta x\leq 1~{}{\rm pc}$) even with $\dot{\sigma}=\dot{\sigma}_{\rm Sn}$. At low resolution ($\delta x=4~{}{\rm pc}$), with consequentially higher numerical diffusion, and thus unphysically strong cooling, thermal runaway was still excited for $\dot{\sigma}=10\dot{\sigma}_{\rm Sn}$. Rapid steady acceleration of the SSD occurs, which is excluded here with lower $\dot{\sigma}$.

In Figure 10(a) we compare evolution of ${\langle e_{\rm B}\rangle}$ for $\dot{\sigma}=0.2\dot{\sigma}_{\rm Sn}$ and $1\dot{\sigma}_{\rm Sn}$ at $\delta x=2~{}{\rm pc}$ and 4 pc, and (b) the fractional volume of the hot gas $f_{V,h}$. The higher rate $\dot{\sigma}$ produces more hot gas and sustains the dynamo, more so for $\delta x=2~{}{\rm pc}$ (blue) than $\delta x=4~{}{\rm pc}$ (brown). For lower $\dot{\sigma}$, ${\langle e_{\rm B}\rangle}$ decays, but growth occurs briefly in the subinterval between the first two vertical lines, while $f_{V,h}$ increases, particularly for $\delta x=2~{}{\rm pc}$ (green). Least-squares fits for ${\overline{\gamma}}$ are shown for this subinterval (green and cyan). Higher $f_{V,h}$ occurs for higher $\dot{\sigma}$ with fits for ${\overline{\gamma}}$ (blue and red) indicated in the second subinterval. Another boost in SSD occurs for $\delta x=4~{}{\rm pc}$ after 270 megayears when $f_{V,h}$ is at an even higher peak. In general for $\delta x=4~{}{\rm pc}$ there is less hot gas than for $\delta x=2~{}{\rm pc}$ at each given SN rate, because of higher numerical diffusion and cooling.

Mean kinetic energy density is not strongly affected by the SN rate (Gent et al., 2021, see their Figure 3[d]), so the change in SSD growth rate would appear to be due to the increased heating of the ISM at higher $\dot{\sigma}$. Figure 10 demonstrates that this is indeed the case. High SSD growth rate ${\overline{\gamma}}$ coincides with high $f_{V,h}$ for hot gas within each model and between models, higher ${\overline{\gamma}}$ occurs in models that support higher $f_{V,h}$.

Anticipating future examination of the LSD, in Model LP5r8sm-S we consider how stratification affects the SSD. To exclude the LSD, we omit large-scale rotation and shear. Inhomogeneities due to non-periodic boundaries vertically might still cause an LSD effect, so any residual mean magnetic field volume and horizontal component average is subtracted at each timestep. In contrast to the periodic domains, the open vertical boundary enables gas flow and magnetic helicity flux to and from the disk.

In Figure 11(a) we plot ${\langle e_{\rm B}\rangle}$ for this model with $\dot{\sigma}=0.5\dot{\sigma}_{\rm Sn}$ alongside models with $\dot{\sigma}=\dot{\sigma}_{\rm Sn}$. Given that ${\overline{\gamma}}$ always increases with $\dot{\sigma}$ in otherwise equivalent models, it is evident that the substantially increased growth rate in model LP5r8sm-S is due to its stratification. We also see from Figure 11(b) that the stratified Model LP5r8sm-S has the largest fractional volume of hot gas, although none of the low resolution periodic models have large $f_{V,h}$ (see also Figure 10]b]).

From the summary statistics in Figure 4(b), we find for $\delta x=4~{}{\rm pc}$ that the absolute growth $\tilde{\Gamma}$ of magnetic energy is fastest in the cold gas, but occurs mostly in the warm gas due to its higher fractional volume. In contrast the stratification supports high ${\overline{\gamma}}$ and high $f_{V,h}$ for the hot gas, depicted in Figures 4–7 by a large orange plus symbol.

As ${\overline{\gamma}}$ increases, the absolute growth rate $\tilde{\Gamma}$ in the hot gas reduces, as shown in Figure 4(b) (left to right), for $\delta x=4~{}{\rm pc}$, While for $\delta x\geq 2~{}{\rm pc}$ it increases. In the cold gas $\tilde{\Gamma}$ increases for $\delta x=4~{}{\rm pc}$ and reduces for higher resolution. At $\delta x=4~{}{\rm pc}$, however, ${\overline{\gamma}}$ does not vary over time as much as for higher resolution. At all stages $\tilde{\Gamma}$ in the warm gas is five to ten times larger in the stratified model than the periodic models.

Vorticity in all phases is much lower (Figure 4) at $\delta x=4~{}{\rm pc}$, but of these it is highest in the stratified model. The values of $\cal{M}$, ${\rm Rm}$, and velocities (Figs. 5 and 7), are also highest in the stratified model. Higher velocities, and vorticity, in all phases are possible with stratification, because the gas can expand into lower pressure diffuse gas away from the midplane.

Consensus based on studies of isothermal MHD turbulence predicts that ${\rm Rm}_{\rm crit}$, the critical Rm above which SSD can be excited, reduces as Pm increases above ${\rm Pm}\gtrsim 0.1$ with an asymptotic limit of order ${\rm Rm}_{\rm crit}\sim 20$ for ${\rm Pm}\gg 1$ (Schekochihin et al., 2002; Schober et al., 2012; Seta et al., 2020). Less directly relevant to our context of ISM turbulence, for low Pm there is a range $0.1>{\rm Pm}>0.01$ where ${\rm Rm}_{\rm crit}$ has a maximum (Schekochihin et al., 2005; Iskakov et al., 2007; Brandenburg, 2011), and then for ${\rm Pm}<0.01$ ${\rm Rm}_{\rm crit}$ again begins to drop (Warnecke et al., 2022). Once excited, the rate of growth of SSD also increases, up to some asymptote, for higher Pm. In these isothermal systems, it has also been demonstrated that as $\cal{M}$ increases it becomes more difficult to excite SSD (${\rm Rm}\gtrsim 60$ for ${\rm Pm}\geq 5$ Haugen et al., 2004) and growth rates are correspondingly reduced (Federrath et al., 2011, 2014). Hence, according to these findings, SSD can be excited in our models, as this threshold value is easily exceeded.

Superficially, these conclusions might appear irrelevant to an SSD in multiphase SN-driven turbulence with respect to the averages across all phases. The response of models with respect to Pm is inconsistent. Higher SN rates with correspondingly higher compressive forcing appear to excite SSD more easily and with larger growth rate. SSD does not appear to follow a dominant eigenmode, establishing a distinct exponential amplification until saturation.

However, detailed analysis reveals that what is inconsistent in this more realistic ISM is the distribution of temperatures, gas densities, and flow properties throughout the thermal phases and over time. Determination of Pm, Rm, Re, and $\cal{M}$ must be specific to each region of the ISM. It is likely that the variation we see by position in these properties occurs also in isothermal simulations, but that the statistics of these remain quite consistent throughout the domain and duration. In a multiphase system, on the other hand, these statistics vary by orders of magnitude as a function of position and time.

The breakdown of the models by phase and overall SSD growth rate shows that the insights from isothermal SSD experiments do apply, but contrary trends, such as increasing growth in SSD with increasing $\cal{M}$, seen in Figure 5(a), are due to the SSD supportive counter trends, such as increasing velocity (with $\cal{M}$). The expected effect of $\cal{M}$ on growth rate is well demonstrated in the anticorrelation going from hot to cold phase. Conversely, in Figure 5(b), Rm seems anticorrelated with SSD growth rate in the hot phase, indicating that SSD is more easily excited as the hot filling fraction and typical temperature of the hot gas increase, even for lower Rm. The expected trend of positive correlation between Rm and SSD growth rate is evident in the differences between phases.

In their two-phase simulations, Seta & Federrath (2022) determined that the magnetic energy in their models grows more slowly than in isothermal models with the same forcing and equivalent mean Mach numbers $\cal{M}$$\sim 1$ for the warm gas and 5 for the cold gas. Their SSD saturates after around 400 megayears for solenoidal forcing and 800 megayears for compressive forcing, for a grid resolution comparable to Model UP0r10sl. For this and all our 1 pc resolution periodic models, the SSD saturates within 40–60 megayears with mean $\cal{M}$$\sim 1$, 0.6 and 0.25 in the cold, warm and hot gas, respectively. However, in all phases $\cal{M}$ spans around two orders of magnitude over a large proportion of the volume, and outliers extend even further. The variance is much larger than in the two-phase simulations. The faster SSD in our models may, therefore, be less sensitive to the mean $\cal{M}$, as growth occurs predominantly in low Mach number, high vorticity regions, particularly in the hot gas.

The SSD persists longest in the cold phase and next longest in the warm phase (Fig. 8[c]). We speculate that an explanation for the growth rates differing between phases in three-phase and not in two-phase simulations is that high vorticity generated in the hot gas transfers to warm then cold regions, but with a time delay that is significantly longer than from the warm to cold phase. As a result, we find faster growth in all phases.

Estimates for the post inflationary seed magnetic field are in the range $10^{-65}$–$10^{-20}\upmu$G (Grasso & Rubinstein, 2001; Subramanian, 2016). With an exponent of ${\overline{\gamma}}=500~{}{\rm Gyr}^{-1}$ this could be amplified to sub-equipartition saturation of around $1\upmu$G within 170 to 590 Myr. Given these results consistently yield saturation of the SSD significantly below equipartition with the kinetic energy density (see also Schober et al., 2015), this does not exclude LSD action being sufficient to amplify even the lowest estimates of post-inflation magnetic field to the values observed in high redshift galaxies.

It is perhaps surprising, and inconsistent with current understanding of dynamo, to see in Figure 6 locations with ${\rm Rm}\leq 1$, well below isothermal predictions of ${\rm Rm}_{\rm crit}$, coincident with growth rate $\tilde{\Gamma}>10^{3}~{}{\rm Gyr}^{-1}$ in cold gas. These are outliers in a highly inhomogeneous multiphase structure, in which there is phase mixing and energy transfer within and between phases. The bulk of the simulation data corresponds to high Rm and high $\tilde{\Gamma}$ with the positive correlation reflected in the increasing log mean from the cold to hot phase.

To examine whether the high apparent growth rates at low Rm can be explained without SSD, let us approximate the effect of a blast wave following Sedov-Taylor (Sedov, 1959; Taylor, 1950) travelling at $150~{}{\rm km~{}s}^{-1}$, as would occur in an SN remnant with radius near 30 pc in a $10^{4}$ K medium of uniform gas number density $1~{}{\rm cm^{-3}}$. From Equation (8) with the normalized quantity $\hat{{{\mbox{\boldmath$B$}}}}={{\mbox{\boldmath$B$}}}\langle{\mu_{0}e_{B}}\rangle^{-1/2}$ we can write

Considering the first term with a shock front traveling at $150~{}{\rm km~{}s}^{-1}$ into warm gas with sound speed around $12~{}{\rm km~{}s}^{-1}$ resolved in our model over three cells, $-\nabla\cdot{{\mbox{\boldmath$u$}}}\approx 138~{}{\rm km~{}s}^{-1}(3~{}{\rm pc})^{-1}\approx 4.6\times 10^{4}~{}{\rm Gyr}^{-1}$. If the local magnetic field is even twice the volume averaged strength, this would yield $\tilde{\Gamma}\gtrsim 10^{5}~{}{\rm Gyr}^{-1}$. The effect of the other terms and nonadiabatic cooling would adjust this estimate up or down, but nevertheless it remains entirely consistent with the rare high values of $\Gamma$ and $\tilde{\Gamma}$ in regions of low Rm in Figure 6.

For a magnetic field perpendicular to the blast wave

reaching a maximum at the shock centre, where the first derivative approaches infinity and the second derivative vanishes, but near zero either side of the shock where the inverse is true. Both can coexist with the compression $\tilde{\Gamma}$ estimated above, explaining these low Rm outliers in the joint histograms.

Conventionally, a single characteristic Rm for a domain is related to ${\overline{\gamma}}$, associated with the single exponential growth rate of the dominant eigenmode of the SSD. The analysis of Grete et al. (2017) found that the Lorentz forces could induce weak nonlocal energy transfers, but that magnetic and kinetic energy transfers in both directions are predominantly highly localised in $k$-space, including for isothermal compressible MHD. Hence, the differences in characteristic length scales between the phases would tend to yield Rm supportive of very different dominant dynamo modes. Thus, in the multiphase ISM the SSD growth rate ${\overline{\gamma}}$ represents a superposition of varying dynamo modes, in which the dominant modes can switch from one phase to another over time due to heating, cooling and phase mixing.

The mean $\Gamma$ of the instantaneous localised growth rates collected while $\gamma(t)$ is highest have similar magnitude to ${\overline{\gamma}}$ as fitted in each simulation, supporting the concept of a superposition. Sufficiently high mean $\Gamma$ or $\tilde{\Gamma}$ correlated with a supercritical mean Rm would dominate the overall growth rates irrespective of the outliers.

We considered the properties in the simulations of velocity, Mach number, fluid and magnetic Reynolds numbers, Prandtl number, kinetic and magnetic helicity, and kinetic energy to identify how they relate to the relative and absolute growth rates of magnetic energy $\Gamma$ and $\tilde{\Gamma}$. By far the clearest correlation to growth rate was in the enstrophy, or norm squared vorticity.

The dependence of growth rate on velocity shows a weaker dependence, which is very similar to that of the Helmholtz decomposed flow displayed in Figure 7, and even less correlation is seen in the kinetic energy density. The latter indicates that it is the magnitude of the flow rather than the kinetic energy density that is important to the SSD efficiency.

We break down the statistics into epochs of various rates of volume averaged growth ${\overline{\gamma}}$. Periods of higher ${\overline{\gamma}}$ have higher velocities and vorticity. From the stronger correlation between vorticity and growth, we conclude that the vorticity of the flow drives the SSD. Despite the compressive structure of the SN forcing there is, as expected, a high fraction of rotational flow present, as identified from the Helmholtz decomposition. The excess of rotational to compressive flow is indicated by the compressive ratio at less than 40% and even below 10% in the hot gas. The efficiency of vortex generation is primarily driven in the ISM by the strong baroclinicity, the cross of orthogonal gradients of pressure and density, at the SN shock interactions (see Käpylä et al., 2018). Typically the magnitude of the velocity and the proportion of enstrophy increase with temperature. Fitting data from Figure 4(a) and for the high resolution data only from Figure 7(c), we find that the local growth rate $\langle\Gamma\rangle\propto\langle|\omega|\rangle$ and $\langle\Gamma\rangle\propto\langle{\rm Re}\rangle$.

A concern might exist that vorticity could be artificially generated, as has been identified at the intersection of adaptive grids of mixed mesh refinement (Plewa & Müller, 2001), or due to the use of locally varying diffusive schemes (Robertson et al., 2010) at lower resolution. Diffusivity is independent of position in our models, except for its use in resolving shocks. Inspection of the energy spectra (as in Gent et al., 2021, and here in Figure 12 of the Appendix) confirms in all cases that energy peaks are far from the Nyquist frequency and with very low energy near theses scales. Hence, any transfer of energy from the flow to the magnetic field must occur at length scales independent of the shock resolution over a few cells.

The study of vortex generation by Käpylä et al. (2018) was a purely hydrodynamical study, although stratified, with similar levels of vortex generation. It had no hyperdiffusion. Given the relatively low magnetic energy here, even after saturation of the SSD, we can also exclude the generation of vorticity being due to the magnetic field.

We have shown that SSD activity increases when the fractional volume of the hot gas increases, and we have shown that the magnetic energy density in the hot gas grows the most rapidly during these epochs. At other times, though, the magnetic field does not grow as quickly in the hot gas as in the warm gas. The vorticity in the hot gas is also particularly high during these rapid bursts of SSD. So, both the volume and the stirring of the hot gas increase, combining to drive strong dynamo action for a limited period.

In Kirchschlager et al. (2022), we use model UP0r10sl to examine dust-processing due to SN blast waves. The progress of individual remnants are followed in detail for isolated explosions in a dense turbulent region or a diffuse turbulent region. The effects on the magnetic field inside each remnant is instructive for the understanding of SSD in hot gas (see their Figure 1). In the case of the explosion in a dense region, the magnetic field is evacuated from within the remnant and is swept along by the the blast wave for up to 1 megayear. In the case of the diffuse region, the surrounding regions of strong magnetic field are pushed away by the blast wave, but inside the remnant the magnetic field grows, with the interior field strength orders of magnitude stronger at 1 megayear than at the time of the explosion. The SSD is enhanced in the hot gas where the interaction of the blast wave with the inhomogeneous density structure excites strong, interacting reverse shocks that can persist to late times in the diffuse medium. In dense regions, on the other hand, the reverse shocks from the ambient ISM are weak relative to the blast wave and then quickly slow in the relatively high gas density remnant interior.

This is consistent with our finding that the SSD is more sensitive to the vorticity than the kinetic energy, as the diffuse remnant contains relatively low kinetic energy. In our stratified model LP5r8sm-S the SSD is more active, with the ISM typically more diffuse. The SSD accelerates when the hot gas develops more solenoidal turbulence through SN remnant interactions and mergers (such as emerging superbubbles in the local and nearby galactic disks Bik et al., 2018; Sofue & Kataoka, 2021). Simulated superbubbles with magnetic fields (e.g., Ferrière et al., 1991; Stil et al., 2009) have tended to apply uniform ambient fields and smooth ambient gas, such that amplification of the magnetic fields are concentrated around the outer shell due to compression, as in the case of the dense region. Superbubbles in turbulent MHD have been modelled (e.g., Korpi et al., 1999; de Avillez & Breitschwerdt, 2005; Breitschwerdt & de Avillez, 2006; Gressel et al., 2008; Gent et al., 2013a, b), but the local structure of the magnetic field not analysed in such detail.

We conclude that the SSD in the SN-driven multiphase turbulence of the ISM:

1. a.

   has an average growth rate linearly correlated with the average vorticity of the flow, which is efficiently generated by baroclinic effects due to SN blast wave interactions;

2. b.

   correspondingly has an average growth rate well correlated with the average magnitude of the rotational flow velocity, as well as the compressible and total velocities and the fluid Reynolds number, and anticorrelates with the compressive ratio;

3. c.

   grows as a superposition of varying dynamo modes in each phase and within different regions of each phase;

4. d.

   grows intermittently, predominantly during periods when the hot gas, with low sonic Mach number $\cal{M}$, has strong enstrophy and large volume filling fraction;

5. e.

   is thus most efficient in hot gas produced when supernovae explode in diffuse regions, where interacting reverse shocks can have high velocities;

6. f.

   is more efficient in a stratified disk, where diffuse gas away from the disk is more easily heated and supports high velocities with low $\cal{M}$;

7. g.

   and already appears to approach asymptotic solutions at values of ${\rm Rm}$, ${\rm Re}$ and ${\rm Pm}$ significantly below the actual values estimated to apply in the ISM, although higher resolution will be required to confirm this point (see Appendix).

Prior to the onset of SSD the magnetic field strength aligns with the gas density, consistent with turbulent mixing, and following saturation of the SSD the field redistributes back into alignment with the gas concentration. The hot gas saturates significantly below equipartition, and even below the the ISM average of around 5% of equipartition with kinetic energy density, while the cold gas saturates at close to the equipartition energy.