AppendixAppendix References

# Growth or decay: universality of the turbulent dynamo saturation

James R. Beattie${}^{\href https://orcid.org/0000-0001-9199-7771\,1,2}$,
Christoph Federrath${}^{\href https://orcid.org/0000-0002-0706-2306\,1,3}$,
Neco Kriel${}^{\href https://orcid.org/0000-0002-3558-3926\,1}$, Philip
Mocz${}^{\href https://orcid.org/0000-0001-6631-2566\,4,5}$, and Amit
Seta${}^{\href https://orcid.org/0000-0001-9708-0286\,1}$
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
2Department of Astronomy and Astrophysics, University of California, Santa
Cruz, 1156 High Street, Santa Cruz, CA 96054
3Australian Research Council Centre of Excellence in All Sky Astrophysics
(ASTRO3D), Canberra, ACT 2611, Australia
4Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane,
Princeton, NJ 08544, USA
5Lawrence Livermore National Laboratory, 7000 East Ave, Livermore, CA, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The small-scale turbulent dynamo (SSD) is likely to be responsible for the
magnetisation of the interstellar medium (ISM) that we observe in the Universe
today. The SSD efficiently converts kinetic energy $E_{\rm kin}$ into magnetic
energy $E_{\rm mag}$, and is often used to explain how an initially weak
magnetic field with $E_{\rm mag}\ll E_{\rm kin}$ is amplified, and then
maintained at a level $E_{\rm mag}\lesssim E_{\rm kin}$. Usually, this process
is studied by initialising a weak seed magnetic field and letting the
turbulence grow it to saturation. However, in this study, using three-
dimensional, non-ideal magnetohydrodynamical turbulence simulations, we show
that the same saturated state can also be achieved if initially $E_{\rm
mag}\gg E_{\rm kin}$ or $E_{\rm mag}\sim E_{\rm kin}$. This is realised
through a two-stage exponential decay (1. a slow backreaction that converts
$E_{\rm mag}$ into $E_{\rm kin}$, and 2. Ohmic dissipation concentrated in
anisotropic current sheets) into the saturated state, for which we provide an
analytical model. This means that even if there are temporary local
enhancements of $E_{\rm mag}$ in the interstellar medium, such that $E_{\rm
mag}>E_{\rm kin}$, e.g., through amplifications such as compressions, over a
long enough time the field will decay into the saturated state set by the SSD,
which is determined by the turbulence and magnetic dissipation. However, we
also provide analytical models for the decay times and utilise wait-time
statistics from compressive supernova events to show that if the magnetic
field is enhanced above the saturated state, it will not have enough time to
decay the field before the next supernova event. Hence, unless there exists a
mechanism for destroying magnetic fields that is not in our non-ideal MHD
models, the amplitudes of interstellar magnetic fields may also be a function
of a complex history of compressions, rather than only the action of a SSD.

###### keywords:

MHD – turbulence – ISM: kinematics and dynamics – ISM: magnetic fields –
dynamo

††pubyear: 2022††pagerange: Growth or decay: universality of the turbulent
dynamo saturation–14

##  1 Introduction

The present day Universe is magnetised, thus beckoning the question: how did
it become so? One answer is the small-scale turbulent dynamo (SSD; or
fluctuation dynamo), which is a mechanism for transforming turbulent kinetic
energy into turbulent magnetic energy until both are statistically stationary
and approximately in energy equipartition – the so-called saturation of the
SSD. In the interstellar medium (ISM) of galaxies, saturation probably
occurred at redshifts $z\approx 25-8$, depending upon the nature of density
fluctuations (Xu & Lazarian, 2016; McKee et al., 2020), consequentially making
the study of present-day ISM magnetic fields the study of the saturated stage
of the SSD. However, understanding both the physics, statistics, and
consequently constructing a predictive model for the saturation of the SSD
remains an active problem in the dynamo community.

In a Markovian fashion, the saturated state of the magnetic field that
develops from the small-scale turbulent dynamo forgets the field that seeded
it (Seta & Federrath, 2020). These initial, primordial fields may have been
incredibly weak (maybe as weak as $10^{-16}\,\rm{G}$ in the gas of the
intergalactic medium), and perhaps were generated from a battery process
(e.g., Biermann, 1950) or through spontaneous magnetic field creation during
the electroweak symmetry breaking epoch in the primordial Universe (e.g.,
Brandenburg & Subramanian, 2005; Subramanian, 2016, 2019). Such a weak,
fluctuating seed field $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$ can be
amplified exponentially fast in time, viz.,
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}\propto\exp\left\\{\gamma
t\right\\}$ by turbulent motions in the plasma, where
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$ is
the volume ($L^{3}\equiv\operatorname{\mathcal{V}}$) integral magnetic energy
and $\gamma\sim v_{\nu}/\ell_{\nu}$ the growth rate that goes with the
reciprocal dynamical time of the turbulence at the viscous scale $\ell_{\nu}$
(e.g., McKee et al., 2020). This is termed the kinematic, exponential growth,
or linear induction equation stage of the turbulent dynamo (e.g., Brandenburg
& Subramanian, 2005), and lasts until
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
becomes large enough to cause a backreaction on the momentum transport through
$(\nabla\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})\times\mathrm{{\mn@boldsymbol{\mathit{b}}}}$,
e.g., when
$\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})\sim\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})$.
In this stage, the induction equation becomes strongly nonlinear because
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ satisfies the momentum equation with
non-negligible
$(\nabla\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})\times\mathrm{{\mn@boldsymbol{\mathit{b}}}}$,
which then contributes to the induction equation and the overall time-
evolution of $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$. This stage is termed the
linear growth (because
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}\propto
t$) or the nonlinear induction equation stage of the dynamo (e.g.,
Schekochihin et al., 2002a; Maron et al., 2004; Cho et al., 2009; Xu &
Lazarian, 2016). Finally, after a sufficient amount of time
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
saturates, such that
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}\sim\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$111Note
that previous works suggested
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}\sim\text{Re}^{-1/2}\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
(Batchelor, 1950). This would result in an exceptionally weak magnetic field
for astrophysical media, such as molecular clouds in the interstellar medium,
which are characterised by Re as high as $\text{Re}\sim 10^{9}$ (Krumholz,
2014). However SSD simulations and now laboratory experiments (Tzeferacos et
al., 2018; Liao et al., 2019) show higher levels of saturation (see §3.2.1 in
Rincon, 2019, for more details).. The integral quantities of the saturation
depend upon the compressibility of the plasma (e.g., the turbulent Mach number
$\mathcal{M}=\left\langle{v^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/c_{s}$,
where $c_{s}$ is the sound speed; see Haugen et al. 2004; Federrath et al.
2011; Seta et al. 2020; Seta & Federrath 2021), the nature of the turbulent
driving source (Federrath et al., 2011; Achikanath Chirakkara et al., 2021),
and the diffusion timescales for the magnetic and velocity fields – the
magnetic Prandtl number, Pm (for finite Pm) (e.g., Schober et al., 2012).

However,
$\left\langle{\mathit{b}^{2}}\right\rangle_{\ell^{3}}\sim\left\langle{\mathit{v}^{2}}\right\rangle_{\ell^{3}}$
need not be true for all $\ell$ in the magnetised plasma, and the saturation
of the SSD is a scale-dependent (saturation looks different on different wave
modes) phenomenon (e.g., Schekochihin et al., 2002a; Maron et al., 2004;
Schober et al., 2015, 2012). Here we summarise the phenomenology championed
most recently by Galishnikova et al. (2022). As the dynamo approaches
saturation, the magnetic field fluctuations on each scale $b_{\ell}^{2}$ reach
energy equipartition with the turbulence on that scale $v_{\ell}^{2}\sim
b_{\ell}^{2}$, successively moving the energy equipartition scale (equivalent
to the shearing scale in Maron et al. 2004) $\ell_{\rm eq}\sim k_{\rm
eq}^{-1}$ from small, viscous-dominated scales $\ell_{\rm eq}\sim\ell_{\nu}$
(where $t_{\ell}=\ell/v_{\ell}$ is short), to the larger and slower eddies on
$\ell\gg\ell_{\nu}$ until a maximal stretching rate is achieved, given by
$t_{\text{max}}^{-1}\sim\big{(}\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/L\big{)}\big{(}\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}\big{)}^{2}$
on $\ell_{\rm eq}$ (Galishnikova et al., 2022). At this point, only
$\ell>\ell_{\rm eq}$ are not suppressed by the magnetic tension, and hence the
final value of
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
is sensitive to where this scale is. Clearly, for supersonic dynamo
experiments, where
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}$
is reduced compared to the simulated subsonic counterparts (see for example,
Federrath et al., 2011; Seta & Federrath, 2021) and $\ell_{\rm eq}$ must be on
smaller scales, which allows for hydrodynamic motions to dominate over a
larger range of $k$. The details of the exact saturation mechanism are still,
however, an active area of research (some ideas and discussion in Rincon,
2019; Seta et al., 2020; Seta & Federrath, 2021).

However, what happens if
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}\gg
1$ to begin with, potentially from phenomena that enhance the magnetic field
beyond the dynamo saturation level? This situation may arise from
amplifications, such as large-scale compressions of the ISM from galaxy-galaxy
interactions and radial flows boosting the magnetic field (Steinwandel et al.,
2020), or even small-scale events such as supernova driven shocks (e.g., Korpi
et al., 1999; Lu et al., 2020; Chevance et al., 2022), or stellar feedback
(e.g., Lancaster et al., 2021; Menon et al., 2022). Such compressions may
enhance a magnetic field on scales where flux-freezing holds (valid above the
scales of ion-neutral damping, $\mathcal{O}(10^{-3}-10^{-2}\rm{pc})$ for
Alfvén modes; Krumholz et al. 2020), such that $B\propto\rho^{\alpha}$, where
$\rho$ is the gas density and $\alpha$ is the enhancement factor that depends
upon the detailed geometry of the compression (Tritsis & Tassis, 2016; Mocz &
Burkhart, 2018; Beattie et al., 2021). In these circumstances, it is not clear
if different, strong magnetic fields (accompanied with a different magnetic
morphology) intrinsically change the $t\rightarrow\infty$ behaviour of the
turbulent magnetic field, or more specifically, if the properties of the
saturated magnetic field are functions of the detailed workings of the
kinematic and non-linear stages in the SSD or indifferent to them.
Understanding these properties and the length and timescales involved in this
process, is the key motivation for this study.

Table 1: Main simulation parameters and derived quantities.

Sim. ID | Re | Rm | $\nu t_{0}/\ell_{0}^{2}$ | $\eta t_{0}/\ell_{0}^{2}$ | $\mathcal{M}$ | $b_{\rm init}$ | $E_{\rm mag,0}/E_{\rm kin}$ | $(E_{\rm mag}/E_{\rm kin})_{\text{sat}}$ | $(\operatorname{\mathcal{M}_{\text{A}}})_{\text{sat}}$ | driven
---|---|---|---|---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
$\text{Pm}=1$
weakPm1 | 500 | 500 | $2\times 10^{-3}$ | $2\times 10^{-3}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{-10}$ | 0.09 $\pm$ 0.03 | 3.46 $\pm$ 0.60 | ✓
satPm1 | 500 | 500 | $2\times 10^{-3}$ | $2\times 10^{-3}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{-2}$ | 0.09 $\pm$ 0.03 | 3.45 $\pm$ 0.64 | ✓
strongPm1 | 500 | 500 | $2\times 10^{-3}$ | $2\times 10^{-3}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{2}$ | 0.09 $\pm$ 0.03 | 3.44 $\pm$ 0.59 | ✓
initbPm1 | 500 | 500 | $2\times 10^{-3}$ | $2\times 10^{-3}$ | 0.5 | $7\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 9$ | $10^{2}$ | 0.09 $\pm$ 0.03 | 3.45 $\pm$ 0.62 | ✓
decayPm1 | 500 | 500 | $2\times 10^{-3}$ | $2\times 10^{-3}$ | 0.0 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{2}$ | - | - | ✗
$\text{Pm}=2$
weakPm2 | 500 | 1000 | $2\times 10^{-3}$ | $1\times 10^{-3}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{-10}$ | 0.23 $\pm$ 0.05 | 2.16 $\pm$ 0.26 | ✓
satPm2 | 500 | 1000 | $2\times 10^{-3}$ | $1\times 10^{-3}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{-2}$ | 0.23 $\pm$ 0.06 | 2.14 $\pm$ 0.27 | ✓
strongPm2 | 500 | 1000 | $2\times 10^{-3}$ | $1\times 10^{-3}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{2}$ | 0.24 $\pm$ 0.09 | 2.12 $\pm$ 0.29 | ✓
initbPm2 | 500 | 1000 | $2\times 10^{-3}$ | $1\times 10^{-3}$ | 0.5 | $7\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 9$ | $10^{2}$ | 0.22 $\pm$ 0.06 | 2.17 $\pm$ 0.28 | ✓
$\text{Pm}=4$
weakPm4 | 500 | 2000 | $2\times 10^{-3}$ | $5\times 10^{-4}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{-10}$ | 0.35 $\pm$ 0.08 | 1.72 $\pm$ 0.19 | ✓
satPm4 | 500 | 2000 | $2\times 10^{-3}$ | $5\times 10^{-4}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{-2}$ | 0.33 $\pm$ 0.08 | 1.74 $\pm$ 0.20 | ✓
strongPm4 | 500 | 2000 | $2\times 10^{-3}$ | $5\times 10^{-4}$ | 0.5 | $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ | $10^{2}$ | 0.33 $\pm$ 0.07 | 1.76 $\pm$ 0.19 | ✓
initbPm4 | 500 | 2000 | $2\times 10^{-3}$ | $5\times 10^{-4}$ | 0.5 | $7\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 9$ | $10^{2}$ | 0.34 $\pm$ 0.08 | 1.74 $\pm$ 0.19 | ✓

  * Notes. Column (1): the unique simulation ID. Column (2): the Reynolds number of the plasma, Equation 5, in the saturated phase of the dynamo. Column (3) the same as column (2) but for the magnetic Reynolds number, Equation 6. Column (4): the coefficient for the kinematic viscosity (see viscous stress tensor in Equation 2) expressed in units of correlation time of the driving $t_{0}$ and driving scale $\ell_{0}$. Column (5): the same as column (4) but for the Ohmic resistivity (Equation 3). Column (6): the turbulent Mach number, $\mathcal{M}=\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/c_{s}$, where $c_{S}$ is the sound speed, in the saturated phase of the dynamo. Column (7): the domain of the parabola for the initial magnetic field. Column (8): the initial magnetic and kinetic energy ratio. Column (9): the magnetic and kinetic energy ratio in the saturated state of the dynamo. Column (10): the same column (9) but for the Alfvén Mach number. Column (11): an indication of whether the turbulence is continuously driven ✓ or initialised and left to decay ✗.

In this study we establish that the dynamo saturation and morphology of the
saturated magnetic field is universal, and completely independent of the first
two stages of the small-scale turbulent dynamo, as hypothesised by Maron et
al. (2004), but not explored systematically. Specifically, the
$t\rightarrow\infty$ integral quantities, energy spectra, and characteristic
scales of the magnetic field do not depend upon the initial conditions of the
magnetic field. This means that the kinematic and non-linear growth stages of
the dynamo are not a set of initial conditions that uniquely define the
saturated state of the small-scale turbulent dynamo and that the properties of
the saturated field in the ISM ought to be completely determined by the
turbulence222Assuming that the magnetic dissipative mechanisms (based on the
microphysics of the gas) are somewhat universal on small scales across the
ISM.. We also establish a mechanism akin to an “anti-dynamo", where, if
$\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}\gg
1$, the magnetic field will undergo a two-stage decay. The first decay stage
is from the magnetic energy being transferred into kinetic energy on diffusion
timescales of the magnetic field, and the second decay stage from Ohmic
diffusion in anisotropic current sheets. Furthermore, our analysis suggests
that the rate at which compressive events in the ISM may grow magnetic fields
beyond a dynamo saturation is much shorter than the timescale for the
amplified magnetic field to decay back into the saturated state. Hence, the
amplitudes of magnetic fields in the ISM may be a function of a complex
history of amplification processes, rather than only the action of a SSD,
unless there exists a mechanism (perhaps fast reconnection; Lazarian &
Vishniac 1999) for destroying magnetic fields that is not captured in our
turbulence box simulations – a well-known problem that spans many
astrophysical applications (see e.g., the excess magnetic flux problem in the
star formation process; Braithwaite 2012; Krumholz & Federrath 2019;
Hennebelle & Inutsuka 2019).

Our study is organised as follows. In Section 2 we discuss the numerical
simulations that we use to probe how the saturation of the SSD responds to
changing the initial ratio between the magnetic and kinetic energy and the
scales that the magnetic energy exist upon. In Section 3 we report upon the
ratio of integral energies. Next in Section 4 we explore the saturation on a
scale-by-scale manner, and then in Section 5 as a function of characteristic
scales of the magnetic and current density field. In Section 6 we provide
analytical models for the magnetic field decay into the saturated state of the
dynamo and the timescales involved in the decay, and discuss the implications
that this has on situations where the magnetic field in the ISM may grow
through reversible processes, such as compressions via supernovae, beyond the
dynamo saturation, followed by decay back to the saturation. Finally, in
Section 7 we summarise the key results of our study.

##  2 Numerical simulations

We use a modified version of the flash code (Fryxell et al., 2000; Dubey et
al., 2008), utilising a second-order conservative MUSCL-Hancock 5-wave
approximate Riemann scheme (Bouchut et al., 2010; Waagan et al., 2011;
Federrath et al., 2021) to solve the 3D, non-ideal, isothermal, compressible
MHD equations with a stochastic acceleration field acting to drive the
turbulence,

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}})$ $\displaystyle=0,$
(1)
$\displaystyle\rho\frac{\partial\mathrm{{\mn@boldsymbol{\mathit{v}}}}}{\partial
t}-\nabla\cdot\left[\frac{1}{4\pi}\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}-\rho\mathrm{{\mn<EMAIL_ADDRESS>50.00008pt\ignorespaces$
$\displaystyle\left.-\left(c_{s}^{2}\rho+\frac{b^{2}}{8\pi}\right)\mathbb{I}-2\nu\rho\mathcal{S}\right]$
$\displaystyle=\rho\mathrm{{\mn@boldsymbol{\mathit{f}}}},$ (2)
$\displaystyle\frac{\partial\mathrm{{\mn@boldsymbol{\mathit{b}}}}}{\partial
t}-\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})-\eta\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}}$
$\displaystyle=0,$ (3)
$\displaystyle\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{b}}}}$
$\displaystyle=0,$ (4)

where $\otimes$ is the tensor product333For vector fields
$\mathrm{{\mn@boldsymbol{\mathit{g}}}},\,\mathrm{{\mn@boldsymbol{\mathit{h}}}}\in\mathbb{R}^{n}$,
$\mathrm{{\mn@boldsymbol{\mathit{g}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{h}}}}$
is a multi-linear map
$\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times n}$ such
that
$\mathrm{{\mn@boldsymbol{\mathit{g}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{h}}}}=g_{i}h_{j}e_{ij}$,
where $\left\\{\mathrm{{\mn@boldsymbol{\mathit{e}}}}\right\\}$ form the set of
orthonormal (normal Euclidean) basis vectors for $\mathbb{R}^{n}$, and we use
Einstein notation, summing over repeated indices. and $\mathbb{I}$, the
identity matrix. We solve the equations on a periodic domain of dimension
$L^{3}\equiv\operatorname{\mathcal{V}}$, with $144^{3}$ grid cells, where
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ is the fluid velocity, $\rho$ is the
gas density, $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$ is the magnetic field,
where
$\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}(t)}\right\rangle_{\operatorname{\mathcal{V}}}=0$,
$c_{s}$ is the sound speed, $\mathcal{S}$ is the strain rate tensor,
$\mathcal{S}=(1/2)(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}+[\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}]^{T})+(1/3)(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})\mathbb{I}$444
Note that for an incompressible fluid,
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}=0$,
$\rho=\rho_{0}=\text{const.}$, and then
$\nabla\cdot\nu\mathcal{S}=\nu\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}$,
as expected., and $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$, the stochastic
turbulent acceleration source term that drives the turbulence. In the ISM,
$\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ could be from, for example, supernova
shocks, internal instabilities in the gas, gravity, galactic-scale shocks and
shear, or ambient pressure from the galactic environment (Brunt et al., 2009;
Elmegreen, 2009; Federrath, 2015; Krumholz & Burkhart, 2016; Grisdale et al.,
2017; Jin et al., 2017; Körtgen et al., 2017; Federrath et al., 2017; Krumholz
et al., 2018; Colling et al., 2018; Schruba et al., 2019; Lu et al., 2020).

###  2.1 Turbulent driving

The forcing term $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ follows an Ornstein-
Uhlenbeck process with finite $e$-fold correlation time,
$t_{0}=\ell_{0}/\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}$.
$\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ is constructed in Fourier space with
energy injected on the peak scale
$|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|=2$ (equivalently,
$\ell_{0}=L/2$) and falls off to zero with a parabolic spectrum within
$1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$. On $\ell_{0}$, we
use $t_{0}$ to control the rms velocity, which we set to
$\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/c_{s}=\mathcal{M}=0.5$
in the kinematic stage of the dynamo, which, up to
$2\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}$
in the velocity distribution, corresponds to an incompressible flow. We inject
energy isotropically and solely into the solenoidal
($\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{f}}}}=0$) mode component of
$\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ (as solenoidal driving gives a higher
dynamo efficiency in comparison to compressive driving, see Federrath et al.,
2011; Achikanath Chirakkara et al., 2021). See Federrath et al. (2008);
Federrath et al. (2009); Federrath et al. (2010); Federrath et al. (2022) for
more details about the turbulent driving. We run each of the experiments from
$t/t_{0}=0$ to $t/t_{0}=1000$, writing the 3D field variables to disk every
$t/t_{0}=0.1$ to ensure we produce a holistic dataset that resolves each of
the dynamo stages in time, and is able to capture well-sampled statistics from
the magnetic field in the large $t$ limit.

###  2.2 Plasma Reynolds numbers

Apart from $\mathcal{M}$, there are three main dimensionless numbers that
parameterise our simulations. The first is the hydrodynamic Reynolds number,

$\displaystyle\text{Re}=\frac{|\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})|}{|\nu\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}|}=\frac{\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}\ell_{0}}{\nu},$
(5)

which informs us of the relative strength for the Reynolds stress
$|\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})|$
compared with the dissipation
$|\nu\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{v}}}}|$ (assuming
incompressibility; see footnote 4) in Equation 2 (the momentum equation). This
number also provides a measure of the width for the range of scales that are
self-similar in the turbulence – part of the non-linear turbulent cascade,
i.e., $\ell_{\nu}\ll\ell\ll\ell_{0}$, where
$\ell_{\nu}\sim\text{Re}^{3/4}\ell_{0}$. By setting $\nu$, for a fixed
$\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}\ell_{0}$,
we are able to control Re for each of our simulations. In this study we use
$\text{Re}=500$ for all of our simulations, as indicated in column (2) of
Table 1.

The second dimensionless parameter in our simulations is the magnetic Reynolds
number,

$\displaystyle\text{Rm}=\frac{|\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})|}{|\eta\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}}|}=\frac{\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}\ell_{0}}{\eta},$
(6)

which is analogous to Re, and compares the induction
$|\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})|$
and dissipation $|\eta\nabla^{2}\mathrm{{\mn@boldsymbol{\mathit{b}}}}|$ terms
in Equation 3 (the induction equation). By setting $\eta$ we control Rm and
vary it between 500 and 2000, as indicated in column (3) of Table 1, ensuring
that we are significantly above the critical Rm for the SSD to take place
($\text{Rm}\sim 100$; Ruzmaikin & Sokolov 1981; Haugen et al. 2004;
Schekochihin et al. 2004a; Federrath et al. 2014).

The final dimensionless parameter in our simulations is the Prandtl number,
which is simply the ratio between the two plasma Reynolds numbers,

$\displaystyle\text{Pm}=\frac{\nu}{\eta}=\frac{\text{Rm}}{\text{Re}}.$ (7)

In units of the correlation time of the forcing555Note
$[\nu]=[\eta]\sim\ell^{2}/t_{\rm diffuse}$, and that
$t_{\eta}/t_{0}=(\ell_{\eta}^{2}/\ell_{0})\eta^{-1}\left\langle{\mathit{v}^{2}}\right\rangle^{1/2}$,
which is $t_{\eta}/t_{0}=\text{Rm}^{-1}$ for $\ell_{\eta}\sim\ell_{0}$, which
is valid for the types of highly-magnetised, decaying plasmas that we will be
dealing with throughout our study., Pm is $\sim$ the ratio between diffusion
timescales for the magnetic $t_{\eta}$ and kinetic $t_{\nu}$ fluctuations,
$\text{Pm}=t_{\eta}/t_{\nu}$. For $\text{Pm}=1$, the diffusion timescales are
equal, and thus it naturally follows that $\ell_{\eta}\sim\ell_{\nu}$ where
$\ell_{\eta}$ is the magnetic dissipation scale. Equivalently, Pm also
provides a measure of the scale separation between $\ell_{\nu}$ and
$\ell_{\eta}$. For Kolmogorov (1941) turbulence in a magnetised plasma with
$\text{Pm}\gg 1$, as is the case for the ISM (see Section 6 and footnote 13
for a more detailed discussion), $\ell_{\nu}/\ell_{\eta}\gg 1$, and
$\text{Pm}^{1/2}\sim\ell_{\nu}/\ell_{\eta}$ (Schekochihin et al. 2002b
proposed this relation, which was recently confirmed by Kriel et al. 2022 and
Brandenburg et al. 2022 using direct numerical simulations). Similarly, for
$\text{Pm}\ll 1$ plasmas, $\ell_{\nu}/\ell_{\eta}\ll 1$, and
$\text{Pm}^{1/4}\sim\ell_{\nu}/\ell_{\eta}\ll 1$ (Moffatt, 1961) as is the
case for liquid metal experiments, stars, and planetary plasmas (Rincon,
2019). In this study we will be focusing on $\text{Pm}\geq 1$ plasmas, varying
$\text{Pm}=1-4$, relevant to the ISM, albeit without being able to venture
very far from $\text{Pm}=1$ due to the limited simulation resolution available
to us.

Figure 1: Time-dependent integral energy ratios (left axis) and
$\operatorname{\mathcal{M}_{\text{A}}}$ (right axis) for the $\text{Rm}=500$,
$\text{Rm}=1000$ and $\text{Rm}=2000$ simulations with
$b_{\text{init}}=1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$
(see Table 1). Each column corresponds to a different Rm ensemble, indicated
by the annotation above the column. Different initial magnetic field strengths
are indicated by the colouring, shown in the legend of the rightmost panel.
Regardless of the initial magnetic field strengths, the turbulent dynamo
produces a universal saturation value for each Rm simulation ensemble, which
we annotate in the bottom right corner of each panel.

###  2.3 Initial conditions

The initial velocity field in our simulations is set to
$|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(x,y,z,t=0)|/c_{s}=0$, with units
$c_{s}=1$, and the density field is initialised to a constant value,
$\rho(x,y,z,t=0)=\rho_{0}$, where the density has units $\rho_{0}=1$. We do
not initialise our simulations with a mean-field component of the magnetic
field,
$\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}}\right\rangle_{\operatorname{\mathcal{V}}}=0$,
instead including only a fluctuating-component of the magnetic field,
$\mathrm{{\mn@boldsymbol{\mathit{b}}}}(x,y,z,t=0)/(c_{s}\rho_{0}^{1/2})$,
which has units $c_{s}\rho_{0}^{1/2}=1$. This is an obvious and useful non-
dimensionalisation of Equation 1–4 for an isothermal fluid, which allows us to
scale to an arbitrary dimensionalisation of an isothermal systems, for
example, any approximately isothermal phase of the interstellar gas (Wolfire
et al., 1995).

To explore the universality of the dynamo saturation, we set the initial
magnetic field $\mathrm{{\mn@boldsymbol{\mathit{b}}}}(t=0)$ to be in one of
four configurations for each Pm experiment. In three of the four
configurations we set $\mathrm{{\mn@boldsymbol{\mathit{b}}}}(t=0)$ with the
same $\mathit{b}_{\text{init}}\equiv
1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$ parabolic field as
the driving momentum field (see, Section 2.1) but with
$E_{\text{mag},0}/E_{\rm kin}=\left\\{10^{-10},10^{-2},10^{2}\right\\}$
(column (8) in Table 1), where

$\displaystyle\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}\equiv\frac{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{8\pi}=\frac{L}{16\pi^{2}}\int_{0}^{\infty}\operatorname{d}\\!{k}\,\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta},$
(8)

is the integral magnetic energy. The integral kinetic energy is likewise
defined as

$\displaystyle\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}\equiv\frac{\rho_{0}\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{2}=\frac{\rho_{0}L}{4\pi}\int_{0}^{\infty}\operatorname{d}\\!{k}\,\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta},$
(9)

where
$\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}$
and
$\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}$
are the 1D shell-integrated power spectra. The first (which we call the weak
experiments) of the three is a classical dynamo experiment, which leads to
$E_{\rm mag}$ evolving through all three of the dynamo stages: kinematic,
nonlinear, saturation. The second (sat experiments) are initialised such that
$E_{\text{mag},0}/E_{\rm kin}\approx(E_{\rm mag}/E_{\rm kin})_{\rm sat}$,
where $(E_{\rm mag}/E_{\rm kin})_{\rm sat}$ is the saturated state of the
energy ratio measured from the weak experiment. Next, the third (strong
experiments) are initialised with roughly four orders of magnitude more
$\left\langle{E_{\rm mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ than is
supported by the SSD. The final configuration is also a strong experiment, but
$\mathrm{{\mn@boldsymbol{\mathit{b}}}}(t=0)$ is initialised on smaller scales
than in the other experiments, also with an initial parabolic profile, but
peaking at $|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|=8$ (equivalently,
$\ell=L/8$) and falling off within $\mathit{b}_{\rm
init}=7\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 9$. We call these
experiments initb and compare them to the other strong experiments in Appendix
A. The final configuration allows us to probe how the strong-field experiments
respond to having the bulk of the $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ initialised on high-$k$
modes, close to scales that ought to be dominated by dissipation.

###  2.4 Pure decay experiment

In addition to those experiments already described, we also include one
experiment, decayPm1, which is run without any turbulent forcing field, viz.,
$|\mathrm{{\mn@boldsymbol{\mathit{f}}}}|=0$. This experiment is strong, and
thus allows us to compare decay statistics with our other strong experiments
that are continuously driven. Now we turn our attention to the integral energy
quantities of the plasmas.

##  3 Integral energy quantities

In Figure 1 we plot the integral energy ratio $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ (left axis) and
$\operatorname{\mathcal{M}_{\text{A}}}$ (right axis) as a function of
correlation times of the turbulence driving, $t/t_{0}$. In each panel we show
a different Rm (labelled in the top right; and hence a different Pm), showing
the most resistive simulations on the left and the least on the right. Each
colour represents a different $E_{\text{mag},0}/E_{\rm kin}$, indicated in the
legend on the rightmost panel (weak, black; sat, red; strong, aqua). This
colour scheme will be maintained throughout the remainder of the study.

For each Rm (panel), regardless of $E_{\text{mag},0}/E_{\rm kin}$, the
integral energies for each of the experiments reach the same saturated value,
which varies between roughly $10\%$ and $30\%$ in $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ (corresponding to
$\operatorname{\mathcal{M}_{\text{A}}}=3.5-1.7$, respectively), labelled on
each plot with the orange line (band showing 1$\sigma$ within the averaging
range) and clearly depending upon Rm. We list all the saturated values of the
energy ratio and $\operatorname{\mathcal{M}_{\text{A}}}$ in Table 1. However,
the journey towards the saturation is quite different between the three
different $E_{\text{mag},0}/E_{\rm kin}$. The weak runs explore the full three
stages of the SSD, whilst the sat experiments are immediately in, and
maintained at the saturated state, at least according to $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$. Of course, the magnetic
field may be being reorganised by the turbulence in $k$-space, but regardless,
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ remains constant at the
saturated level. For the first $10t_{0}$ the strong experiments remain strong,
and in fact $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ increases by roughly an order
of magnitude (discussed further in the following section). But these dynamics
are short-lived, and beyond $10t_{0}$ $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ is clearly decaying
exponentially into the universal $t\rightarrow\infty$ state. The timescale of
the dissipation, which we will explore in detail in Section 6, is clearly a
function of Rm (or Pm), with the most resistive experiments taking shorter
$t/t_{0}$ paths towards the saturation. To explore the universality of the
saturation in greater detail we now turn to the scale-by-scale representation
of the same plots in $k$-space.

Figure 2: Top: the temporal evolution (black, $t/t_{0}\sim 0$; aqua,
$t/t_{0}\sim 500$) of the magnetic energy as a function of $k$-modes $E_{\rm
mag}(k)$, normalised by the kinetic energy $E_{\rm kin}(k)$ for the classical
small-scale dynamo experiment, (weak). The red horizontal line indicates
energy equipartition. Bottom: The same as the top panel, except for the
strong-field (strong) decaying dynamo experiment. In the decaying regime, the
low-$k$ modes lose energy slowly (only settling after $t\gtrsim 100t_{0}$).
However, simultaneous to the decay, the high-$k$ modes are being amplified (or
maintained for $\text{Pm}=1$, on roughly the same timescale as the low-$k$
mode decay) by the turbulence. When the low-$k$ modes have decayed, and the
high-$k$ modes have been amplified beyond energy equipartition and the
saturated state of the small-scale dynamo is reached. For all experiments, the
ratios of the spectra exhibit an extended self-similar structure, $E_{\rm
mag}(k)/E_{\rm kin}(k)\propto k^{2}$ (shown in bottom-right plot). We also
show $E_{\rm mag}(k)/E_{\rm kin}(k)\propto k^{3/2}/k^{-5/3}=k^{19/6}$ (bottom-
middle panel) for Kazantsev (1968) magnetic and Kolmogorov (1941) velocity
spectrum. We show a similar plot but for just the $E_{\rm kin}(k)$ in Figure
13 and just $E_{\rm mag}(k)$ in Figure 14.

##  4 Energy spectra

###  4.1 The two different journeys towards saturation

In Figure 2 we show the ratio between the time-dependent 1D magnetic $E_{\rm
mag}(k)=\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}/(8\pi)$
and kinetic $E_{\rm
kin}(k)=\rho_{0}\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}/2$
energy spectrum for the weak experiment (top panel) and strong experiment
(bottom panel). Each spectrum is coloured by $t/t_{0}$, varying from black
$t/t_{0}\sim 0$ to light aqua $t/t_{0}=500$. In each column we show each Rm-Pm
combination, which are all in the saturated state by $500t_{0}$ (see Figure 1)
– the maximum correlation time we plot in this figure. We annotate the energy
equipartition $E_{\rm mag}(k)=E_{\rm kin}(k)$ with the red-dashed line in both
panels, and directly plot the time-evolution of $k_{\text{eq}}:E_{\rm
mag}(k_{\text{eq}})=E_{\rm kin}(k_{\text{eq}})$ (see Appendix B for details on
how we define this scale, which is consistent with the saturation
phenomenology we described in Section 1) for the weak (top) and strong
(bottom) panels in Figure 3, coloured by Rm and Pm.

The top panels of Figure 2 correspond to the classical SSD experiment, which
is where we will begin our analysis. Firstly, within a few $t_{0}$, $E_{\rm
mag}(k)/E_{\rm kin}(k)$ is quickly organised into a self-similar state, peaked
on high-$k$ modes – scales close to (if not at) the resistive scale
$k_{\eta}$. This ratio is maintained through the whole kinematic stage (the
black curves) and saturated state (light aqua curves). We will discuss this
morphology in more detail when focusing on the strong experiments.
$\operatorname{d}\\!{E_{\rm mag}(k)}/\operatorname{d}\\!{(t/t_{0})}$ can be
roughly read off by eye666Note that all spectra are sampled at the same rate,
$\sim 10/t_{0}$, hence “large” differences (see the kinematic stage in the
top-right $\text{Pm}=4$ panel) between any successive spectra correspond to
“large” $\operatorname{d}\\!{E_{\rm mag}(k)}/\operatorname{d}\\!{(t/t_{0})}$,
and likewise for small differences (see low-$k$ modes decaying in the bottom-
right $\text{Pm}=4$ panel)., and seems to be constant across all $k$ in the
kinematic stage (the turbulence grows all modes evenly once the field has been
reorganised). $k_{\text{eq}}$ (the $k$ scale intersecting the red-dashed line)
in Figure 3, shows that as the experiments (top panel) approach saturation,
$k_{\text{eq}}$ moves towards larger scales (smaller $k$), as eddies on
successively larger scales become responsible for amplifying the magnetic
field (Galishnikova et al., 2022). However, as found previously in, e.g.,
Maron et al. (2004), an exact scale-by-scale energy equipartition is not
realised, and in the saturated stage, the magnetic field is able to be
maintained by the turbulence in superequipartition on scales where $E_{\rm
mag}(k>k_{\text{eq}})>E_{\rm kin}(k>k_{\text{eq}})$. This means that the
timescale for the coupling between the kinetic and magnetic energy $t_{\rm
couple}\sim\big{(}\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/L\big{)}^{-1}\big{(}\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}/\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}\big{)}^{-2}$
is shorter than the magnetic energy diffusion timescale
$t_{\eta}\sim\eta^{-1}\ell^{2}$ on these scales (assuming that these scales
are dominated by diffusion) $t_{\rm couple}<t_{\eta}$, i.e., there are
multiple across-field (not necessarily local) coupling events, feeding and
growing the magnetic energy, per events that are able to dissipate it via
Ohmic dissipation (or other means), consistent with the transfer function
analysis performed in Galishnikova et al. (2022) (they frame this phenomena as
the injection energy doing work against the Lorentz force). This is
exacerbated as Pm increases, which increases $t_{\eta}$, and shifts
$k_{\text{eq}}$ to lower $k$-modes. Now we turn our attention to the bottom
panel of Figure 2.

In the bottom panels of Figure 2 we observe a different pathway to the same
saturated $E_{\rm mag}(k)$ state, similarly to what Maron et al. (2004)
discussed in §5.4 of their study. Unlike the weak experiment,
$\operatorname{d}\\!{E_{\rm mag}(k)}/\operatorname{d}\\!{(t/t_{0})}$ is very
different on large and small scales in the strong experiments. Due to the
initialisation of the $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$-field, the
$E_{\rm mag}(k)$ is concentrated on the largest scales for $t_{0}\sim 0$. As
is demonstrated from the top panel in the weak experiment, these modes are not
able to be maintained by the SSD, and start decaying slowly until $E_{\rm
mag}(k)\sim E_{\rm kin}(k)$, and then quickly into the final saturated state
where $E_{\rm mag}(k)<10^{-1}E_{\rm kin}(k)$ (but which depend upon Rm, even
at these low modes). This final fast decay stage can be seen in the integral
energy plots, Figure 1, where $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ crashes shortly before being
maintained in the saturated stage. In contrast, as the low-$k$ magnetic modes
decay, $k_{\text{eq}}$ moves to smaller scales (the opposite as the weak
experiments; bottom panel in Figure 3) and the high-$k$ modes grow into the
self-similar peaked $E_{\rm mag}(k)/E_{\rm kin}(k)$ structure that we observed
in the saturated state of the weak experiment. Because the low-$k$ modes begin
their decay slowly, the growth of the high-$k$ modes boost the integral energy
by some amount for small $t\lesssim 10t_{0}$, as we saw in Figure 1 (at
$t\lesssim 10t_{0}$ the $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ initially grows by an order
of magnitude). But even with these modes growing, the amount of energy being
lost through the decay of low-$k$ modes surpasses the growth and gives rise to
the exponential decay that we observed in Figure 1. To summarise, the magnetic
field initially undergoes a slow decay, but then after some characteristic
time undergoes a fast decay into the saturated state. We show the separate
energy spectra (not the ratio) in Appendix E.

For both the strong and weak experiments, $E_{\rm mag}(k)/E_{\rm kin}(k)$
shows a self-similar structure that extends all the way from $k_{0}$ to
$k_{\rm max}=\text{max}\left\\{E_{\rm mag}(k)/E_{\rm kin}(k)\right\\}$, which
is reminiscent of an extended self-similarity – fractal structure in the
turbulence encapsulating much more than just the scales that exist in the
cascade (inertial range for incompressible turbulence, Benzi et al.,
1993)777Note that extended self-similarity is classically invoked for
measuring extended power-law scalings in velocity structure functions when
each order is normalised by the $3^{\rm rd}$-order velocity structure
function. This means, to strengthen this analogy, one may seek to construct
higher-order statistics of the ratios between the magnetic and kinetic energy,
which we leave for future studies to pursue.. In the bottom-central panel, we
plot $E_{\rm mag}(k)/E_{\rm kin}(k)\propto k^{19/6}=k^{3/2}/k^{-5/3}$, for a
Kazantsev (1968) magnetic field and Kolmogorov (1941) velocity field, but this
is not preferred, and instead a $E_{\rm mag}(k)/E_{\rm kin}(k)\propto k^{2}$
describes the data well over a broad range of $k$, regardless of Rm. However,
the measured saturation as a function of scale is the integral of this ratio
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}\propto\int_{0}^{k}\operatorname{d}\\!{k^{\prime}}\,E_{\rm
mag}(k^{\prime})/E_{\rm kin}(k^{\prime})\propto k^{3}$, adding up all $k$
contributing to the saturation on each scale. This shows that the turbulence
becomes preferentially magnetised on small scales, highlighting the small-
scale (compared to kinetic energies) nature of magnetic field energies.

We expect that this narrative does not change in general with Rm (for a fixed
Re), however the saturated state that is found by the magnetic field, either
through growing from below or decaying from above, will change with increasing
Pm (e.g., $k_{\text{eq}}$ shifts towards lower-$k$ scales, resulting in more
and more scales where $E_{\rm mag}(k)>E_{\rm kin}(k)$, increasing
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$). Next we look at the average
energy spectra in the saturated state for all Rm explored in this study.

Figure 3: The energy equipartition scale, $k_{\rm eq}:E_{\rm mag}(k)=E_{\rm
kin}(k)$ as a function of $t/t_{0}$, coloured by Pm. Top: the classical dynamo
experiment, weak, growing into the saturated stage. Bottom: the experiment
initialised with a strong magnetic field strong, decaying into the saturated
stage.

###  4.2 Saturated energy spectra

We show the kinetic (dashed) and magnetic (solid) energy spectra averaged in
the saturated regime in Figure 4, with the same colouring scheme as in Figure
1, for each Rm (and hence Pm), increasing from left-to-right in each of the
panels, annotated in the top-right of each panel. We label the Boldyrev (2006)
(dynamically aligned turbulence888Or Kraichnan (1965) turbulence, i.e.,
Kolmogorov (1941) with an irreducible magnetic field.; $k^{-3/2}$) and
Kolmogorov (1941) (homogeneous, isotropic turbulence; $k^{-5/3}$) kinetic
energy spectra scalings in the top-right of each panel

As we have shown in the previous section, but now more clearly for different
Rm, the overall shape of the $E_{\rm mag}$ spectra does not depend upon
whether the saturated state is reached from below via the three stages of the
SSD, or from above in the decaying of low-$k$ modes and growth of high-$k$
modes – that is, the saturated state of $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$, for a given Pm, is
universal, and set by the turbulence. The only deviation that can be seen
between the spectra in each panel is at the very high-$k$ modes, which is
likely due to numerical effects, hence we ignore these differences.

Figure 4: Time-averaged energy spectra in the saturated dynamo state, coloured
and organised by $\text{Rm}-\text{Pm}$ as in Figure 1, with different
linestyles for $E_{\rm kin}(k)$ and $E_{\rm mag}(k)$, as indicated by the
legend in the middle panel. Magnetic turbulent wavenumbers directly computed
from the $E_{\rm mag}(k)$ are annotated in each panel, where $k_{\text{cor}}$
is the correlation wavenumber (Equation 10; the scale of the largest
correlated turbulent structures), $k_{\text{u}}$ is the micro or rms scale
(Equation 12; the scale of the smallest correlated turbulent structures) and
$k_{\text{eq}}$ is the energy equipartition scale (Equation 29; the lowest
wavenumber where $E_{\rm mag}(k)=E_{\rm kin}(k)$). The Kazantsev (1968)
($k^{3/2}$; top-left), Kolmogorov (1941) and Boldyrev (2006) ($k^{-5/3}$ and
$k^{-3/2}$, respectively; top-right) spectra are annotated in black.

##  5 Characteristic magnetic field scales

In this section we define six characteristic scales of the magnetic field that
probe either the correlation (and hence turbulence cascade) or the morphology
of the intermittent dissipative structures in the magnetic field, namely
current sheets. We start with the scales directly computed from the magnetic
spectra.

###  5.1 Probing the structure of the magnetic energy

The first scale of interest is the integral scale of the magnetic field,
$k_{\text{cor}}$, which is computed directly from the 1D energy spectrum,

$\displaystyle
k_{\text{cor}}=\left(\int_{0}^{\infty}\operatorname{d}\\!{k}\,k^{-1}\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}\middle/\int_{0}^{\infty}\operatorname{d}\\!{k}\,\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}\right)^{-1},$
(10)

and probes the characteristic size of the largest field correlation. In the
kinematic stage $k_{\text{cor}}\sim\text{Rm}^{1/2}\sim k_{\eta}$ (Galishnikova
et al., 2022; Kriel et al., 2022). The next scale we define is the peak scale
of the spectrum,

$\displaystyle
k_{\text{peak}}=\text{argmax}_{k}\left[\mathcal{I}\left\\{\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}\right\\}\right],$
(11)

where $\mathcal{I}\left\\{\ldots\right\\}$ is a cubic interpolation operator
applied to the spectrum, allowing us to compute $k_{\text{peak}}$ in a model-
free and bin-free fashion. In the kinematic regime, $k_{\text{peak}}\sim
k_{\eta}$ (Kriel et al., 2022). Also note that because $k_{\text{peak}}\sim
k_{\eta}$, $k_{\text{peak}}\sim k_{\text{cor}}$ (see Appendix C for a more
detailed comparison of $k_{\text{peak}}$ and $k_{\text{cor}}$), this tells us
that correlations are being destroyed on $k<k_{\eta}$, all through the
Kazantsev interval ($k^{3/2}$, for which fields are folded; Schekochihin et
al. 2004a) of the spectrum, and piling up into folds around $k_{\eta}$. As the
dynamo approaches the saturated stage, $k_{\text{peak}}$ shifts back to lower
$k$ modes, $k_{\text{peak}}\ll k_{\eta}$ (possibly at the resistive dynamical
timescale, Schekochihin et al. 2002a, which we comment upon in more detail
shortly) consistent with previous theory (Schekochihin et al., 2002a; Xu &
Lazarian, 2016; McKee et al., 2020) and simulations (Seta et al., 2020; Seta &
Federrath, 2021), and correlations are able to develop $k<k_{\eta}$. The third
is the microscale of the magnetic field, $k_{\text{u}}$ (a scale analogous to
the Taylor microscale or the rms scale of the field),

$\displaystyle k_{\text{u}}=$
$\displaystyle\left(\frac{\left\langle{|\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}|^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}\right)^{1/2},$
(12) $\displaystyle=$
$\displaystyle\left(\int_{0}^{\infty}\operatorname{d}\\!{k}\,k^{2}\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}\middle/\int_{0}^{\infty}\operatorname{d}\\!{k}\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})|^{2}}\right\rangle_{\theta}\right)^{1/2}.$
(13)

This scale gives us a characteristic size of the magnetic field gradients,
i.e., structures in the magnetic field that are not smoothed out by
dissipation. In Kolmogorov (1941) turbulence, the energy cascade is defined on
the interval $k_{\text{cor}}>k>k_{\text{u}}$.

First, we annotate the correlation and microscale onto Figure 4. The
separation between the scales defines the range of $k$ within the turbulence
cascade. This range is limited, only over a few $k$ modes, but this is the
nature of doing turbulence studies in the presence of limited grid resolution,
which in turn limit the values for Re and Rm and the separation between
injection and dissipation in the plasma. As Rm increases, we see that
$k_{\text{eq}}$ (shown in Figure 3) moves closer to $k_{\text{cor}}$. We
hypothesise that in the $\text{Rm}\rightarrow\infty$ limit,
$k_{\text{eq}}\rightarrow k_{\text{cor}}$, in turn making $(E_{\rm mag}/E_{\rm
kin})_{\text{sat}}$ independent of Pm. This is because, following the
phenomenology presented in Galishnikova et al. (2022), $k_{\text{eq}}$
describes the largest scales that are being significantly fed by the kinetic
energy modes, and naturally, they ought to support a correlated structure.
Therefore, if $k_{\text{eq}}<k_{\text{cor}}$ are significantly coupled to the
kinetic energy modes, $k_{\text{eq}}\sim k_{\text{cor}}$ as the turbulence
correlates and grows them.

Now we explore these scales throughout the whole temporal evolution of the
simulations. In Figure 5 we show $k_{\text{cor}}$, $k_{\text{u}}$ and
$k_{\text{peak}}$ as a function of $t/t_{0}$, coloured in the same fashion as
Figure 1, and with different linestyles for different Rm experiments. For
fixed Rm, the saturated scales are exactly the same, as we showed previously
in Figure 4. However, similarly to the integral statistics, the journey
towards saturation is very different for the different initial magnetic field
strengths, similar to what we saw in the time-dependent energy spectra, Figure
2. Firstly, for the weak experiments (black) all three scales move from
low-$k$ to high-$k$ in the kinematic regime, $0\leq t/t_{0}\leq 20$. This
represents the whole spectra shifting to high-$k$, where $k_{\text{peak}}\sim
k_{\text{cor}}\sim k_{\eta}$ (Schekochihin et al., 2004b; Xu & Lazarian, 2016;
Kriel et al., 2022; Galishnikova et al., 2022). Xu & Lazarian (2016) predicts
that $k_{\text{peak}}$ starts to move to lower $k$-modes, which we see happen
for $t/t_{0}\gtrsim 20$, into the saturation. The extent of the change between
the scales in the kinematic and saturated regime is larger with increasing Rm.
By comparing $\operatorname{d}\\!{k_{\text{peak}}}$ with
$\operatorname{d}\\!{(}t/t_{0})$ it is clear that
$|\operatorname{d}\\!{k_{\text{peak}}}/\operatorname{d}\\!{(}t/t_{0})|\ll 1$,
and decreases as Rm gets smaller. This suggests that
$\operatorname{d}\\!{k_{\text{peak}}}/\operatorname{d}\\!{(}t/t_{0})$ changes
on diffusive timescales ($t_{\eta}/t_{0}\sim\text{Rm}^{-1})$ rather than
dynamical timescales ($t_{0}$) but we leave the further investigation of this
for future studies.

In contrast to the weak simulations, the strong simulations (aqua) begin very
large-scale (dominated by low-$k$ modes; $k_{\text{cor}}\sim
k_{\text{peak}}\sim 1/L$), and peaked $(k_{\text{u}}\sim 2/L)$. After a
critical $t/t_{0}$, which we will discuss and model in Section 6, the scales
respond to the decaying field and over a short timescale, before finding the
saturation beyond $t/t_{0}\gtrsim 200$. We can see from Figure 2, that this is
roughly at the $t/t_{0}$ where the high-$k$ modes are growing through the
dynamo action (operating on dynamical timescales), once the low-$k$ modes have
all but decayed. Naturally, once $k_{\text{u}}$ moves to high-$k$, the regular
turbulent cascade in the magnetic field (see Figure 8, row BB in Grete et al.
2017) can proceed, and the Ohmic diffusion that we set on the small scales can
destroy the magnetic field. However, it is well known that spontaneous (self-
organised), intermittent sheets of current form in turbulent MHD plasmas
(Makwana et al., 2015; Zhdankin et al., 2015; Chernoglazov et al., 2021),
which occupy a small volume in the plasma, but account for most of the
dissipation (Zhdankin et al., 2013, 2015; Zhdankin et al., 2016; Richard et
al., 2022). Next we study the characteristic scales of these structures using
the tools developed in Schekochihin et al. (2004b).

Figure 5: The characteristic magnetic energy spectra wavenumbers, correlation
$k_{\text{cor}}$ (Equation 10), peak scale $k_{\text{peak}}$ (Equation 11) and
micro $k_{\text{u}}$ (Equation 12), top-to-bottom, respectively, as a function
of $t/t_{0}$. Colours and linestyles are the same as in Figure 1. Smoothing
splines are shown overlayed to the data to reveal the general trends in the
evolution. Characteristic decay timescales for the strong experiments are
shown with vertical lines, which are discussed in detail in Section 6. Figure
6: The same as Figure 5 but for the characteristic scales of the morphology
for the dissipative structures in the magnetised turbulence. From top-to-
bottom, the characteristic length $k_{\parallel}$ (Equation 14), width
$k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
(Equation 15), thickness
$k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
(Equation 16), respectively, of the current sheets.

###  5.2 Probing the magnetic field morphology

The next three scales, which are explored in detail in Schekochihin et al.
(2004b); Seta & Federrath (2021) and Galishnikova et al. (2022), we define are
more intimately connected with the dissipation, morphology and intermittency
of the magnetic field. All of them are fourth-order quantities, hence describe
non-Gaussian structures in the magnetised plasma. The dissipative structures
in MHD turbulence are self-organised into three-dimensional anisotropic sheets
of current (Zhdankin et al., 2013, 2015; Zhdankin et al., 2016), and hence we
consider the following scales as a proxy for measuring the characteristic
size-scales of those structures over the course of the time-evolution in our
simulations. The first scale we define is the length of the current sheet
along $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$. It is the characteristic
fluctuation scale of the magnetic tension force which acts to restore the
magnetic field line curvature, $k_{\parallel}$,

$\displaystyle
k_{\parallel}=\left(\frac{\left\langle{|\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})|^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}^{4}}\right\rangle_{\operatorname{\mathcal{V}}}}\right)^{1/2}.$
(14)

This provides a size scale for the structures along the local magnetic field,
and ought to be associated with the characteristic length that structures
become decorrelated along the field lines. An analogous, but across-field
scale is
$k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$,

$\displaystyle
k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}=4\pi\left(\frac{\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}|^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}^{4}}\right\rangle_{\operatorname{\mathcal{V}}}}\right)^{1/2},$
(15)

where
$\mathrm{{\mn@boldsymbol{\mathit{j}}}}=(\nabla\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})/4\pi$
is the current density, which defines the characteristic size of magnetic
field gradients perpendicular to the magnetic field, providing the
characteristic width of the intermittent structures in the magnetic field. The
final scale is akin to the thickness of magnetic field structures,
$k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$,

$\displaystyle
k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}=4\pi\left(\frac{\left\langle{|\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}|^{2}}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathrm{{\mn@boldsymbol{\mathit{b}}}}^{4}}\right\rangle_{\operatorname{\mathcal{V}}}}\right)^{1/2},$
(16)

which is the characteristic size scale that
$\mathrm{{\mn@boldsymbol{\mathit{b}}}}$ changes signs on, in the orthogonal
direction to $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$, i.e., the characteristic
size scale of the magnetic field reversals (Schekochihin et al., 2004b;
Galishnikova et al., 2022).
$k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}\sim
k_{\eta}$ is a natural organisation of the magnetic field in a shearing
velocity flow, akin to the small-scales in a turbulent plasma (see the
illustrative Figure 13 in Rincon, 2019). Because we are able to define the
three scales of the intermittent structures, we are able to also probe the
three-dimensional geometry, e.g., $k_{\parallel}\lesssim
k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}<k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
is a sheet,
$k_{\parallel}<k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}\sim
k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
is a filament or tube, and $k_{\parallel}\sim
k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}\sim
k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
is a cubic structure. We use the scales defined in Eqs. (14)–(16) to
characterise the morphological differences for the intermittent structures in
our experiments.

We plot each of these scales as a function of $t/t_{0}$ for each of our
simulations in Figure 6, following the same conventions as Figure 5. Once
again, regardless of the initial $\mathrm{{\mn@boldsymbol{\mathit{b}}}}$, the
statistics of the intermittent dissipative structures find the same values in
the saturation, following the same general path through $t/t_{0}$ as the
magnetic field scales we discussed in the previous section. This has
implications for how the dissipation proceeds in our experiments.

For the weak experiments (black), the ordering
$k_{\parallel}<k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}<k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
(a sheet) is consistent with previous findings from Schekochihin et al.
(2004b) and (Galishnikova et al., 2022). This suggests that sheet-like
intermittent structures are ubiquitous throughout the entire SSD process.
However, in the kinematic regime, where magnetic fields are undergoing
exponential growth, the current sheets become smaller in all dimensions. To
see how this effects the overall geometry we plot the length-to-width
($\ell_{\parallel}/\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}=k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}/k_{\parallel}$)
and thickness-to-width
($\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}/\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}=k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}/k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$)
aspect ratios in Figure 7, following the same conventions as the previous
figures, but now instead plotting $\ell\sim k^{-1}$ for each of the scales – a
more natural representation of the aspect ratios. Interestingly, all of the
weak and sat (red) experiments reveal universal aspect ratios, regardless of
Rm, i.e., $\ell_{\parallel}\approx
2\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$,
and
$\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}\approx
0.5\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$,
i.e.,
$\ell_{\parallel}:\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}:\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}\approx
4:2:1$. This means
$\ell_{\parallel}\propto\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}\propto\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
(which should be also $\propto\text{Rm}^{-1/2}$ based on the results from
Galishnikova et al. (2022), in our Pm regime). To summarise, even though the
absolute size-scale of these dissipative structures are changing, for the weak
and sat experiments, the aspect ratio is remaining constant. However, the same
is not true for the strong experiments.

The strong experiments have very different aspect ratios than the other
simulations during the decay of the magnetic field. At early times, well
before saturation, they are elongated along the field by up to roughly a
factor of 4,
$\ell_{\parallel}/\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}\approx
4$, greater than the other simulations. In contrast
$\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}/\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}>1$,
which gives rise to the ordering
$\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}\gtrsim\ell_{\parallel}>\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
(a thin sheet stretched between field reversals). This is coincident with
$\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
being located on the $kL/(2\pi)=1-2$ modes, as shown in Figure 6, which is a
consequence of the magnetic field being dominated by the decay of these modes.
This follows because
$\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
is a proxy for both the field-reversal and dissipation scale, and with only
low-$k$ modes, field reversals and dissipation can only happen on these large
scales. At later times, both
$\ell_{\parallel}/\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
and
$\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}/\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
shrink until finding the same value as the weak and sat experiments. This is
quite remarkable, and highlights how the magnetic field morphology, and even
the characteristics of dissipative structures are being completely controlled
by the turbulence.

Figure 7: The characteristic aspect ratio of the dissipative structures,
current sheets, as a function of $t/t_{0}$, utilising the same conventions as
Figure 8. Top: the characteristic length-to-width ratio,
$\ell_{\parallel}/\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}=k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}/k_{\parallel}$.
Bottom: the characteristic thickness-to-width ratio,
$\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}/\ell_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}=k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}/k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$.

##  6 Implications

We have demonstrated that the saturated state of the SSD is universal
regardless of the initial magnetic field configuration (Appendix A) or
strength. Now let us consider the implications of this result. Consider a
scale-height ($\sim$100$\,\rm{pc}$; Li et al. 2014; Karlsson et al. 2013;
Krumholz & Ting 2018; Beattie et al. 2022) volume of diffuse atomic
interstellar gas (or warm neutral medium; WNM) in the disc of a galaxy in the
present day, with rms velocities such that $\mathcal{M}\sim 1$ (Gaensler et
al., 2011), stirred by detonating supernovae. Through SSD action, the magnetic
field has been amplified into the saturated regime long ago (Xu & Lazarian,
2016; McKee et al., 2020), however, as each shockwave from a supernova
permeates through the medium, the low-$k$ mode field is also amplified through
compressions of the gas, such that $\rho/\rho_{0}\propto(b/b_{0})^{\alpha}$,
where $\alpha$ changes with the magnetic field geometry (Molina et al., 2012;
Tritsis & Tassis, 2016; Mocz & Burkhart, 2018; Beattie et al., 2021) and
$b_{0}$ is the background, saturated field999Note that there is a secondary
effect where solenoidal modes develop in the post-shocked gas, which in turn
facilitate a secondary SSD effect, as shown in Hu et al. (2022). However, it
is not clear if this effect would be very strong if the entire medium was
already in the saturated dynamo state, so for now, we just focus on the
compressively-amplified magnetic field.. In this simple, idealised thought
experiment, three situations may come to fruition:

  1. 1.

if the timescale for dissipating $(t_{\rm decay})$ the compressively-amplified
magnetic field is longer than the timescale between compressions $(t_{\rm
delay})$, i.e., $t_{\rm decay}\gg t_{\rm delay}$, then the magnetic field of
the medium will continuously grow beyond the dynamo saturation101010Of course
the field cannot grow indefinitely via this mechanism, and nor can
$\rho/\rho_{0}$, which would violate mass conservation but let us ignore these
details for the sake of simplicity., and the SSD will not be the mechanism for
maintaining the magnetic field,

  2. 2.

if however $t_{\rm decay}\ll t_{\rm delay}$, then the magnetic field will be
able to decay to the saturation level set by the SSD between compressions,

  3. 3.

if $t_{\rm decay}\sim t_{\rm delay}$, then the magnetic field will be
sustained at $\rho/\rho_{0}\propto(b/b_{0})^{\alpha}$, where $b_{0}$ is the
background magnetic field set by the saturation of the SSD.

Using the stochastic stellar population synthesis code, slug (Krumholz et al.,
2015), Chevance et al. (2022) finds $t_{\rm delay}\sim 10\,\rm{Myr}$ (shown in
their Figure 4), i.e., every $\sim 10\,\rm{Myr}$ a shock permeates through the
medium providing a $\rho/\rho_{0}\propto(b/b_{0})^{\alpha}$ boost to the
magnetic field. Therefore, our task is to determine $t_{\rm decay}$ directly
from our simulations to see which of the three scenarios may be realised.

Figure 8: The integral magnetic field energy normalised by initial magnetic
field (aqua) as a function of the correlation time of the driving for the
strong experiments. The magnetic field undergoes a two-stage decay into the
saturated state of the turbulent dynamo, captured by the Ohmic $+$ current
sheets $+$ slow backreaction model shown in blue, Section 6.1.2. First there
is a slow transfer of magnetic energy into the kinetic energy, and then Ohmic
diffusion in anisotropic current sheets provide a faster decay into the
saturated state. Models for just Ohmic dissipation in current sheets (yellow;
Section 6.1.1) and isotropic Ohmic dissipation (purple; Equation 21) retain
too much magnetic energy. We show a comparison of a purely decaying simulation
decayPm1 (brown-red) to compare with the simulations that decay into the
saturated regime.

###  6.1 Magnetic decay into a saturated state

In this section we provide three models for the decaying magnetic field, all
of which utilise an exact corollary of the induction equation in $k$-space
from Schekochihin et al. (2002a), modified to describe the time-evolution of
$\left\langle{E_{\rm mag}}\right\rangle_{\operatorname{\mathcal{V}}}$. The
first model describes decay as being controlled by Ohmic dissipation (or Joule
heating) at the microscale of the magnetic field (Equation 12). The second is
the same model but modified to include a geometric factor that captures that
the dissipation ought to proceed in the anisotropic current sheets we
described in our previous analysis in Section 5.2. The third has the same
properties as the second but has an additional term that represents the flow
of $E_{\rm mag}$ into $E_{\rm kin}$.

####  6.1.1 Ohmic dissipation at the microscale

As highlighted by Schekochihin et al. (2002a), an exact solution for the
evolution of the integral magnetic energy evolution from the induction
equation is

$\displaystyle\frac{\operatorname{d}\\!{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}}{\operatorname{d}\\!{t}}=\underbrace{8\pi\gamma(t)\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}_{\text{growth}}-\overbrace{8\pi\eta
k^{2}_{\rm u}(t)\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}^{\text{decay}},$ (17)

where,

$\displaystyle\gamma(t)=\frac{\left\langle{(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}):(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})}\right\rangle_{\operatorname{\mathcal{V}}}}{\left\langle{\mathit{b}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}},$
(18)

which simply describes the magnetic energy in terms of growth,
$8\pi\gamma(t)\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$, where $\gamma(t)$ is the
double-dot product of the
$\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}$
and $\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ tensors111111A map
$\mathbb{R}^{n\times n}\times\mathbb{R}^{n\times n}\rightarrow\mathbb{R}$,
defined in Einstein notation
$\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}:\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}}\equiv(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}})_{ij}(\nabla\otimes\mathrm{{\mn@boldsymbol{\mathit{v}}}})_{ij}$.,
and decay $-8\pi\eta k^{2}_{\rm u}(t)\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$. When
$\gamma(t)=\gamma=\text{const.}>0$, it is easy to see that the magnetic field
will undergo exponential growth if $|\gamma|>|8\pi\eta k^{2}_{\rm u}(t)|$. For
a plasma that is undergoing pure decay in units of the correlation scale of
the forcing $t_{0}$, e.g., the $\sim$ timescale of the velocity structures at
the outer scale of the turbulence, the model becomes

$\displaystyle\frac{\operatorname{d}\\!{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}}{\operatorname{d}\\!{t/t_{0}}}$
$\displaystyle=-\frac{8\pi\eta\ell_{0}k^{2}_{\rm
u}(t/t_{0})}{\left\langle{\mathit{v}^{2}}\right\rangle_{\operatorname{\mathcal{V}}}^{1/2}}\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}},$ (19)

which can naturally be expressed in terms of the magnetic Reynolds number,

$\displaystyle\frac{\operatorname{d}\\!{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}}{\operatorname{d}\\!{t/t_{0}}}$
$\displaystyle=-\frac{8\pi}{\text{Rm}}\left(\frac{k_{\rm
u}(t/t_{0})}{k_{0}}\right)^{2}\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}},$ (20)

and solutions exist of the form

$\displaystyle\frac{\left\langle{E_{\rm
mag}(t/t_{0})}\right\rangle_{\operatorname{\mathcal{V}}}}{E_{\rm mag,0}}$
$\displaystyle=\exp\left\\{-\frac{8\pi}{\text{Rm}}\int_{0}^{t/t_{0}}\operatorname{d}\\!{\left(\frac{t^{\prime}}{t_{0}}\right)}\left[\frac{k_{\rm
u}(t^{\prime}/t_{0})}{k_{0}}\right]^{2}\right\\},$ (21)

which result in a time-dependent exponential critically relying upon the
integral of the isotropic scale $k_{\rm u}(t/t_{0})$. We have made all the
measurements that we require to construct this model: we directly measured
$k_{\rm u}(t/t_{0})$ for each of the simulations in the middle panel of Figure
5 and show the numerical integral
$\int_{0}^{t/t_{0}}\operatorname{d}\\!{\left(t^{\prime}/t_{0}\right)}k^{2}_{\rm
u}(t^{\prime}/t_{0})k_{0}^{-2}$ in Figure 12, in Appendix D. We plot the
Equation 21 in Figure 8, in purple, which significantly over-predicts the
amount of $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ present in the plasma (aqua),
and under-predicts the turnover. This means dissipation is too weak at the
start of the decay process, and too strong at the end. Let us take the next
steps into addressing these shortcomings of the model.

We know Ohmic dissipation need not be isotropic, and in fact it ought to be
concentrated within the current sheets (e.g., Zhdankin et al., 2013) explored
in the last section. We therefore define a geometric factor due to the
anisotropic nature of the current sheets by taking a product of the geometric
ratios,

$\displaystyle\Omega(t/t_{0})=\frac{k_{\parallel}(t/t_{0})}{k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}(t/t_{0})}\frac{k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}(t/t_{0})}{k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}(t/t_{0})}=\frac{k_{\parallel}(t/t_{0})}{k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}(t/t_{0})}$
(22)

which is a dimensionless measure for the aspect ratio of a characteristic
current sheet as a function of the correlation time of the forcing. Note from
Figure 7, $\Omega(t/t_{0})<1\;\forall(t/t_{0})$, acting to reduce the
dissipation compared to the isotropic assumption. Placing this geometric
factor in Equation 20, such that the model becomes
$\partial_{t/t_{0}}\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}=-8\pi\Omega(t/t_{0})\text{Rm}^{-1}k_{\rm
u}^{2}(t/t_{0})k_{0}^{-2}\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$, the Ohmic dissipation in
current sheets model is,

$\displaystyle\frac{\left\langle{E_{\rm
mag}(t/t_{0})}\right\rangle_{\operatorname{\mathcal{V}}}}{E_{\rm mag,0}}=$
$\displaystyle\exp\left\\{-\frac{8\pi}{\text{Rm}}\int_{0}^{t/t_{0}}\operatorname{d}\\!{\left(\frac{t^{\prime}}{t_{0}}\right)}\right.$
$\displaystyle\left.\hskip
60.00009pt\ignorespaces\Omega(t^{\prime}/t_{0})\left[\frac{k_{\rm
u}(t^{\prime}/t_{0})}{k_{0}}\right]^{2}\right\\},$ (23)

We perform now the same numerical integration as before, but with
$\Omega(t/t_{0})$ in the integrand, directly constructed from the scales
computed in Section 5.2 (see Appendix D for more details on this integral as a
function of $t/t_{0}$). We plot the Section 6.1.1 in Figure 8, in light
orange. Now the placement of the turnover is reasonable, suggesting that the
Ohmic dissipation in the current sheets plays a role in the later stages of
the decay. However, we now have even more $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ than previously since the
dissipation is concentrated in just current sheets in this model, and we still
are yet to capture the correct slope of the decay at early times. Hence, now
we return to Equation 17 to understand what we may have overlooked.

####  6.1.2 Ohmic dissipation at the microscale and backreaction

As we saw in the previous section, clearly the decay is faster than the purely
Ohmic diffusion case. We assumed that $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ was being lost solely through
dissipation. However, we show in Appendix E, specifically Figure 13, that
while the low-$k$ magnetic field is decaying, it is also transferring energy
into the velocity field via the Lorentz force. This process can be captured in
the Equation 17 model through the first $\gamma(t)$ term. In the regular
dynamo, $\gamma(t)>0$ as energy transfers from the velocity to the magnetic
field. However, if $\gamma(t)<0$, this would capture the opposite effect.
Physically, this may be from reconnection events in the turbulence that
annihilate the low-$k$ field and convert $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ to $\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$. For Sweet (1958)-Parker
(1957)-like reconnection events, in-flow velocities responsibly for squeezing
together the X-point field topology are set by the diffusion timescale of the
plasma. Because the decay is a low-$k$ mode phenomena, where the diffusion is
the slowest, the enhancement in $\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ from the reconnection will
also be slow. Another viable mechanism where $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ may grow $\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ is simply through a strong
backreaction from
$(\nabla\times\mathrm{{\mn@boldsymbol{\mathit{b}}}})\times\mathrm{{\mn@boldsymbol{\mathit{b}}}}\sim\nabla\cdot(\mathrm{{\mn@boldsymbol{\mathit{b}}}}\otimes\mathrm{{\mn@boldsymbol{\mathit{b}}}}+\mathit{b}^{2}\mathbb{I})$
on the momentum, critically relying on gradient in the
$\mathrm{{\mn@boldsymbol{\mathit{b}}}}$-field. Naturally, these gradients are
destroyed on the diffusion timescale, and therefore, like the reconnection,
must only slowly enhance the kinetic energy. Both cases we have described let
us consider the simple case where $\gamma(t)=\gamma=-\eta k_{0}^{2}$, i.e.,
the rate is time-independent and evolving over the (slow) diffusion timescale
of the low-$k$ mode motions. Hence, in units of correlation time, Equation 17
becomes,

$\displaystyle\frac{\operatorname{d}\\!{\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}}}{\operatorname{d}\\!{t/t_{0}}}$
$\displaystyle=-\frac{8\pi}{\text{Rm}}\left(1+\Omega(t/t_{0})\left[\frac{k_{\rm
u}(t^{\prime}/t_{0})}{k_{0}}\right]^{2}\right)\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}},$ (24)

and therefore,

$\displaystyle\frac{\left\langle{E_{\rm
mag}(t/t_{0})}\right\rangle_{\operatorname{\mathcal{V}}}}{E_{\rm mag,0}}=$
$\displaystyle\exp\left\\{-\frac{8\pi}{\text{Rm}}\left(\frac{t}{t_{0}}\right.\right.$
$\displaystyle\left.\left.+\int_{0}^{t/t_{0}}\operatorname{d}\\!{\left(\frac{t^{\prime}}{t_{0}}\right)}\Omega(t^{\prime}/t_{0})\left[\frac{k_{\rm
u}(t^{\prime}/t_{0})}{k_{0}}\right]^{2}\right)\right\\}.$ (25)

We show this model in blue in Figure 8. We find, now that we include the
reduction of $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ by a slow backreaction it
fits perfectly without any free parameters. We interpret this in the following
manner. The decaying system has two stages, (just as the regular dynamo has
two growth stages: exponential, and linear growth stages) but both are
exponential. The first stage is dominated by a slow (diffusion timescale)
backreaction that destroys the large-scale field (Figure 2). This could
possibly through reconnection events that are diffusion-limited121212Limited
by the negligible diffusion at low-$k$ that determines the in-flow velocities
for Sweet (1958)-Parker (1957) reconnection. or just the Lorentz force pushing
against the momentum as it smooths out any $\nabla\mathit{b}^{2}/(8\pi)$ and
reorganises the field. However, once the $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}$ approaches the saturation,
coincident with $k_{\text{u}}$ moving to small-scales (Figure 5), and the
field develops more energy in the high-$k$ modes (Figure 2), Ohmic dissipation
in anisotropic current sheets dominate the dissipation until the saturation is
reached. The latter stage is faster (because it happens at higher-$k$ modes
set by $k_{\text{u}}$) than the former stage, and the Ohmic dissipation stage
is only present when the turbulence is present, i.e., the purely decaying
experiment (dark red in Figure 8) follows the slow backreaction decay. This
indicates that the turbulence controls the final stage of dissipation by
correlating the magnetic field on small scales, in turn shifting
$k_{\text{u}}>1$, where dissipation happens over shorter timescales.

We now consider the characteristic or $e$-fold time of the decay, where
$\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}=(1/e)E_{\text{mag},0}\approx
0.37E_{\text{mag},0}$. Because the Ohmic dissipation only dominates the
dynamics at late-$t/t_{0}$, we consider just the slow backreaction term in
Section 6.1.2, hence, $t_{\eta}/t_{0}=\text{Rm}/(8\pi)$, which we annotate in
Figure 8 with the grey vertical lines for each blue curve, for each Rm, and
for Figure 5, Figure 6, Figure 7 and Figure 12. Noticeably, in all three
figures, $t_{\eta}/t_{0}$ is the timescale where the characteristics scales of
the high-$k$, intermittent, and dissipative scales, $k_{\text{u}}$,
$k_{\parallel}$,
$k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
and
$k_{\mathrm{{\mn@boldsymbol{\mathit{b}}}}\times\mathrm{{\mn@boldsymbol{\mathit{j}}}}}$
react to the decay by shifting to higher-$k$ modes in a dramatic fashion. The
scales of the low-$k$ correlations, i.e., $k_{\text{cor}}$, $k_{\text{peak}}$
lag to respond to $t_{\eta}/t_{0}$, even though this is roughly when the
low-$k$ modes in the magnetic field have decayed significantly.
Mathematically, $k_{\text{cor}}$ and $k_{\text{peak}}$, which we know are
almost identical scales from Section 5, are not sensitive to the high-$k$
modes like the other, higher-order moments of the spectra that we compute for
the intermittent and dissipative scales. They do not respond to the decay for
another few $e$-fold times, which corresponds to the low-$k$ $E_{\rm mag}$
modes moving below $E_{\rm kin}$ energy equipartition, and the exact
transition timescale where the
$\Omega(t/t_{0})k_{\text{u}}^{2}(t/t_{0})k_{0}^{-2}$ integral dominates the
decay model, Section 6.1.2.

####  6.1.3 Interstellar magnetic field decay into saturation versus growth
timescales

We found that
$t_{\eta}=(8\pi)^{-1}\text{Rm}\,t_{0}=(8\pi)^{-1}\text{Rm}\,\ell_{0}(c_{s}\mathcal{M})^{-1}$
is the characteristic timescale for the slow backreaction decay stage in our
simulations. Now, with an Rm estimate we are able to complete our idealised
model for the $\mathcal{M}=1$, $c_{s}\sim 10\,\text{km}\,\text{s}^{-1}$,
$(100\,\rm{pc})^{3}$ volume of the WNM. We use $\text{Rm}=10^{10}$ – a value
that represents somewhere in-between the important ISM diffusion
processes131313 Note that one can in principle choose a number of different Rm
that reflect the different mechanisms for magnetic field diffusion. For Ohmic
diffusion, which may be dominant at sub-AU scales, Ferrière (2020) suggests
$\text{Rm}\sim 10^{18}$ at $\ell_{0}\sim 100\,\rm{pc}$ (see Table 2) and using
Spitzer resistvity values, Schekochihin et al. (2004b) finds $\text{Rm}\sim
10^{14}$. However, considering the strong ionisation in the gas at $\lesssim
1\,\rm{pc}$, hydromagnetic fluctuations may be damped by ion-neutral damping.
Following Krumholz & Burkhart (2016), we estimate $\eta\sim
10^{23}-10^{22}\,\text{cm}^{2}\,\text{s}^{-1}$ for the dissipation associated
with ion-neutral damping on smaller than molecular cloud ($10\,\text{pc}$)
scales. Hence $\text{Rm}=c_{s}\mathcal{M}\ell_{0}/\eta\sim 3\times 10^{3}$ for
$\ell_{0}=100\,\text{pc}$, $c_{s}=10\,\text{km}\,\text{s}^{-1}$ and
$\mathcal{M}=1$. For the purposes of our calculation in this section, using
either of these Rm values does not change the conclusions.. This means,

$\displaystyle t_{\eta}=4\times
10^{9}\left(\frac{\text{Rm}}{10^{10}}\right)\left(\frac{\ell_{0}}{100\,\text{pc}}\right)\left(\frac{c_{s}}{10\,\text{km}\,\text{s}^{-1}}\right)^{-1}\left(\frac{\mathcal{M}}{1.0}\right)^{-1}\,\text{Myr},$
(26) 141414Roughly 3$\times 10^{5}$ times the age of the Universe.

and using $t_{\rm delay}=10\,\rm{Myr}$ we find $t_{\eta}/t_{\rm
delay}\sim\mathcal{O}(10^{8})$, i.e., $t_{\eta}\gg t_{\rm delay}$. We
interpret this as there are $\mathcal{O}(10^{8})$ irreversible amplification
events that compress the medium per unit of time required to decay a single
amplification event. In this idealistic model, this means the magnetic field
grows beyond the dynamo saturation, and the saturated magnetic field from the
SSD no longer becomes the magnetic field that one measures in the interstellar
medium – instead it is a history of all of the reversible amplification events
that are unable to decay quickly. However, in our non-ideal MHD models we may
be missing key phenomena that allow the magnetic field to decay on much
shorter timescales. Whatever this phenomena may be, such as stochastic fast
turbulent reconnection (e.g., Lazarian & Vishniac, 1999), we do not observe it
in our non-ideal decaying strong experiments. We also note that our
calculations throughout this study are based on isothermal, non-ideal
magnetohydrodynamic turbulence, which may approximate some regions of the ISM,
but certainly do not capture the detailed thermal physics, nor the fact that
the ISM is embedded in a galactic gravitational potential (and is also self-
gravitating) with both large-scale and small-scale magnetic field components.
We leave the details of how these phenomena may effect the universality of the
SSD saturated state and the magnetic field decay from superequipartition for
future studies.

##  7 Summary and key findings

Using an ensemble of non-ideal, three-dimensional magnetohydrodynamic
turbulence simulations we study the universality of the small-scale turbulent
dynamo by changing the initial energy of the magnetic field and the
initialisation of the magnetic field in $k$-space. We find that the saturation
and the characteristic scales of the magnetic field are attracted to the same
values, regardless of the initial magnetic field configuration, and even if
the magnetic field is allowed to decay into the saturation. This suggests that
the structure and energy in the magnetic field is completely controlled by the
turbulence and the material dissipation properties of the plasma. We provide a
thorough analysis of the situation where the magnetic field decays into the
saturated dynamo state, and identify two decay stages. Early in the decay, is
it controlled by a flow of magnetic to kinetic energy at the Ohmic diffusion
timescale, which removes energy from the low-$k$ modes. Closer to the
saturation, the decay is dominated by Ohmic diffusion concentrated in
anisotropic current sheets at the microscale of the magnetic field (Equation
12). We list the key findings of this analysis below:

  * •

We establish that the saturated state of the small-scale turbulent dynamo is
independent of the initial conditions of the magnetic field amplitude and
structure, over a large range of highly-relevant statistics, such as the
integral energies (see Figure 1), magnetic spectrum (see Figure 2 and Figure
4), characteristic magnetic scales (see Figure 5) and morphology of the
dissipative structures (see Figure 6 and Figure 7). Hence, the kinematic and
non-linear stages of the classical dynamo are not a unique set of conditions
that determine the final state of the magnetic field, and for a given Pm, it
is instead the turbulence that sets this state for any different magnetic
field configuration.

  * •

We provide a detailed analysis of simulations that decay into the saturated
state of the dynamo in Section 6. We find that as the magnetic field undergoes
a two-stage decay into the final saturated state (see Figure 8). In the first
stage, there is a backreaction between the velocity and magnetic field which
operates at the diffusion timescale of the magnetic field, and feeds the
kinetic energy, growing all modes (see Figure 13). The second stage is
dominated by Ohmic diffusion at the microscale in anisotropic current sheets,
which we characterise in our simulations in Section 5.2. We provide an
analytical model for this two-stage decay, Section 6.1.2, derived directly
from the time-evolution equation for the integral magnetic energy Equation 17.

  * •

In Section 6 we consider the case where compressions from supernova shocks
lead to amplification of the magnetic field above a background dynamo
saturation level in a $(100\,\rm{pc})^{3}$ volume of the warm neutral medium.
We compare the characteristic decay time it would take to get back to the
saturation from our model with the timescale between supernova events and show
that there will be thousands (the exact numbers depend upon the Rm
prescription) of events per decay time. This suggests, unless a faster
mechanism for annihilating magnetic fields is present in the ISM but not
present in our simulations, the magnetic fields we measure may be complex
histories of all of the compressions that have permeated through the gas, and
not from the SSD.

## Acknowledgements

We thank Roland Crocker and Mark R. Krumholz for the discussions regarding
supernova wait time distributions, and Siyao Xu, Jim Stone, Bart Ripperda and
Lachlan Lancaster for the general discussions regarding this work.

J. R. B. acknowledges financial support from the Australian National
University, via the Deakin PhD and Dean’s Higher Degree Research (theoretical
physics) Scholarships and the Australian Government via the Australian
Government Research Training Program Fee-Offset Scholarship and the Australian
Capital Territory Government funded Fulbright scholarship. C. F. acknowledges
funding provided by the Australian Research Council (Future Fellowship
FT180100495), and the Australia-Germany Joint Research Cooperation Scheme (UA-
DAAD). We further acknowledge high-performance computing resources provided by
the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants
pr32lo, pn73fi, and GCS Large-scale project 22542), and the Australian
National Computational Infrastructure (grant ek9) in the framework of the
National Computational Merit Allocation Scheme and the ANU Merit Allocation
Scheme. P. M. acknowledges this work was in part performed under the auspices
of the U.S. Department of Energy by Lawrence Livermore National Laboratory
under contract DE-AC52-07NA27344. Lawrence Livermore National Security, LLC

The simulation software, flash, was in part developed by the Flash Centre for
Computational Science at the Department of Physics and Astronomy of the
University of Rochester. We utilise Federrath et al. (2022) for the turbulent
forcing function, $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$. Data analysis and
visualisation software used in this study: C++ (Stroustrup, 2013), numpy
(Oliphant, 2006; Harris et al., 2020), matplotlib (Hunter, 2007), cython
(Behnel et al., 2011), visit (Childs et al., 2012), scipy (Virtanen et al.,
2020), scikit-image (van der Walt et al., 2014).

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Achikanath Chirakkara et al. (2021) Achikanath Chirakkara R., Federrath C., Trivedi P., Banerjee R., 2021, Phys. Rev. Lett., 126, 091103
  * Batchelor (1950) Batchelor G. K., 1950, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 201, 405
  * Beattie et al. (2021) Beattie J. R., Mocz P., Federrath C., Klessen R. S., 2021, MNRAS, 504, 4354
  * Beattie et al. (2022) Beattie J. R., Krumholz M. R., Federrath C., Sampson M., Crocker R. M., 2022, arXiv e-prints, p. arXiv:2203.13952
  * Behnel et al. (2011) Behnel S., Bradshaw R., Citro C., Dalcin L., Seljebotn D. S., Smith K., 2011, Computing in Science & Engineering, 13, 31
  * Benzi et al. (1993) Benzi R., Ciliberto S., Tripiccione R., Baudet C., Massaioli F., Succi S., 1993, Phys. Rev. E, 48, R29
  * Biermann (1950) Biermann L., 1950, Zeitschrift Naturforschung Teil A, 5, 65
  * Boldyrev (2006) Boldyrev S., 2006, Phys. Rev. Lett., 96, 115002
  * Bouchut et al. (2010) Bouchut F., Klingenberg C., Waagan K., 2010, Numerische Mathematik, 115, 647
  * Braithwaite (2012) Braithwaite J., 2012, MNRAS, 422, 619
  * Brandenburg & Subramanian (2005) Brandenburg A., Subramanian K., 2005, Phys. Rep., 417, 1
  * Brandenburg et al. (2022) Brandenburg A., Rogachevskii I., Schober J., 2022, arXiv e-prints, p. arXiv:2209.08717
  * Brunt et al. (2009) Brunt C. M., Heyer M. H., Mac Low M. M., 2009, A&A, 504, 883
  * Chernoglazov et al. (2021) Chernoglazov A., Ripperda B., Philippov A., 2021, ApJ, 923, L13
  * Chevance et al. (2022) Chevance M., et al., 2022, MNRAS, 509, 272
  * Childs et al. (2012) Childs H., et al., 2012, in , High Performance Visualization–Enabling Extreme-Scale Scientific Insight. Taylor & Francis, pp 357–372
  * Cho et al. (2009) Cho J., Vishniac E. T., Beresnyak A., Lazarian A., Ryu D., 2009, ApJ, 693, 1449
  * Colling et al. (2018) Colling C., Hennebelle P., Geen S., Iffrig O., Bournaud F., 2018, A&A, 620, A21
  * Dubey et al. (2008) Dubey A., et al., 2008, in Pogorelov N. V., Audit E., Zank G. P., eds, Astronomical Society of the Pacific Conference Series Vol. 385, Numerical Modeling of Space Plasma Flows. p. 145
  * Elmegreen (2009) Elmegreen B. G., 2009, in Andersen J., Nordströara m B., Bland -Hawthorn J., eds, IAU Symposium Vol. 254, The Galaxy Disk in Cosmological Context. pp 289–300 (arXiv:0810.5406), doi:10.1017/S1743921308027713
  * Federrath (2015) Federrath C., 2015, MNRAS, 450, 4035
  * Federrath et al. (2008) Federrath C., Klessen R. S., Schmidt W., 2008, ApJ, 688, L79
  * Federrath et al. (2009) Federrath C., Klessen R. S., Schmidt W., 2009, ApJ, 692, 364
  * Federrath et al. (2010) Federrath C., Roman-Duval J., Klessen R., Schmidt W., Mac Low M. M., 2010, A&A, 512
  * Federrath et al. (2011) Federrath C., Chabrier G., Schober J., Banerjee R., Klessen R. S., Schleicher D. R. G., 2011, Phys. Rev. Lett., 107, 114504
  * Federrath et al. (2014) Federrath C., Schober J., Bovino S., Schleicher D. R. G., 2014, ApJ, 797, L19
  * Federrath et al. (2017) Federrath C., et al., 2017, in Crocker R. M., Longmore S. N., Bicknell G. V., eds, IAU Symposium Vol. 322, The Multi-Messenger Astrophysics of the Galactic Centre. pp 123–128 (arXiv:1609.08726), doi:10.1017/S1743921316012357
  * Federrath et al. (2021) Federrath C., Klessen R. S., Iapichino L., Beattie J. R., 2021, Nature Astronomy
  * Federrath et al. (2022) Federrath C., Roman-Duval J., Klessen R. S., Schmidt W., Mac Low M. M., 2022, TG: Turbulence Generator, Astrophysics Source Code Library, record ascl:2204.001 (ascl:2204.001)
  * Ferrière (2020) Ferrière K., 2020, Plasma Physics and Controlled Fusion, 62, 014014
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Gaensler et al. (2011) Gaensler B. M., et al., 2011, Nature, 478, 214
  * Galishnikova et al. (2022) Galishnikova A. K., Kunz M. W., Schekochihin A. A., 2022, arXiv e-prints, p. arXiv:2201.07757
  * Grete et al. (2017) Grete P., O’Shea B. W., Beckwith K., Schmidt W., Christlieb A., 2017, Physics of Plasmas, 24, 092311
  * Grisdale et al. (2017) Grisdale K., Agertz O., Romeo A. B., Renaud F., Read J. I., 2017, MNRAS, 466, 1093
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Haugen et al. (2004) Haugen N. E. L., Brandenburg A., Mee A. J., 2004, MNRAS, 353, 947
  * Hennebelle & Inutsuka (2019) Hennebelle P., Inutsuka S.-i., 2019, Frontiers in Astronomy and Space Sciences, 6, 5
  * Hu et al. (2022) Hu Y., Xu S., Stone J. M., Lazarian A., 2022, arXiv e-prints, p. arXiv:2207.06941
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Jin et al. (2017) Jin K., Salim D. M., Federrath C., Tasker E. J., Habe A., Kainulainen J. T., 2017, MNRAS, 469, 383
  * Karlsson et al. (2013) Karlsson T., Bromm V., Bland-Hawthorn J., 2013, Reviews of Modern Physics, 85, 809
  * Kazantsev (1968) Kazantsev A. P., 1968, Soviet Journal of Experimental and Theoretical Physics, 26, 1031
  * Kolmogorov (1941) Kolmogorov A. N., 1941, Doklady Akademii Nauk Sssr, 30, 301
  * Korpi et al. (1999) Korpi M. J., Brandenburg A., Shukurov A., Tuominen I., Nordlund Å., 1999, ApJ, 514, L99
  * Körtgen et al. (2017) Körtgen B., Federrath C., Banerjee R., 2017, MNRAS, 472, 2496
  * Kraichnan (1965) Kraichnan R. H., 1965, The Physics of Fluids, 8, 1385
  * Kriel et al. (2022) Kriel N., Beattie J. R., Seta A., Federrath C., 2022, MNRAS, 513, 2457
  * Krumholz (2014) Krumholz M. R., 2014, Physics Reports, 539, 49
  * Krumholz & Burkhart (2016) Krumholz M. R., Burkhart B., 2016, MNRAS, 458, 1671
  * Krumholz & Federrath (2019) Krumholz M. R., Federrath C., 2019, Frontiers in Astronomy and Space Sciences, 6, 7
  * Krumholz & Ting (2018) Krumholz M. R., Ting Y.-S., 2018, MNRAS, 475, 2236
  * Krumholz et al. (2015) Krumholz M. R., Fumagalli M., da Silva R. L., Rendahl T., Parra J., 2015, MNRAS, 452, 1447
  * Krumholz et al. (2018) Krumholz M. R., Burkhart B., Forbes J. C., Crocker R. M., 2018, MNRAS, 477, 2716
  * Krumholz et al. (2020) Krumholz M. R., Crocker R. M., Xu S., Lazarian A., Rosevear M. T., Bedwell-Wilson J., 2020, MNRAS, 493, 2817
  * Lancaster et al. (2021) Lancaster L., Ostriker E. C., Kim J.-G., Kim C.-G., 2021, ApJ, 922, L3
  * Lazarian & Vishniac (1999) Lazarian A., Vishniac E. T., 1999, ApJ, 517, 700
  * Li et al. (2014) Li H. B., Goodman A., Sridharan T. K., Houde M., Li Z. Y., Novak G., Tang K. S., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 101 (arXiv:1404.2024), doi:10.2458/azu_uapress_9780816531240-ch005
  * Liao et al. (2019) Liao A. S., et al., 2019, Physics of Plasmas, 26, 032306
  * Lu et al. (2020) Lu Z.-J., Pelkonen V.-M., Padoan P., Pan L., Haugbølle T., Nordlund Å., 2020, arXiv e-prints, p. arXiv:2007.09518
  * Makwana et al. (2015) Makwana K. D., Zhdankin V., Li H., Daughton W., Cattaneo F., 2015, Physics of Plasmas, 22, 042902
  * Maron et al. (2004) Maron J., Cowley S., McWilliams J., 2004, ApJ, 603, 569
  * McKee et al. (2020) McKee C. F., Stacy A., Li P. S., 2020, MNRAS, 496, 5528
  * Menon et al. (2022) Menon S. H., Federrath C., Krumholz M. R., 2022, arXiv e-prints, p. arXiv:2206.14190
  * Mocz & Burkhart (2018) Mocz P., Burkhart B., 2018, MNRAS, 480, 3916
  * Moffatt (1961) Moffatt K., 1961, Journal of Fluid Mechanics, 11, 625–635
  * Molina et al. (2012) Molina F. Z., Glover S. C. O., Federrath C., Klessen R. S., 2012, MNRAS, 423, 2680
  * Oliphant (2006) Oliphant T., 2006, NumPy: A guide to NumPy, USA: Trelgol Publishing, http://www.numpy.org/
  * Parker (1957) Parker E. N., 1957, J. Geophys. Res., 62, 509
  * Richard et al. (2022) Richard T., Lesaffre P., Falgarone E., Lehmann A., 2022, arXiv e-prints, p. arXiv:2206.03045
  * Rincon (2019) Rincon F., 2019, Journal of Plasma Physics, 85, 205850401
  * Ruzmaikin & Sokolov (1981) Ruzmaikin A. A., Sokolov D. D., 1981, Pisma v Astronomicheskii Zhurnal, 7, 701
  * Schekochihin et al. (2002a) Schekochihin A. A., Cowley S. C., Hammett G. W., Maron J. L., McWilliams J. C., 2002a, New Journal of Physics, 4, 84
  * Schekochihin et al. (2002b) Schekochihin A. A., Boldyrev S. A., Kulsrud R. M., 2002b, ApJ, 567, 828
  * Schekochihin et al. (2004a) Schekochihin A. A., Cowley S. C., Maron J. L., McWilliams J. C., 2004a, Physical review letters, 92, 054502
  * Schekochihin et al. (2004b) Schekochihin A. A., Cowley S. C., Taylor S. F., Maron J. L., McWilliams J. C., 2004b, The Astrophysical Journal, 612, 276
  * Schober et al. (2012) Schober J., Schleicher D., Federrath C., Klessen R., Banerjee R., 2012, Phys. Rev. E, 85, 026303
  * Schober et al. (2015) Schober J., Schleicher D. R. G., Federrath C., Bovino S., Klessen R. S., 2015, Phys. Rev. E, 92, 023010
  * Schruba et al. (2019) Schruba A., Kruijssen J. M. D., Leroy A. K., 2019, The Astrophysical Journal, 883, 2
  * Seta & Federrath (2020) Seta A., Federrath C., 2020, MNRAS, 499, 2076
  * Seta & Federrath (2021) Seta A., Federrath C., 2021, Physical Review Fluids, 6, 103701
  * Seta et al. (2020) Seta A., Bushby P. J., Shukurov A., Wood T. S., 2020, Physical Review Fluids, 5, 043702
  * Steinwandel et al. (2020) Steinwandel U. P., Dolag K., Lesch H., Burkert A., 2020, arXiv e-prints, p. arXiv:2012.12905
  * Stroustrup (2013) Stroustrup B., 2013, The C++ Programming Language, 4th edn. Addison-Wesley Professional
  * Subramanian (2016) Subramanian K., 2016, Reports on Progress in Physics, 79, 076901
  * Subramanian (2019) Subramanian K., 2019, Galaxies, 7, 47
  * Sweet (1958) Sweet P. A., 1958, in Lehnert B., ed., Vol. 6, Electromagnetic Phenomena in Cosmical Physics. p. 123
  * Tritsis & Tassis (2016) Tritsis A., Tassis K., 2016, MNRAS, 462, 3602
  * Tzeferacos et al. (2018) Tzeferacos P., et al., 2018, Nature Communications, 9, 591
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Waagan et al. (2011) Waagan K., Federrath C., Klingenberg C., 2011, Journal of Computational Physics, 230, 3331
  * 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
  * Xu & Lazarian (2016) Xu S., Lazarian A., 2016, ApJ, 833, 215
  * Zhdankin et al. (2013) Zhdankin V., Uzdensky D. A., Perez J. C., Boldyrev S., 2013, ApJ, 771, 124
  * Zhdankin et al. (2015) Zhdankin V., Uzdensky D. A., Boldyrev S., 2015, ApJ, 811, 6
  * Zhdankin et al. (2016) Zhdankin V., Boldyrev S., Uzdensky D. A., 2016, Physics of Plasmas, 23, 055705
  * van der Walt et al. (2014) van der Walt S., et al., 2014, PeerJ, 2, e453

##  Appendix A Influence of changing initial magnetic field structure

In the main study we explored the effect of changing different $E_{\rm
mag,0}/E_{\rm kin}$ ratios, and found that the saturated state of the magnetic
field was universal. Now we also how that regardless of the initial structure
of the magnetic field (the initial filter $\mathit{b}_{\rm init}$ convolved
with the $k$-space initialisation), that the magnetic field is reorganised
into the same attractor saturated state.

We show in Figure 9 the same energy ratio plot as in Figure 1 but for two sets
of experiments where $\mathit{b}_{\rm
init}=1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}/2\pi|\leq 3$ (aqua; the
simulations from the main study) and $\mathit{b}_{\rm
init}=7\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}/2\pi|\leq 9$ (purple).
Energy ratios are scaled by Rm to separate the simulations at $t/t_{0}\gg 1$,
and show that regardless of the $\mathit{b}_{\rm init}$, $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ and
$\operatorname{\mathcal{M}_{\text{A}}}$ remains the same for each different
Rm.

Figure 9: $\left\langle{E_{\rm
mag}}\right\rangle_{\operatorname{\mathcal{V}}}/\left\langle{E_{\rm
kin}}\right\rangle_{\operatorname{\mathcal{V}}}$ (left) and
$\operatorname{\mathcal{M}_{\text{A}}}$ (right) as a function of $t/t_{0}$,
scaled by Rm to clearly show the distinct behaviour of the energy ratios as
$t/t_{0}\gg 1$. The strong experiments from the main study are shown in aqua
and experiments with $\mathit{b}_{\rm
init}=7\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}/2\pi|\leq 9$ shown in
purple.

##  Appendix B Energy equipartition scale calculation

We define the energy equipartition mode $k_{\text{eq}}(t/t_{0})$ in Figure 3
as

$\displaystyle\left\\{k(t/t_{0})|\mathcal{I}\left\\{E_{\rm
mag}(k,t/t_{0})\right\\}-\mathcal{I}\left\\{E_{\rm kin}(k,t/t_{0})\right\\}=0\
\right\\},$ (27)

where $\mathcal{I}\left\\{E_{\rm mag}(k,t/t_{0})\right\\}$ and
$\mathcal{I}\left\\{E_{\rm kin}(k,t/t_{0})\right\\}$ are the interpolated
energy spectra. For each $t/t_{0}$ we have a spectrum of ordered $N$
$k(t/t_{0})$ that satisfy Equation 27,

$\displaystyle\left\\{k_{\text{eq},i}\right\\}=\left\\{k_{\text{eq},1},k_{\text{eq},2},k_{\text{eq},3},\ldots
k_{\text{eq},N}\right\\},$ (28)

due to fluctuations in modes deep in the numerical dissipation regime (see
bottom panel of Figure 2). Hence, to be consistent with the relevant
equipartition scale in e.g., (Galishnikova et al., 2022), we take the root at
the lowest-$k$ modes,

$\displaystyle k_{\text{eq}}\equiv\min\left\\{k_{\text{eq},i}\right\\},$ (29)

for each $t/t_{0}$. Naturally, this provides us with a length scale where the
plasma transitions from being dominated by $E_{\rm kin}(k)$ to being dominated
by $E_{\rm mag}(k)$ modes in the plasma, as we show directly in Figure 4.
$k_{\rm eq}$ appears later in the time-evolution of the simulation, and based
on Figure 3 one can see that this changes for different Rm, so we begin
plotting the scale upon the first appearance of it in the simulations.

##  Appendix C Magnetic spectra scale correlations

Schekochihin et al. (2004b) and Galishnikova et al. (2022) assume that
$k_{\text{cor}}$ (Equation 10), the correlation of scale of $E_{\rm mag}(k)$
is proportional to peak energy scale, $k_{\text{peak}}$ (Equation 11), in the
kinematic stage of the SSD. Qualitatively, Figure 5 shows that this seems true
in not only the kinematic stage, but the nonlinear and saturated stage of the
dynamo, and also in the decaying stage of the “anti-dynamo". Here we show two
representative plots in Figure 10 to quantitatively confirm that this is
indeed the case, no matter what the state of the magnetic field.

In the top panel of Figure 10 we show the 2D histogram of
$k_{\text{peak}}$-$k_{\text{cor}}$ for the whole $t/t_{0}$ in the
$\text{Pm}=4$ weak simulation and in bottom we show the corresponding strong
simulation. The colouring indicates whether the data for the scales is before
the saturation (red) or after (purple). As we showed in Figure 5, both of the
scales move to lower-$k$ in the saturated state compared to the kinematic
stage in the weak calculations, and the opposite for the strong calculations.
The opacity in each hexagonal bin corresponds to the amount of data in that
bin, which in turn corresponds to the amount of time spent at that
$(k_{\text{cor}},k_{\text{peak}})$ value. This shows that most of the time,
whether the magnetic field is growing, decaying, or in the saturated state,
$k_{\text{peak}}=k_{\text{cor}}$. During the evolution between the stages this
never deviates by a factor of $\gtrsim 2$, as we show with the $1:2$ and $2:1$
lines.

We make the same plot for $k_{\text{u}}$-$k_{\text{cor}}$ in Figure 11. It
shows that for both the weak and strong experiments, $k_{\text{u}}\approx
2k_{\text{u}}$ through the kinematic phase (top panel) and decay phase (bottom
panel), $t<t_{\rm sat}$. As $t\rightarrow t_{\rm sat}$, $k_{\text{u}}$ and
$k_{\text{cor}}$ move to smaller wavenumbers for the weak experiments, and
higher for the strong experiments, all whilst maintaining a similar
$k_{\text{u}}\approx 2k_{\text{u}}$ relation. The net result is that there is
separation between the $k_{\text{cor}}$ and $k_{\text{u}}$ that is maintained
throughout the entire dynamo and anti-dynamo process.

Figure 10: The 2D distribution of $k_{\text{peak}}$ as a function of
$k_{\text{cor}}$ for the $\text{Pm}=4$ weak (top) and strong (bottom)
simulations, coloured by red for $t/t_{0}$ that are less than $t_{\rm
sat}/t_{0}$, where $t_{\rm sat}$ is the time it takes to get into the
saturated state of the small-scale dynamo and purple for $t/t_{0}$ in the
saturated state. The transparency for each bin is proportional to the
$t/t_{0}$ spent at that $k_{\text{peak}}$/$k_{\text{cor}}$ value. Both
simulations exhibit very close to $1:1$ scaling between the scales, regardless
of the the dynamo stage. Figure 11: The same as Figure 10, but for
$k_{\text{u}}$ (Equation 12) as a function of $k_{\text{cor}}$.

##  Appendix D Integrated microscale of the magnetic field

The analytical decay models in Section 6 all rely upon integrals over the
microscale of the magnetic field $k_{\text{u}}$,
$\int_{0}^{t/t_{0}}\operatorname{d}\\!{\left(t^{\prime}/t_{0}\right)}k^{2}_{\rm
u}(t^{\prime}/t_{0})k_{0}^{-2}$ and
$\int_{0}^{t/t_{0}}\operatorname{d}\\!{\left(t^{\prime}/t_{0}\right)}\Omega(t^{\prime}/t_{0})k^{2}_{\rm
u}(t^{\prime}/t_{0})k_{0}^{-2}$. In Figure 12 we show both the integrals (top
and bottom, respectively) as a function of $t/t_{0}$, with the same
conventions as the rest of the figures in the main study.

All integrals exhibit the same power-laws in the saturated state (and other
states) of the dynamo, e.g.,

$\displaystyle\int_{0}^{t/t_{0}}\operatorname{d}\\!{\left(\frac{t^{\prime}}{t_{0}}\right)}\Omega(t^{\prime}/t_{0})\left[\frac{k_{\text{u}}(t^{\prime}/t_{0})}{k_{0}}\right]^{2}\propto(t/t_{0})^{\alpha},\;t/t_{0}\gg
1,$ (30)

for some power $\alpha$. Furthermore, the strong experiments (aqua) exhibit a
transition between two power laws with the same index but different
proportionality constants. Including $\Omega(t/t_{0})$ in the integrand has a
two-fold effect: (1) the overall dissipation is reduced, as discussed in
Section 6, and (2) the transition between the two power laws in the strong
experiments shifts to larger $t/t_{0}$ – exactly the desired behaviour that we
required to change the decay transition $t/t_{0}$.

Figure 12: Top: the integrated square of $k_{\text{u}}/k_{0}$, Equation 12,
for the magnetic field as a function of the correlation time for the turbulent
forcing. Bottom: the same as the top but with the anisotropic geometric
factor, Equation 22, for the dissipation in the integrand.

##  Appendix E Energy spectra

In Section 4 we explored the time-evolution of the $E_{\rm mag}(k)/E_{\rm
kin}(k)$ ratio and the saturated $E_{\rm mag}(k)$ and $E_{\rm kin}(k)$
spectras, however neither of these plots clearly showed details of the time-
evolution of $E_{\rm kin}(k)$ or $E_{\rm mag}(k)$ separately. In Figure 13 we
show $E_{\rm kin}(k,t)$, using the same panel configuration as in Figure 2,
and likewise for $E_{\rm mag}(k,t)$ in Figure 14.

Notably, for the strong experiments (bottom panel) in Figure 13, as $t/t_{0}$
moves towards saturation, $E_{\rm kin}(k,t)$ grows on all $k$-modes,
indicating that as $E_{\rm mag}$ decays, there is a conversion into $E_{\rm
kin}$ (possibly via very slow reconnection events or the Lorentz force slowly
smoothing out the gradients in the strong magnetic field). The classical
dynamo experiments (top panel) show a shallowing kinetic energy spectrum as
the transition from the kinematic (black) regime, to the saturated (aqua)
regime, opposite to what is expected to happen based on the scale-by-scale
equipartition between the magnetic tension and strain (Galishnikova et al.,
2022). Now we turn our attention to Figure 14. As we discuss throughout the
main text, $E_{\rm mag}(k,t)$ is initially dominated by low-$k$ modes, which
decay as the high-$k$ modes that are coupled to the turbulence grow. This
creates a simultaneous process of growth and decay at in different parts of
the energy spectra.

Figure 13: The same as Figure 2 but for the kinetic energy spectra, showing
that the saturated kinetic energy spectra reduces in energy around the modes
close to the equipartition scale for the weak experiments (top row) and grows
on all scales during the magnetic decay in the strong spectra (bottom panel).
This highlights the duality between decaying and growing into the saturated
state of the small-scale dynamo. Figure 14: The same as Figure 13 but for the
magnetic energy spectra.