# Transition from small-scale to large-scale dynamo in a supernova-driven,
multiphase medium

Frederick A. Gent Astroinformatics, Department of Computer Science, Aalto
University, P.O. Box 15400, FI-00076 Espoo, Finland  School of Mathematics,
Statistics and Physics, Newcastle University, NE1 7RU, UK  Mordecai-Mark Mac
Low Department of Astrophysics, American Museum of Natural History, New York,
NY 10024, USA  Maarit J. Korpi-Lagg Astroinformatics, Department of Computer
Science, Aalto University, P.O. Box 15400, FI-00076 Espoo, Finland  Max Planck
Institute for Solar System Research, D-37707 Göttingen, Germany  Nordic
Institute for Theoretical Physics, SE-106 91 Stockholm, Sweden

###### Abstract

Magnetic fields are now widely recognised as critical at many scales to
galactic dynamics and structure, including multiphase pressure balance, dust
processing, and star formation. Using imposed magnetic fields cannot reliably
model the interstellar medium’s (ISM) dynamical structure nor phase
interactions. Dynamos must be modelled. ISM models exist of turbulent magnetic
fields using small-scale dynamo (SSD). Others model the large-scale dynamo
(LSD) organising magnetic fields at scale of the disc or spiral arms.
Separately, neither can fully describe the galactic magnetic field dynamics
nor topology. We model the LSD and SSD together at sufficient resolution to
use explicit Lagrangian resistivity and viscosity. The galactic SSD saturates
within 20 Myr. We show that the SSD is quite insensitive to the presence of an
LSD and is even stronger in the presence of a large-scale shear flow. The LSD
grows more slowly in the presence of SSD, saturating after 5 Gyr vs. 1–2 Gyr
in studies where the SSD is weak or absent. The LSD primarily grows in warm
gas in the galactic midplane. Saturation of the LSD occurs due to
$\alpha$-quenching near the midplane as the growing mean field produces a
magnetic $\alpha$ that opposes the kinetic $\alpha$. The magnetic energy in
our models of the LSD shows slightly sublinear response to increasing
resolution, indicating that we are converging towards the physical solution at
1 pc resolution. Clustering supernovae in OB associations increases the growth
rates for both the SSD and the LSD, compared to a horizontally uniform
supernova distribution.

dynamo — magnetohydrodynamics (MHD) — ISM: supernova remnants — ISM: magnetic
fields — turbulence

††journal: ApJ††software: Pencil Code (Brandenburg & Dobler, 2002; Brandenburg
et al., 2020)

##  1 Introduction

Magnetic fields have a major impact on our observation and interpretation of
astrophysical entities, their structure and dynamics. The intensity and
polarization of emission that we measure to observe the universe are affected
by the location, orientation, and strength of the magnetic field both in the
Solar neighbourhood (Tsouros et al., 2023; Hutschenreuter et al., 2023) and
the more distant universe where these emissions originate or propagate (Beck
et al., 2003). Dickey et al. (2022) compare observations of Faraday rotation
measures of diffuse Galactic synchrotron emission against extragalactic
sources. They show that measurements are highly sensitive to the local
structure of the magnetic field by comparing correlations between the surveys
across the sky (see also, Girichidis, 2021). Measurement of this magnetic
field relies on inferences, based on existing models. Simulations of the
magnetic field can be used to test such models and inferences (e.g., Betti et
al., 2019), particulary if dynamo models can generate realistic turbulent and
large-scale structures.

Magnetic fields influence the dynamics, morphology, and evolution of galaxies.
Presence of a turbulent magnetic field improves by 30–50% (Kirchschlager et
al., 2023) the dust survival rates against supernova (SN) blast waves, a
critical factor in the recycling of the interstellar medium (ISM) and
subsequent star formation processes. The large-scale magnetic field also
affects the gas scale height in galactic discs (Gressel, 2008; Hill et al.,
2012; Gent, 2012) and the multiphase structure of the ISM (Hill et al., 2012;
Evirgen et al., 2017, 2019). Observations of external spiral galaxies show
magnetic fields aligned along the disc and parallel to the spiral arms. The
total field is often stronger in the spiral arms (Fletcher et al., 2011;
Adebahr et al., 2013; Beck, 2015a), although occasionally in the interarm
regions (Beck & Hoernes, 1996; Harnett et al., 2004).

However, some inconsistencies between measurements tracing differing
properties of the magnetic field raise questions about its strength and
structure as a function of galactic radius (Beck, 2015b), in different thermal
phases, or specific environments such as star-forming regions (Borlaff et al.,
2021). We seek to explore numerically the extent differences in the action of
dynamos at various scales and in different thermal phases account for the
observations.

Simulations of the ISM in disc galaxies at kiloparsec scales, which include
magnetic fields by imposing a background (de Avillez & Breitschwerdt, 2005;
Dobbs & Price, 2008; Kim & Ostriker, 2015a), or an initially strong (Hill et
al., 2012; Iffrig & Hennebelle, 2017) horizontal field, evolve SN driven
turbulence to compare the multiphase structure of the ISM with and without
magnetic fields. However, to produce magnetic fields with realistic structure
and energetics, the large-scale dynamo (LSD) must be modelled explicitly.

Figure 1:  Structure of high resolution ($\delta x=1$ pc) Model H2 with
properties given in Table 1. Figure shows (left to right) gas number density,
magnetic energy density, and azimuthal magnetic field during the nonlinear
stage of the dynamo. The lower offset planes show horizontal midplane slices.
The online version of this article includes an animation of this Figure for
the full run time of the model.

Some global models of galactic dynamos that span much larger scales (e.g.,
Otmianowska-Mazur et al., 1997; Hanasz et al., 2009; Kotarba et al., 2010;
Pakmor & Springel, 2013; Butsky et al., 2017) do not resolve the turbulent
scales of the local magnetic field so that the small-scale dynamo (SSD) is
absent or only present at kiloparsec scale (e.g., Rieder & Teyssier, 2016;
Martin-Alvarez et al., 2018).

Models in spiral galaxies with SN-driven turbulence on length scales of at
most a few kiloparsecs (e.g., Korpi et al., 1999; Gressel et al., 2008; Gent
et al., 2013b; Bendre, 2018) show the LSD saturates on time scales of 1–2 Gyr
near equipartion with the ISM kinetic energy density, depending on the rates
of galactic rotation and SN explosions. However, the underlying turbulent
structure of the magnetic field produced by the SSD is either absent or under-
resolved.

The parameters and properties of the SSD in SN-driven turbulence are studied
by Balsara et al. (2004); Balsara & Kim (2005); Gent et al. (2021) and Gent et
al. (2023). Gent et al. (2021) demonstrate that SSD models numerically
converge at cell size $\lesssim 1$ pc and that the SSD is present if the
microphysical resistivity $\eta\lesssim 5\cdot 10^{-3}$ kpc km s-1.

In this work, we model the combination of LSD and SSD in a stratified galactic
disc containing SN-driven turbulence, including high enough grid resolution to
resolve the SSD. We compare to solutions with the LSD suppressed to identify
the effect the LSD has on the SSD; we consider the effect of clustering of SN
in OB associations; and we examine the effect of resolution on the solution of
the LSD. In Section 2 we explain the numerical setup, model parameters and
notation conventions. In Section 3 we describe how the saturated SSD is
affected by the LSD, the effect of SN clustering, and the effect of
resolution. In Section 4 we discuss how the results relate to previous
understanding of the galactic magnetic field. We provide supplementary
illustrative material in Appendix A and outline in Appendix B how the versions
of the Pencil Code adapted during the months and years over which the models
were evolved. Data from the models presented in this paper are available in
the American Museum of Natural History Digital Library (Gent et al., 2023).

##  2 Model description

Table 1:  Models included in analysis of LSD & SSD interactions.  Model | ${\delta{{\mbox{\boldmath$x$}}}}$ | $\Omega$ | Dynamo | $\nu_{6},\chi_{6},\eta_{6}$
---|---|---|---|---
| [pc] | [$\Omega_{\rm Sn}$] |  | [kpc5 km s-1]
H2 | 1.0 | 2 | LSD+SSD | 1 (-15)
L2 | 4.0 | 2 | LSD+SSD | 5 (-12)
L2-cl | 4.0 | 2 | LSD+SSD | 5 (-12)
L0-B0 | 4.0 | 0 | SSD | 5 (-12)
L2-B0 | 4.0 | 2 | SSD | 5 (-12)

###  2.1 Model parameters

We use the shearing sheet approximation to model a cylindrical wedge segment
of a galactic disc using Cartesian coordinates $(x,y,z)$ aligned with
galactocentric radius, azimuth, and the normal to the plane of the disc. We
use dimensions of $L=$(0.512, 1.024, 3.072) kpc, centered vertically on the
midplane. Grid cells in our models are cubical with size $\delta x=4~{}{\rm
pc}$, aside from one high-resolution model with $\delta x=1~{}{\rm pc}$. The
SSD has been shown to numerically converge at these resolutions (Gent et al.,
2021). Model L2 has been initialized with a random field of a few nanogauss,
while all other models start with a random seed field of a few picogauss. This
difference evidently has no effect on the growth rates, once an eigenmode of
the SSD has been established. The density and magnetic field structure of the
high-resolution model are illustrated in Figure 1.

The set of models included are listed in Table 1. Models L2 and H2 apply our
standard set of parameters, differing only by having low or high resolution.
Models otherwise differ from Model L2 as follows. SNe occur in OB associations
in Model L2-cl. To isolate the properties of the SSD from the LSD in this
stratified system, we numerically exclude the mean field by subtracting at
every time step the averages of each magnetic field component in each
horizontal plane for Model L0-B0 without rotation and Model L2-B0 with the
same differential rotation as Model L2.

The radial domain size is sufficient to contain the typical scale of
superbubbles that evolve (Tenorio-Tagle & Bodenheimer, 1988; Mac Low & McCray,
1988; Norman & Ikeuchi, 1989; Ferrière et al., 1991) from multiple SNe, while
the azimuthal domain size of 1 kpc supports multiple independent eddies of the
largest turbulent structures. The vertical scale is sufficient to support two
scale heights of the warm gas, above and below the plane. Open vertical
boundaries allow flows in and out of the domain, but our grid has insufficient
vertical extent to fully model a galactic fountain. Hence, to sustain the
simulation against gas losses due to vertical advection out of the domain, we
continuously maintain the total gas mass by compensating total gains or losses
each time step with a density adjustment in every cell proportional to the
local gas density. The natural rate of gas outflows yields mass loss of around
10% per gigayear, so the instantaneous adjustments are very small, and
relative gradients in momentum are conserved locally in space and time.

The initial condition for all models has a vertical distribution of gas
density and temperature in hydrodynamic and thermodynamic equilibrium. This is
solved numerically in a one-dimensional vertical column of the ISM to ensure
that the disc is in numerical equilibrium. The total mass in the initial state
is approximately equivalent to the combined mass of the cold and warm neutral
medium and the ionized medium, as summarised in Ferrière (2001) Equation (2).

We choose the resistivity $\eta=8\cdot 10^{-4}$ kpc km s-1 and the viscosity
$\nu=4\cdot 10^{-3}$ kpc km s-1 to have high enough values to exceed the
numerical diffusivity at the grid scale for the resolution of these models
(see Gent et al., 2021, for details of how this was demonstrated). We omit
explicit thermal conductivity, which acts at time and length scales below our
model resolution.

###  2.2 Dynamical equations

We use the Pencil Code (Brandenburg & Dobler, 2002; Brandenburg et al., 2020)
to model SN-driven turbulence as described previously in Gent (2012) and Gent
et al. (2020).

Variations of the system of non-ideal, compressible, non-adiabatic,
magnetohydrodynamics, Equations (1)–(4), have been used to simulate SN-driven
turbulence in Gent et al. (2013b, a); Käpylä et al. (2018); Gent et al. (2021)
and Gent et al. (2023). Equations in Gent et al. (2021) and Gent et al. (2023)
have no terms with ${{\mbox{\boldmath$\Phi$}}},\,{{\mbox{\boldmath$\Omega$}}}$
or $S$. The other three studies omit hyperdiffusion and have explicit
diffusivities proportional to sound speed, including thermal diffusivity
$\chi$. Hydrodynamic models of Gent et al. (2013b) and Käpylä et al. (2018)
exclude Equation (3), Ohmic heating and Lorentz force. Gent et al. (2013a)
include shock resistivity for Equation (3).

In this study we solve

$\displaystyle\frac{\text{D}\rho}{\text{D}t}$ $\displaystyle=$
$\displaystyle-\rho\nabla\cdot{{\mbox{\boldmath$u$}}}+\nabla\cdot\zeta_{D}\nabla\rho,$
(1) $\displaystyle\rho\frac{\text{D}{{\mbox{\boldmath$u$}}}}{\text{D}t}$
$\displaystyle=$ $\displaystyle{E_{\rm m}\dot{\sigma}}h^{-1}-\rho c_{\rm
s}^{2}\nabla\left({s}/{c_{\rm
p}}+\ln\rho\right)-\rho\nabla{{\mbox{\boldmath$\Phi$}}}$ (2) $\displaystyle-$
$\displaystyle\rho
Su_{x}{{\mbox{\boldmath$\hat{y}$}}}-2\rho{{\mbox{\boldmath$\Omega$}}}\times{{\mbox{\boldmath$u$}}}+\mu_{0}^{-1}\nabla\times{{\mbox{\boldmath$B$}}}\times{{\mbox{\boldmath$B$}}}$
$\displaystyle+$
$\displaystyle\nabla\cdot\left(2\rho\nu{\textbf{{W}}}\right)+\rho\nabla\left(\zeta_{\nu}\nabla\cdot{{\mbox{\boldmath$u$}}}\right)$
$\displaystyle+$
$\displaystyle\nabla\cdot\left(2\rho\nu_{6}{\textbf{{W}}}^{(5)}\right)-{{\mbox{\boldmath$u$}}}{{\mbox{\boldmath$\nabla$}}}\cdot\left(\zeta_{D}{{\mbox{\boldmath$\nabla$}}}\rho\right),$
$\displaystyle\frac{\partial{{\mbox{\boldmath$A$}}}}{\partial t}$
$\displaystyle=$
$\displaystyle{{\mbox{\boldmath$u$}}}\times{{\mbox{\boldmath$B$}}}-SA_{y}{{\mbox{\boldmath$\hat{y}$}}}-Sx\frac{\partial{{\mbox{\boldmath$A$}}}}{\partial
y} $ (3) $\displaystyle+$
$\displaystyle\eta\nabla^{2}{{\mbox{\boldmath$A$}}}+\eta_{6}\nabla^{6}{{\mbox{\boldmath$A$}}},$
$\displaystyle\rho T\frac{\text{D}s}{\text{D}t}$ $\displaystyle=$
$\displaystyle E_{\rm
th}\dot{\sigma}h^{-1}+\rho\Gamma-\rho^{2}\Lambda+\eta\mu_{0}^{-1}|\nabla\times{{\mbox{\boldmath$B$}}}|^{2}$
(4) $\displaystyle+$ $\displaystyle
2\rho\nu\left|{\textbf{{W}}}\right|^{2}+\rho~{}\zeta_{\nu}\left(\nabla\cdot{{\mbox{\boldmath$u$}}}\right)^{2}$
$\displaystyle+$ $\displaystyle\nabla\cdot\left(\zeta_{\chi}\rho T\nabla
s\right)+\rho T\chi_{6}\nabla^{6}s$ $\displaystyle-$ $\displaystyle
c_{v}~{}T\left(\zeta_{D}\nabla^{2}\rho+\nabla\zeta_{D}\cdot\nabla\rho\right).$

The ideal gas equation of state closes the system, assuming an adiabatic index
(ratio of specific heats) $c_{\rm p}/c_{v}=5/3$. Treating the ISM as a
monatomic, fully ionized plasma we apply a mean molecular weight of 0.531.
Most variables take their usual meanings; a list of notations used is included
in Appendix B.

###  2.3 Gravity and shear

We include an external gravity field with stellar and dark halo contributions
following Kuijken & Gilmore (1989), such that the vertical gravitational
acceleration

$g(z)=\frac{\partial\Phi_{z}}{\partial
z}=-\frac{a_{s}z}{\sqrt{z_{s}^{2}+z^{2}}}-\frac{a_{h}z}{z_{h}},$ (5)

with $a_{s}=1.4~{}{\rm km~{}s}^{-1}~{}{\rm Myr}^{-1}$, $z_{s}=0.2~{}{\rm
kpc}$, $a_{h}=0.5~{}{\rm km~{}s}^{-1}~{}{\rm Myr}^{-1}$ and $z_{h}=1~{}{\rm
kpc}$.

We do not follow the self-gravity of the gas. Taking the mean midplane warm
gas density of 0.1 cm-3 and a mean sound speed of $12~{}{\rm km~{}s}^{-1}$ in
the warm gas, $\lambda_{\rm J}\gg 1~{}{\rm kpc}$ and more so for the diffuse
hot gas. As we do not allow the cold gas to cool below $\sim 100$ K at maximal
number density of a few tens cm-3, star-forming structures cannot form. In
reality, star-forming molecular clouds indeed have a filling fraction of only
around 5% (e. g. Ferrière, 2001) and so can be neglected for the purposes of
modelling the dynamo.

In Equation (2) ${{\mbox{\boldmath$\Omega$}}}=(0,0,\Omega)$ with angular
momentum $\Omega\propto\Omega_{\rm Sn}$ and the Solar neighborhood value
$\Omega_{\rm Sn}=30~{}{\rm km~{}s}^{-1}~{}{\rm kpc}^{-1}$, assuming an orbital
speed for the Sun of $255~{}{\rm km~{}s}^{-1}$ at galactocentric radius
$8.5~{}{\rm kpc}$. The advective derivative,

$\frac{\text{D}~{}}{\text{D}t}=\frac{\partial~{}}{\partial
t}+({{\mbox{\boldmath$U$}}}+{{\mbox{\boldmath$u$}}})\cdot\nabla,$ (6)

includes transport by an imposed shear flow
${{\mbox{\boldmath$U$}}}=(0,Sx,0)$, in which the rate of shear $S=-\Omega$
assuming a flat rotation curve, and a perturbation velocity $u$ that
represents the local deviation from the overall rotational velocity $U$.

###  2.4 Cooling and heating

The uniform far ultraviolet heating follows the functional form given by
Wolfire et al. (1995) with temperature dependence approximated as

$\Gamma_{\rm UV}(T)=\frac{\Gamma_{0}}{2}\left(1+\tanh\left[\frac{2\cdot
10^{4}\,{\rm K}-T}{2000\,{\rm K}}\right]\right),$ (7)

with $\Gamma_{0}=0.0147$ erg g-1 s-1. However, we apply an enhanced heating
rate $\Gamma=3.5\,\Gamma_{\rm UV}$ in order to support the thickness of the
disc in the absence of ionization and cosmic ray energy (Hill et al., 2018).
The cooling, applies a piecewise power-law temperature dependence
$\Lambda(T)=\Lambda_{k}T^{\beta_{k}}$ within the range $T_{k}<T<T_{k+1}$ based
on Wolfire et al. (1995) for cold and warm gas phases and Sarazin & White
(1987) for the hot phase, as parameterized previously (see e.g, Sánchez-
Salcedo et al., 2002; Gressel et al., 2008) and see Gent et al. (2013b, their
Table 1), but with $\Lambda=0$ at $T<90$ K. This cooling law includes
thermally unstable branches at high and low temperatures as expected for the
ISM.

###  2.5 Artificial diffusivities

Terms in the dynamical equations that contain
$\zeta_{D}=2{\cal{C}},~{}\zeta_{\nu}=5{\cal{C}}$ and $\zeta_{\chi}=2{\cal{C}}$
resolve shock discontinuities with artificial diffusion of mass, momentum, and
thermal energy, respectively. They depend on the local convergence
${\cal{C}}=-\nabla\cdot{{\mbox{\boldmath$u$}}}|_{-\rm ve}$, effectively
quadratic in Equation (2), as described in detail in Gent et al. (2020).
Equations (2) and (4) include momentum and energy conserving corrections for
the artificial mass diffusion $\zeta_{D}$ applied in Equation (1). Following
Gent et al. (2021), resistive shock diffusion is omitted from Equation (3),
where hyperdiffusion is sufficient to resolve grid-scale instabilities without
suppressing the important physical process of field compression in SN remnant
shells.

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}$. Hyperdiffusivity
coefficients are $\nu_{6}=\chi_{6}=\eta_{6}=$ $5\times 10^{-12}$ and $1\times
10^{-15}$ for resolution at 4 pc and 1 pc, respectively.

Due to continual code development in response to the performance and evolution
of the simulations over a period exceeding two years in real time, some run
parameters and the version of the code being executed vary. The versions and
key adjustments to the Pencil Code applied to each model at the times
indicated are tabulated in Appendix B.

###  2.6 Supernova explosions

For each SN explosion thermal energy $E_{\rm th}=10^{51}~{}{\rm erg}$ is
instantaneously added to a source region with spherical Gaussian distribution
of radial scale $r_{0}$. For models with $\delta x=4~{}{\rm pc}$, $r_{0}\geq
20~{}{\rm pc}$, while with $\delta{{\mbox{\boldmath$x$}}}=1~{}{\rm pc}$,
$r_{0}\geq 12~{}{\rm pc}$ in order to adequately approximate a sphere on a
Cartesian grid. Where the ambient ISM is very diffuse the radius is adjusted
such that the total mass in the ambient ISM is approximately 180 ${\rm
M}_{\odot}$ to avoid excessive temperatures negatively impacting the timestep.

In dense regions at $\delta x=4~{}{\rm pc}$ up to 5% of the energy is instead
added via momentum injection, $E_{\rm m}$ to compensate for excessive cooling
(see Kim & Ostriker, 2015b; Gent et al., 2020). This applies an axially
symmetric diverging flow, $\Delta{{\mbox{\boldmath$u$}}}$, with speed profile
following the same Gaussian distribution as the thermal energy. The pre-
existing ambient ISM is not modified at the explosion site, so the injected
thermal energy, and velocity field when used, is added to the existing ambient
turbulent ISM.

In most cases no ejecta mass is added, as this is negligible relative to the
ambient ISM mass, but in cases where the temperature exceeds about $10^{8}$ K,
some mass is included, with the same Gaussian profile as the energy, to cool
the gas and maintain the timestep (Olovsson et al., 2005; Joung & Mac Low,
2006). Cumulatively, this added ejecta is $\sim 2M_{\odot}$ per SN and
$\lesssim 40M_{\odot}$ for any single remnant.

SN Type II and I are distinguished only by their Poisson rates of
$\dot{\sigma}_{\rm II}=15~{}{\rm kpc}^{-2}~{}{\rm Myr}^{-1}$ and
$\dot{\sigma}_{\rm I}=2\,~{}{\rm kpc}^{-2}~{}{\rm Myr}^{-1}$, respectively.
The model SN rate is therefore $\dot{\sigma}=\dot{\sigma}_{\rm
II}+\dot{\sigma}_{\rm I}=0.5\dot{\sigma}_{\rm Sn}$ (van den Bergh & Tammann,
1991; Mannucci et al., 2005). The rate half that of the Solar neighbourhood is
adopted to reduce mass losses via galactic wind, while maintaining a
multiphase structure relevant to disc galaxies.

Both types of SNe are uniform randomly distributed horizontally, with Gaussian
distributions in $z$ of scale height $h_{\rm II}=90~{}{\rm pc}$ and $h_{\rm
I}=325~{}{\rm pc}$, respectively. The SN rate $\dot{\sigma}$ applies the
galactic area rate of SNe independently of the vertical size of the domain
$h$, which is included in Equations (2) and (4) to express the average volume
rate. The number of SNe does not change significantly with domain height,
provided it is $\gtrsim 6h_{\rm I}$, so $h=3.072$ kpc is not a model
parameter.

Breathing and sloshing oscillations of the disc (Walters & Cox, 2001) grow
castastrophic if SNe are located independent of the movement of the disc.
Therefore, we track the vertical centre of the Type II SN probability
distribution to the most dense layer of the gas $n_{\rm peak}={\rm
max}(\langle n\rangle_{xy})$. The Type II SN scale height also varies in
response to the thickness of the disc. It is assumed that $n_{\rm peak}$ would
tend to be 1 cm-3. The scale height of the probability distribution is
adjusted to $h_{\rm II}/n_{\rm peak}$. The model total mass is conserved, so
$h_{\rm II}/n_{\rm peak}<(>)1$ cm-3 suggests a higher (lower) scale height for
the gas. The stellar population is not modelled, so this assumes a simple
correlation between star formation distribution and gas disc position and
thickness.

In Model L2-cl, 70% of Type II SNe are also clustered within associations of
radius 40 pc. We assume a Solar neighbourhood OB area rate $\lesssim
0.5~{}{\rm kpc}^{-2}~{}{\rm Myr}^{-1}$. With the Model Type II SN area rate
$\simeq 15~{}{\rm kpc}^{-2}~{}{\rm Myr}^{-1}$ the average cluster includes
about 22 SNe.

When each cluster is initiated the time to the next cluster is determined
using a Poisson rate. The random location of a Type II SN remnant identifies
the centre of the new cluster. There is a 70% probability that a remnant will
occur within the cluster. The cluster has a spherical normal probability
distribution of radial scale 40 pc. Otherwise an SN is located using the
parameters standard to all models.

For the small area of our model it is reasonable to assume that there are no
concurrent clusters, and that consecutive clusters follow immediately.

###  2.7 Averaging and normalization conventions

Angular brackets without a subscript indicate the quantity is averaged over
the entire volume, while the subscript $xy$ limits the average to horizontal
planes, and the subscript $\ell$ indicates that the quantity has been Gaussian
smoothed over a kernel of length $\ell$ prior to averaging. An overbar
indicates averaging over a domain and an interval of time, which may vary, as
specified in the text.

The velocity deviation from the shear flow $u$, solved for in Equation (2),
excludes the imposed galactic shear velocity $U$. However, the system includes
a stratified disc and so is susceptible to the anisotropic kinematic
$\alpha$-effect (AKA), dependent on the rate of shear (see, e.g., Brandenburg
& Rekowski, 2001; Käpylä et al., 2018). Hence, the deviation velocity $u$
includes mean vertical winds and horizontal vortical flows. These are excluded
to determine the turbulent velocity ${{\mbox{\boldmath$u$}}}^{\prime}$, which
is most relevant to the SSD. The turbulent velocities are derived from the
horizontally averaged velocities as

$u_{j}^{\prime}=u_{j}-\langle u_{j}\rangle_{xy}.$ (8)

Similarly the kinetic energy

${\overline{e_{K}}}=\frac{1}{2}\overline{\langle\rho u^{2}\rangle}$ (9)

also includes mean flows, so to isolate the contribution from the turbulent
kinetic energy from the horizontally averaged quantities we use

$\langle e_{K}^{\prime}\rangle_{xy}=\langle
e_{K}\rangle_{xy}-\frac{1}{2}\langle\rho\rangle_{xy}\sum_{j=1}^{3}\langle
u_{j}^{2}\rangle_{xy},$ (10)

To present the solutions for the magnetic energy growth with respect to the
kinetic energy density these are normalised by time- and volume-averaged
quantities. In the case of ${\overline{e_{K}}}$ or
${\overline{e_{K}^{\prime}}}$ only the time interval after the magnetic field
has passed its minimum energy is averaged to exclude initial transients in the
kinetic quantities.

###  2.8 Mean electromotive force

For the LSD the large-scale or mean-field evolution version of Equation (3)
can be derived in terms of the mean magnetic field
${\langle{{\mbox{\boldmath$B$}}}\rangle}$ as

$\frac{\partial{\langle{{\mbox{\boldmath$B$}}}\rangle}}{\partial
t}=\nabla\times({\langle{{\mbox{\boldmath$u$}}}\rangle}\times{\langle{{\mbox{\boldmath$B$}}}\rangle}+{{{\mbox{\boldmath$\cal{E}$}}}}-\eta\nabla\times{\langle{{\mbox{\boldmath$B$}}}\rangle})$
(11)

in which $\cal{E}$ contributes the electromotive force (EMF) from averaging
turbulent fluctuations
$\langle{{\mbox{\boldmath$u$}}}^{\prime}\times{{\mbox{\boldmath$B$}}}^{\prime}\rangle$.
Applying a second order correlation approximation to an inhomogeneous,
anisotropic turbulence as found in the ISM of a spiral galaxy, a general
expression for the EMF has the form

$\cal{E}$ $\displaystyle=$
$\displaystyle\tensorsym\alpha\cdot{\langle{{\mbox{\boldmath$B$}}}\rangle}+{{\mbox{\boldmath$\gamma$}}}\times{\langle{{\mbox{\boldmath$B$}}}\rangle}-\tensorsym\beta\cdot(\nabla\times{\langle{{\mbox{\boldmath$B$}}}\rangle})$
(12) $\displaystyle-$
$\displaystyle{{\mbox{\boldmath$\delta$}}}\times(\nabla\times{\langle{{\mbox{\boldmath$B$}}}\rangle})-\tensorsym\kappa\cdot(\nabla{\langle{{\mbox{\boldmath$B$}}}\rangle})^{(s)},$

in which tensors $\tensorsym\alpha$ and $\tensorsym\beta$ are second order and
$\tensorsym\kappa$ third order, and each term represents a physical process
acting on the dynamo. The $\tensorsym\alpha$ tensor applies effects from
small-scale helicity, $\tensorsym\beta$ turbulent diffusivity, $\delta$
turbulent pumping, $\gamma$ shear or rotation current effects, and
$\tensorsym\kappa$ includes residual effects, dependent on the symmetric part
of the magnetic gradient tensor
$\nabla{\langle{{\mbox{\boldmath$B$}}}\rangle})^{(s)}$. A comprehensive study
of each of these effects will be required to explain fully the LSD, but that
is beyond the scope of the current investigation.

However, some indication of the EMF properties may be extracted from examining
a simplified analog of the $\tensorsym\alpha$ term. Under the assumption of
isotropic, homogeneous turbulence $\tensorsym\alpha$ can be reduced to a
scalar $\alpha=\alpha_{K}+\alpha_{M}$. The kinetic helicity contributes to the
LSD as

$\alpha_{K}\approx-\frac{1}{3}\tau\langle{{\mbox{\boldmath$u$}}}^{\prime}\cdot{{\mbox{\boldmath$\omega$}}}^{\prime}\rangle,$
(13)

where $\tau$ is the correlation time of the turbulence and
${{\mbox{\boldmath$\omega^{\prime}$}}}=\nabla\times{{\mbox{\boldmath$u^{\prime}$}}}$
is its vorticity. Due to the conservation of magnetic helicity, for the LSD to
exist some small-scale helicity flux is required. The magnetic helicity
contributes to the LSD as

$\alpha_{M}\approx\frac{1}{3}\tau\langle{{\mbox{\boldmath$B$}}}^{\prime}\cdot{{\mbox{\boldmath$J$}}}^{\prime}\rho^{-1}\rangle.$
(14)

While these simplified expressions are insufficient to describe the full range
of turbulent transport motions driving the LSD, they are relatively easy to
obtain and offer some indication of how the LSD changes over time and differs
between models. The Lorentz force acts against the flow as the dynamo
approaches saturation, where opposite sign of $\alpha_{M}$ can lead to
$\alpha$-quenching (Krause & Raedler, 1980) of the LSD.

##  3 Results

###  3.1 Small-scale dynamo

Figure 2:  Volume-averaged magnetic energy density for models L2-B0 (_blue_)
with $\Omega=2\Omega_{\rm Sn}$ and L0-B0 (_green_) with $\Omega=0$, in which
LSD is suppressed by continuously removing the horizontally averaged magnetic
field components to ensure only SSD is present. Exponential fits for the
growth rates measured during the dotted intervals are listed in the legend.
Normalisation is by the time-averaged turbulent kinetic energy
${\overline{e_{K}^{\prime}}}$. The horizontal dotted line indicates 2%
equipartition.

To isolate the properties of the SSD in the stratified and differentially
rotating ISM we exclude the presence of a mean field in models L0-B0 and L2-B0
by continuously subtracting the averages of each component of the magnetic
field on each horizontal plane. The energy density of the turbulent magnetic
field generated by the SSD is plotted in Figure 2. The saturation of the
magnetic field for Model L0-B0 occurs at about 2% of equipartition with the
time-averaged turbulent kinetic energy density ${\overline{e_{K}^{\prime}}}$
of Equation (10). The SSD in Model L0-B0, without shear, grows with
exponential growth rate $\gamma=24.6~{}{\rm Gyr}^{-1}$, while Model L2-B0
grows slightly faster, with $\gamma=26.7~{}{\rm Gyr}^{-1}$, and saturates at a
higher level. Large-scale shear and rotation, therefore, partially enhance the
SSD.

(a)(b)

Figure 3:  _(a)_ Magnetic energy spectra at specified times in pure SSD Models
L0-B0 and L2-B0, each normalised by its maximum power. The spectra for Model
L0-B0 are shifted $10^{-2}$ to differentiate from Model L2-B0. _(b)_ Kinetic
energy spectra with matching line styles, times and axis shift. Listed
wavenumbers $k$, in multiples of $k_{1}=2\uppi{L_{x}}^{-1}$, locate the
maximal magnetic power at each time.

Samples of the spectral properties of the magnetic energy and kinetic energy
of both pure SSD models are displayed in Figure 3. Spectra have been smoothed
with a Gaussian kernel of length 5 Myr and 1.5 $k_{1}$ in time and wavelength.
The magnetic energy spectra, which grow many orders of magnitude over time,
are normalised by maximum energy density to ease comparison between times. In
both models the peak wavenumber of the magnetic energy during the kinematic
stage of the SSD is at $17k_{1}$, with $k_{1}=2\uppi{L_{x}}^{-1}\simeq
12.3~{}{\rm kpc}^{-1}$. This approaches the resistive scale for these models,
as expected for the SSD.

The statistically steady kinetic spectra are derived only from the velocity
$u$, but are multiplied by the conserved volume-averaged mean gas density to
approximate the kinetic energy density. The spectral properties of the flow
appear quite insensitive to the large-scale rotation and shear, but in the
saturated state the shear appears to excite a peak magnetic energy with
slightly larger scales than without.

(a)(b)

Figure 4:  Time evolution of magnetic energy for SSD models _(a)_ L0-B0 and
_(b)_ L2-B0 at sampled scales of $k$, as noted in the legends in multiples of
$k_{1}=2\uppi{L_{x}}^{-1}$.

(a)(c)(b)(d)

Figure 5:  Time-latitude diagrams of the horizontally averaged magnetic field
_(a)_ $\langle B_{x}\rangle_{xy}$ and _(b)_ $\langle B_{y}\rangle_{xy}$, for
the fiducial model L2 in microgauss and, respectively, _(c)_ and _(d)_
normalised by the time-dependent maxima to highlight the evolving structure
and scales of the LSD.

To examine the evolution of the magnetic field over time, we plot the power in
selected wavenumbers of the magnetic spectra in Figure 4. Both models follow
an eigenmode growing at the same rate at all wavelengths. The ranking of the
power in each wavenumber $k$ is quite rigid throughout for both models. These
results show that the SSD behaviour and structure is not affected by the
large-scale rotation and shear, other than to increase the growth rate and
saturation energy of the SSD.

###  3.2 Large-scale dynamo

Having examined the characteristics of the SSD in a rotating shear flow with
Model L2-B0, we now apply the same parameters, but do not suppress the growth
of the large scale field in Model L2. Figure 5 displays the evolution of the
mean magnetic field as horizontal averages of the magnetic field components
$\langle B_{x}\rangle_{xy}$ and $\langle B_{y}\rangle_{xy}$ up to saturation
of the LSD. Relative to the horizontal components, $B_{z}$ has stronger local
fluctuations driven by the vertical flows. However, due to periodic boundary
conditions on $A_{x}$ and $A_{y}$, its horizontal averages must vanish.

The azimuthal shear favours a steep pitch angle of the mean magnetic field,
which ultimately primarily aligns azimuthally. The late time value of $\langle
B_{y}\rangle_{xy}$ is strongly positive and an order of magnitude stronger
than mean $\langle B_{x}\rangle_{xy}$. In the period beyond 3.9 Gyr, after
which the mean field remains quite steady, along the midplane
$\overline{\langle B_{x}\rangle_{xy}}=-0.17\pm 0.20\,\upmu$G, while
$\overline{\langle B_{y}\rangle_{xy}}=9.16\pm 3.56\,\upmu$G, with errors given
by the standard deviation.

In Figures 5 _(c)_ and _(d)_ , horizontal averages of $B_{x}$ and $B_{y}$ are
shown normalised by their time-dependent maximal values, so that their
structure at early times can be easily compared to the saturated state. It is
evident that the mean of both components during the kinematic stage of the SSD
fluctuates on short time and spatial scales. More coherent and persistent
structures emerge after 1 Gyr, when the SSD has already saturated. At first
these are typically of opposite sign to the eventual state, with more positive
$B_{x}$ and negative $B_{y}$. Although the sign reversals persist in $B_{x}$,
these fluctuations remain much weaker than in $B_{y}$. The vertical topology
of the field remains steady after about 3.75 Gyr, although the strength of the
field continues to grow as seen in panel _(b)_.

Figure 6:  Fiducial LSD Model L2: energy density of total magnetic field
$\langle B^{2}\rangle(8\uppi)^{-1}$ (_blue_), horizonally averaged mean-field
${\langle B\rangle_{xy}}^{2}(8\uppi)^{-1}$ (_green_) and Gaussian smoothed
mean-field ${\langle B\rangle_{\ell}}^{2}(8\uppi)^{-1}$ (_red_), with
$\ell=50~{}{\rm pc}$. (Factors of $(8\uppi)^{-1}$ are dropped in the legend
for readability in this and following figures showing magnetic energy
evolution.) Normalisation is by the time-averaged turbulent kinetic energy
density ${\overline{e_{K}^{\prime}}}$ (_blue dotted_), while the total kinetic
energy is given by ${\overline{e_{K}}}$ (_red dotted_). Exponential fits for
the growth rates are listed in the legend, applying to intervals $365~{}{\rm
Myr}<t<800~{}{\rm Myr}$, $0.9~{}{\rm Gyr}<t<1.2~{}{\rm Gyr}$ and $2.5~{}{\rm
Gyr}<t<4.8~{}{\rm Gyr}$, respectively.

All profiles of Figure 5 exhibit oscillatory structure. These correspond to
the breathing and sloshing modes predicted by Walters & Cox (2001), due to
gravitational interaction with gas outflows. (This is also shown in gas number
density in Appendix A.1 for a zoomed-in interval.) The period of 100–150 Myr
supports the model of Walters & Cox (2001) well for Solar neighbourhood
parameters.

In Figure 6 we plot the growth of the total magnetic field (blue) for Model
L2. The mean field is also plotted as horizontal averages $\langle
B\rangle_{xy}$ and Gaussian smoothed $\langle B\rangle_{\ell}$, with a
smoothing length $\ell=50~{}{\rm pc}$, as described previously (Gent et al.,
2013b). Fits for the exponential growth rate $\gamma$ are listed in the legend
for intervals spanning the SSD, a transitional stage, and the LSD. The growth
rate of the total field during the SSD matches that of Model L0-B0 at
$24.6~{}{\rm Gyr}^{-1}$ and lies about 10% below that of Model L2-B0, with
otherwise equal parameters to this model. The kinematic SSD therefore appears
to be mildly inhibited by the action of the LSD.

During the SSD, the mean field grows almost identically with the SSD, and
counter to the finding of Gent et al. (2013b) is near insensitive to the
definition of the mean field. This suggests the growth of the mean field
during this phase is the result of winding up by the large scale shear of the
growing random field. It is reflected in the small scale structure of the
horizontal averages during this period as shown in Figure 5.

Before the total field reaches 2%${\overline{e_{K}^{\prime}}}$, at which the
SSD saturates for Model L0-B0 and L2-B0 (Figure 2), the mean-field growth
slows down. During this transition the definition of the mean field does
substantially affect the measured growth rate, which agrees with the analysis
of Gent et al. (2013b), who omitted the kinematic phase and considered only
this transition. Later, as the LSD becomes dominant, the definition makes
little difference to evolution of the mean-field energy, although it affects
very much its topology (e.g., see Figure 7.6 Shukurov & Subramanian, 2021).

(a)(b)

Figure 7:  _(a)_ Magnetic energy spectra from snapshots at different times
from fiducial Model L2. The spectra are normalised by maximum power in each
snapshot. _(b)_ Kinetic energy spectra with line styles and colour matching
_(a)_ for the same snapshots. Spectra have been smoothed with a Gaussian
kernel of length 5 Myr and 1.5 $k_{1}$ in time and wavenumber. Snapshot time
and wavenumber $k$ in multiples of $k_{1}=2\uppi{L_{x}}^{-1}$ of the maximal
magnetic power in the legends are common to of both panels.

The SSD does not saturate, but continues to grow at a smaller rate
$\gamma=4.2~{}{\rm Gyr}^{-1}$ above 2%${\overline{e_{K}^{\prime}}}$. The
strength of the growth of the fluctuating field is unexpected. The mean field
remains at least two orders of magnitude weaker than the turbulent field.
Tangling must therefore be capable of generating turbulent field at rates
substantially larger than the source. The mean field grows at a rate $\gamma$
between 7 and $14~{}{\rm Gyr}^{-1}$, but this slows down as its energy
approaches 2%${\overline{e_{K}^{\prime}}}$.

Figure 8:  Time evolution of magnetic energy for LSD model L2 over sampled
scales of $k$, as noted in the legends in multiples of
$k_{1}=2\uppi{L_{x}}^{-1}$.

During the kinematic stage of the SSD, the spectral structure of model L2, as
displayed for early snapshots in Figure 7, is equivalent to those of Models
L0-B0 and L2-B0 with peak energy at $k=17k_{1}$. After saturation of the SSD
and transition to the LSD, the magnetic energy shifts to larger scales
($k\rightarrow k_{1}$). The dominance of the largest scale increases over time
until the LSD saturates around 5 Gyr.

The kinetic energy varies strongly over time due to the sporadic and
heterogeneous nature of the SN forcing, but the dominant energy persistently
remains at large scales $k<5k_{1}$. As the LSD approaches saturation there is
some reduction in kinetic energy at all scales, but large scale flows are
particularly suppressed, with peak power tending to shift to $k\gtrsim
5k_{1}$.

(a)(b)(c)

Figure 9:  Magnetic energy spectra during epochs of the SSD as indicated in
each panel (L2 is shifted 365 Myr, as discussed in Sect. 3.3). Models listed
in the legends have power law fits for comparison with Kazantsev $k^{3/2}$
scaling. Spectra have been smoothed with a Gaussian filter of length 5 Myr in
time and 1.5 $k_{1}$ in wavelength.

In Figure 8 the evolution of the energy in wavenumber bins with $k\geq
17k_{1}$ is insensitive to the LSD and retains the same order as for the SSD
in Figure 4. In contrast, wavenumber bins with $k<17k_{1}$, which had lower
energies during the SSD, gain energy, quickly becoming dominant over the
larger wavenumber bins with $k>17k_{1}$ following saturation of the SSD. Thus,
the LSD can be identified as acting primarily on scales $k\lesssim 17k_{1}$.
However, in the nonlinear stage of the LSD, even the smallest scales of the
magnetic field increase in energy, tenfold compared to the pure SSD models,
with $k=60k_{1}$ exceeding $10^{-16}$ erg cm-3.

The energy growth rate in wavenumber bin $k\leq 1.3k_{1}$ (Figure 8) separates
first from the SSD eigenmode, followed by energy binned at successively
increasing wavenumber $k<17k_{1}$. This mirrors contrasting growth in Figure 6
between $\langle e_{B}\rangle_{xy}$, consisting of only the longest length
scale, and $\langle e_{B}\rangle_{\ell}$, comprising a range of intermediate
scales. The growth of $\langle e_{B}\rangle_{xy}$ slows down earliest, but it
subsequently grows more rapidly than $\langle e_{B}\rangle_{\ell}$.

In Figure 9 we plot the magnetic energy spectra together for Models L2, L0-B0,
and L2-B0, at three characteristic stages of the SSD. To compare equivalent
stages of the evolution of the SSD between the three models the times selected
for Model L2 are shifted by 365 Myr. In the early kinematic stage the Models
L0-B0 and L2-B0 have a power law approximating $k^{3/2}$ (Kazantsev, 1968) at
small wavenumbers, consistent with the expected Kazantsev theory.

Even at these early times the power law is somewhat more shallow for Model
L2-B0, due to the differential rotation, and yet more so for L2, due to the
inclusion of LSD, increasing energy at larger scales. Throughout the kinematic
stage, however, all have the peak energy at $k\simeq 17k_{1}$. All models show
a shift in energy toward lower $k$ as the SSD saturates, with the power law
$\rightarrow k^{1}$ scaling at low $k$, but only for Model L2 with LSD does
the power law $\rightarrow k^{0}$ at late times.

(a)(b)(c)

Figure 10:  Energy spectra for _(a)_ Model L0-B0, _(b)_ L2-B0 and _(c)_ L2.
Power law fits are shown for times indicated (not time shifted). Each spectrum
is normalised to its own maximum, but shifted by a small factor to
differentiate each time, while adequately comparing spectral changes that
differ by a factor of $10^{10}$ over time.

Figure 10 displays evolving spectral structure for each model from the
kinematic stage to saturation of the SSD. Even in the presence of the LSD,
evolution is predominantly self-similar until the SSD saturates, with scaling
$\simeq k^{5/4}$ for Model L2 between 0.08 and 0.78 Gyr in Figure 10 _(c)_ ,
and $\simeq k^{3/2}$ for the equivalent stages in _(a)_ and _(b)_. The
steepening resistive tail in Figure 10 _(c)_ after the SSD saturates indicates
that LSD continues to pile energy into larger scales within the intertial
range, which is not so for L0-B0  _(a)_ or L2-B0  _(b)_.

###  3.3 Numerical convergence

(a)(b)(c)

Figure 11:  _(a)_ Volume-averaged magnetic energy density for all models
listed in the legend, normalised by the time-averaged turbulent kinetic energy
density ${\overline{e_{K}^{\prime}}}$. _(b)_ Volume-averaged turbulent kinetic
energy density for the LSD models, and _(c)_ volume-averaged root-mean-square
turbulent velocity $u^{\prime}$.

We next study the numerical convergence of the combined dynamos. We first
consider the convergence behavior of the SSD. Figure 11 _(a)_ shows the
extreme sensitivity of the growth rate of the SSD to numerical resolution
described by Gent et al. (2021) and Gent et al. (2023), with the high-
resolution Model H2 saturating within 10–20 Myr, in contrast to around 0.8–1.2
Gyr for L2/L2-B0. The SN rate in models L2 and H2 are 2.5 times higher than
modelled in Gent et al. (2023) and the density is stratified, facilitating
even higher velocity and sound speed. These models suggest that the saturation
time for the SSD in the SN-driven ISM likely has a converged limit of only a
few megayears.

In order to study the numerical convergence of the combination of the LSD and
SSD after saturation of the SSD, the large variation in time of SSD saturation
needs to be taken into account, even though the actual saturation energy of
the magnetic field has converged well at these resolutions. For comparisons of
Model L2 to Model H2 during the period that the LSD dominates field growth,
Model H2 can be shifted by 970 Myr to align the LSD solutions after the SSD
has saturated. (Similarly, shifting the dynamo solution for Model L2
illustrated in Figure 11(a) by 365 Myr aligns it with the SSD solutions for
Models L0-B0 and L2-B0 for the comparisons in Section 3.2 above.)

Figure 12:  Model H2: magnetic energy density normalised by time-averaged
turbulent kinetic energy ${\overline{e_{K}^{\prime}}}$. Total magnetic energy
$\langle B^{2}\rangle(8\uppi)^{-1}$ (_blue_) and horizontally-averaged mean-
field energy ${\langle B\rangle_{xy}}^{2}(8\uppi)^{-1}$ (_green_) have fits
for exponential growth rate $\gamma$ as listed in the legend, spanning
$12~{}{\rm Myr}<t<19~{}{\rm Myr}$, $25~{}{\rm Myr}<t<33~{}{\rm Myr}$ and
$59~{}{\rm Myr}<t<0.2~{}{\rm Gyr}$, respectively.

The turbulent velocity of model H2 (Fig. 11 _(c)_) does slightly exceed that
of L2, but the main driver of the high rate of SSD at high resolution is the
resolution of sufficiently small scale turbulent structures. The turbulent
kinetic energy density (Fig. 11 _(b)_) is well converged at the resolution of
these models, as is the saturation strength of the SSD.

In Figure 12 we show the evolution of the total and horizontally averaged mean
magnetic energy of Model H2. For comparison with Model L2 from here and Figure
6 we align intervals relevant to the LSD in Figure 13. In the kinematic stage
of the SSD the growth rate $\gamma=1300~{}{\rm Gyr}^{-1}$ for H2 vs
$24.6~{}{\rm Gyr}^{-1}$ for L2, quite consistent with our previous results for
SSD in the multiphase ISM. Thereafter, the LSD appears to grow somewhat more
rapidly at high resolution than low.

(a)(b)

Figure 13:  Volume averaged energy density _(a)_ of total magnetic energy
$\langle B^{2}\rangle(8\uppi)^{-1}$ and horizontally averaged mean-field
energy ${\langle B\rangle_{xy}}^{2}(8\uppi)^{-1}$, normalised by time-averaged
turbulent kinetic energy density ${\overline{e_{K}^{\prime}}}$ for Models L2
and H2 (with the latter shifted 0.97 Gyr to align saturation of their
respective SSDs). Exponential fits of growth rates, as listed in the legend,
apply to the interval $1.08~{}{\rm Gyr}<t<1.5~{}{\rm Gyr}$. _(b)_ Ratio of
horizontally-averaged mean-field energy to turbulent magnetic field energy for
same models and intervals as _(a)_.

Figure 13 _(a)_ compares the total and horizontally-averaged magnetic energy
of Model L2 to the time-shifted solution of Model H2. The growth rate at four
times higher linear resolution is less than a factor of 1.5 higher. The
amplitude of the magnetic energy, 500 Myr following saturation of the SSD,
differs by a factor of 3.8.

Figure 13 _(b)_ shows the ratios obtained from the values in panel _(a)_ ,
which are $\lesssim 0.1$ and depend more on the evolution of the mean-field
energy than the resolution of the SSD. Estimates from observation place the
ratio of mean to turbulent magnetic energy at 0.1–1 (e.g., Tabatabaei et al.,
2008; Houde et al., 2013; Haverkorn, 2015; Beck, 2015a). In the later stages
the proportion of turbulent energy is even higher for Model L2 than H2.
However, this is explained by the increased efficiency of the LSD in Model H2
which remains well below saturation, so the system is not yet in a statistical
steady state.

Finally, Figure 14 compares magnetic energy spectra from the two models
evolving over time. During the SSD the energy peak is quite consistent around
$k\simeq 17k_{1}$.

This Figure shows the continued growth of the large scale (low $k$) field in
Model H2 and its saturation in model L2. In the high-resolution model, the
largest available scale has the most energy, unlike in the lower resolution
model at late times where the spectrum flattens. The amplitude of the lowest
$k$ power at the final time in Model H2 is 3.9 times that of Model L2.

The examined quantities all show substantial but, at least marginally,
sublinear dependence on numerical resolution. This suggests that our models
are converging towards a physical solution, though higher resolution remains
desirable.

(a)(b)

Figure 14:  Energy spectra from _(a)_ low-resolution Model L2 and _(b)_ high-
resolution Model H2. Power law fits are shown at times indicated (970 Myr
earlier for Model H2 to account for its fast SSD).

In Appendix A.2 comparisons of spectra from equivalent times and time-latitude
diagrams of the magnetic fields during this corresponding interval illustrate
the convergence in the structure and topology of the magnetic field. At this
early stage of the LSD the azimuthal field is predominantly negative and has
similar scale height independent of resolution, although it later evolves, at
least in the low-resolution models, to be positive and symmetric.
Computationally, the LSD solution for H2 requires a long enough integration
time that it has not yet been determined whether the similarity in the LSD
solution to the low-resolution model extends to the later stages of the
dynamo.

###  3.4 Effect of clustering

(a)(b)

Figure 15:  Volume-averaged total kinetic energy density
${\overline{e_{K}^{\prime}}}$ _(a)_ for all models as listed in the legend and
_(b)_ volume-averaged root-mean-square total velocity $u$ for the same models.

Figure 11 _(a)_ shows that Model L2-cl, with OB clustering included, has
growth of magnetic energy driven by the SSD over four times faster than Model
L2, without clustering. (Note that, although they start with somewhat
different seed fields, the growth rates are unaffected by the difference.)
Clustering increases the volume filling fraction of hot superbubbles in the
midplane. The increased effective sound speed supports rapid growth of the
SSD, twice that of Model L2, and 4/3 faster after the SSD has saturated.

In Figure 15(a) we see that OB clustering supports higher total kinetic
energy, but that this primarily results from stronger large-scale motions as
the turbulent contribution to the kinetic energy, plotted separately in Figure
11(b), is actually weaker than for Model L2. By contrast the turbulent
velocity in Figure 11 _(c)_ has a higher magnitude (Fig. 15 _(b)_) when OB
clustering is included, because the high turbulent velocities are in regions
of low density and thus contribute little to the kinetic energy. Gent et al.
(2023) have demonstrated that it is the flow velocity, rather than the kinetic
energy, that is critical to the SSD growth rate in the multiphase ISM, while
saturation strength is related to the kinetic energy and is thus independent
of the inclusion of clustering.

For Model L2-cl dips in total kinetic energy can be traced to the blowout of
gas occurring in the northern latitudes at these times. We include in Appendix
A.1 time-latitude diagrams of the gas density and azimuthal magnetic field to
contrast the topology of the gas and resultant field between Models L2 and
L2-cl.

By shifting the solution displayed in Figure 11 for L2-cl by 420 Myr the stage
immediately following the saturation of the SSD aligns with Model L2. Taking
this comparative interval of the LSD for Model L2-cl as that fitted for Model
L2 in Figure 13 yields $\gamma=4~{}{\rm Gyr}^{-1}$ vs 3 for the total magnetic
energy and 11.1 vs 7.8 for the mean magnetic energy, respectively (See also
Appendix A.1).

###  3.5 The $\alpha$-effect

We approximate the proxy for $\alpha=\alpha_{K}+\alpha_{M}$ in Eqs. (13) and
(14) from the time series of horizontal averages with

$\langle{{\mbox{\boldmath$u$}}}^{\prime}\cdot{{\mbox{\boldmath$\omega$}}}^{\prime}\rangle_{xy}=\langle{{\mbox{\boldmath$u$}}}\cdot{{\mbox{\boldmath$\omega$}}}\rangle_{xy}-\sum_{k=1}^{3}\langle
u_{k}\rangle_{xy}\langle\omega_{k}\rangle_{xy}$ (15)

and

$\left<\frac{{{\mbox{\boldmath$B$}}}^{\prime}\cdot{{\mbox{\boldmath$J$}}}^{\prime}}{\rho}\right>_{xy}=\frac{\langle{{\mbox{\boldmath$B$}}}\cdot{{\mbox{\boldmath$J$}}}\rangle_{xy}}{\langle\rho\rangle_{xy}}-\sum_{k=1}^{3}\frac{\langle
B_{k}\rangle_{xy}\langle J_{k}\rangle_{xy}}{\langle\rho\rangle_{xy}},$ (16)

where the first expression on the right hand side is used, because
$\langle({{\mbox{\boldmath$B$}}}\cdot{{\mbox{\boldmath$J$}}})\rho^{-1}\rangle_{xy}$
is unavailable. In Figure 16 this $\alpha$ is illustrated as time-latitude
diagrams for Models L2 and H2, during their entire evolution, and also during
a period of Model L2 comparable to Model H2 to assist comparison.

We assume a correlation time of $\tau=10~{}{\rm Myr}$, common for both models
and both $\alpha_{K}$ and $\alpha_{M}$. This is consistent with an application
of the test-field method (Bendre, 2018) and yields rms values for $\alpha$ in
Models L2, H2 and L2-cl of 15, 30, and $89~{}{\rm km~{}s}^{-1}$, which
correspond reasonably to calculations of $u^{\prime}$ included in Figure
11(c). It is slightly higher than the $5~{}{\rm Myr}$ determined by Hollins et
al. (2017), the $1.9~{}{\rm Myr}$ of Brandenburg et al. (2013), or the order
of magnitude estimate 3–7 Myr of (Chamandy & Shukurov, 2020). Whether $\tau$
differs between magnetic and kinetic $\alpha$, or overall, a factor of 1/5
would not qualitatively alter the dynamics.

L2L2H2(a)(b)(c)

Figure 16:  Time-latitude diagrams of horizontally averaged $\alpha$-effect
for Model L2  _(a)_ throughout and _(b)_ during a restricted interval
comparable with _(c)_ Model H2Ȧ linear colour-scale is applied in the range
$-5~{}{\rm km~{}s}^{-1}<\alpha<5~{}{\rm km~{}s}^{-1}$ and symmetric
logarithmic colour-scale for larger magnitudes.

Considering the turbulent $\alpha$ displayed for Model L2 in Figure 16(a) the
northern latitudes are predominantly positive, while the southern latitudes
are negative, which is consistent with the action of the Coriolis force on
each side of the midplane. We provide a breakdown of each $\alpha$
contribution in Appendix A.3. This shows that the evolution of $\alpha$ is
dominated by the contribution of kinetic $\alpha_{K}$ in the first 2 Gyr. The
negative sign in the northern latitudes near 2 Gyr is due to current helicity,
as are the strong high latitude signatures after 2.5 Gyr. However, the high
altitude $\alpha$-effect where the magnetic field is quite weak likely has
little effect on the LSD.

The reversal in sign of the mean field after 2 Gyr corresponds to damping of
$\alpha$ in the dynamo active warm gas near the midplane, but also with growth
of magnetic $\alpha_{M}$ of opposite sign to $\alpha_{K}$ in this region (see
Appendix A.3). The $\alpha$-quenching near the midplane, due to reducing
$\alpha_{K}$ and growing magnetic $\alpha$ of opposite sign saturates the LSD.

Comparing Models L2 and H2 in Figures 16(b) and (c) during corresponding 700
Myr intervals, their $\alpha$ topologies are qualitatively similar. At high
resolution, $\alpha$ is typically stronger, and also has reversals in sign,
which are more ubiquitous and on smaller spatial scales and time scales. In
both models helicity of opposite sign crosses the midplane and propagates away
from the disc. The small-scale structure of the helicity occupies the
outflows, while large-scale same-sign helicity populates the inflows. Thus,
the small-scale helicity is sporadically removed from the disc allowing the
LSD to grow.

##  4 Discussion and conclusions

We have verified that the large-scale dynamo (LSD) in disc galaxies is robust
in the presence of SSD at high magnetic Reynolds (Rm) and saturates in
approximate equipartition with the kinetic energy density. However, the LSD
growth rate appears to be impeded by the presence of a vigorous small-scale
dynamo (SSD), in comparison with the LSD in cases where the SSD is weaker or
absent altogether (Korpi et al., 1999; Gressel et al., 2008; Gent et al.,
2013b; Bendre, 2018). In contrast, the SSD appears quite insensitive to the
presence of an LSD. The presence of large-scale rotation and shear does,
however, inject more vorticity, increasing the growth rate of the SSD (Gent et
al., 2023). The increased growth rates apply primarily to larger scales.

In the kinematic stage of the SSD, all scales in the magnetic field grow self
-similarly. Upon saturation and in the absence of SSD, the hierarchy of the
magnetic field spectrum is preserved. The power peaks at a wavenumber $k\simeq
17k_{1}$ independent of resolution (or effective Rm), which remains the case
in the nonlinear phase of the SSD. When the LSD is present, the self similar
growth of the magnetic energy spectra continues in the same manner for all
$k\gtrsim 17k_{1}$, through to saturation of the SSD and beyond. This
identifies the turbulent field of the SN driven ISM as occupying structures
organised on length scales $\ell\lesssim 30~{}{\rm pc}$.

Unlike the case lacking LSD, the small scales of the magnetic field do
continue to grow after the SSD has saturated, but much more slowly than the
large scales, which might be explained by tangling of the growing large-scale
field. For $k\lesssim 17k_{1}$, before the SSD saturates the self-similar
growth ceases and growth at the largest scales slows down. After the LSD
becomes dominant, the larger scales overtake the energy at $k\simeq 17k_{1}$
and the growth at these smaller scales slows even further, until eventually
the largest scale contains the peak energy and grows faster than any other
scale. This indicates that a horizontal domain of 1 kpc is insufficient to
identify the peak energy scale of the mean magnetic field in a disc galaxy
with an active LSD.

###  4.1 Mean-field decomposition

When considering mean-field dynamo theory it can be helpful to decompose the
field into the mean and fluctuating fields subject to a large separation of
scales. In practice, however, it is instead common in simulations,
experiments, and observations to find a continuum of power across scales
without any clear minimum in the spectrum between large and small scales.
Horizontal averages have an arbitrary dependence on the choice of domain,
while kernel averaging can rely on an arbitrary smoothing scale. Following
Gent et al. (2013b); Hollins et al. (2017) we adopted a Gaussian smoothing
scale of $\ell=50~{}{\rm pc}$ to separate the mean and random magnetic fields
for $\langle B\rangle_{\ell}$ in Figure 6.

Our analysis of the temporal evolution of the spectrum provides a physically
motivated criterion for separating the mean field from the turbulent field.
This depends on their sensitivity, or lack thereof, to the action of the LSD.
Truncating the spectrum of the magnetic field at $k\geq 17k_{1}$ allocates to
the mean-field only those scales sensitive to the LSD, while masking the
spectrum at $k<17k_{1}$ allocates to the random field those scales insensitive
to the LSD.

###  4.2 OB clustering

In a horizontal domain of 1 kpc we find that simulations including clustering
of OB stars in associations can rapidly develop thermal runaway (Kim &
Ostriker, 2015b) on either side of the midplane. The expanded fractional
volume of hot gas is beneficial to the SSD (Gent et al., 2023), but does not
favour the LSD, as can be seen by comparing Figure 17(c) to Figure 5(d). With
OB clustering, the field is confined to the warm gas at the midplane, while
without OB clustering the field is substantial at large heights. Figure 11(a)
shows that the LSD actually grows more efficiently in the thinner, more dense
disc of the OB-clustered model than without. Both observations strongly
support the hypothesis that the mean-field dynamo grows mainly in the warm
rather than the hot gas.

###  4.3 Quenching of the $\alpha$-effect

We have used kinetic and current helicities as proxies for measuring the
collective inductive action of SN turbulence, called the $\alpha$-effect. A
more thorough analysis, yielding also other important turbulent transport
coefficients, would require the usage of a suitable test field method (see
Käpylä et al., 2022), but due to significantly increased computational costs,
such analysis is not performed here.

As Model L2 approaches equipartition in Figure 16(a) we observe
$\alpha$-quenching near the midplane, with magnetic helicity of opposite sign
suppressing the dynamo. While the LSD is active, small-scale opposite sign
helicity is removed from the disc by the outflows. At large $|z|$, the
$\alpha$-effect remains non-zero, with magnetic and kinetic helicity typically
of the same sign. However, both magnetic and kinetic energy density are
extremely weak relative to the disc.

It is possible that the reversal of the midplane field at around 2 Gyr is an
indication of the nonlinear magnetic buoyancy instability. The coincident
change of sign of $\alpha$ in the northern latitudes, which is a magnetic
$\alpha$-effect, is a signature of this behaviour as identified by Machida et
al. (2013); Tharakkal et al. (2023), and Qasi et al. (2023). Subsequent
reversals do not occur, as the scale height of the magnetic field expands,
reducing the vulnerability of the system to magnetic buoyancy instability.

F.A.G. and M.J.K.-L. acknowledge support from the Academy of Finland ReSoLVE
Centre of Excellence (grant 307411), the Ministry of Education and Culture
Global Programme USA Pilot 9758121 and the ERC under the EU’s Horizon 2020
research and innovation programme (Project UniSDyn, grant 818665) and generous
computational resources from CSC – IT Center for Science, Finland, under Grand
Challenge GDYNS Project 2001062. F.A.G. benefited from in-depth discussions at
“Towards a Comprehensive Model of the Galactic Magentic Field” at Nordita 2023
supported by NordForsk. M-M.M.L. was partly supported by US NSF grant
AST18-15461, and thanks the Inst. für Theoretische Astrophysik der Uni.
Heidelberg for hospitality.

## References

  * Adebahr et al. (2013) Adebahr, B., Krause, M., Klein, U., et al. 2013, A&A, 555, A23
  * Balsara & Kim (2005) Balsara, D. S., & Kim, J. 2005, ApJ, 634, 390
  * Balsara et al. (2004) Balsara, D. S., Kim, J., Mac Low, M.-M., & Mathews, G. J. 2004, ApJ, 617, 339
  * Beck (2015a) Beck, R. 2015a, A&A Rev., 24, 4
  * Beck (2015b) —. 2015b, A&A, 578, A93
  * Beck & Hoernes (1996) Beck, R., & Hoernes, P. 1996, Nature, 379, 47
  * Beck et al. (2003) Beck, R., Shukurov, A., Sokoloff, D., & Wielebinski, R. 2003, A&A, 411, 99
  * Bendre (2018) Bendre, A. B. 2018, PhD thesis, Universität Potsdam, http://nbn-resolving.de/urn:nbn:de:kobv:517-opus4-407517
  * Betti et al. (2019) Betti, S. K., Hill, A. S., Mao, S. A., et al. 2019, ApJ, 871, 215
  * Borlaff et al. (2021) Borlaff, A. S., Lopez-Rodriguez, E., Beck, R., et al. 2021, ApJ, 921, 128
  * Brandenburg & Dobler (2002) Brandenburg, A., & Dobler, W. 2002, Computer Physics Communications, 147, 471
  * Brandenburg et al. (2013) Brandenburg, A., Gressel, O., Käpylä, P. J., et al. 2013, ApJ, 762, 127
  * Brandenburg & Rekowski (2001) Brandenburg, A., & Rekowski, B. V. 2001, A&A, 379, 1153
  * Brandenburg & Sarson (2002) Brandenburg, A., & Sarson, G. R. 2002, Phys. Rev. Lett., 88, 055003
  * Brandenburg et al. (2020) Brandenburg, A., Johansen, A., Bourdin, P. A., et al. 2020, arXiv e-prints, arXiv:2009.08231
  * Butsky et al. (2017) Butsky, I., Zrake, J., Kim, J.-h., Yang, H.-I., & Abel, T. 2017, ApJ, 843, 113
  * Chamandy & Shukurov (2020) Chamandy, L., & Shukurov, A. 2020, Galaxies, 8, 56
  * de Avillez & Breitschwerdt (2005) de Avillez, M. A., & Breitschwerdt, D. 2005, A&A, 436, 585
  * Dickey et al. (2022) Dickey, J. M., West, J., Thomson, A. J. M., et al. 2022, ApJ, 940, 75
  * Dobbs & Price (2008) Dobbs, C. L., & Price, D. J. 2008, MNRAS, 383, 497
  * Evirgen et al. (2017) Evirgen, C. C., Gent, F. A., Shukurov, A., Fletcher, A., & Bushby, P. 2017, MNRAS, 464, L105
  * Evirgen et al. (2019) Evirgen, C. C., Gent, F. A., Shukurov, A., Fletcher, A., & Bushby, P. J. 2019, MNRAS, 488, 5065
  * Ferrière (2001) Ferrière, K. M. 2001, Reviews of Modern Physics, 73, 1031
  * Ferrière et al. (1991) Ferrière, K. M., Mac Low, M.-M., & Zweibel, E. G. 1991, ApJ, 375, 239
  * Fletcher et al. (2011) Fletcher, A., Beck, R., Shukurov, A., Berkhuijsen, E. M., & Horellou, C. 2011, MNRAS, 412, 2396
  * Gent (2012) Gent, F. A. 2012, PhD thesis, Newcastle University School of Mathematics and Statistics, http://hdl.handle.net/10443/1755
  * Gent et al. (2020) Gent, F. A., Mac Low, M.-M., Käpylä, M. J., Sarson, G. R., & Hollins, J. F. 2020, Geophysical and Astrophysical Fluid Dynamics, 114, 77
  * Gent et al. (2021) Gent, F. A., Mac Low, M.-M., Käpylä, M. J., & Singh, N. K. 2021, ApJ, 910, L15
  * Gent et al. (2023) Gent, F. A., Mac Low, M.-M., & Korpi-Lagg, M. J. 2023, Data for Transition from small-scale to large-scale dynamo in a supernova-driven, multiphase medium, v1.0, American Museum of Natural History, doi:10.5531/sd.astro.6.  https://doi.org/10.5531/sd.astro.6
  * Gent et al. (2023) Gent, F. A., Mac Low, M.-M., Korpi-Lagg, M. J., & Singh, N. K. 2023, ApJ, 943, 176
  * Gent et al. (2013a) Gent, F. A., Shukurov, A., Fletcher, A., Sarson, G. R., & Mantere, M. J. 2013a, MNRAS, 432, 1396
  * Gent et al. (2013b) Gent, F. A., Shukurov, A., Sarson, G. R., Fletcher, A., & Mantere, M. J. 2013b, MNRAS, 430, L40
  * Girichidis (2021) Girichidis, P. 2021, MNRAS, 507, 5641
  * Gressel (2008) Gressel, O. 2008, PhD thesis, Astrophysikalisches Institut Potsdam
  * Gressel et al. (2008) Gressel, O., Ziegler, U., Elstner, D., & Rüdiger, G. 2008, Astronomische Nachrichten, 329, 619
  * Hanasz et al. (2009) Hanasz, M., Wóltański, D., & Kowalik, K. 2009, ApJ, 706, L155
  * Harnett et al. (2004) Harnett, J., Ehle, M., Fletcher, A., et al. 2004, A&A, 421, 571
  * Haugen & Brandenburg (2004) Haugen, N. E. L., & Brandenburg, A. 2004, Phys. Rev. E, 70, 036408
  * Haverkorn (2015) Haverkorn, M. 2015, in Astrophysics and Space Science Library, Vol. 407, Magnetic Fields in Diffuse Media, ed. A. Lazarian, E. M. de Gouveia Dal Pino, & C. Melioli, 483
  * Hill et al. (2012) Hill, A. S., Joung, M. R., Mac Low, M. M., et al. 2012, ApJ, 750, 104
  * Hill et al. (2018) Hill, A. S., Mac Low, M.-M., Gatto, A., & Ibáñez-Mejía, J. C. 2018, ApJ, 862, 55
  * Hollins et al. (2017) Hollins, J. F., Sarson, G. R., Shukurov, A., Fletcher, A., & Gent, F. A. 2017, ApJ, 850, 4
  * Houde et al. (2013) Houde, M., Fletcher, A., Beck, R., et al. 2013, ApJ, 766, 49
  * Hutschenreuter et al. (2023) Hutschenreuter, S., Haverkorn, M., Frank, P., Raycheva, N. C., & Enßlin, T. A. 2023, arXiv e-prints, arXiv:2304.12350
  * Iffrig & Hennebelle (2017) Iffrig, O., & Hennebelle, P. 2017, A&A, 604, A70
  * Joung & Mac Low (2006) Joung, M. K. R., & Mac Low, M.-M. 2006, ApJ, 653, 1266
  * Käpylä et al. (2018) Käpylä, M. J., Gent, F. A., Väisälä, M. S., & Sarson, G. R. 2018, A&A, 611, A15
  * Käpylä et al. (2022) Käpylä, M. J., Rheinhardt, M., & Brandenburg, A. 2022, ApJ, 932, 8
  * Kazantsev (1968) Kazantsev, A. P. 1968, Soviet Journal of Experimental and Theoretical Physics, 26, 1031
  * Kim & Ostriker (2015a) Kim, C.-G., & Ostriker, E. C. 2015a, ApJ, 815, 67
  * Kim & Ostriker (2015b) —. 2015b, ApJ, 802, 99
  * Kirchschlager et al. (2023) Kirchschlager, F., Mattsson, L., & Gent, F. A. 2023, submitted Nat. Astron., doi:10.21203/rs.3.rs-2405487/v1
  * Korpi et al. (1999) Korpi, M. J., Brandenburg, A., Shukurov, A., & Tuominen, I. 1999, A&A, 350, 230
  * Kotarba et al. (2010) Kotarba, H., Karl, S. J., Naab, T., et al. 2010, ApJ, 716, 1438
  * Krause & Raedler (1980) Krause, F., & Raedler, K. H. 1980, Mean-field magnetohydrodynamics and dynamo theory
  * Kuijken & Gilmore (1989) Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571
  * Mac Low & McCray (1988) Mac Low, M.-M., & McCray, R. 1988, ApJ, 324, 776
  * Machida et al. (2013) Machida, M., Nakamura, K. E., Kudoh, T., et al. 2013, ApJ, 764, 81
  * Mannucci et al. (2005) Mannucci, F., Della Valle, M., Panagia, N., et al. 2005, A&A, 433, 807
  * Martin-Alvarez et al. (2018) Martin-Alvarez, S., Devriendt, J., Slyz, A., & Teyssier, R. 2018, MNRAS, 479, 3343
  * Norman & Ikeuchi (1989) Norman, C. A., & Ikeuchi, S. 1989, ApJ, 345, 372
  * Olovsson et al. (2005) Olovsson, L., Simonsson, K., & Unosson, M. 2005, International Journal for Numerical Methods in Engineering, 63, 1436
  * Otmianowska-Mazur et al. (1997) Otmianowska-Mazur, K., von Linden, S., Lesch, H., & Skupniewicz, G. 1997, A&A, 323, 56
  * Pakmor & Springel (2013) Pakmor, R., & Springel, V. 2013, MNRAS, 432, 176
  * Qasi et al. (2023) Qasi, Y., Shukurov, A., Tharakkal, D., Gent, F. A., & Bendre, A. B. 2023, in preparation
  * Rieder & Teyssier (2016) Rieder, M., & Teyssier, R. 2016, MNRAS, 457, 1722
  * Sánchez-Salcedo et al. (2002) Sánchez-Salcedo, F. J., Vázquez-Semadeni, E., & Gazol, A. 2002, ApJ, 577, 768
  * Sarazin & White (1987) Sarazin, C. L., & White, III, R. E. 1987, ApJ, 320, 32
  * Shukurov & Subramanian (2021) Shukurov, A., & Subramanian, K. 2021, Astrophysical Magnetic Fields: From Galaxies to the Early Universe, Cambridge Astrophysics (Cambridge University Press), doi:10.1017/9781139046657
  * Tabatabaei et al. (2008) Tabatabaei, F. S., Krause, M., Fletcher, A., & Beck, R. 2008, A&A, 490, 1005
  * Tenorio-Tagle & Bodenheimer (1988) Tenorio-Tagle, G., & Bodenheimer, P. 1988, ARA&A, 26, 145
  * Tharakkal et al. (2023) Tharakkal, D., Shukurov, A., Gent, F. A., Sarson, G. R., & Snodin, A. 2023, arXiv e-prints, arXiv:2305.03318
  * Tsouros et al. (2023) Tsouros, A., Edenhofer, G., Enßlin, T., Mastorakis, M., & Pavlidou, V. 2023, arXiv e-prints, arXiv:2303.10099
  * van den Bergh & Tammann (1991) van den Bergh, S., & Tammann, G. A. 1991, ARA&A, 29, 363
  * Walters & Cox (2001) Walters, M. A., & Cox, D. P. 2001, ApJ, 549, 353
  * Wolfire et al. (1995) Wolfire, M. G., Hollenbach, D., McKee, C. F., Tielens, A. G. G. M., & Bakes, E. L. O. 1995, ApJ, 443, 152

##  Appendix A Supplementary material

###  A.1 Vertical density structure and effect of clustering

In Figure 17 the time-latitude diagrams of gas density are shown for Model L2
without and Model L2-cl with OB clustering. The effect of clustering is to
more efficiently evacuate the gas from the lower halo. The secondary effect on
the topology of the magnetic field shown in panel _(c)_ is to confine the mean
field to the disc. Compare panel _(c)_ to Figure 5 for Model L2.

L2L2-clL2-cl(a)(b)(c)

Figure 17:  Time-latitude diagrams of horizontally averaged gas number density
for _(a)_ Model L2 and _(b)_ L2-cl and _(c)_ $B_{y}$ for Model L2-cl.

Figure 18 shows SSD and LSD growth rates for Model L2-cl, for comparison with
the growth rates in Figure 6 for Model L2. The latter 300 Myr interval shifted
0.42 Myr corresponds to the LSD stage fitted for Model L2 in Figure 13.

Figure 18:  Model L2-cl: magnetic energy density normalised by time-averaged
turbulent kinetic energy ${\overline{e_{K}^{\prime}}}$. Total magnetic energy
$\langle B^{2}\rangle(8\uppi)^{-1}$ (_blue_) and horizontally-averaged mean-
field energy ${\langle B\rangle_{xy}}^{2}(8\uppi)^{-1}$ (_green_) have fits
for exponential growth rate $\gamma$ as listed in the legend, spanning
$0.2~{}{\rm Gyr}<t<0.48~{}{\rm Gyr}$ and $0.48~{}{\rm Gyr}<t<0.78~{}{\rm
Gyr}$, respectively.

###  A.2 Convergence

L2H2L2H2

Figure 19:  Time-latitude diagrams of horizontally averaged $B_{x}$ and
$B_{y}$ immediately after the SSD has initially saturated for Models L2 and
H2, as indicated by the labels in each panel. Note the different starting time
for the two models, chosen to capture the saturated SSD in each case.

A further indication of the level of convergence of the solution (see Sect.
3.3)is a comparison of the topology of the mean magnetic field produced
between Models H2 and L2. In Figure 19 the time-latitude diagrams for the
horizontally-averaged magnetic field components are plotted for both models.
During this stage of the LSD, the field strength, orientation, and topologies
are qualitatively similar prior to the later reversal in the sign of the
azimuthal field in Model L2. Model H2 has not yet progressed far enough to
determine whether such a field reversal occurs at high resolution.

(a)(b)(c)

Figure 20:  _(a)_ Magnetic energy spectra normalised by maximum power in each
snapshot. Wavenumber $k$ corresponding to maximal power for each snapshot are
listed in the legend. _(b)_ Kinetic energy spectra. The snapshot time in Gyr
is listed in the legend. Both panels show the matching snapshots for matching
line style and colour. Spectra have been smoothed with a Gaussian kernel of
length 5 Myr and 1.5 $k_{1}$ in time and wavelength. _(c)_ Time evolution of
magnetic energy over sampled scales of $k$, as noted in the legends in
multiples of $k_{1}=2\uppi{L_{x}}^{-1}$.

In Figure 20 from Model H2 magnetic (a) and kinetic (b) energy spectra are
shown for comparison to the equivalent spectra from Model L2 in Figure 7. The
peak wavenumber for the magnetic spectra during the kinematic stage of the SSD
are at $k\simeq 17k_{1}$ for L2 but vary $16k_{1}\leq k\leq 26k_{1}$ at high
resolution contributing to the greater efficiency of the SSD. Although both
models fluctuate over time, the peak energy for the kinetic spectra at low $k$
are almost an order of magnitude higher for H2 than for L2, which has
previously been shown for SSD in Gent et al. (2021). The inertial range for
both models follows a consistent $k$-scaling, although extending much further
at high resolution, as would be expected. Critically, lower resolution affects
not only the energetics of the small-scale turbulence, but also the systemic
flows. Nevertheless, as can be seen from Figure 19 these differences do not
appear to affect the horizontally-averaged structure of the large-scale field.

Figure 20(c) shows how different scales $k$ within the magnetic field grow
over time for Model H2, to compare with that of Model L2 in Figure 8(b).
Saturation of the SSD occurs at about $10^{-17}~{}{\rm erg}~{}{\rm cm^{-3}}$
for all $k\gtrsim 17k_{1}$ in both models, after which the hierarchy between
these scales is conserved in the nonlinear stage of the SSD. For $k\lesssim
17k_{1}$ the energy continues to grow, and the hierarchy between these scales
are eventually inverted relative to their hierarchy during the kinematic stage
of the SSD. This spectral evolution is not only qualitatively consistent, but
quantitatively compatible independent of resolution subsequent to the
kinematic stage of the SSD.

Gent et al. (2021) demonstrate that SSD solutions converge for resolution
$\delta{{\mbox{\boldmath$x$}}}\lesssim 1~{}{\rm pc}$. The magnetic energy
spectra plotted in Figure 20 _(a)_ show that the peak energy during the SSD
varies around slightly smaller scale, $17k_{1}\lesssim k\lesssim 27k_{1}$
rather than more consistently near $17k_{1}$. Comparing the kinetic energy
spectra between Figure 7 _(b)_ and Figure 20 _(b)_ it is evident that tail of
the inertial range occurs $k$ just below $17k_{1}$ for Model L2, but continues
as far as $25k_{1}$ for Model H2. This extends the scales at which kinetic
energy can be transferred to the magnetic field, but in addition the growth of
the magnetic field is more rapid at higher resolution at all scales up to
$k=k_{1}$ (Figure 4 of Gent et al., 2021). The primary driver of this enhanced
SSD is the increased vorticity, lower sound speed, and higher fractional
volume of the hot gas available at higher resolution (Gent et al., 2023).

L2L2

Figure 21:  Time-latitude diagrams of horizontally averaged
$\alpha_{K}$-effect _(a)_ and $\alpha_{M}$-effect _(b)_ for Model L2. A linear
colour-scale is applied $-5~{}{\rm km~{}s}^{-1}<\alpha<5~{}{\rm km~{}s}^{-1}$
and symmetric logarithmic colour-scale otherwise.

Figure 20 _(c)_ shows how the magnetic energy evolves at sampled wavenumbers
$k$. At the onset of the transition from SSD to LSD, the growth rates at the
largest scales of the magnetic field are the first to slow down and separate
from the eigenmode. From this point on scales larger than $k=17k_{1}$ grow
more rapidly, while scales smaller than $17k_{1}$ grow more slowly or
stagnate. Up to 200 Myr, $k=17k_{1}$ remains the scale with the highest
energy, but is thereafter surpassed.

H2H2

Figure 22:  Time-latitude diagrams of horizontally averaged
$\alpha_{K}$-effect _(a)_ and $\alpha_{M}$-effect _(b)_ for Model H2. A linear
colour-scale is applied $-5~{}{\rm km~{}s}^{-1}<\alpha<5~{}{\rm km~{}s}^{-1}$
and symmetric logarithmic colour-scale otherwise.

###  A.3 $\alpha$ coefficients

In Figures 21 and 22 for Models L2 and H2, respectively, we show the separate
contributions of kinetic $\alpha_{K}$ and magnetic $\alpha_{M}$. For Model L2
the systemic sign of the $\alpha_{K}$ either side of the midplane is permeated
by small-scale reversals, especially in the SN active region nearest the
midplane. The strength of the kinetic helicity and the reversals dissipate
over time. The magnetic $\alpha_{M}$ is negligible almost to 2 Gyr.
Subsequently it has large-scale structure away from the midplane, typically
with same sign to the kinetic $\alpha_{K}$. These correspond to the
intersections between outflows and inflows and are dynamically weak, occuring
in diffuse gas with weak magnetic energy.

For Model H2 we see that $\alpha_{M}$ at the midplane is initially negligible.
When it becomes dynamically important after about 400 Myr it has the same sign
as that typical of $\alpha_{K}$, thus enhancing the LSD. In the outflows,
however, we see that $\alpha_{K}$ tends to be of opposite sign to the large-
scale $\alpha$, so also has a positive effect on the LSD by supporting the
removal of small-scale magnetic helicity.

Table 2:  Meanings of variables Symbol |  Denoting |  Units/Definition
---|---|---
$U$ |  galactic shear velocity |  [cm s-1]
$\Omega=-S$ |  angular velocity, rate of shear |  [s-1]
$u$ |  deviation of gas velocity from shear |  [cm s-1]
${{\mbox{\boldmath$u$}}}^{\prime}$ |  gas turbulent velocity |  [cm s-1]
$\rho$ |  gas density |  [g $~{}{\rm cm^{-3}}$]
$t$ |  time |  [s]
$T$ |  gas temperature |  [K]
$s$ |  specific entropy |  [cm2 s-2 K-1]
$A$ |  magnetic vector potential |  [G cm]
$B$ |  magnetic field |  [G]
$J$ |  current density |  [erg cm-4 G-1]
$\dfrac{\text{D}~{}}{\text{D}t}$ |  advective (material) derivative |  $\dfrac{\partial~{}}{\partial t}+({{\mbox{\boldmath$U$}}}+{{\mbox{\boldmath$u$}}})\cdot\nabla$
$\zeta_{D},\zeta_{\nu},\zeta_{\chi}$ |  shock diffusion coefficients |  $\propto{\cal C}=-\left(\nabla\cdot{{\mbox{\boldmath$u$}}}|_{-\rm ve}\right)$
$\nu,\eta$ |  viscosity, resistivity coefficients |  [cm2 s-1]
$\nu_{6},\chi_{6},\eta_{6}$ |  hyperdiffusion coefficients |  [cm6 s-1]
$\mu_{0}$ |  vacuum magnetic permeability |  1
$c_{\rm s}$ |  sound speed |  [cm s-1]
$c_{\rm p}$ |  specific heat at constant pressure |  [cm2 s-2 K-1]
$c_{v}$ |  specific heat at constant volume |  [cm2 s-2 K-1]
$\dot{\sigma}$ |  SN explosion rate |  [cm-2 s-1]
$\dot{\sigma}_{\rm Sn}$ |  Solar neighbourhood SN rate |  [cm-2 s-1]
$\Gamma$, $\Gamma_{\rm UV}$ |  far ultraviolet heating function |  [cm2 s-3]
$\Lambda$ |  radiative cooling |  [cm5 s-3 g-1]
$E_{\rm th}$ |  SN explosion energy |  1051 [erg]
$E_{\rm m}$ |  SN explosion momentum injection |  $\sum_{\rm Sn}\rho\Delta{{\mbox{\boldmath$u$}}}$ [g cm s-1]
$h$ |  vertical domain size |  $\leq$3.072
$h_{\rm II}\,,h_{\rm I}$ |  scale height for SN Type II, I |  0.09, 0.325 kpc
$e_{B}$ |  magnetic energy density |  [erg cm-3]
$e_{K}$ |  total kinetic energy density |  [erg cm-3]
$e_{K}^{\prime}$ |  turbulent kinetic energy density |  [erg cm-3]
$\gamma$ |  volume-averaged ${\langle e_{B}\rangle}$ growth rate |  [s-1]
$k_{1}$ |  reference spectral wavenumber |  $\sim 12.3$ kpc-1
$\alpha$ |  inductive term in EMF |  [cm s-1]
$\nabla$ |  gradient vector |  e.g., $\left(\dfrac{\partial~{}}{\partial x},\dfrac{\partial~{}}{\partial y},\dfrac{\partial~{}}{\partial z}\right)$
$\nabla^{2n}$ |  2nd or 6th order Lapacian |  e.g., $\dfrac{\partial^{2n}~{}}{\partial x^{2n}}+\dfrac{\partial^{2n}~{}}{\partial y^{2n}}+\dfrac{\partial^{2n}~{}}{\partial z^{2n}}$, $n=1,3$
W |  traceless rate of strain tensor |  ${\mathsf{W}}_{ij}=\dfrac{1}{2}\left(\dfrac{\partial u_{i}}{\partial x_{j}}+\dfrac{\partial u_{j}}{\partial x_{i}}-\dfrac{2}{3}\delta_{ij}\nabla\cdot{{\mbox{\boldmath$u$}}}\right)$
$|\textbf{{W}}|^{2}$ |  contraction of W |  $|\textbf{{W}}|^{2}={\mathsf{W}}_{ij}{\mathsf{W}}_{ij}$
$\textbf{{W}}^{(5)}$ |  5th order rate of strain tensor |  ${\mathsf{W}}_{ij}^{(5)}=\dfrac{1}{2}\left(\dfrac{\partial^{5}u_{j}}{\partial x_{i}^{5}}+\dfrac{\partial^{4}}{\partial x_{i}^{4}}\left(\dfrac{\partial u_{i}}{\partial x_{j}}\right)-\dfrac{1}{3}\dfrac{\partial^{4}}{\partial x_{i}^{4}}\left(\nabla\cdot{{\mbox{\boldmath$u$}}}\right)\right)$

##  Appendix B Code and parameter evolution during run time

Table 3:  Tracking parameter changes over duration of simulations: For the models listed the parameters applied commencing at time $t$ are shown: $h_{\rm adj}$ indicates the scale height of the SN vertical distribution varies with the scale height of the gas; $\Gamma_{\rm UV}$ is the photoelectric heating factor relative to the the Habing rate $\Gamma_{0}$. $T_{\rm max}$ is the maximum temperature allowed with $r_{0}$ of the SN explosion, and $T_{\rm ratio}\times T_{\rm max}$ is the maximum allowed within 2.25$r_{0}$ of the explosion. $n_{\rm cool}$ indicates whether cooling mass is used to reduce the temperature, otherwise the SN site is rejected and an alternative selected. $n_{\rm max}$ is the maximum gas number density at the grid location of an explosion, to avoid locations where SN evolution is not well resolved at ${\delta x}=4~{}{\rm pc}$ $n_{\rm ratio}$ is the maximum ratio between maximum and minimum gas density within $r_{0}$ of an explosion, to avoid excessive viscous forces.  Model | $t$ | $h_{\rm adj}$ | $\Gamma_{\rm UV}$ | $T_{\rm max}$ | $n_{\rm max}$ | $n_{\rm ratio}$ | $T_{\rm ratio}$ | $E_{m}$ | $n_{\rm cool}$ | $\nu_{6}\,\chi_{6}\,\eta_{6}$
---|---|---|---|---|---|---|---|---|---|---
| [Gyr] |  |  [$\Gamma_{0}$] |  [$10^{6}$ K] |  [cm-3] |  [$n_{\rm max}{n_{\rm min}}^{-1}$] |  [${T_{\rm max}}^{-1}$] |  |  |  [kpc-5 km s-1]
H2 | 0 | F | 1 | 25 | 50 | 2e4 | 8e2 | T, 7.5% | F | 1e-15
H2 | 0.02 | F | 1 | 25 | 150 | 3e4 | 9e2 | T, 7.5% | F | 1e-15
H2 | 0.028 | T | 2.5 | 25 | 150 | 3e4 | 9e2 | T, 7.5% | F | 1e-15
H2 | 0.037 | T | 3.5 | 25 | 100 | 1e4 | 7 | T, 7.5% | T | 1e-15
H2 | 0.465 | T | 3.5 | 25 | 100 | 1e4 | 7 | T, 7.5% | T | 1.1e-15
H2 | 0.652 | T | 3.5 | 7.5 | 100 | 1e4 | 12 | F | F | 1.2e-15
H2 | 0.672 | T | 3.5 | 50 | 100 | 1e4 | 10 | F | F | 1.2e-15
L2, L2-cl, |  |  |  |  |  |  |  |  |  |
L0-B0, L2-B0 | 0 | F | 3.5 | 12.5 | 50 | 1e4 | 7 | T, 7.5% | T | 5e-12
L2 | 0.365 | F | 3.5 | 12.5 | 50 | 1e4 | 7 | T, 7.5% | F | 5e-12
L2 | 0.461 | F | 3.5 | 12.5 | 50 | 1e4 | 7 | T, 7.5% | T | 5e-12
L2 | 3.272 | T | 3.5 | 25 | 100 | 1e4 | 7 | T, 7.5% | T | 5.2e-12
L2 | 3.435 | T | 3.5 | 25 | 100 | 1e4 | 7 | T, 7.5% | T | 5.25e-12
L2-cl | 0.1 | T | 3.5 | 25 | 50 | 1e4 | 4 | T, 7.5% | T | 5e-12
L2-cl | 0.128 | T | 3.5 | 20 | 50 | 1e4 | 5 | T, 7.5% | T | 5e-12
L2-cl | 0.349 | T | 3.5 | 25 | 50 | 1e4 | 4 | T, 7.5% | T | 5e-12
L2-cl | 0.375 | T | 3.5 | 12.5 | 50 | 1e4 | 7 | T, 7.5% | T | 5e-12

A summary of the symbols used in the manuscript are listed in Table 2. It is
in the nature of such simulations, which require long integration time,
spanning several months and longer in real time, that new code adaptations to
stability issues or improvements to algorithms are discovered and applied
while the simulations are progressing. Such is the expense in terms of
computation and storage that it is not feasible to restart the simulations
with the new implementation and repeat through to saturation of the LSD. This
is particularly the case for Model H2, which has used tens of millions of
computational billing units.

Prior to the suite of simulations presented in this study, we explored a
variety of lower resolution runs with larger and smaller domains, and a
variety of far ultraviolet heating rates and vertical distributions of SNe, so
that the parameters chosen during this suite of experiments has been fairly
robust and consistent. Parameter adjustments are listed for each Model in
Table 3. Nevertheless, some updates, as detailed in Table 4, needed to be
applied during the course of the experiment. When these updates have been
committed to address our specific issues, we have also acquired any updates
applied by other developers between our commits, which are not particular to
SN-driven turbulence, but nevertheless may affect code performance.

Table 4:  Tracking code version over duration of simulations: date, repository
hash code at github.com/pencil-code, model, and final column summary of
update. The start time in Gyr when each version applies is listed for each
model.
†: Model L2 originally initiated with $\dot{\sigma}=\dot{\sigma}_{\rm Sn}$,
inducing a higher SSD growth rate. Beyond 0.3 Gyr continued with
$\dot{\sigma}=0.5\dot{\sigma}_{\rm Sn}$. The period 0–0.365 Gyr runs
retrospectively with higher initial field strength to replace the original
data.
**: Update directly affects the included models. The cooling ejecta is absent
only for a small span of three models: H2 $t<37$ Myr; L2 $365<t<461$ Myr and
L2-cl $t<208$ Myr, and will have little effect on the properties of the
turbulence.
General changes between versions could affect the timestep or precise realisations of the inherently chaotic calculations. The effect does not alter the statistical properties of the simulations.  Date | hash | H2 | L2 | L2-cl | L0-B0 | L2-B0 |
---|---|---|---|---|---|---|---
2021-01-22 | 8b05b3884 | 0 |  |  |  |  | Machine specific
2021-04-13 | 163e3e15d | 0.025 |  | 0 |  |  | Citations update
2021-10-14 | f7bb7be7a |  |  0† |  |  |  | Test solutions
2021-06-01 | 92794628c | 0.033 | 0.365 | 0.1 |  |  | Complex numbers
2021-06-03 | b482858d2 |  | 0.461 | 0.208 |  |  | Added cooling ejecta**
2021-06-04 | 533ad9dad |  | 0.584 | 0.334 |  |  | Add no mass in dense gas
2021-06-15 | 9426d1a5f |  | 1.017 | 0.513 |  |  |
2021-06-02 | f7bb7be7a | 0.037 | 1.237 |  | 0 | 0 | Test fields initialized
2022-10-18 | 6b8c45449 | 0.192 |  |  |  |  | Streamlined SN location
2022-10-19 | 524ed6ca2 | 0.196 |  |  |  |  | Omit redundant MPI call
2022-11-08 | 81cc4c1d0 | 0.205 |  |  |  |  |
2022-12-12 | 2b0563676 | 0.216 |  |  |  |  | Updated SN list length
2023-02-15 | c4b34885f | 0.368 |  |  |  |  |  Use log $T$ if calculated
2023-04-11 | 86e11c9c6 | 0.651 |  |  |  |  | Added diagnostics
end time |  | 0.672 | 5.500 | 1.181 | 2.195 | 2.008 |

Here, we seek to explain what issues have arisen in our experiments, how these
have been addressed and apply to each model, and list which version of the
code has been used during each experiment. The primary causes of code
instability arise from the location of specific SN remnants, which are
difficult to stably resolve. This is most likely when applied to a high
density region, in which momentum injection is required, but may also include
neigbourhoods of low density and therefore high temperatures and high viscous
stresses. Usually restarting from a recent snapshot with an updated random
seed for SN locations avoids this location and advances the simulation with a
statistically equivalent realisation. High density sites that do not cause
stability crises continue to be included. This is more common at low
resolution. At high resolution it is easier to resolve more challenging SN
remnant structures, although occasionally a seed reset is needed. This does
not substantially affect the comparability of the simulations over time or
between models.

One update to the code that was included during the experiments is to permit
the use of additional SN mass ejecta, where the temperature exceeds a limit of
75 MK to cool the site to less than this temperature and therefore maintain a
longer timestep. Prior to this change these locations would be avoided and an
alternate site without such high temperatures arising randomly selected, so
that remnants in diffuse locations may have been less frequent. This change
should have a minor effect on the statistics of the dynamo, as these locations
are in any case selected rarely and higher temperatures continue to arise,
when remnants expand into adjacent diffuse ISM, rapidly heating the gas
outside the initial SN remnant radius.

Another occasional cause of code instability is too low hyperdiffusivity
coefficients. The values for $\nu_{6}$, $\chi_{6}$ and $\eta_{6}$ reported in
Section 2.5 were adopted following experiments in periodic unstratified
simulations, where they adequately resolved grid-scale instabilities during
the linear and nonlinear stages of the SSD. As our simulations approached the
nonlinear stage of the LSD, we discovered some grid-scale instabilities that
could not be attributed to remnant location, as they occurred far away from
any explosions. These were solved by increasing the hyperdiffusivity
coefficients modestly. It is unlikely that these affect the statistics of the
dynamo, because they still apply only at the very largest wavenumbers in the
model, where the energies are very weak compared to the energies at even the
tail end of our inertial range, let alone those in the LSD range of
wavenumbers, as is evident in the later spectra of Model L2, which are
displayed in Figure 7 and Model H2, in Figure 20(a). In Model H2 these
coefficients were increased to $1.2\times 10^{-15}$ after 652 Myr. In this
study, $5\times 10^{-12}$ was applied throughout for
${\delta{{\mbox{\boldmath$x$}}}}=4~{}{\rm pc}$, except as listed in Table 3 at
late times for Model L2.