Small-Scale Dynamo in Supernova-Driven Interstellar Turbulence

Magnetic fields grow quickly even at early cosmological times, suggesting the action of a small-scale dynamo (SSD) in the interstellar medium (ISM) of galaxies. Many studies have focused on idealized turbulent driving of the SSD. We here simulate more realistic supernova-driven turbulence to determine whether it can drive an SSD. We vary the physical resistivity (and implicitly magnetic Reynolds number), as well as the numerical resolution and supernova rate to delineate the regime in which an SSD occurs. We find convergence for SSD growth rate with resolution approaching sub-parsec scale for a given supernova distribution {\sigma}. Despite higher Mach numbers and negative impact of compressibility expected on SSD, growth rates increase for {\sigma}={\sigma}_{\rm sn} versus 0.2{\sigma}_{\rm sn}, with {\sigma}_{\rm sn} the solar neighbourhood rate. Across the modelled range of 0.5 to 4 parsec resolution we find that for sufficiently low resistivity the SSD saturates consistently at about 5% of energy equipartion, independent of growth rate and low Prandtl numbers in our experiments. As the grid becomes coarser, the minimum resolved physical resistivity increases. The trend suggests that numerical resistivity suppresses SSD for grid spacing much exceeding 4 pc.


INTRODUCTION
Corresponding author: Maarit Käpylä Email: frederick.gent@aalto. fi, mordecai@amnh.org, maarit.kapyla@aalto.fi, nishant@iucaa.in This letter addresses the necessary conditions and properties of a small-scale dynamo (SSD) in the interstellar medium (ISM). SSD acts at small eddy scales of the turbulence, thus driving magnetic field growth at correspondingly short turnover times. The fastest growing SSD modes are far smaller than the large-scale dynamo (LSD) modes generating the magnetic field structures organised at the systemic scale of the galactic disk. Hence, simulations capable of capturing LSD alongside the faster growing modes of SSD are computationally challenging. However, the interaction of SSD modes with the LSD likely fundamentally alters the evolution and structure of the magnetic field.
Many simulations of supernova-(SN)-driven turbulence with realistic vertical stratification (e.g., de Avillez & Breitschwerdt 2005;Piontek & Ostriker 2007;Hill et al. 2012;Hennebelle & Iffrig 2014) have no mechanisms for inducing LSD, such as rotation and shear. The effect of strong ordered magnetic fields is modelled by initial imposition of a background, typically uniform, magnetic field. If an imposed field is sufficiently close to equipartition, it would excessively influence the results. Some large-scale models do seek to include LSD (e.g., Korpi et al. 1999;Gressel et al. 2008;Hanasz et al. 2009;Wang & Abel 2009;Pakmor et al. 2017;Steinwandel et al. 2019;Gressel & Elstner 2020), but most show no SSD. Steinwandel et al. (2019) do find an SSD, but no LSD. Gent et al. (2013a); Evirgen et al. (2017) do appear to include an SSD. To confirm this and determine its effect on LSD, we must understand the properties of the SSD.
Any magnetic noise produced by tangling will also grow exponentially due to an LSD if present. This noise plays an important role in quenching the LSD due to the magnetic α-effect, causing saturation of large-scale fields. We need to discriminate this effect from an SSD.
Previous experiments (e.g., Balsara et al. 2004;Mac Low et al. 2005) examined the SN-driven SSD using ideal magnetohydrodynamics (MHD) with a limited set of resolution and parameters that did not allow demonstration of convergence of the solutions nor dependence on magnetic Reynolds (Rm) or Prandtl (Pm) numbers. They included a weak imposed uniform field; we shall show that the amplification of their field is a result of SSD action and not just tangling of the field.
In this letter we compare the SSD to tangling in an idealized simulation (Sect. 2), and then determine its presence in simulations of SN-driven turbulence in isolation from any drivers of an LSD (Sect. 3). Our numerical implementations use the Pencil Code 1 . A broad resolution and parameter study allows us to identify the critical resistivity and resolution for excitation of an SSD, and follow the growth to saturation (Sect. 4). This provides objective criteria with which to determine the presence of SSD in simulations (such as Gent et al. 2013a;Gressel & Elstner 2020;Steinwandel et al. 2019). Finally, we conclude in Sect. 5.

DISENTANGLING THE DYNAMO
To illustrate differences between tangling and SSD we adopt a simplified model.
Nonhelical random forcing is applied at wavenumber k f /k 1 = 8 to 256 3 zone, 2π-periodic, isothermal boxes. The lowest wavenumber in the domain is k 1 = 1 and the largest is the Nyquist frequency k/k 1 = 128. The imposed uniform field has e B 6 · 10 −22 e K , where e K is the timeaveraged kinetic energy density. Two simulations are distinguished by use of dimensionless resistivity η = 10 −4 and η = 2 · 10 −3 , and viscosity ν = 5 · 10 −3 . These yield Rm = 150, with Pm = 50, exciting SSD, and Rm = 7.4 with Pm = 2.5, inhibiting the dynamo so that amplification is limited to tangling of the imposed field. Figure 1(a) shows the SSD growing exponentially in just over 400 eddy turnover times; see Zeldovich et al. (1983) for SSD properties and excitation conditions. Tangling induces only linear growth (see inset), saturating just above the imposed field energy within 50 turnover times.
Power spectra for the magnetic energy of the evolving SSD and tangling model are plotted in Figures 1(b) and 1(c), respectively, alongside a late kinetic energy spectrum. Magnetic energy spectra are compensated by Kazantsev's k 3/2 scaling (Schekochihin et al. 2002;Bhat & Subramanian 2014), and the kinetic energy by the Kolmogorov's k −5/3 spectrum. The SSD magnetic spectrum (Figure 1b) experiences only a slight shift of its peak to smaller wavenumber during its evolution. The forcing scale k f is strongly reflected in kinetic energy and tan-gling magnetic energy, but does not affect the SSD magnetic spectra. The Kazantsev range extends larger than k f for SSD, while confined below k f with tangling. Thus, in the SSD, kinetic energy along the Kolmogorov range transfers to the magnetic field at these scales, inducing an inverse Kazantsev range at scales below k f , while tangling transfers energy only at scales between k f and the scale of the imposed field. Dissipation controls scales below the Kazantsev range.

SN TURBULENCE MODEL DESIGN
In our SN-driven turbulence models, we exclude large-scale magnetic field dynamics by neglecting rotation, shear, and stratification. Our simulation domain is a periodic cube of length 256 pc and zone size of δx = 0.5, 1, 2 or 4 pc. A random 1 nG seed field excludes tangling of an imposed field as the source of any magnetic amplification.
We solve the system of non-ideal compressible MHD equations ρT with the ideal gas equation of state closing the system. Most variables take their usual meanings. Terms containing ζ D , ζ ν and ζ η resolve shock discontinuities with artificial diffusion of mass, momentum, and energy proportional to shock strength (see Gent et al. 2020, for details). Unlike Gent et al. (2013a) shock diffusion is not applied to equation (4), to avoid excessive magnetic dissipation in shocks, which actually enhance the field by compression. Terms containing ν 3 , χ 3 and η 3 apply sixth-order hyperdiffusion to resolve grid-scale instabilities (see, e.g., Brandenburg & Sarson 2002;Haugen & Brandenburg 2004). In Sect. 2, we only solve equations (2) and (4); set ρ = 1; do not apply shock-dependent diffusion nor hyperdiffusion; and take B = ∇ × A + B imposed .
SNe are exploded at random positions at a Poisson rateσ normalized by the solar neighborhood valueσ sn 50 kpc −3 Myr −1 . Explosions inject E th = 10 51 erg of thermal energy, except in dense regions, where a proportion is kinetic E kin (see Gent et al. 2020). For comparability, explosions in all models with the sameσ occur at the same times and places. Non-adiabatic heating Γ and cooling Λ(T ) are included (Gent et al. 2013b) following Wolfire et al. (1995) and Sarazin & White (1987).
The subset of experiments presented in this letter have viscosity ν = 0, and hyperviscosity ν 3 set optimally for each δx to ensure the flow is well resolved. We benchmark the magnetic field evolution by setting η = 0, using only hyperresistivity η 3 . We determine how low a physical resistivity η can be resolved by varying it from 10 −5 to 10 −3 kpc km s −1 (units assumed henceforth). Unlike our earlier experiments (Gent et al. 2013b(Gent et al. ,a, 2020, we do not include thermal diffusivity, χ, as the artificial diffusivities are adequate to ensure numerical stability and physical effects of thermal conductivity can be expected to be relevant only at the unresolved or marginally resolved Field length defined by Begelman & McKee (1990, named after George Field, not the magnetic field).
A commonly used expression for Rm is where = 2π/k f is the forcing scale. However, the inhomogeneity of the turbulence and changing multiphase composition of the ISM render defining a single forcing scale or rms velocity unreliable. We, therefore, consider computing Rm at each grid point directly from the ratio of advective to diffusive terms in the induction equation Using averaged values (eq. [5]), Rm = 1.5 · 10 4 , for η = 10 −4 , = 50 pc and u rms = 30 km s −1 .
With the local definition (eq. [6]) in the model with δx = 0.5 pc and η = 10 −4 at 20 Myr we find Rm ∈ (1.0 · 10 −4 , 5.8 · 10 3 ), and mean value Rm = 4.1 ± 5.7. Thus, the mean Rm appears unhelpful: in the ISM some phases and regions may exceed the critical Rm for an SSD while others do not. Microscopic resistivity is temperature dependent (Cohen et al. 1950). Even considering only turbulent resistivity to have relevance, typical length scales and rms velocities differ between phases. Therefore, we use the explicitly chosen resistivity η to discriminate between models rather than Rm. We have set ν = 0 and ν 3 = η 3 sufficient to numerically resolve the grid scale. The magnetic Prandtl number Pm = Rm/Re, and like Rm we have Pm ∈ (4.2 · 10 −6 , 4.8 · 10 4 ), due to inclusion of shock capturing viscosity and differences in definition of hyperviscosity and hyperresistivity. Hence, some part of each simulation will be characterised by high Pm, but Pm < 1. This is a regime less conducive to exciting the SSD than the high Pm regime typical of the ISM (Haugen et al. 2004a). In separate experiments with ν > η (high Pm), which we shall describe in future work, we do in fact see in-creased growth rates for lower η, as expected (Schekochihin et al. 2007).

RESOLUTION AND RESISTIVITY
Figure 2(a) shows improving convergence of dynamo growth rate for η = 10 −4 as resolution increases, although even at 0.5 pc resolution full convergence does not appear to have been reached. Saturation at around 5% of e K appears to be a well-converged result. At η = 10 −3 , there is false convergence (Fryxell et al. 1991) of solutions with similar magnetic energy decay at δx = 2 and 4 pc. However, higher resolution solutions converge to rapid magnetic amplification occurring between 40 and 60 Myr.
For δx ≥ 2 pc growth rates vary sporadically. Growth is accelerated in Figure 2 Other times show decay or lower growth rates, depending on δx or η. The irregular SN forcing and varying fractional volume of different phases of the ISM account for these alternating regimes, but we defer analysis of the details to a future work. Figure 3 shows that profiles at η = 10 −5 are indistinguishable from η = 0, and so numerical resistivity still controls dynamics, even at δx = 0.5 pc. At η = 10 −4 the dynamo diverges from η = 0, showing that physical resistivity is dynamically dominant, at least in part of the  domain. At all δx in our study physical resistivity η = 10 −3 is clearly well resolved and determines the dynamics. Models with commonσ have the same schedule and location of SN, but the timestep at low resolution is larger, such that actual timing and explosive environment can differ between models. Increased statistical noise is therefore evident, particularly in panels (e) and (f). Growth rates are sporadic within models, but are consistent between models. For example, between 8 and 15 Myr the strength of growth (or decay) in Figure 3(a, b) reduces as η increases. Accelerated growth near 19 Myr for η = 10 −3 coincides with even higher acceleration for η ≤ 10 −4 , consistent with theory. The low η models at 0.5 pc saturate already by 20 Myr, but for 1 pc there is another epoch of lower growth rate up to 40 Myr and then a subsequent acceleration resulting in saturation for all models within 60 Myr, including for η = 10 −3 .
In Figure 3(c, d) at low resolution withσ = 0.2σ sn profiles for η = 10 −4 are not distinct from η = 0, and there is no dynamo at all for δx = 4 pc with well resolved η = 10 −3 . At 100 Myr there is a period of accelerated dynamo at δx = 2 pc, but for η = 10 −3 this is not sustained and subsequently diffuses away.
SN-driven turbulence does not have a welldefined forcing scale, unlike the simplified model in Section 2, because of the heterogeneous ISM structure and random explosions. The forcing scale will be distributed at scales greater than about 60 pc (Hollins et al. 2017 , Table 3), or k 17 kpc −1 . Figure 4 shows the evolving compensated spectra for δx = 0.5 pc between 9 and 32 Myr. We compare spectra for η = 10 −4 , (a, b), spanning the period of dynamo growth and saturation, and for η = 10 −3 , (c, d) during the period that the field persists near seed strength. The kinetic spectra over time match between models, until saturation of the dynamo at 32 Myr (b) diminishes its energy relative to (d). The Kolmogorov range extends at all times to k > 20 kpc −1 and mostly to k > 40 kpc −1 .
The compensated magnetic energy spectra in Figures 4(a) and (c) have ranges conforming to the Kazantsev inverse cascade. For η = 10 −4 this range extends to k 20 kpc −1 , above the forcing scale and consistent with SSD (Figure 1b), until it contracts upon saturation to k < 10 kpc −1 , consistent with no dynamo (Figure 1c). In Figure 4(c) the Kazantsev range occurs at k 10 except for the period of 19-22 Myr, which corresponds to a short growth spurt in Figure 2(b). The Kazantsev range does not extend to as high k as the Kolmogorov range, a consequence of the low Pm in these SN-driven models. In the high Pm ISM, where transfer from kinetic energy can occur at every wavenumber in the Kolmogorov spectrum, an SSD is even easier to excite. Figure 5(a) -(f) shows spectra at t = 19.5 Myr when an SSD is active for the higher resolution models, and at t = 100 Myr for the lower resolution models, which shows SSD for the δx = 2 pc model, but not for the 4 pc model. (Hydrodynamic parameters are fixed at each resolution.) The kinetic energy spectrum Figure 4. Compensated magnetic (a, c) and kinetic (b, d) power spectra for δx = 0.5 pc at times given in megayears by the legends. Resistivity is η = 10 −4 (a, b) or η = 10 −3 (c, d). Compensation is by the Kazantsev spectrum k 3/2 (a, c) or the Kolmogorov spectrum k −5/3 (b, d).
is convergent for δx = 0.5 pc and δx = 1 pc, except for a viscous cutoff at lower k for δx = 1. The kinetic energy at all k is reduced for lower resolution, not just at the viscous dissipation scale, but even at the largest scales. To fully resolve the kinetic energy for scales larger than k = 24 kpc −1 , that is 40 pc, therefore requires δx 1 pc.
The kinetic energy spectra display a bottleneck effect (Falkovich 1994;Haugen et al. 2003), an energy cascade less efficient than k −5/3 leading to an accumulation of power and then rapid dissipation at high k. This bottleneck shifts to lower k as δx increases, but always at higher k than the Kazantsev range in the corresponding magnetic spectrum, due to the low Pm.
Perhaps surprisingly, there is little difference in the kinetic energy spectrum, Figure 5(h), between simulations withσ = 0.2σ sn andσ sn , despite five times as much energy being applied to the forcing in the latter case. For η = 10 −3 there is more energy near the bottleneck in the higḣ σ case (cyan dotted line), so more energy can be tranferred to the SSD at the smallest scales. The comparison is obscured for η = 10 −4 , because the highσ magnetic field is already 3 or was evidence of an SSD in SN-driven ISM turbulence and not just caused by tangling of their imposed field. Through the most extensive resolution and parameter study to date, we show that SN-driven turbulence easily excites an SSD even at SN rates well below the Galactic value. Our conclusion is supported by noting that the resistivity of the ISM is far smaller than we can resolve numerically, so the ISM is far more susceptible to dynamo action than our models. Our models with δx = 0.5 and 1 pc with η = 10 −3 kpc km s −1 show that a seed field of less than 1 nG can be amplified to saturation at microgauss levels within about 10 Myr (Figure 2). Unlike isothermal turbulence high Mach numbers do not necessarily impede SSD in the ISM with increasing SN forcing. We further show that simulations with insufficient resolution can appear to converge to a false solution lacking dynamo activity (Figure 2b). This can occur because these simulations are not scale independent. The SN energy input and the physically motivated ISM cooling processes impose length and time scales that must be adequately resolved. To reach true convergence requires resolution of 1 pc or better. In our models, resolutions δx ≥ 2 pc not only give an incorrect dynamo solution, but also exhibit significant kinetic energy losses due to excess dissipation. This also affects the energy spectrum at the largest scales, so should be considered when interpreting results using adaptive mesh refinement. We do, however find the wellconverged result at all resolutions that when an SSD is excited it saturates at about 5% of the energy equipartition level. At low resolution this is a lower bound, because the timeaveraged kinetic energy density is understated, due to dissipative losses. This might account for the discrepancy between the energy density of the mean magnetic field and the random magnetic field in Gent et al. (2013a), due to the LSD being well resolved while the SSD remained underresolved.
We find that the conventional approach from dynamo theory of categorising the turbulence according to Rm based on a forcing scale , random velocity u and resistivity η is inadequate for such a complicated system. We will show in upcoming work that the ISM appears to have multiple regimes with thermal phases occupying changing fractional volumes (e.g. Gatto et al. 2015) and hosting SSD instabilities with different thresholds and growth rates. Explaining the interaction of these SSDs will require more sophisticated statistical and perhaps toplogical techniques, in advance of being able to address how such an SSD interacts with an LSD.
With grid resolution δx ≥ 2 pc, an SSD in SN-driven turbulence simulations can be excluded with η 10 −3 kpc km s −1 . In Gressel et al. (2008) and Gressel & Elstner (2020), δx is 8.3 and 6.7 pc, respectively and η 6.5 · 10 −3 kpc km s −1 , which appears to exclude an SSD. Gent et al. (2013a) applied η 8 · 10 −4 kpc km s −1 with δx = 4 pc, which would support SSD with SN rates similar to the solar neighbourhood. The latter obtain an LSD with galactic angular momentum Ω = Ω sn , where Ω sn = 25 km s −1 kpc −1 is the rate in the solar neighbourhood. The former require Ω ≥ 4Ω sn to excite an LSD. The value of Rm at the largest scales would be 7.5 times higher for the latter model, which alone can explain the LSD at lower Ω. Early efforts by Korpi et al. (1999) were unable to detect even LSD. With a resolution even less than Gressel et al. (2008) Ω Ω sn would be required to excite LSD, and SSD would be ruled out.