The Small-scale Dynamo in a Multiphase Supernova-driven Medium

We use the Pencil Code (Brandenburg & Dobler 2002; The Pencil Code Collaboration et al. 2021) 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

$$\begin{eqnarray}&&\displaystyle \frac{D\rho }{{Dt}}=-\rho {\rm{\nabla }}\cdot {\boldsymbol{u}}+{\rm{\nabla }}\cdot {\zeta }_{D}{\rm{\nabla }}\rho ,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\begin{array}{c}\rho \frac{D{\boldsymbol{u}}}{{Dt}}={\rm{\nabla }}{E}_{\mathrm{kin}}\sigma -\rho \,{c}_{s}^{2}{\rm{\nabla }}\left(s/{c}_{p}+\mathrm{ln}\rho \right)+{\boldsymbol{j}}\times {\boldsymbol{B}}\\ \,+\,{\rm{\nabla }}\cdot \left(2\rho \nu {\mathsf{W}}\right)+\rho \,{\rm{\nabla }}\left({\zeta }_{\nu }{\rm{\nabla }}\cdot {\boldsymbol{u}}\right)\\ \,+{\rm{\nabla }}\cdot \left(2\rho {\nu }_{6}{{\mathsf{W}}}^{(6)}\right)-{\boldsymbol{u}}{\boldsymbol{\nabla }}\cdot \left({\zeta }_{D}{\boldsymbol{\nabla }}\rho \right),\end{array}\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{A}}}{\partial t}={\boldsymbol{u}}\times {\boldsymbol{B}}+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{A}}+{\eta }_{6}{{\rm{\nabla }}}^{6}{\boldsymbol{A}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\begin{array}{c}\rho T\frac{{Ds}}{{Dt}}={E}_{\mathrm{th}}\dot{\sigma }+\rho {{\rm{\Gamma }}}_{\mathrm{UV}}-{\rho }^{2}{\rm{\Lambda }}+\eta {\mu }_{0}{{\boldsymbol{j}}}^{2}\\ \,+\,2\rho \nu {\left|{\mathsf{W}}\right|}^{2}+\rho \,{\zeta }_{\nu }{\left({\rm{\nabla }}\cdot {\boldsymbol{u}}\right)}^{2}\\ \,+\,{\rm{\nabla }}\cdot \left({\zeta }_{\chi }\rho T{\rm{\nabla }}s\right)+\rho T{\chi }_{6}{{\rm{\nabla }}}^{6}s\\ \,-\,{c}_{v}\,T\left({\zeta }_{D}{{\rm{\nabla }}}^{2}\rho +{\rm{\nabla }}{\zeta }_{D}\cdot {\rm{\nabla }}\rho \right),\end{array}\end{eqnarray} \tag{ 4 }$$

with the ideal gas equation of state closing the system, assuming an adiabatic index (ratio of specific heats) c _p /c _v = 5/3. Most variables take their usual meanings; a list of notations used is given in Table 2.

Terms containing ν₆, χ₆, and η₆ 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 ≃1 at the scale of the zone size δx. Terms including ζ_D , ζ_ν , and ζ_χ 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 ζ_D applying in Equation (1). Following Gent et al. (2021), resistive shock diffusion is omitted from Equation (3).

SNe at rate $\dot{\sigma }$ inject E_th = 10⁵¹ erg thermal energy. In dense regions up to 5% is instead injected as kinetic energy E_kin (see Kim & Ostriker 2015; Gent et al. 2020). Nonadiabatic heating Γ_UV and cooling Λ(T) (as detailed in Gent et al. 2013a) follow Wolfire et al. (1995) and Sarazin & White (1987). Some symbols used here and below are listed in Table 2 to ease reference.

We simulate a nonstratified, nonrotating domain with initial gas number density n = 1 cm⁻³ covering (256 pc)³ and periodic across all boundaries. Resolution along each edge spans 0.5−4 pc, corresponding to grid sizes of 64³−512³ zones. Each model SN rate is given relative to the estimated rate in the solar neighborhood of the Milky Way, ${\dot{\sigma }}_{\mathrm{Sn}}=50\,\,{\mathrm{Myr}}^{-1}\,{\mathrm{kpc}}^{-3}$. SNe occur uniformly and 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 ≃ 1 cm⁻³ and $\dot{\sigma }\simeq {\dot{\sigma }}_{\mathrm{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 two to four 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 cm⁻³ 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 ≃ 1 cm⁻³ at the midplane. With this we examine the effect that vertical expansion of SN remnants or superbubbles has on the SSD. Here Γ_UV is amplified to 3.5× 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 pc along the disk and ±1.534 kpc perpendicular to it. SNe are randomly located uniformly in the horizontal plane and normally in the vertical direction  with scale heights of 90 and 325 pc for Type I and Type II SNe, respectively (Ferrière 2001). This model otherwise is as described in Gent et al. (2013a, 2013b).

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

Table 1. Models Included

Model	Resolution	$\dot{\sigma }$	η	ν	η3, ν3
	(pc)	(${\dot{\sigma }}_{\mathrm{Sn}})$	(kpc km s−1)	(kpc km s−1)
UP0r10sl	0.5	0.2	1 (−3)	0.0	1.6 (−16)
HP0r10sl	1.0	0.2	1 (−3)	0.0	3.5 (−15)
HP1r10sl	1.0	0.2	1 (−3)	1 (−3)	2.5 (−15)
HP2r5sl	1.0	0.2	5 (−4)	1 (−3)	2.5 (−15)
HP5r2sl	1.0	0.2	2 (−4)	1 (−3)	2.5 (−15)
HP5r8sm	1.0	0.5	8 (−4)	4 (−3)	2 (−15)
MP0r10sl	2.0	0.2	1 (−3)	0.0	8.25 (−14)
MP2r5sl	2.0	0.2	5 (−4)	1 (−3)	8.25 (−14)
MP0r10sh	2.0	1.0	1 (−3)	0.0	8.25 (−14)
MP0r5sh	2.0	1.0	5 (−4)	0.0	8.25 (−14)
LP0r10sl	4.0	0.2	1 (−3)	0.0	2 (−12)
LP2r5sl	4.0	0.2	5 (−4)	1 (−3)	2 (−12)
LP5r8sm-S	4.0	0.5	8 (−4)	4 (−3)	1 (−11)
LP0r10sh	4.0	1.0	1 (−3)	0.0	2 (−12)
LP0r5sh	4.0	1.0	5 (−4)	0.0	2 (−12)

Note. Prefixes "U," "H," "M," and "L" refer, respectively, to ultrahigh, high, medium, and low resolution. Numbers after "P" indicate the nominal magnetic Prandtl number ν/η. The "r" numbers indicate the resistivity coefficient η in units of 10⁻⁴ kpc km s⁻¹. The SN rates "sh," "sm," and "sl," denote 1, 0.5, and 0.2 times the solar neighborhood rate, respectively. "S" denotes the stratified model. All models use coefficients 2, 5, and 2 for shock handling terms with ζ_D , ζ_ν , and ζ_χ , respectively.

Download table as:  ASCIITypeset image

Angular brackets indicate that 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, where the intervals may vary, as explained in the text.

In the case of $\bar{{e}_{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 〈e _B 〉 is then normalized by $\bar{{e}_{K}}$ to ease model comparison.

The erratic growth of the volume-averaged magnetic energy shown in Figure 1 and discussed in Section 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

$$\begin{eqnarray}&&\langle {e}_{B}\rangle (t)=\langle {e}_{B}\rangle ({t}_{0})\exp [\bar{\gamma }(t-{t}_{0})].\end{eqnarray} \tag{ 5 }$$

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 behavior 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 Equation (3) to identify instantaneous changes to the magnetic energy. Hyperdiffusion, being purely numerical, is neglected. Taking the curl of Equation (3) and contracting it with ${\boldsymbol{B}}{\mu }_{0}^{-1}$, we obtain an equation for the change of the magnetic energy density

$$\begin{eqnarray}&&\frac{\partial {e}_{B}}{\partial t}={\boldsymbol{B}}\cdot {\rm{\nabla }}\times \left({\boldsymbol{u}}\times {\boldsymbol{B}}+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{A}}\right){\mu }_{0}^{-1}.\end{eqnarray} \tag{ 6 }$$

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

$$\begin{eqnarray}&&{\rm{\Gamma }}({\boldsymbol{x}},t)=\frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}\times \left({\boldsymbol{u}}\times {\boldsymbol{B}}+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{A}}\right)}{{e}_{B}{\mu }_{0}}.\end{eqnarray} \tag{ 7 }$$

Here Γ( 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 Γ in relation to various physical properties will help determine how the varying growth depends on the multiphase structure of the ISM.

The growth rate Γ 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

$$\begin{eqnarray}&&\tilde{{\rm{\Gamma }}}({\boldsymbol{x}},t)=\frac{{\boldsymbol{B}}\cdot {\rm{\nabla }}\times \left({\boldsymbol{u}}\times {\boldsymbol{B}}+\eta {{\rm{\nabla }}}^{2}{\boldsymbol{A}}\right)}{\langle {e}_{B}\rangle (t){\mu }_{0}},\end{eqnarray} \tag{ 8 }$$

replacing e _B in the denominator with the time-dependent volume-averaged 〈e _B 〉(t). Rescaling by time-dependent 〈e _B 〉 assists comparison between all stages of the SSD. Both Γ and $\tilde{{\rm{\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,

$$\begin{eqnarray}&&\mathrm{Rm}({\boldsymbol{x}},t)=\displaystyle \frac{| {\rm{\nabla }}\times {\boldsymbol{u}}\times {\boldsymbol{B}}| }{| \eta {\mu }_{0}{\rm{\nabla }}\times {\boldsymbol{j}}| },\end{eqnarray} \tag{ 9 }$$

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

Table 2. Meanings of Variables

Symbol	Denoting	Units/Definition
e B	magnetic energy density	erg cm−3
$\bar{{e}_{K}}$	time-averaged kinetic energy density	erg cm−3
$\overline{\gamma }$	volume-averaged 〈e B 〉 growth rate	Gyr−1
Γ	relative e B growth rate	Gyr−1
$\tilde{{\rm{\Gamma }}}$	absolute e B growth rate	Gyr−1
$\dot{\sigma }$	SN explosion rate	kpc−3 Myr−1
$\displaystyle \frac{D}{{Dt}}$	material derivative	$\displaystyle \frac{\partial }{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}$
∇	gradient vector	e.g., $\left(\displaystyle \frac{\partial }{\partial x},\displaystyle \frac{\partial }{\partial y},\displaystyle \frac{\partial }{\partial z}\right)$
ρ	gas density	g cm−3
u	gas velocity	km s−1
t	time	Myr
s	specific entropy	erg g−1 K−1
T	gas temperature	K
A	magnetic vector potential	G cm
B	magnetic field	G
j	current density	erg cm−4 G−1
W	traceless rate of strain tensor	${{\mathsf{W}}}_{{ij}}=\displaystyle \frac{1}{2}\left(\displaystyle \frac{\partial {u}_{i}}{\partial {x}_{j}}+\displaystyle \frac{\partial {u}_{j}}{\partial {x}_{i}}-\displaystyle \frac{2}{3}{\delta }_{{ij}}{\rm{\nabla }}\cdot {\boldsymbol{u}}\right)$
∣W∣2	contraction of W	∣W∣2 = W ij W ij
W (6)	sixth-order rate of strain tensor	${{\mathsf{W}}}_{{ij}}^{(6)}=\displaystyle \frac{1}{2}\left(\displaystyle \frac{{\partial }^{5}{u}_{j}}{\partial {x}_{i}^{5}}+\displaystyle \frac{{\partial }^{4}}{\partial {x}_{i}^{4}}\left(\displaystyle \frac{\partial {u}_{i}}{\partial {x}_{j}}\right)-\displaystyle \frac{1}{3}\displaystyle \frac{{\partial }^{4}}{\partial {x}_{i}^{4}}\left({\rm{\nabla }}\cdot {\boldsymbol{u}}\right)\right)$
ζD , ζν , ζχ	shock diffusion coefficients	$\propto {\left(-{\rm{\nabla }}\cdot {\boldsymbol{u}}\right)}_{+}^{2}$
ν, η	viscosity, resistivity coefficients	kpc km s−1
ν6, χ6, η6	hyperdiffusion coefficients	kpc5 km s−1
ΓUV	UV heating	erg g−1 s−1
Λ	radiative cooling	erg cm3 g−2 s−1
Ekin+Eth	SN explosion energy	1051erg
μ0	vacuum magnetic permeability	1
c s	sound speed	km s−1
c p	specific heat at constant pressure	erg g−1 K−1
c v	specific heat at constant volume	erg g−1 K−1

Download table as:  ASCIITypeset image

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

$$\begin{eqnarray}&&{f}_{V,i}=\displaystyle \frac{{V}_{i}}{V},\end{eqnarray} \tag{ 10 }$$

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 Γ 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 γ(t) as the rate of change in 〈e _B 〉 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 γ(t) is lowest, median, and highest. ⁶

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 $\mathrm{log}{\rm{\Gamma }}$ or $\mathrm{log}\tilde{{\rm{\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{{\rm{\Gamma }}}$ and (twice the) enstrophy or the norm squared vorticity ∣ω∣² 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.

Zoom In Zoom Out Reset image size

The most striking correlation we find is between relative growth rate Γ and the norm squared vorticity ∣ω∣², which we show in Figure 4(a). We find that the growth rate $\mathrm{log}{\rm{\Gamma }}$ is strongly proportional to $\mathrm{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 γ(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 Γ and γ are correlated with increased vorticity in the hot gas, consistent with results from models using more idealized turbulence (Federrath et al. 2011; Achikanath Chirakkara et al. 2021).

Zoom In Zoom Out Reset image size

To confirm whether this correlation is reflected in the absolute gains in magnetic energy, we plot the same summary for $\tilde{{\rm{\Gamma }}}$ in Figure 4(b). At high resolution the same trends are preserved as for Γ, with the changes in $\tilde{{\rm{\Gamma }}}$ most associated with changes in the characteristics of the hot gas. However, for low resolution $\tilde{{\rm{\Gamma }}}$ is actually anticorrelated 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{{\rm{\Gamma }}}$, consistent with results from isothermal modeling, except at δ x = 2 pc for the upper bin of γ(t) (right), when vorticity and fractional volume increase 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 ${ \mathcal M }$ against relative growth rate Γ for each bin of volume-averaged magnetic energy growth rate. Overall the expected trend of Γ reducing as ${ \mathcal 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 ${ \mathcal M }$ and high relative growth rate Γ and cold gas the inverse. What is, perhaps, unexpected is to see that within each phase there is a trend of increasing Γ with increasing ${ \mathcal M }$. This can be explained by the countereffect of increased velocities, which drive higher Mach numbers for a given sound speed, but also higher Rm. So increased Γ with ${ \mathcal M }$ should have complementary trends in Rm.

Zoom In Zoom Out Reset image size

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 Equation (9). With strongly fluctuating characteristic velocities and length scales applying between the continuously shifting phases, as well as 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, respectively, there are regions in all phases where Rm > 10³, and in the hot phase Rm > 10⁴. 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.

Zoom In Zoom Out Reset image size

Results for Rm versus Γ are shown in Figure 5(b). Velocities increase with temperature, increasing Rm and supporting faster SSD, an overall trend as would be expected for Γ. Within each phase and within each γ(t) bin, however, Γ is weakly anticorrelated with Rm, and across each bin the patterns are largely unaffected for cold and warm gas, suggesting that higher velocities may be associated with higher Mach numbers. Results are similar for Re versus Γ as shown in Figure 5(c) for samples from models that include explicit viscosity ν. The higher-resolution models (δ 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 underresolved density gradients. The comparison between the 1 and 0.5 pc models without viscosity suggests that the 1 pc model is well converged. For the hot gas Rm, Re, and Γ all increase for the upper growth rate bin, suggesting that 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₁ + u₂ into purely solenoidal u₁ = ∇ × F and compressive flows u₂ = ∇Φ, in which F and Φ are, respectively, a vector potential and scalar potential pair from which orthogonal flows u₁ and u₂ can be derived.

Zoom In Zoom Out Reset image size

Both flows appear strongly correlated with Γ, 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 Γ 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 γ(t) even as flow strength varies. In the hot gas Γ clearly increases when ∣∇ × F ∣ increases between low (left) and high (right) γ(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

$$\begin{eqnarray}&&{R}_{\mathrm{cs}}=\displaystyle \frac{\langle | {\rm{\nabla }}\cdot {\boldsymbol{u}}{| }^{2}\rangle }{\langle | {\rm{\nabla }}\cdot {\boldsymbol{u}}{| }^{2}\rangle +\langle | {\rm{\nabla }}\times {\boldsymbol{u}}{| }^{2}\rangle }.\end{eqnarray} \tag{ 11 }$$

This shows that the ratio of solenoidal energy to total compressive and solenoidal energy (1 − R_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 cause 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 ν and resistivity η. For a fixed η, the SSD did not depend on Pm, and the energy at which the dynamo saturated was consistently around 5% of equipartition. For a fixed ν, on the other hand, SSD growth rate did increase with Pm for Pm ≤ 5 in the parameter space considered. The energy at saturation increased with Pm, remaining steady for Pm ≳ 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 η. Where the viscosity, ν, and the resistivity, η, are explicitly defined, we can nominally express the magnetic Prandtl number, Pm = ν/η, which is what we shall use here. Alternatively, as a function of position Pm can also vary owing 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 Pm_eff in various models to examine the relevance of the explicit diffusivities to the dynamics (see Figure 12).

In Figure 8(a) we show the volume-averaged magnetic energy density 〈e _B 〉 for models with Pm ∈ [1, 5] with resolution δ x = 1 pc. Around 20 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 Myr < t < 20 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 Pm = 1 at 41 Myr < t < 44 Myr, where the growth rate exceeds that for Pm = 5 earlier. Under most circumstances we would expect $\overline{\gamma }$ to be lower for the larger η, 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.

Zoom In Zoom Out Reset image size

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 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 Pm ∈ [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 Myr. The resistive time step differs between the models; small changes in time step 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 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 Pm = 5. The value of f_V,h then rises again at 40 Myr, particularly for 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 〈e _B 〉_T and plot these in Figure 8(c) for the models with Pm = 1 and 5. These are normalized by $\bar{{e}_{K}}$. Horizontal dotted lines of matching color for each phase indicate the time-averaged kinetic energy density by phase ${\bar{{e}_{K}}}_{T}/\bar{{e}_{K}}$. These lines, plotted for Pm = 5, are very similar to those for Pm = 1.

Usually 〈e _B 〉_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{\bar{{e}_{K}}}_{h}$ at 40 Myr for Pm = 5 and 55 Myr for Pm = 1. In contrast, the saturation energy in the cold and warm phases is affected by Pm. In the nonlinear phase the magnetic energy in the hot phase 〈e _B 〉_h decays to around $0.01{\bar{{e}_{K}}}_{h}$, while the magnetic energy in the cold phase 〈e _B 〉_c grows, above $\bar{{e}_{K}}$ _c in the case of Pm = 5.

Is the latter consistent with compression in the statistical steady MHD state? Under compression, magnetic field strength relates to density ∣B∣ ∝ ρ^α , where α ∈ [0, 1] for compression along field lines up to disk 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 rms strength of the magnetic field normalized by cold gas density $\langle {b}_{\mathrm{rms}}{\rangle }_{C}=\langle {B}^{2}{\rangle }_{C}^{1/2}\langle \rho {\rangle }_{C}^{-\alpha }$ (right axis), taking α = 1/2. In the interval 50 Myr < t < 60 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 ∼65 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 < α < 1 in the analysis does not qualitatively alter the development of b.

Zoom In Zoom Out Reset image size

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\leqslant \dot{\sigma }\leqslant 40{\dot{\sigma }}_{\mathrm{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 (δ x ≤ 1 pc) even with $\dot{\sigma }={\dot{\sigma }}_{\mathrm{Sn}}$. At low resolution (δ x = 4 pc), with consequentially higher numerical diffusion, and thus unphysically strong cooling, thermal runaway was still excited for $\dot{\sigma }=10{\dot{\sigma }}_{\mathrm{Sn}}$. Rapid steady acceleration of the SSD occurs, which is excluded here with lower $\dot{\sigma }$.

In Figure 10(a) we compare evolution of 〈e _B 〉 for $\dot{\sigma }=0.2{\dot{\sigma }}_{\mathrm{Sn}}$ and $1{\dot{\sigma }}_{\mathrm{Sn}}$ at δ x = 2 and 4 pc, and in Figure 10(b) for these models we compare 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 δ x = 2 pc (blue) than δ x = 4 pc (brown). For lower $\dot{\sigma }$, 〈e _B 〉 decays, but growth occurs briefly in the subinterval between the first two vertical lines, while f_V,h increases, particularly for δ x = 2 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 δ x = 4 pc after 270 Myr, when f_V,h is at an even higher peak. In general, for δ x = 4 pc there is less hot gas than for δ x = 2 pc at each given SN rate because of higher numerical diffusion and cooling.

Zoom In Zoom Out Reset image size

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 nonperiodic boundaries vertically might still cause an LSD effect, so any residual mean magnetic field volume and horizontal component average is subtracted at each time step. 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 〈e _B 〉 for this model with $\dot{\sigma }=0.5{\dot{\sigma }}_{\mathrm{Sn}}$ alongside models with $\dot{\sigma }={\dot{\sigma }}_{\mathrm{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)).

Zoom In Zoom Out Reset image size

From the summary statistics in Figure 4(b), we find for δ x = 4 pc that the absolute growth $\tilde{{\rm{\Gamma }}}$ of magnetic energy is fastest in the cold gas but occurs mostly in the warm gas owing 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 sign.

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

Vorticity in all phases is much lower (Figure 4) at δ x = 4 pc, but of these it is highest in the stratified model. The values of ${ \mathcal M }$, Rm, and velocities (Figures 5 and 7) are also highest in the stratified model. Higher velocities, as well as 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_crit, the critical Rm above which SSD can be excited, reduces as Pm increases above Pm ≳ 0.1 with an asymptotic limit of order Rm_crit ∼ 20 for Pm ≫ 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 > Pm > 0.01 where Rm_crit has a maximum (Schekochihin et al. 2005; Iskakov et al. 2007; Brandenburg 2011), and then for Pm < 0.01 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 ${ \mathcal M }$ increases it becomes more difficult to excite SSD (Rm ≳ 60 for Pm ≥ 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 ${ \mathcal 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 ${ \mathcal M }$, seen in Figure 5(a), are due to the SSD supportive countertrends, such as increasing velocity (with ${ \mathcal M }$). The expected effect of ${ \mathcal 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 ${ \mathcal M }$ ∼ 1 for the warm gas and 5 for the cold gas. Their SSD saturates after around 400 Myr for solenoidal forcing and 800 Myr 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 Myr with mean ${ \mathcal M }$ ∼ 1, 0.6, and 0.25 in the cold, warm, and hot gas, respectively. However, in all phases ${ \mathcal 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 ${ \mathcal 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 (Figure 8(c)). We speculate that an explanation for the growth rates differing between phases in three-phase simulations and not in two-phase simulations is that high vorticity generated in the hot gas transfers to warm and 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 range from 10⁻⁶⁵ to 10⁻²⁰ μG (Grasso & Rubinstein 2001; Subramanian 2016). With an exponent of $\overline{\gamma }=500\,{\mathrm{Gyr}}^{-1}$ this could be amplified to sub-equipartition saturation of around 1 μG within 170–590 Myr. Given that 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 the current understanding of dynamo, to see in Figure 6 locations with Rm ≤ 1, well below isothermal predictions of Rm_crit, coincident with growth rate $\tilde{{\rm{\Gamma }}}\gt {10}^{3}\,{\mathrm{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{{\rm{\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 (Taylor 1950; Sedov 1959) traveling at 150 kms⁻¹, as would occur in an SN remnant with radius near 30 pc in a 10⁴ K medium of uniform gas number density 1 cm⁻³. From Equation (8) with the normalized quantity $\hat{{\boldsymbol{B}}}={\boldsymbol{B}}\langle {\mu }_{0}{e}_{B}{\rangle }^{-1/2}$ we can write

$$\begin{eqnarray}\begin{array}{rcl}\tilde{{\rm{\Gamma }}} & = & -{\hat{B}}^{2}{\rm{\nabla }}\cdot {\boldsymbol{u}}+\hat{{\boldsymbol{B}}}\cdot (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}){\boldsymbol{u}}\\ & & -\,\hat{{\boldsymbol{B}}}\cdot ({\boldsymbol{u}}\cdot {\rm{\nabla }})\hat{{\boldsymbol{B}}}+\eta \hat{{\boldsymbol{B}}}\cdot {{\rm{\nabla }}}^{2}\hat{{\boldsymbol{B}}}.\end{array}\end{eqnarray} \tag{ 12 }$$

Considering the first term with a shock front traveling at 150 kms⁻¹ into warm gas with sound speed around 12 kms⁻¹ resolved in our model over three cells, −∇ · u ≈ 138 kms⁻¹(3 pc)⁻¹ ≈ 4.6 × 10⁴ Gyr⁻¹. If the local magnetic field is even twice the volume-averaged strength, this would yield $\tilde{{\rm{\Gamma }}}\gtrsim {10}^{5}\,{\mathrm{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 Γ and $\tilde{{\rm{\Gamma }}}$ in regions of low Rm in Figure 6.

For a magnetic field perpendicular to the blast wave

$$\begin{eqnarray}&&\mathrm{Rm}\sim \displaystyle \frac{| {\partial }_{\perp }{u}_{\perp }{B}_{\parallel }| }{| \eta {\partial }_{\perp }^{2}{B}_{\parallel }| },\end{eqnarray} \tag{ 13 }$$

reaching a maximum at the shock center, where the first derivative approaches infinity and the second derivative vanishes, but near zero on either side of the shock, where the inverse is true. Both can coexist with the compression $\tilde{{\rm{\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 localized 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 owing to heating, cooling, and phase mixing.

The mean Γ of the instantaneous localized growth rates collected while γ(t) is highest have similar magnitude to $\overline{\gamma }$ as fitted in each simulation, supporting the concept of a superposition. Sufficiently high mean Γ or $\tilde{{\rm{\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 Γ and $\tilde{{\rm{\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 〈Γ〉 ∝ 〈∣ω∣〉 and $\langle {\rm{\Gamma }}\rangle \propto \langle \mathrm{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 these 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 is 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 are 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 blast wave for up to 1 Myr. 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 Myr 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 is concentrated around the outer shell owing to compression, as in the case of the dense region. Superbubbles in turbulent MHD have been modeled (e.g., Korpi et al. 1999; de Avillez & Breitschwerdt 2005; Breitschwerdt & de Avillez 2006; Gressel et al. 2008; Gent et al. 2013a, 2013b), but the local structure of the magnetic field is not analyzed in such detail.

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

• 1.  
  It 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.  
  It 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.  
  It grows as a superposition of varying dynamo modes in each phase and within different regions of each phase.
• 4.  
  It grows intermittently, predominantly during periods when the hot gas, with low sonic Mach number ${ \mathcal M }$, has strong enstrophy and large volume filling fraction.
• 5.  
  It is thus most efficient in hot gas produced when SNe explode in diffuse regions, where interacting reverse shocks can have high velocities.
• 6.  
  It is more efficient in a stratified disk, where diffuse gas away from the disk is more easily heated and supports high velocities with low ${ \mathcal M }$.
• 7.  
  It already appears to approach asymptotic solutions at values of Rm, Re, and Pm significantly below the actual values estimated to apply in the ISM, although higher resolution will be required to confirm this point (see the 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 ISM average of around 5% of equipartition with kinetic energy density, while the cold gas saturates at close to the equipartition energy.