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

Magnetic fields are 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 the low explicit Lagrangian resistivity required. 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.


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. 2017Evirgen et al. , 2019)).Observations of external spiral galaxies often show magnetic fields aligned along the disc and parallel to the spiral arms, or otherwise ordered over kpc scales (e.g., Beck et al. 2019, their Table 4).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.
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. 2008b;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  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 a 52 second animation of this Figure, demonstrating the rapid transition to turbulence, the saturation of the turbulent magnetic energy within about 30 Myr and the much slower organisation into large scales of the azimuthal component of the magnetic field up to 652 Myr.The breathing and sloshing modes of the disc are also apparent.The vertical striping evident in the magnetic energy occurs from inflows through the boundary which has a vertical magnetic field condition.A high resolution version is included in the data repository (Gent et al. 2023).
produced by the SSD is either absent or underresolved.
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 ≲ 1 pc and that the SSD is present if the microphysical magnetic diffusivity η ≲ 5 • 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, some insights into the growth rates and drivers of the LSD, and the effect of SN clustering.In Section 4 we consider the effect of resolution and the robustness of the numerical setup.In Section 5 we discuss how the results relate to previous understanding of the galactic magnetic field and in Section 6 we summarise our main conclusions.We provide sup- plementary illustrative material in Appendix A, including a reference table of notation employed in this manuscript, and outline in Appendix B what versions of the Pencil Code were used during the months and years over which the models were evolved.Data from the models presented in this paper are available in a repository at the American Museum of Natural History Digital Library (Gent et al. 2023).

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 δx = 4 pc, aside from one high-resolution model with δx = 1 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 nonadiabatic hydrostatic equilibrium.This is solved numerically in a one-dimensional vertical column of the ISM.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 magnetic diffusivity η = 8 • 10 −4 kpc km s −1 and the viscosity ν = 4 • 10 −3 kpc km s −1 to be low enough to support a vigorous SSD, although numerical diffusivity remains dominant in our low resolution models (see Gent et al. 2021, for details of how this was demonstrated).Thus, the nominal magnetic Prandtl number in our models is Pm = 5.We omit explicit thermal conductivity, which acts at time and length scales below our model resolution.

Dynamical equations
We use the Pencil Code (Brandenburg & Dobler 2002;Pencil Code Collaboration et al. 2021) to model SN-driven turbulence as described previously in Gent (2012) and Gent et al. (2020).
In this study we solve The ideal gas equation of state closes the system, assuming an adiabatic index (ratio of specific heats) c 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 tabulated within Appendix B.

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 with a s = 1.4 km s −1 Myr −1 , z s = 0.2 kpc, a h = 0.5 km s −1 Myr −1 and z h = 1 kpc.In the local shearing periodic box approximation (e.g., Goldreich & Lynden-Bell 1965;Brandenburg et al. 1995;Yang et al. 2009) the angular velocity Ω is treated as constant.Coriolis 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 km s −1 in the warm gas, the Jeans length λ J ≫ 1 kpc and more so for the diffuse hot gas.As we do not allow the cold gas to cool below ∼ 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) Ω = (0, 0, Ω) with angular momentum Ω ∝ Ω Sn and the Solar neighborhood value Ω Sn = 30 km s −1 kpc −1 , assuming an orbital speed for the Sun of 255 km s −1 at galactocentric radius 8.5 kpc.The advective derivative, includes transport by an imposed shear flow U = (0, Sx, 0), in which the rate of shear S = −Ω assuming a flat rotation curve, and a perturbation velocity u that represents the local deviation from the overall rotational velocity U .

Cooling and heating
The uniform far ultraviolet heating follows the functional form given by Wolfire et al. (1995) with temperature dependence approximated as with Γ 0 = 0.0147 erg g −1 s −1 .However, we apply an enhanced heating rate Γ = 3.5 Γ 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 Λ(T ) = Λ k T β 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. 2008b) and see Gent et al. (2013b, their Table 1), but with Λ = 0 at T < 90 K.This cooling law includes thermally unstable branches at high and low temperatures as expected for the ISM.

Artificial diffusivities
Terms in the dynamical equations that contain ζ D = 2C, ζ ν = 5C and ζ χ = 2C resolve shock discontinuities with artificial diffusion of mass, momentum, and thermal energy, respectively.They depend on the local convergence C = −∇ • u| −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 ζ D applied in Equation (1).Following Gent et al. (2021), resistive shock diffusion is omitted from Equation (3), where hyperdiffusion is sufficient to resolve gridscale instabilities without suppressing the important physical process of field compression in SN remnant shells.
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.

Supernova explosions
For each SN explosion thermal energy E th = 10 51 erg is instantaneously added to a source region with spherical Gaussian distribution of radial scale r 0 .For models with δx = 4 pc, r 0 ≥ 20 pc, while with δx = 1 pc, r 0 ≥ 12 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 M ⊙ to avoid excessive temperatures negatively impacting the timestep.
In dense regions at δx = 4 pc up to 5% of the energy is instead added via momentum injection, E m to compensate for excessive cooling (see Kim & Ostriker 2015b;Gent et al. 2020).This applies an axially symmetric diverging flow, ∆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 ∼ 2M ⊙ per SN and ≲ 40M ⊙ for any single remnant.
SN Type II and I are distinguished only by their Poisson rates of σII = 15 kpc −2 Myr −1 and σI = 2 kpc −2 Myr −1 , respectively.The model SN rate is therefore σ = σII + σI = 0.5 σ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 II = 90 pc and h I = 325 pc, respectively.The SN rate σ 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 ≳ 6h 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 peak = max(⟨n⟩ xy )/(1 cm −3 ).The Type II SN scale height also varies in response to the thickness of the disc.It is assumed that n peak would tend to be 1.The scale height of the probability distribution is adjusted to h II /n peak .The model total mass is conserved, so n peak < (>)1 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 ≲ 0.5 kpc −2 Myr −1 .With the Model Type II SN area rate ≃ 15 kpc −2 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.

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

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 ℓ indicates that the quantity has been Gaussian smoothed over a kernel of length ℓ 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 kinetic α-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 u ′ , which is most relevant to the SSD.The turbulent velocities are derived from the horizontally averaged velocities as Similarly the kinetic energy also includes mean flows, so to isolate the contribution from the turbulent kinetic energy from the horizontally averaged quantities we use 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 e K or e ′ K only the time interval after the magnetic field has passed its minimum energy is averaged to exclude initial transients in the kinetic quantities.

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 ⟨B⟩ as in which E contributes the electromotive force (EMF) from averaging turbulent fluctuations 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 in which tensors α and β are second order and κ third order, and each term represents a physical process acting on the dynamo.The α tensor applies effects from small-scale helicity, β turbulent diffusivity, δ turbulent pumping, γ shear or rotation current effects, and κ includes residual effects, dependent on the symmetric part of the magnetic gradient tensor ∇⟨B⟩) (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 α term.The action of the Coriolis force on the sheared, vertically stratified disc will tend to inject helicity into the systemic vertical flows of opposite sign on either side of the midplane.The dynamo amplifies a magnetic field with opposing small-scale helicity, which will quench the dynamo, if it cannot be removed.(Brandenburg & Subramanian 2005;Shukurov & Subramanian 2021) Under the assumption of isotropic, homogeneous turbulence α can be reduced to a scalar The kinetic helicity contributes to the LSD as where τ is the correlation time of the turbulence and ω ′ = ∇ × u ′ is its vorticity (Moffatt 1978;Krause & Rädler 1980).Due to the conservation of magnetic helicity, for the LSD to exist some small-scale helicity flux is required (Pouquet et al. 1976;Subramanian & Brandenburg 2004).The magnetic helicity contributes to the LSD as 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 α M can lead to αquenching (Krause & Rädler 1980) of the LSD.

Small-scale dynamo
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.
Over a range of models and regimes the SSD is constrained to saturate at or below equipartition with the kinetic energy density, but has so far been found to actually saturate well below Volume-averaged magnetic energy density for models L2-B0 (blue) with Ω = 2Ω Sn and L0-B0 (green) with Ω = 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 e ′ K .The horizontal dotted line indicates 2% equipartition.
equipartition.Archontis et al. (2007) obtain an intermittently superequipartition solution with a smooth flow, but it reduces to ∼ 10% in the turbulent regime.Federrath et al. (2011) report saturation levels of only 0.8%-5%of equipartition for isothermal forced turbulence at sonic Mach number M = 2, reaching the higher value with solenoidal forcing.This drops to 2% for M ≥ 10.As high as 60% is possible only with incompressible flows and solenoidal forcing.Schober et al. (2015) derive a semi-analytic model prediction of 1.3% to 43.8% at Pm > 1, with the lower end applying for highly compressible flows.
The saturation of the magnetic field for Model L0-B0 occurs at about 2% of equipartition with the time-averaged turbulent kinetic energy density e ′ K of Equation ( 10).The saturation ratio is therefore a signature of the compressibility at high Rm (see also, e.g., Schekochihin et al. 2002;Gent et al. 2021) and is likely to evolve rapidly in the ISM to a minimum of 1-2% equipartition.The SSD in Model L0-B0, without shear, grows with exponential growth rate γ = 24.6Gyr −1 , while Model L2-B0 grows slightly faster, with γ = 26.7 Gyr −1 , and saturates at a higher level.Large-scale shear and rotation, therefore, slightly enhance the SSD.

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 3 displays the evolution of the mean magnetic field as horizontal averages of the magnetic field components ⟨B x ⟩ xy and ⟨B y ⟩ 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 ⟨B y ⟩ xy is strongly positive and an order of magnitude stronger than mean ⟨B x ⟩ xy .In the period beyond 3.9 Gyr, after which the mean field remains quite steady, along the midplane ⟨B x ⟩ xy = −0.17±0.20 µG, while ⟨B y ⟩ xy = 9.16 ± 3.56 µG, with errors given by the standard deviation.
In Figures 3(c) and (d), horizontal averages of B x and B y are shown normalised by their timedependent 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).
All profiles of Figure 3 exhibit oscillatory structure.These correspond to the breathing and sloshing modes predicted by Walters & Cox (2001), due to the gravitational interaction of the stellar disc with the gas outflows.(This is also shown in gas number density in Appendix A.2 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 4(a) we plot the growth of the total magnetic energy (blue) for Model L2.The mean (red ), with ℓ = 50 pc.(Factors of (8π) −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 e ′ K (blue dotted ), while the total kinetic energy is given by e K (red dotted ).Exponential fits for the growth rates are listed in the legends, applying to intervals in (a) 365 Myr < t < 800 Myr, 0.9 Gyr < t < 1.2 Gyr and 2.5 Gyr < t < 4.8 Gyr, respectively, and for (b) spanning 12 Myr < t < 19 Myr and 0.11 Gyr < t < 0.35 Gyr, respectively.field energy is also plotted as horizontal averages ⟨B⟩ xy and Gaussian smoothed ⟨B⟩ ℓ , with a smoothing length ℓ = 50 pc, as described previously (Gent et al. 2013b).Fits for the exponential growth rate γ are listed in the legend for intervals spanning the SSD, a transitional stage, and the LSD.
The growth rate of the total magnetic energy during the SSD matches that of Model L0-B0 at 24.6 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 kinematic stage, the mean field magnetic energy 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.
Before the total magnetic energy reaches 2% e ′ K , the level at which the SSD saturates for Model L0-B0 and L2-B0 (Figure 2), the meanfield energy growth slows down.The SSD energy does not saturate at this point, but continues to grow at a smaller rate γ = 4.2 Gyr −1 above 2% e ′ K .The strength of the growth of the fluctuating field energy is unexpected.The mean field energy remains at least two orders of magnitude weaker than the turbulent field energy.Tangling must therefore be capable of generating turbulent field at rates substantially larger than the source.The mean field energy grows at a rate γ between 7 and 14 Gyr −1 , but this slows down to 1 Gyr −1 as it approaches 2% e ′ K .
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).
In Figure 4(b) we show the evolution of magnetic energy of Model H2.For comparison with Model L2 the intervals relevant to the kinematic SSD have γ = 1307 Gyr −1 , while the transition phase has γ = 5.6 Gyr −1 .(The high resolution model H2 has not yet completed the transition to the LSD.)The extremely rapid kinematic growth at high resolution is quite consistent with our previous results for SSD in the multiphase ISM (Gent et al. 2021(Gent et al. , 2023)).The near self-similar growth rates during the kinematic SSD within each model suggest 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.This is reflected in the small scale structure of the horizontal averages during this period as shown in Figure 3(c)-(d).
In modelling of SSD and LSD evolution at high Rm with an α 2 dynamo with helical forcing, Subramanian & Brandenburg (2014) and Bhat et al. (2016) find no LSD during the kinematic phase of the SSD, with single eignemode growth at all scales.Only in the nonlinear stage of the SSD do the growth rates at large scales exceed growth at small scales.Even when shear is included Bhat et al. (2019, their figures 1 and 2) show no distinction between the growth rates at each scale until the SSD becomes nonlinear.With the same presentation (see our Appendix) we could infer the same conclusion.However, this would be incorrect since for Model L2 we Figure 5. Magnetic energy spectra during epochs of the SSD as indicated in each panel (L2 is shifted 365 Myr, as discussed in Sect.4.1).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 centered about t in time and 1.5 k 1 about k.
do find the LSD enhances the growth rate at large scales throughout the kinematic SSD.
This enhancement is evident by inspection of Figure 5.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.An expanded examination of the kinematic stage in Bhat et al. (2019, figure 10) might therefore show a similar trend.
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 ≃ 17k 1 .In common with Haugen et al. (2004); Cho & Ryu (2009); Bhat & Subramanian (2013) and Eyink et al. (2013), all models show a shift in energy toward lower k as the SSD saturates, with the power law → k 1 scaling at low k.Only for Model L2 with LSD does the power law approach k 0 at late times.

The α-effect
We approximate the proxy for α = α K + α M in Eqs. ( 14) and (15) from the time series of horizontal averages with and where the first expression on the right hand side is used because ⟨(B • ∇ × B)ρ −1 ⟩ xy is unavailable, and in the second ϵ ijk denotes the Levi-Cevita symbol and summation over repeated indices i, j, k is assumed.In Figure 6 this α 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 τ = 10 Myr in Equation ( 14) to calculate both α K and α M for both models.This is consistent with an application of the test-field method (Bendre 2018) and yields rms values for α in Models L2, H2 and L2cl of 15, 30, and 89 km s −1 , which correspond reasonably to calculations of u ′ included in Figure 7(c).This value of τ is slightly higher than the 5 Myr determined by Hollins et al. (2017), the 1.9 Myr of Brandenburg et al. (2013), or the order of magnitude estimate of 3-7 Myr of (Chamandy & Shukurov 2020).Whether τ differs between magnetic and kinetic α, or overall, a factor of 1/5 would not qualitatively alter the dynamics.
Considering the turbulent α displayed for Model L2 in Figure 6(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 Gent et al.
α contribution in Appendix A.4.This shows that the evolution of α is dominated by the contribution of kinetic α 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 α-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 α in the dynamo active warm gas near the midplane, but also with growth of magnetic α M of opposite sign to α K in this region (see Appendix A.4).The αquenching near the midplane, due to reducing α K and growing magnetic α of opposite sign saturates the LSD.
Comparing Models L2 and H2 in Figures 6(b) and (c) during corresponding 700 Myr intervals, their α topologies are qualitatively similar.At high resolution, α 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, which are evident as smallscale turbulent structures that appear first at the midplane and propagate quickly towards the halo, while large-scale, same-sign helicity populates the inflows, appearing at the boundary and propagating more slowly towards the midplane.Thus, the small-scale helicity at sufficiently high Rm (Mitra et al. 2010) is sporadically removed from the disc allowing the LSD to grow and for sufficiently high helicity fluxes to saturate even at levels above equipartition (Chamandy et al. 2014).

Effect of clustering
Figure 7(a) plots total magnetic energy density evolution for all models.It 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 turbulent magnetic field, twice that of Model L2, and 4/3 faster after the SSD has saturated.
The turbulent kinetic energy density in Figure 7(b) is similar with and without clustering, while the turbulent velocity in Figure 7(c) is much higher due to clustering.When OB clustering is included, the highest turbulent velocities are in hot 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.

Numerical convergence
We next study the numerical convergence of the combined dynamos.We first consider the convergence behavior of the SSD. Figure 7(a) shows the extreme sensitivity of the growth rate of the SSD to numerical resolution (as reported Gent et al. 2021Gent et al. , 2023)).In the high-resolution Model H2 an eigenmode is established at 7 Myr, saturating 26 Myr later.By contrast, the lowresolution models L2 and L2-B0 take up to 40 Myr to establish the eigenmode, then 900 and 500 Myr, respectively, to saturate.
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.In Gent et al. (2023) the high and ultra high resolution SSD established an eigenmode at 40 Myr and then saturated within 17 Myr.More realistic FUV heat- ing allows the warm gas to cool, which increases the fractional volume of hot gas.It is likely therefore, given higher SN rates that the duration of the SSD would reduce further to a few megayears or even shorter.
The kinematic growth rate γ B differs by a factor of about 50 between model H2 and L2.Theory predicts γ B ∝ Re with some additional weak dependence on ln(Rm) (Kleeorin et al. 1986;Subramanian 1997) and inverse dependence on M (Moss & Shukurov 1996;Federrath 2016).
Without explicit diffusivities, using a second order solver, Teyssier et al. (2006) find that their effective Re scales inversely with the square of grid resolution.For high resolution Gent et al. (2021) show the Lagrangian diffusivities used here dominate, while at low resolution numerical diffusivities dominate.Their direct comparison to second order solutions of Balsara et al. (2004) demonstrates that numerical diffusivity reduces substantially at sixth order.The ratio of numerical Re between H2 and L2 could therefore be as high as 4 6 = 4096, which would more than account for a factor of 50 in γ B .
ISM turbulence is multi-scale due to random explosions and variation in remnant location and size.It is strongly affected by nonadiabatic cooling.To test the applicability of theoretical predictions a dedicated parameter study of SSD at resolution δx ≤ 1 pc is required.
In order to study the numerical convergence of the LSD solution after saturation of the SSD, the large variation in time of SSD saturation needs to be taken into account.The saturation energy of the magnetic field is anyway well converged.For the period that the LSD dominates field growth, Model H2 is shifted by 970 Myr to align the LSD solutions.(Similarly, shifting the dynamo solution for Model L2 illustrated in Figure 7 The turbulent velocity of model H2 (Fig. 7(c)) does slightly exceed that of L2, but the main driver of the high rate of SSD at high resolution is that the filamentary dense structures occupy a lower fractional volume, reducing net cooling and increasing the fractional volume of hot diffuse turbulent structures.The total turbulent kinetic energy density (Fig. 7(b)) is well converged at the resolution of these models, as is the saturation strength of the SSD.Volume averaged energy density (a) of total magnetic energy ⟨B 2 ⟩(8π) −1 and horizontally averaged mean-field energy ⟨B⟩ xy 2 (8π) −1 , normalised by time-averaged turbulent kinetic energy density e ′ K 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 Gyr < t < 1.5 Gyr.(b) Ratio of horizontallyaveraged mean-field energy to turbulent magnetic field energy for same models and intervals as (a).
Figure 8(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.
By shifting the solution displayed in Figure 7 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 8 yields γ = 4 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.2).
Figure 8(b) shows the ratios obtained from the values in panel (a), which are ≲ 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;Beck et al. 2019).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.
In Appendix A.3 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.

Hydrodynamics
In Figure 9(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 7(b), is actually weaker than for Model L2.By contrast the turbulent velocity in Figure 7(c) has a higher magnitude (Fig. 9 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.2 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.

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 ℓ = 50 pc to separate the mean and random magnetic fields for ⟨B⟩ ℓ in Figure 4.
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 ≥ 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.

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.Asymmetry arises in SN-driven turbulence due to sporadic random clustering of explosions on one side of the midplane.This is enforced by OB-clustering, which ensures an extended period when subsequent explosions will reoccur on that side.This is a physical behaviour, but across the span of the disc would tend to average out, leaving a pocklike imprint over its gas surface.
The expanded fractional volume of hot gas, in this particular instance above the midplane for Model L2-cl, is beneficial to the SSD (Gent et al. 2023).However, it does not favour the LSD, as can be seen by comparing the spread of the mean magnetic field to large |z| in Figure 3(c), while it is confined to the midplane with the Gent et al.
warm gas for Model L2-cl (see Appendix A.2). Figure 7(a), and comparison of their mean fields in Appendix A.2, show that after saturation of the SSD the mean field actually grows more efficiently in the thinner, more dense disc of the OB-clustered model L2-cl than L2 without.Gent et al. (2023) show sporadic rapid growth in the diffuse hot gas during the kinematic SSD, but upon saturating the turbulent field becomes increasingly correlated with the gas density.Both absence of mean field in the hot gas and increased LSD growth rate associated with a more uniform warm phase support the hypothesis that the mean-field dynamo grows mainly in the warm rather than the hot gas, but could also indicate correlation of the mean field with mean gas density.Evirgen et al. (2017) demonstrate that the turbulent field is proportionate to gas density in both warm and hot phases, but the mean field preferentially resides in the warm gas.

Quenching of the α-effect
We have used kinetic and current helicities as proxies for measuring the collective inductive action of SN turbulence, called the α-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.Singular value decomposition (Racine et al. 2011;Bendre et al. 2020) may also be applied retrospectively, although the presence of SSD in our models needs to be taken into account.
As Model L2 approaches equipartition in Figure 6(a) we observe α-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 α-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 α in the northern latitudes, which is a magnetic α-effect, is a signature of this behaviour as identified by Machida et al. (2013); Tharakkal et al. (2023), andQazi 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.

CONCLUSIONS
We have verified that the large-scale dynamo (LSD) in disc galaxies is robust in the presence of SSD at high magnetic Reynolds number Rm. (This varies strongly across the domain, reaching values as high as ∼ 10 5 in hot regions (see Gent et al. 2023).)The LSD 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. 2008b;Gent et al. 2013b;Bendre 2018).
Chamandy & Singh (2017, see their Figure 2) seek to relate the growth of the mean galactic field to various ratios of small-scale magnetic to turbulent kinetic energy density.Our normalised growth rate for model L2 would be ∼ 3.5 vs 4.5 (Gent et al. 2013b), which are slightly higher values than the upper estimate from Chamandy & Singh (2017) of ∼ 2.5 for Ω = 2Ω Sn .The increased mean field growth with reduced ratios is however consistent with their model.
In contrast, the SSD appears quite insensitive to the presence of an LSD.Given that the peak in the power spectrum of Model L2 remains at a scale below the range of forcing scales provided by SNe, we can be confident the magnetic energy growth is attibutable to SSD rather than tangling of the mean-field (Gent et al. 2021, their Figure 1).Tangling injects peak energy at the scale of the forcing.The presence of largescale rotation and shear does, however, inject more vorticity, increasing the growth rate of the SSD (Singh et al. 2017;Gent et al. 2023).The increased growth rates apply primarily to larger scales.
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 ≲ 17k 1 , before the SSD saturates the selfsimilar growth ceases and growth at the largest scales slows down.After the LSD becomes dominant, the larger scales overtake the energy at k ≃ 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.
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 ≃ 17k 1 , or 30 pc scale, independent of resolution (or effective Rm).This remains the case in the nonlinear phase of the pure SSD, even though the portion of total magnetic energy at scales larger than 30 pc increases at late times.When the LSD is present, the self similar growth of the magnetic energy spectra continues in the same manner for all k ≳ 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 ℓ ≲ 30 pc.
We have shown that the SSD in the ISM can be expected to saturate at a few percent of equipartition within a few megayears, depending on the SN rate and disc mass.During the kinematic stage of the SSD the mean field is also amplified at slightly above the self-similar growth rate of the SSD up to ∼ 10 −5 -10 −4 of equipartition.The LSD in the presence of high Rm SSD is slower relative to LSD at lower Rm, although its dependence across a range of turbulent Rm remains to be explored.The mean field in spiral galaxies is likely to require more than 2 Gyr to saturate, depending on SN rate, rate of shear and scale height, which also remains to be explored, and its saturation strength depends on the sum of turbulent kinetic energy and energy from the galactic shear.The properties of the saturated magnetic field, including its mean-field fraction, topology, and sensitivity to parameters are pending further investigation.
The total magnetic field already approaches equipartition with the turbulent kinetic energy within 1 Gyr.The turbulent magnetic field continues to grow in the presence of the LSD at relative rates below that of the LSD, but can increment magnetic energy orders of magnitude larger than the mean field itself.The alpha effect appears to support the LSD by transporting helicity on small scales away from the midplane, and the magnetic alpha, which can acquire magnitude comparable with the kinetic alpha within 400 Myr, initially has the same sign and supports the LSD.In the later stages, as the mean field continues to grow, the magnetic alpha reverses sign and quenches the LSD.A more robust analysis of the turbulent transport processes is required.

Gent et al.
We thank the anonymous referee for a thorough, constructive, and intelligent review of the work, which has inspired substantial improvements to the quality of the presentation.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 Nord-Forsk.M-M.M.L. was partly supported by US NSF grants AST18-15461 and AST23-07950, and thanks the Inst.für Theoretische Astrophysik der Uni.Heidelberg for hospitality.Samples of the spectral properties of the magnetic energy and kinetic energy of both pure SSD models are displayed in Figure 10(a) and (b).Spectra have been smoothed with a Gaussian kernel of length 5 Myr and 1.5 k 1 in time and wavelength, with smoothing both forward and backward in time.The magnetic energy spectra, which grow many orders of magnitude over time, are normalised by maximum energy density to ease comparison between times.The equivalent spectra for Model L2 are plotted in panels (c) and (d).In all models the peak wavenumber of the magnetic energy during the kinematic stage of the SSD is at 17k 1 , with k 1 = 2πL x −1 ≃ 12.3 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.After saturation of the SSD and transition to the LSD, the magnetic energy shifts to larger scales (k → 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 .Later, as the mean field becomes significant beyond 3 Gyr, energy is lost at large scales with the spectrum at k ≲ 5k 1 becoming more shallow.
To examine the evolution of the magnetic field energy spectra over time, we plot the power in selected wavenumbers of the pure SSD models in Figure 11(a) and (b).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 quite insensitive to the large-scale rotation and shear, other than to increase the growth rate and saturation energy of the SSD.
In Figure 11(c) the evolution of the energy in wavenumber bins with k ≥ 17k 1 is insensitive to the LSD and retains the same ranking as for the SSD in panels (a) and (b).In contrast, wavenumber bins with k < 17k 1 , which had lower energies during the SSD, first detach from the eigenmode shortly before 1 Gyr, slowing down.However, when the SSD saturates shortly after 1 Gyr these wavenumbers continue to gain energy, eventually becoming dominant over the larger wavenumber bins with k > 17k 1 .During this transition up to about 2 Gyr the ranking within k < 17k 1 becomes inverted, while at smaller scales the ranking from the kinematic stage is conserved.Thus, the LSD can be identified as acting primarily on scales k ≲ 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 even k = 60k 1 exceeding 10 −16 erg cm −3 .
Determining how much the pivot point of k = 17k 1 , corresponding to a scale of 30 pc, depends on parameters such as supernova rate, galactic shear, the mass of the disc, or the Reynolds numbers would require a parameter sweep beyond the scope of this study.We hypothesise that this qualitative behaviour will provide a useful condition to identify the separation be-     tween mean and fluctuation scales, given that spectra from simulations typically do not exhibit minima that could identify separation between mean and fluctuation scales.
The inversion of rank in energy for k < 17k 1 mirrors contrasting growth in Figure 4    earliest, but it subsequently grows more rapidly than ⟨e B ⟩ ℓ .
Figure 12 displays evolving spectral structure for Models L0-B0, L2-B0 and L2 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 ≃ k 5/4 for Model L2 between 0.08 and 0.78 Gyr in Figure 12(c), and ≃ k 3/2 for the equivalent stages in (a) and (b).The steepening resistive tail in Figure 12(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).Given the mean field in the form of horizontal averages is removed from Model L2-B0 and yet the shear adds energy preferentially to scales at k < 17k 1 , this indicates the mean field is not fully expressed in this form.

A.2. Vertical density structure and effect of clustering
In Figure 13    Figure 14 shows SSD and LSD growth rates for Model L2-cl, for comparison with the growth rates in Figure 4 for Model L2.The latter 300 Myr interval when shifted 0.42 Myr corresponds to the LSD stage fitted for Model L2 in Figure 8.

A.3. Convergence
A further indication of the level of convergence of the solution (see Sect. 4.1) is a comparison of the topology of the mean magnetic field produced between Models H2 and L2.In Figure 15 the time-latitude diagrams for the horizontallyaveraged 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.
In Figure 16 from Model H2 magnetic (a) and kinetic (b) energy spectra are shown for comparison to the equivalent spectra from Model L2 in Figure 10.The peak wavenumber for the magnetic spectra during the kinematic stage of the SSD are at k ≃ 17k 1 for L2 but vary 16k 1 ≤ k ≤ 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 15 these differences do not appear to affect the horizontallyaveraged structure of the large-scale field.Figure 16(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 11(b).Saturation of the SSD occurs at about 10 −17 erg cm −3 for all k ≳ 17k 1 in both models, after which the hierarchy between these scales is conserved in the nonlinear stage of the SSD.For k ≲ 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 δx ≲ 1 pc.The magnetic energy spectra plotted in Figure 16(a) show that the peak energy during the SSD varies around slightly smaller scale, 17k 1 ≲ k ≲ 27k 1 rather than more consistently near 17k 1 .Comparing the kinetic energy spectra between Figure 10(b) and Figure 16(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).

A.4. α coefficients
In Figures 17 and 18 for Models L2 and H2, respectively, we show the separate contributions of kinetic α K and magnetic α M .For Model L2 the systemic sign of the α 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 α M is negligible almost to 2 Gyr.Subsequently it has large-scale structure away from the midplane, typically with same sign to the kinetic α 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 α 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 α K , thus enhancing the LSD.This is remarkably consis-tent with the analytic prediction by Gopalakrishnan & Subramanian (2023) of the timescale over which α M would become dynamically effective.At this stage it has the same sign as α K and may enhance the LSD.For Model L2 this is delayed to almost 2 Gyr, and soon exhibits sign opposite to α K .In the outflows we see that α K exhibits fluctuating sign on short time and spatial scales, which may also have a positive effect on the LSD by supporting the removal of small-scale kinetic helicity of opposite sign to the mean field.

B. CODE AND PARAMETER EVOLUTION DURING RUN TIME
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 Gent et al.   5)  5th order rate of strain tensor ϵ ijk Levi-Cevita symbol, i, j, k ∈ (1, 2, 3) +(−)1 if permutation of 1,2,3 is even(odd) or 0 if an index repeats.
Table 3. Tracking parameter changes over duration of simulations: For the models listed the parameters applied commencing at time t are shown: h adj indicates the scale height of the SN vertical distribution varies with the scale height of the gas; Γ UV is the photoelectric heating factor relative to the the Habing rate Γ 0 .T max is the maximum temperature allowed within r 0 , the SN radial scale as explained in Section 2.6 of the SN explosion, and T ratio × T max is the maximum allowed within 2.25r 0 of the explosion.n cool indicates whether cooling mass is used to reduce the temperature, otherwise the SN site is rejected and an alternative selected.n 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 δx = 4 pc n ratio is the maximum ratio between maximum and minimum gas density within r 0 of an explosion, to avoid excessive viscous forces.are not particular to SN-driven turbulence, but nevertheless may affect code performance.
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 ν 6 , χ 6 and η 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 10 and Model H2, in Figure 16(a).In Model H2 these coefficients were increased to 1.2 × 10 −15 after 652 Myr.In this study, 5 × 10 −12 was applied throughout for δx = 4 pc, except as listed in Table 3 at late times for Model L2.
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 σ = σSn , inducing a higher SSD growth rate.Beyond 0.3 Gyr continued with σ = 0.5 σSn .The period 0-0.365Gyr 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.

Figure 1 .
Figure1.Structure of high resolution (δx = 1 pc) Model H2 with properties given in Table1.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 a 52 second animation of this Figure, demonstrating the rapid transition to turbulence, the saturation of the turbulent magnetic energy within about 30 Myr and the much slower organisation into large scales of the azimuthal component of the magnetic field up to 652 Myr.The breathing and sloshing modes of the disc are also apparent.The vertical striping evident in the magnetic energy occurs from inflows through the boundary which has a vertical magnetic field condition.A high resolution version is included in the data repository (Gent et al. 2023).
Figure 2.Volume-averaged magnetic energy density for models L2-B0 (blue) with Ω = 2Ω Sn and L0-B0 (green) with Ω = 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 e ′ K .The horizontal dotted line indicates 2% equipartition.
Figure 3. Time-latitude diagrams of the horizontally averaged magnetic field (a) ⟨B x ⟩ xy and (b) ⟨B y ⟩ 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.

Figure 6 .
Figure 6.Time-latitude diagrams of horizontally averaged α-effect from Equation (13) for Model L2 (a) throughout and (b) during a restricted interval comparable with (c) Model H2.A linear colourscale is applied in the range −5 km s −1 < α < 5 km s −1 and symmetric logarithmic colour-scale for larger magnitudes.
Figure 7.(a) Volume-averaged magnetic energy density for all models listed in the legend, normalised by the time-averaged turbulent kinetic energy density e ′ K .(b) Volume-averaged turbulent kinetic energy density for the LSD models, and (c) volume-averaged root-mean-square turbulent velocity u ′ .

Figure 9 .
Figure 9. Volume-averaged total kinetic energy density e ′ K (a) for all models as listed in the legend and (b) volume-averaged root-mean-square total velocity u for the same models.
Spectral comparison of SSD and LSD

Figure 10 .
Figure10.Magnetic energy spectra (a) in pure SSD Models L0-B0 and L2-B0 and (c) the fiducial Model L2, each normalised by its maximum power.The spectra for Model L0-B0 are shifted 10 −2 to differentiate from Model L2-B0.Kinetic energy spectra (b) and (d) with matching models, line styles, times (and axis shift) to panels (a) and (c), respectively.Listed wavenumbers k, in multiples of k 1 = 2πL x −1 , locate the maximal magnetic power at each time.Spectra have been smoothed with a Gaussian kernel of length 5 Myr about time t and 1.5 k 1 about wavenumber k.

Figure 11 .
Figure 11.Time evolution of magnetic energy for SSD models (a) L0-B0, (b) L2-B0 and (c) L2 at sampled scales of k, as noted in the legends in multiples of k 1 = 2πL x −1 .
between ⟨e B ⟩ xy , consisting of only the longest length scale, and ⟨e B ⟩ ℓ , comprising a range of intermediate scales.The growth of ⟨e B ⟩ xy slows down Gent et al.

Figure 12 .
Figure 12.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.Spectra have been smoothed with a Gaussian kernel of length 5 Myr about time t and 1.5 k 1 about wavenumber k.
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 3 for Model L2.

Figure 13 .
Figure 13.Time-latitude diagrams of horizontally averaged gas number density for (a) Model L2 and (b) L2-cl and (c) ⟨B y ⟩ xy for Model L2-cl.

Figure 15 .
Figure 15.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.
Figure16.(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πL x −1 .

Table 1 .
Models included in analysis of LSD & SSD interactions.Model names include H (L) for high (low) resolution, then the rotation velocity in Solar neighborhood units.OB clustering is labelled "cl"."B0" denotes zero mean magnetic field.
Balbus 2003)orces are included to first order, while centrifugal forces are assumed to balance radial gravity, and thus neglected (see Eq. [76] ofBalbus 2003).