# Dynamo effect in unstirred self-gravitating turbulence

Axel Brandenburg1,2,3,4 and Evangelia Ntormousi5,6
1Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes
Alfvéns väg 12, SE-10691 Stockholm, Sweden
2The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-10691 Stockholm, Sweden
3McWilliams Center for Cosmology and Department of Physics, Carnegie Mellon
University, 5000 Forbes Ave, Pittsburgh, PA 15213, USA
4School of Natural Sciences and Medicine, Ilia State University, 3-5
Cholokashvili Avenue, 0194 Tbilisi, Georgia
5Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy
6Institute of Astrophysics, Foundation for Research and Technology-Hellas,
Vasilika Vouton, GR-70013 Heraklion, Greece<EMAIL_ADDRESS>
###### Abstract

In many astrophysical environments, self-gravity can generate kinetic energy,
which, in principle, is available for driving dynamo action. Using direct
numerical simulations, we show that in unstirred self-gravitating subsonic
turbulence with helicity and a magnetic Prandtl number of unity, there is a
critical magnetic Reynolds number of about 25 above which the work done
against the Lorentz force exceeds the Ohmic dissipation. The collapse itself
drives predominantly irrotational motions that cannot be responsible for
dynamo action. We find that, with a weak magnetic field, one-third of the work
done by the gravitational force goes into compressional heating and the
remaining two-thirds go first into kinetic energy of the turbulence before a
fraction of it is converted further into magnetic and finally thermal
energies. Close to the collapse, however, these fractions change toward 1/4
and 3/4 for compressional heating and kinetic energy, respectively. When the
magnetic field is strong, the compressional heating fraction is unchanged. Out
of the remaining kinetic energy, one quarter goes directly into magnetic
energy via work against the Lorentz force. The fraction of vortical motions
diminishes in favor of compressive motions that are almost exclusively driven
by the Jeans instability. For an initially uniform magnetic field, field
amplification at scales larger than those of the initial turbulence are driven
by tangling.

###### keywords:

dynamo — MHD — turbulence — ISM: general

##  1 Introduction

Dynamo action describes the conversion of kinetic energy into magnetic
(Moffatt, 1978). This can also happen under non-stationary conditions, for
example in decaying turbulence, where the kinetic energy tends to decay in
power-law fashion, so the growth of the magnetic field is no longer
exponential in time, as it would be in stationary turbulence (Brandenburg et
al., 2019; Sur, 2019). This type of unsteady energy conversion is expected to
play a role in many astrophysical settings where the magnetic Reynolds number
is high enough.

In the interstellar medium (ISM), as well as on cosmological scales, self-
gravity can be the dominant driver of turbulence (Field et al., 2008; Klessen
& Hennebelle, 2010). In the context of galaxy growth, the gas accretion flows
typically predicted by numerical simulations around galactic disks (Dekel et
al., 2009; Nelson et al., 2015, for some examples) convert gravitational
potential energy into kinetic energy. In the ISM of galaxies, the
gravitational instability of the disk itself has been proposed as the main
source of turbulence (e.g., Bournaud et al., 2010; Krumholz & Burkhart, 2016),
with a contribution that can be as strong as that of supernova (SN) feedback
(Krumholz et al., 2018). Finally, gravitational accretion is believed to be
the main driver of turbulence within molecular clouds (Ibáñez-Mejía et al.,
2017).

While the exact fraction of potential energy that goes into turbulence is
unknown – and probably depends strongly on the environment – it is clear that,
in general, self-gravity can be an important source of kinetic energy. This
kinetic energy can, in turn, be converted into magnetic energy through dynamo
action. Earlier work has shown that in gravitationally unstable flows, the
magnetic energy increases during the linear phase of the collapse, and that
the magnetic energy declines during the nonlinear phase of the dynamo (Sur et
al., 2010, 2012; Xu & Lazarian, 2020).

The kinematic phase of a turbulent dynamo within a collapsing cloud was
considered by Federrath et al. (2011b), who reported exponential growth of the
magnetic field with a Kazantsev spectrum, as was previously found for forced
turbulence. Their dynamo growth rate exceeded the collapse rate, provided
their resolution criterion of 30 grid cells per Jeans length is obeyed at any
position in space and any point in time. As they demonstrated, at lower
resolution, the magnetic field is only amplified because of flux freezing.

Magnetic fields can also be produced by tangling of a large-scale seed. This
type of growth can still occur in two dimensions, where true dynamo action is
impossible according to the Cowling antidynamo theorem; see Hide & Palmer
(1982) for a generalized antidynamo theorem relevant to compressible flows.

An important goal of this work is to characterize dynamo action in unstirred
decaying turbulence, where gravity provides an energy source that can
eventually revert the decay of the turbulence. We approach the question with
direct numerical simulations of systems with different turbulent initial
conditions, and in various degrees of gravitational instability. Our turbulent
initial conditions are almost exclusively subsonic with a Mach number of 0.2.
One motivation behind this choice is that the critical magnetic Reynolds
number for dynamo action increases by about a factor of two when the flow is
supersonic (Haugen et al., 2004b). The low initial Mach number also allows us
to focus on any dynamo action triggered by the collapse-produced turbulence,
rather than by the decaying turbulence from the initial conditions. In fact,
we will show that, as the models evolve, the collapse itself produces
turbulence that eventually dominates the initial flow. However, since
turbulence in many astrophysical environments, such as molecular clouds, is
supersonic (Schneider et al., 2013), we will also present one run with an
initial Mach number of two. Furthermore, given that the sonic Mach number is
scale-dependent (Federrath et al., 2021), our early subsonic phase might still
be applicable to correspondingly small scales.

In driven turbulence, dynamo action can be adequately characterized by the
growth rate, evaluated as the time derivative of the root-mean-square (rms)
magnetic field. In stationary conditions, this quantity stays reasonably
constant with time. However, in the non-stationary conditions that we study
here, namely decaying turbulence and turbulence generated by gravitational
collapse, the magnetic field no longer grows exponentially, and the growth
rate cannot be used as a dynamo criterion. Therefore, in this study we explore
new, more general dynamo criteria that allow for non-stationary turbulence
conditions. We decided to base our dynamo criterion on the work against the
Lorentz force, where the magnetic curvature force plays the dominant role.
When this work exceeds the Joule dissipation, it might be a dynamo. This
definition of a possible dynamo agrees with the standard definition of a
positive growth rate when the flow is steady, but, unlike any dynamo criterion
proposed so far, it can easily be applied to unsteady flows. It does not,
however, distinguish dynamos in three dimensions from just temporary
amplification through tangling and compression, as can be seen, for example,
in two dimensions. To exclude the effects of two-dimensional (2-D) compression
or tangling, we propose splitting the Lorentz work term into two
contributions, of which one is absent in 2-D. This leads to an additional
criterion that must be satisfied for dynamo action.

We use high-resolution numerical simulations with fixed kinematic viscosity
and magnetic diffusivity to be able to define the threshold for dynamo action.
Note also that, unlike codes with adaptive mesh refinement, where the accuracy
of the solution varies in space (see, for example, Federrath et al., 2010), we
resolve all regions in space equally well. There should therefore be no doubt
that our velocity spectra and other diagnostics are representative of the
domain as a whole.

In this paper, we first define our model (Sect. 2). We then present the
results for the energy spectra and the energy conversion rates, as well as
characteristic wavenumbers and dynamo excitation conditions for weak initial
magnetic fields (Sect. 3). We then compare our results with those for strong
initial magnetic fields (Sect. 4), and conclude in Sect. 5.

##  2 The model

###  2.1 Governing equations

We consider an isothermal gas with sound speed $c_{\rm s}$ in a cubic periodic
domain of size $L^{3}$, so the smallest wavenumber is $k_{1}=2\pi/L$. The
pressure is given by $p=\rho c_{\rm s}^{2}$, where $\rho$ is the density. The
governing equations are (Passot et al., 1995)

$\nabla^{2}\Phi=4\pi G\left(\rho-\rho_{0}\right),$ (1) ${{\rm
D}{}\bm{u}\over{\rm D}{}t}=-{\bm{\nabla}}\left(c_{\rm
s}^{2}\ln\rho+\Phi\right)+\frac{1}{\rho}\left(\bm{J}\times\bm{B}+{\bm{\nabla}}\cdot
2\rho\nu\mbox{\boldmath${\sf S}$}{}\right),$ (2) ${{\rm D}{}\ln\rho\over{\rm
D}{}t}=-{\bm{\nabla}}\cdot\bm{u},$ (3) ${\partial\bm{A}\over\partial
t}=\bm{u}\times\bm{B}-\eta\mu_{0}\bm{J},$ (4)

where $\Phi$ is the gravitational potential, $G$ is Newton’s constant,
$\rho_{0}$ is the spatially averaged density, which is constant in time
because of mass conservation, $\bm{u}$ is the velocity,
$\bm{J}={\bm{\nabla}}\times\bm{B}/\mu_{0}$ is the current density, $\mu_{0}$
is the vacuum permeability, $\bm{B}={\bm{\nabla}}\times\bm{A}$ is the magnetic
field in terms of the magnetic vector potential, ${\sf
S}_{ij}={\textstyle{1\over 2}}(u_{i,j}+u_{j,i})-{\textstyle{1\over
3}}\delta_{ij}{\bm{\nabla}}\cdot\bm{u}$ are the components of the rate-of-
strain tensor with commas denoting partial derivatives, and $\nu$ is the
kinematic viscosity.

The occurrence of the constant $\rho_{0}$ in Eq. (1) comes from a change of
coordinates to a comoving reference frame that is following the global
expansion of the background medium (see Alecian & Leorat, 1988). This is a
consequence of working with an infinite (unbounded) medium. Such a medium can
be stationary, but not static. This was also explained by Falco et al. (2013),
who clarified why the famous Jeans swindle (Binney & Tremaine, 2008) actually
works. In this accelerated frame, Eqs. (1)–(4) describe the departure of
collapsing structures from the background flow.

Linearizing Eqs. (1)–(3) around $\rho=\rho_{0}$ and $\bm{u}=0$, and assuming
the perturbations to be proportional to $e^{{\rm i}\bm{k}\cdot\bm{x}+\sigma
t}$ yields the dispersion relation $\sigma^{2}=\sigma_{\rm J}^{2}-c_{\rm
s}^{2}k^{2}$, where $\sigma_{\rm J}^{2}=4\pi G\rho_{0}$ with $\sigma_{\rm J}$
being the gravitational or Jeans growth rate (Jeans, 1902). The pressureless
free-fall time is $t_{\rm ff}=\sqrt{3/8}\,\pi/\sigma_{\rm J}\approx
1.92/\sigma_{\rm J}$ (Shu, 1992). The Jeans wavenumber is $k_{\rm
J}=\sigma_{\rm J}/c_{\rm s}$, and the Jeans length is then $\lambda_{\rm
J}=2\pi/k_{\rm J}$ (see, e.g., Bonazzola et al., 1987; Truelove et al., 1997).
According to the classical Jeans criterion for gravitational instability, an
interstellar gas cloud will collapse if its free-fall time is shorter than the
sound crossing time in its interior, or, more specifically, $t_{\rm ff}c_{\rm
s}k_{1}<1.92$.

###  2.2 Diagnostic quantities

####  2.2.1 Energetics

Throughout this paper, we use periodic boundary conditions, so all surface
integrals vanish and no mass is lost. In the following, volume averages are
denoted by angle brackets. It is instructive to inspect the evolution
equations of mean potential, kinetic, and magnetic energy densities, ${\cal
E}_{\rm P}=-\langle({\bm{\nabla}}\Phi)^{2}\rangle/8\pi G$, ${\cal E}_{\rm
K}=\langle\rho_{0}\bm{u}^{2}\rangle/2$, and ${\cal E}_{\rm
M}=\langle\bm{B}^{2}\rangle/2\mu_{0}$, respectively (Banerjee & Kritsuk,
2018). They are given by

$\frac{{\rm d}{}{\cal E}_{\rm P}}{{\rm d}{}t}=-W_{\rm J},$ (5) $\frac{{\rm
d}{}{\cal E}_{\rm K}}{{\rm d}{}t}=W_{\rm P}+W_{\rm J}+W_{\rm L}-Q_{\rm K},$
(6) $\frac{{\rm d}{}{\cal E}_{\rm M}}{{\rm d}{}t}=-W_{\rm L}-Q_{\rm M},$ (7)

where $W_{\rm P}=-\langle\bm{u}\cdot{\bm{\nabla}}p\rangle=\langle
p{\bm{\nabla}}\cdot\bm{u}\rangle$ is the work done by the pressure force,
$W_{\rm J}=-\langle\rho\bm{u}\cdot{\bm{\nabla}}\Phi\rangle$ is the work done
by the gravity term, $W_{\rm L}=\langle\bm{u}\cdot(\bm{J}\times\bm{B})\rangle$
is the work done by the Lorentz force, and $Q_{\rm K}=\langle
2\rho\nu\mbox{\boldmath${\sf S}$}{}^{2}\rangle$ and $Q_{\rm
M}=\langle\mu_{0}\eta\bm{J}^{2}\rangle$ are the viscous and Joule dissipation
terms. The thermal energy density is sourced by the terms $-W_{\rm P}+Q_{\rm
K}+Q_{\rm M}$, but with the employed isothermal equation of state, the thermal
energy density is not evolved.

The work done by the gravity term ($W_{\rm J}>0$) leads to a decrease of the
potential energy density and to an increase in the kinetic energy density.
During the collapse, the virial parameter $\alpha_{\rm vir}=2{\cal E}_{\rm
K}/|{\cal E}_{\rm P}|$ is expected to be around unity, but this expectation
can be different at large Mach numbers (Ntormousi & Hennebelle, 2019) and for
strong magnetic fields; see Federrath & Klessen (2012), who also emphasize the
difference between periodic setups, such as ours, and isolated spheres.

We have defined all work terms such that they enter with a plus sign in Eq.
(6), i.e., they lead to an increase in the kinetic energy density if they are
positive, and thus to a loss in some other energy reservoir. As noted above, a
positive $W_{\rm J}$ term leads to a loss of potential energy. Likewise, a
positive $W_{\rm P}$ term leads to a loss in thermal energy. Gravitational
collapse however, leads to compressional heating and $W_{\rm P}$ is therefore
negative. Furthermore, dynamo action leads to a growth in magnetic energy if
$W_{\rm L}$ is negative.

The $W_{\rm L}$ term can be split into three constituents: $W_{\rm L}^{\rm
c}=-\langle\bm{u}\cdot{\bm{\nabla}}\bm{B}^{2}/2\mu_{0}\rangle$, $W_{\rm
L}^{\|}=\langle\bm{u}\cdot(\bm{B}\cdot{\bm{\nabla}}\bm{B}/\mu_{0})_{\|}\rangle$,
and $W_{\rm
L}^{\perp}=\langle\bm{u}\cdot(\bm{B}\cdot{\bm{\nabla}}\bm{B}/\mu_{0})_{\perp}\rangle$.
Here, $-{\bm{\nabla}}\bm{B}^{2}/2\mu_{0}$ is the magnetic pressure
contribution of $\bm{J}\times\bm{B}$, and
$(\bm{B}\cdot{\bm{\nabla}}\bm{B}/\mu_{0})_{\|}$ and
$(\bm{B}\cdot{\bm{\nabla}}\bm{B}/\mu_{0})_{\perp}$ are the stretching terms
along and perpendicular to the magnetic field. The last two forces are also
referred to as tension and curvature forces; see Nordlund et al. (1992) for
their contributions to a convective dynamo.

In the following, we also decompose $W_{\rm L}$ by writing it as $W_{\rm
L}=-\langle\bm{J}\cdot(\bm{u}\times\bm{B})\rangle$ and expanding the curl to
get

$-\langle\bm{J}\cdot(\bm{u}\times\bm{B})\rangle=\langle
J_{i}u_{j}(A_{i,j}-A_{j,i}\rangle\equiv W_{\rm L}^{\rm 2D}+W_{\rm L}^{\rm
3D}.$ (8)

Here, we make use of the fact that the Weyl gauge has been used in Eq. (4). In
two dimensions, the magnetic field can be represented as
$\bm{B}={\bm{\nabla}}\times A_{z}\hat{\mbox{\boldmath$z$}}{}$, with its $x$
and $y$ components lying in the $xy$ plane. Then the term $W_{\rm L}^{\rm
3D}=-\langle J_{i}u_{j}A_{j,i}\rangle$ vanishes in 2-D. Thus, we can identify
$W_{\rm L}^{\rm 3D}$ with a contribution that characterizes the 3-D nature of
the system and can therefore be a proxy for dynamo action, provided $W_{\rm
L}^{\rm 3D}$ is large enough.

To characterize the flow of energy, it is convenient to define the fractions
$\epsilon_{\rm J}^{\rm P}\equiv-W_{\rm P}/W_{\rm J}$, $\epsilon_{\rm J}^{\rm
L}\equiv-W_{\rm L}/W_{\rm J}$, $\epsilon_{\rm J}^{\rm+K}\equiv\dot{\cal
E}_{\rm K}/W_{\rm J}$, and $\epsilon_{\rm J}^{\rm-K}=Q_{\rm K}/W_{\rm J}$.
Likewise, we define the fractions $\epsilon_{\rm L}^{\rm+M}\equiv\dot{\cal
E}_{\rm M}/(-W_{\rm L})$, and $\epsilon_{\rm J}^{\rm-M}=Q_{\rm M}/(-W_{\rm
L})$. To characterize the growth or decay of the magnetic field, we define the
nondimensional ratio $\epsilon_{\rm M}^{\Delta}=(-W_{\rm L}-Q_{\rm M})/Q_{\rm
M}$. A related quantity is the pseudo (or instantaneous) growth rate of
magnetic energy, $\gamma=(-W_{\rm L}-Q_{\rm M})/{\cal E}_{\rm M}$, which can
be divided into the contributions $\gamma_{\rm
c}+\gamma_{\|}+\gamma_{\perp}=\gamma$ from compression and stretching parallel
and perpendicular to the magnetic field, where $\gamma_{\perp}=(-W_{\rm
L}^{\perp}-Q_{\rm M})/E_{\rm M}$ will play the most important role, and
$\gamma_{\rm c}=-W_{\rm L}^{\|}/E_{\rm M}$ and $\gamma_{\|}=-W_{\rm L}^{\rm
c}/E_{\rm M}$ contribute either later or in the presence of strong initial
magnetic fields. Likewise, we define $\gamma_{\rm 2D}=-(W_{\rm L}^{\rm
2D}+Q_{\rm M})/E_{\rm M}$ and $\gamma_{\rm 3D}=-W_{\rm L}^{\rm 3D}/E_{\rm M}$,
so that $\gamma_{\rm 2D}+\gamma_{\rm 3D}=\gamma$.

####  2.2.2 Characteristic wavenumbers

To characterize the compressive and solenoidal flow components, it is
convenient to compute the rms velocity divergence,
$({\bm{\nabla}}\cdot\bm{u})_{\rm
rms}=\left\langle({\bm{\nabla}}\cdot\bm{u})^{2}\right\rangle^{1/2}$, and the
rms vorticity, $\omega_{\rm rms}=\langle\bm{\omega}^{2}\rangle^{1/2}$, where
$\bm{\omega}={\bm{\nabla}}\times\bm{u}$, and to define

$k_{\,{\bm{\nabla}}\cdot\bm{u}}=({\bm{\nabla}}\cdot\bm{u})_{\rm rms}/u_{\rm
rms},$ (9) $k_{\bm{\omega}}=\omega_{\rm rms}/u_{\rm rms},$ (10)

which have the dimension of a wavenumber. Since the flow is helical, we can
also define the wavenumber

$k_{\bm{\omega}\cdot\bm{u}}=|\langle\bm{\omega}\cdot\bm{u}\rangle|/u_{\rm
rms}^{2},$ (11)

which characterizes the typical wavenumber where helicity plays a role. Large
values of $k_{\,{\bm{\nabla}}\cdot\bm{u}}$, $k_{\bm{\omega}}$, and
$k_{\bm{\omega}\cdot\bm{u}}$ imply strong flow divergences or compressions,
strong vortices, and strong swirls, respectively. To characterize the flow
compression from the gravitational collapse, we also define

$k_{p{\bm{\nabla}}\cdot\bm{u}}=-\langle
p{\bm{\nabla}}\cdot\bm{u}\rangle/p_{0}u_{\rm rms}\quad\mbox{(when
$k_{p{\bm{\nabla}}\cdot\bm{u}}>0$)},$ (12)

where $p_{0}=\rho_{0}c_{\rm s}^{2}$ has been introduced for brevity. The
relevance of $k_{p{\bm{\nabla}}\cdot\bm{u}}$ to the collapse phenomenon is
motivated by the fact that a strong flow compression or flow convergence
(${\bm{\nabla}}\cdot\bm{u}<0$) correlates with pressure (defined in the
beginning of Sect. 2.1). Indeed, it turns out that
$k_{p{\bm{\nabla}}\cdot\bm{u}}$ is very small prior to collapse, but it
approaches $k_{{\bm{\nabla}}\cdot\bm{u}}$ close to the collapse.

####  2.2.3 Spectra

We define the kinetic and magnetic energy spectra, $E_{\rm K}(k,t)$ and
$E_{\rm M}(k,t)$, respectively. They are normalized such that $\int E_{\rm
K}(k,t)\,{\rm d}{}k={\cal E}_{\rm K}$ and $\int E_{\rm M}(k,t)\,{\rm
d}{}k={\cal E}_{\rm M}$. It can be advantageous to express them as wavenumber-
dependent Reynolds and Lundquist numbers by defining a velocity and a magnetic
field,

$u_{k}(t)=\sqrt{2kE_{\rm K}(k,t)/\rho_{0}},\quad
B_{k}(t)=\sqrt{2\mu_{0}kE_{\rm M}(k,t)},$ (13)

respectively. We then define

$\mbox{\rm Re}_{k}(t)=u_{k}(t)/\nu k\quad\mbox{and}\quad\mbox{\rm
Lu}_{k}(t)=B_{k}(t)/(\sqrt{\mu_{0}\rho_{0}}\,\eta k).$ (14)

A Kolmogorov-type spectrum with $E_{\rm K}(k)\propto k^{-5/3}$ corresponds
then to $u_{k}\propto k^{-1/3}$ and $\mbox{\rm Re}_{k}\propto k^{-4/3}$. In
the following, we also quote these values at $k=k_{\rm f}$, where the initial
kinetic energy spectrum peaks. We also define the Reynolds number $\mbox{\rm
Re}_{t}$ based on the actual rms velocity, with $\mbox{\rm Re}_{t_{*}}$
denoting the value at the time $t_{*}$, when the exponentially growing gas
motions from the Jeans instability begin to dominate over the initial
turbulence. In the present normalization, the values of $\mbox{\rm Re}_{k_{\rm
f}}$ and $\mbox{\rm Re}_{t}$ are close to the Taylor microscale Reynolds
number (Tennekes & Lumley, 1972), which is universally defined as $\mbox{\rm
Re}_{\lambda}=v^{\prime}\lambda_{\rm Tay}/\nu$. Here, $v^{\prime}=u_{\rm
rms}/\sqrt{3}$ is the one-dimensional rms velocity and $\lambda_{\rm
Tay}=\sqrt{15\nu\rho_{0}/Q_{\rm K}}\,v^{\prime}$ is the Taylor microscale.

Kinetic and magnetic helicity spectra, $H_{\rm K}(k,t)$ and $H_{\rm M}(k,t)$,
are normalized such that $\int H_{\rm K}(k,t)\,{\rm
d}{}k=\langle\bm{\omega}\cdot\bm{u}\rangle$ and $\int H_{\rm M}(k,t)\,{\rm
d}{}k=\langle\bm{A}\cdot\bm{B}\rangle$ are the mean kinetic and magnetic
helicities. They obey the realizability conditions $|H_{\rm K}(k)|\leq
2kE_{\rm K}(k)$ and $|H_{\rm M}(k)|\leq(2/k)E_{\rm M}(k)$ (Moffatt, 1978). It
is then convenient to plot the relative helicities given by the corresponding
ratios $H_{\rm K}(k)/2kE_{\rm K}(k)$ and $kH_{\rm M}(k)/2E_{\rm M}(k)$,
respectively.

We also plot the enstrophy and logarithmic density spectra, $E_{\omega}(k)$
and $E_{\ln\rho}(k)$, respectively. They are normalized such that $\int
E_{\omega}(k)\,{\rm d}{}k=\langle\bm{\omega}^{2}\rangle/2$ and $\int
E_{\ln\rho}(k)\,{\rm d}{}k=\langle(\ln\rho)^{2}\rangle$. Here,
$E_{\omega}(k)/k^{2}$ corresponds to the kinetic energy spectrum of the
vortical part of the velocity, and $E_{\rho}(k)$ reflects its irrotational
part. Finally, the potential energy spectrum is normalized such that $\int
E_{\rm P}(k)\,{\rm d}{}k={\cal E}_{\rm P}$.

###  2.3 Units and parameters

In the plots shown in this paper, we express velocities in units of $c_{\rm
s}$, lengths in units of $k_{1}^{-1}$ (defined in the beginning of Sect. 2.1),
density in units of $\rho_{0}$, and the magnetic field in units of
$\sqrt{\mu_{0}\rho_{0}}c_{\rm s}$. In practice, we do this by choosing in our
simulations $c_{\rm s}=k_{1}=\rho_{0}=\mu_{0}=1$. In the following, to remind
the readers, we often include relevant combinations of $c_{\rm s}$ and $k_{1}$
when specifying the numerical values, but in other cases, especially in the
table entries, we simply omit them to avoid lengthy notation.

Since we do not invoke cooling or any other processes that depend on
dimensions, our simulations can be scaled to any arbitrary system by choosing
physical values for $c_{\rm s}$, $k_{1}$, and $\rho_{0}$. To illustrate the
normalization used in the simulations, let us consider, as an example, the
case $c_{\rm s}=1\,{\rm km\,s}^{-1}$ and $k_{1}=1\,{\rm pc}^{-1}$. Then, our
time unit is $(c_{\rm s}k_{1})^{-1}=0.98\,{\rm Myr}$, so we can think of our
normalized time as $1\,{\rm Myr}$. Considering a typical density of
$\rho_{0}=10^{-21}\,{\rm g}\,{\rm cm}^{-3}$ for the dense regions of the ISM,
we have $\sigma_{\rm J}^{-1}=1.1\,{\rm Myr}$, or $t_{\rm ff}=2.1\,{\rm Myr}$.
Then, the corresponding Jeans wavenumber is $k_{\rm J}=0.9\,{\rm pc}^{-1}$, so
the Jeans length is $\lambda_{\rm J}=7\,{\rm pc}$. For $k_{1}=1\,{\rm
pc}^{-1}$, the side length of the computational domain is $6.28\,{\rm
pc}\approx 0.9\lambda_{\rm J}$. The corresponding normalized (nondimensional)
quantities are then $\sigma_{\rm J}^{-1}c_{\rm s}k_{1}=1.1$ and $\sigma_{\rm
J}/c_{\rm s}k_{1}=k_{\rm J}/k_{1}=0.9$. All work terms are given in units of
$\rho_{0}c_{\rm s}^{3}k_{1}$, which corresponds to $1\,{\rm g}\,{\rm
cm}^{-3}(\,{\rm km}\,{\rm s}^{-1})^{3}/\,{\rm pc}=10^{-11}\,{\rm erg}\,{\rm
cm}^{-3}\,{\rm Myr}^{-1}$ or $0.0024\,L_{\odot}\,{\rm pc}^{-3}$. In this work,
we choose two values for $\sigma_{\rm J}/c_{\rm s}k_{1}$, 2 and 5, which means
that our computational domain is two or five Jeans lengths long, and our mesh
of $2048^{3}$ cells resolves the Jeans length initially with 1024 or 410
points, respectively. As the collapse proceeds, the maximum density increases
and the nominal Jeans length decreases. To stay within the resolution
criterion of Federrath et al. (2011b) of 30 mesh points per Jeans length, we
can only trust the time before the maximum density has exceeded the initial
value by a factor of $(1024/30)^{2}\approx 1200$ and $(410/30)^{2}\approx 200$
for the cases with $\sigma_{\rm J}/c_{\rm s}k_{1}=2$ and 5, respectively.

As mentioned earlier, we focus on subsonic turbulence, where dynamo action is
most easily obtained (Haugen et al., 2004b; Federrath et al., 2011a). In the
context of molecular cloud contraction, this choice puts them in the regime of
low-mass pre-stellar cores. While molecular clouds are supersonically
turbulent on large scales, low-mass pre-stellar cores are subsonic (Larson,
1981; Myers et al., 1991; André et al., 2007). Such weak motions could
originate from the decay of larger-scale turbulence (e.g., Hennebelle &
Falgarone, 2012).

We choose the amplitude of the initial velocity field such that the initial
Mach number $\mbox{\rm Ma}=u_{\rm rms}/c_{\rm s}$ is around 0.1. The turnover
time is given by $\tau=(u_{\rm rms}k_{\rm f})^{-1}$. We are interested in the
cases where the turnover time is comparable to or less than the free-fall time
scale $\sigma_{\rm J}^{-1}$, where $\sigma_{\rm J}$ must be larger than unity
(in units of $c_{\rm s}k_{1}$) for the Jeans instability to be excited. Given
that $\mbox{\rm Ma}\approx 0.1$, this automatically implies that $k_{\rm
f}/k_{1}\geq 10$, provided that $\sigma_{\rm J}$ is not much larger than
unity. We focus on the case with $\sigma_{\rm J}=5c_{\rm s}k_{1}$, but we have
also experimented with smaller values of two and even 1.1. However, when
$\sigma_{\rm J}$ is that small, it takes a long time for the instability to
develop and by that time the initial turbulence would have decayed too much.

The magnetic field strength can also be specified in terms of $v_{\rm A}^{\rm
rms}=B_{\rm rms}/\sqrt{\mu_{0}\rho_{0}}$. In the second part of the paper, we
consider values of $v_{\rm A}^{\rm rms}/c_{\rm s}$ in the range 0.04–0.4. In
the first part of the paper, however, we are interested in the kinematic
regime and therefore consider values of around $10^{-18}$. In all cases, we
adopt a magnetic Prandtl number of unity, i.e., $\nu/\eta=1$, so the Reynolds
number is always equal to the magnetic Reynolds number.

###  2.4 Initial conditions

As initial conditions, we assume $\rho=\rho_{0}$, so there is no density
perturbation. However, we assume that the velocity and magnetic fields have a
random distribution with a $k^{4}$ spectrum below a given wavenumber $k_{\rm
f}$ and a $k^{-5/3}$ spectrum above $k_{\rm f}$. We assume the initial
velocity to be maximally helical at all wavenumbers, but take the magnetic
field to be nonhelical. This then leads to perturbations in the system that
trigger the Jeans instability. Again, with a few exceptions, we deliberately
choose a very weak initial magnetic field so as to see the possibility of a
kinematic dynamo at early times. A dynamo effect in decaying turbulence has
been found in an earlier study (Brandenburg et al., 2019), where one saw a
significant temporal increase of the magnetic field over several orders of
magnitude when the initial field is sufficiently weak, but no significant
increase was found for fields that start off in near-equipartition with the
turbulence.

Figure 1:  Mach numbers for (a) the rms velocity and (b) the rms Alfvén speed.
The red and black dashed lines in (b) have been scaled up by $10^{16}$ to make
them visible. The vertical dashed-dotted lines denote the times $t_{*}$ when
the Mach number has recovered to the original value of about 0.2. (For Run S,
the initial value of 2 was not recovered before the collapse.) The dotted
lines correspond to $\exp(\sigma t)$ with $\sigma=0.9\sqrt{24}$ and
$\sqrt{3}$, respectively. The upper abscissa gives time as $tc_{\rm s}k_{1}$
for the cases with $\sigma_{\rm J}/c_{\rm s}k_{1}=5$.

Figure 2:  $\omega_{z}$ (upper left), $B_{z}$ (upper right),
${\bm{\nabla}}\cdot\bm{u}$ (lower left), and $\ln\rho$ (lower right) near the
end of the run. Note the close correlation of the magnetic field with the
vorticity and their concentration toward regions of strong flow convergence
(${\rm div}\,{}\bm{u}<0$) and high density.

###  2.5 Numerical simulations

We use the Pencil Code (Pencil Code Collaboration et al., 2021), which employs
sixth order accurate derivatives in space and a third order accurate time
stepping scheme. Selfgravity was implemented by Johansen et al. (2007) for
modeling planetesimal formation and employs Fourier transformation. That same
module is also being used for studying dust formation in the ISM (Mattsson &
Hedvall, 2022).

Figure 3:  (a) $\mbox{\rm Re}(k)$ and (b) $\mbox{\rm Lu}(k)$ for Run B at six
different times, indicated by line types and color.  Table 1:  Energy flux
ratios and Reynolds numbers for the runs discussed in the paper. Here, $k_{\rm
J}$ and $k_{\rm f}$ are in units of $k_{1}$.

Run $k_{\rm J}$ $k_{\rm f}$ $\mbox{\rm Re}_{k_{\rm f}}$ $\mbox{\rm
Re}_{t_{*}}$ $\mbox{\rm Lu}_{t_{*}}$ $\epsilon_{\rm J}^{\rm P}$ $\epsilon_{\rm
J}^{\rm+K}$ $\epsilon_{\rm J}^{\rm-K}$ $\epsilon_{\rm J}^{\rm L}$
$\epsilon_{\rm L}^{\rm+M}$ $\epsilon_{\rm J}^{\rm-M}$ $\epsilon_{\rm
M}^{\Delta}$ $\gamma$ $N^{3}$ O1 2 10 500 1000 $2.3\times 10^{-13}$ 0.30 0.68
0.02 0.000 0.28 0.72 0.40 $~{}~{}1.39$ $2048^{3}$ O2 5 2 100 1000 $3.9\times
10^{-14}$ 0.32 0.67 0.01 0.000 0.49 0.51 0.96 $~{}~{}0.83$ $1024^{3}$ A 5 10
500 1000 $1.0\times 10^{-13}$ 0.33 0.64 0.03 0.000 0.35 0.65 0.70 $~{}~{}3.0$
$2048^{3}$ B 5 10 100 200 $9.9\times 10^{-15}$ 0.34 0.63 0.03 0.000 0.21 0.79
0.26 $~{}~{}0.44$ $2048^{3}$ b 5 10 100 200 $9.9\times 10^{-15}$ 0.31 0.66
0.03 0.000 0.20 0.76 0.32 $~{}~{}0.54$ $1024^{3}$ C 5 10 20 40 $9.6\times
10^{-16}$ 0.34 0.63 0.03 0.000 0.06 0.94 0.06 $~{}~{}0.08$ $2048^{3}$ D 5 10 5
10 $1.1\times 10^{-16}$ 0.31 0.66 0.03 0.000 $-0.82$ 1.82 $-0.48$ $-0.44$
$1024^{3}$ E 5 10 1 2 $5.6\times 10^{-18}$ 0.25 0.73 0.02 0.000 $-11.5$ 12.5
$-0.96$ $-1.13$ $1024^{3}$ S 5 10 500 1300 $8.7\times 10^{-14}$ 0.31 0.60 0.09
0.000 0.34 0.66 0.79 $~{}~{}2.47$ $1024^{3}$ M1 5 10 500 1000 $1.3\times
10^{3}$ 0.36 0.46 0.13 0.05 $-17.6$ 18.6 $-0.94$ $-0.38$ $2048^{3}$ M2 5 10
500 1000 $6.4\times 10^{2}$ 0.31 0.67 0.01 0.01 $-0.97$ 1.97 $-0.45$ $-0.16$
$2048^{3}$ M3 5 10 500 1000 $4.1\times 10^{2}$ 0.32 0.65 0.01 0.02 0.03 0.97
0.17 $~{}~{}0.11$ $2048^{3}$ M4 5 10 100 200 $9.8\times 10^{0}$ 0.33 0.65 0.02
0.00 0.20 0.80 0.25 $~{}~{}0.42$ $2048^{3}$ I1 5 10 500 1000 $1.6\times
10^{3}$ 0.29 0.68 0.01 0.02 0.67 0.33 2.63 $~{}~{}0.05$ $2048^{3}$ I2 5 10 500
1000 $7.5\times 10^{2}$ 0.31 0.67 0.01 0.01 0.20 0.80 0.26 $~{}~{}0.03$
$2048^{3}$ I3 5 10 500 1000 $4.3\times 10^{2}$ 0.32 0.65 0.01 0.02 0.32 0.68
0.60 $~{}~{}0.33$ $2048^{3}$

Many of the diagnostic quantities are calculated during run time, including
spectra and slices. Most of the secondary data that are used for the plots are
publicly available; see the code and data availability statement at the end of
the paper.

##  3 Results for weak magnetic fields

###  3.1 Visualizations and spectra

In Figure 1 we show the evolution of the Mach number and the rms Alfvén speed
normalized to the sound speed for $\sigma_{\rm J}=2c_{\rm s}k_{1}$ (Run O1)
and $5c_{\rm s}k_{1}$ (Runs A–E), so $t_{\rm ff}c_{\rm s}k_{1}=0.96$ and
$0.38$, respectively. In both cases, an exponential growth of Ma commences at
some time. For Run O1, the growth rate agrees with that expected from the
dispersion relation, i.e., $\sigma/c_{\rm s}k_{1}=\sqrt{3}$, but for Run B,
the actual value is 10% smaller than the theoretically expected value,
$\sigma/c_{\rm s}k_{1}=\sqrt{24}$, which could be related to the finite
viscosity. We define the moment when the rms velocity has recovered to its
initial value (denoted by the horizontal line for $u_{\rm rms}(0)/c_{\rm
s}=0.2$) as $t_{*}$. Those characteristic times ($t_{*}c_{\rm s}k_{1}\approx
1.5$ and 4.9 for $\sigma_{\rm J}/c_{\rm s}k_{1}=5$ and 2, respectively) are
denoted by vertical dashed-dotted lines in the corresponding colors. Those
times correspond approximately to the moment when the negative potential
energy density begins to exceed the kinetic energy density, i.e., when the
virial parameter $\alpha_{\rm vir}$ drops below two; see Appendix A for a
demonstration. In the supersonic case of Run S, the collapse is found to occur
earlier. This is mainly a consequence of the stronger initial perturbations
(see, e.g. Mac Low & Klessen, 2004, for a review).

The growth of the magnetic field is quantified by the ratio $v_{\rm A}/c_{\rm
s}$; see Appendix B for a comparison with other measures such as $B_{\rm rms}$
and the ratio $|\bm{B}|/\rho^{2/3}$. Given that the velocity is decaying
during the first part of the evolution, we cannot expect an exponential growth
of the magnetic field. During the second part, when the velocity is
exponentially increasing, the magnetic field does not show a corresponding
increase. In the following, we tentatively associate the slow growth of the
magnetic field during the decay phase of the velocity field with a dynamo, and
the second part, which is dominated by the Jeans instability, with just
compressional amplification. With these observations in mind, we continue
using the term dynamo, but leave it for further analysis to establish more
rigorous and convincing criteria.

Figure 4:  (a) Work terms $W_{\rm J}$ (red) and $-W_{\rm P}$ (orange), as well
as $-W_{\rm P}+\dot{\cal E}_{\rm K}$ (black dashed). Also shown is the kinetic
energy dissipation, $Q_{\rm K}$. (b) Work terms $-W_{\rm L}$ (red) and $Q_{\rm
M}$ (blue), as well as $Q_{\rm M}+\dot{\cal E}_{\rm M}$ (black dashed). We
recall that all work terms are in units of $\rho_{0}c_{\rm s}^{3}k_{1}$.

For the rest of the paper, we focus on the case $\sigma_{\rm J}=5c_{\rm
s}k_{1}$, so we expect growth for $k<k_{\rm J}\equiv 5k_{1}$. We summarize our
runs in Table LABEL:Tsummary. In Figure 2 we show a visualization of the $z$
components of the vorticity and magnetic field for Run B, as well as its flow
divergence and the logarithmic density near the end of the run at $tc_{\rm
s}k_{1}=2.02=5.24\,t_{\rm ff}$, shortly before further compression can no
longer be resolved. The fact that this time is much longer than unity is due
to the periodicity of the solution to the Poisson equation (Federrath &
Klessen, 2012; Lane et al., 2022). We see that $\omega_{z}$ and $B_{z}$ show
strong concentrations toward regions where the density also increases. Note
that the regions of negative flow divergence are more strongly concentrated
that those of positive flow divergence, but both ${\bm{\nabla}}\cdot\bm{u}$
and $\ln\rho$ are dominated by a spatially smooth component that, unlike
$\omega_{z}$ and $B_{z}$, lack small scales. Very weak small-scale
perturbations can be seen in the visualizations of ${\bm{\nabla}}\cdot\bm{u}$,
but not in those of $\ln\rho$. Gradients of $\ln\rho$ do, however, show marked
small-scale structure. In our runs, the density contours are more often
aligned with the magnetic field vector than being perpendicular to it. This is
demonstrated in Appendix B for Runs M1 and I1, where we see that $\bm{B}$ is
mostly perpendicular to ${\bm{\nabla}}\ln\rho$.

To see whether at any scale, the Reynolds number is large enough for dynamo
action, we employ the $k$-dependent Reynolds numbers. In Figure 3, we show
spectra of $\mbox{\rm Re}(k)$ and $\mbox{\rm Lu}(k)$ at different times. We
clearly see that the spectra display instability for $k<k_{\rm J}=5\,k_{1}$,
and only at late times ($tc_{\rm s}k_{1}>1.7$), somewhat smaller scales begin
to grow as well. However, the effect on the rest of the spectrum is
surprisingly weak. There is a small increase of $\mbox{\rm Lu}(k)$ at all
wavenumbers, but there is no visible effect from the Jeans instability itself,
except for the time close to the end of the simulation where one sees an
increase of both $\mbox{\rm Re}(k)$ and $\mbox{\rm Lu}(k)$ at the highest
wavenumbers, indicating that more energy is now being channeled through the
turbulent cascade. We also see that the value of $\mbox{\rm Re}(k)$ near the
wavenumber $k_{\rm f}$, where the initial kinetic energy spectrum peaks, is
around $\mbox{\rm Re}_{k_{\rm f}}=100$. According to previous studies, this
value is high enough for dynamo action. Our more detailed studies below
confirm that this is indeed the case. When changing $\nu$ and $\eta$, the peak
values of $\mbox{\rm Re}(k)$ and $\mbox{\rm Lu}(k)$ change correspondingly. In
Appendix C, we show that $\mbox{\rm Re}_{k_{\rm f}}=100$ is close (within a
factor of two) to the Taylor microscale Reynolds number.

###  3.2 Work terms

Next, we consider the evolution of the various work terms; see Figure 4. The
work done by the gravity force,
$-\langle\rho\bm{u}\cdot{\bm{\nabla}}\Phi\rangle\equiv W_{\rm J}>0$, leads to
flow compression, $\langle p{\bm{\nabla}}\cdot\bm{u}\rangle\equiv W_{\rm
P}<0$, and an increase of the kinetic energy density, so

$-\langle\rho\bm{u}\cdot{\bm{\nabla}}\Phi\rangle=-\langle
p{\bm{\nabla}}\cdot\bm{u}\rangle+\frac{{\rm d}{}{\cal E}_{\rm K}}{{\rm
d}{}t}...,$ (15)

where the ellipsis denotes the sum of two additional, subdominant terms:
$Q_{\rm K}$ and $-W_{\rm L}$. The dominant balance in the kinetic energy
evolution is given by the Jeans work, which is found to be balanced to 1/3 by
the pressure work and to 2/3 by the growth in kinetic energy, i.e.,

$-W_{\rm P}\approx\frac{1}{3}W_{\rm J},\quad\dot{\cal E}_{\rm
K}\approx\frac{2}{3}W_{\rm J}.$ (16)

The latter can be integrated to give ${\cal E}_{\rm K}\approx(2/3)\int W_{\rm
J}\,{\rm d}{}t$. Likewise, integrating Eq. (5) gives $-{\cal E}_{\rm
P}\approx\int W_{\rm J}\,{\rm d}{}t$, which implies $\alpha_{\rm vir}=2{\cal
E}_{\rm K}/|{\cal E}_{\rm P}|\approx 4/3$. Its value would be unity, if only
half of $W_{\rm J}$ went into the growth of kinetic energy, but this is not
the case. It is important to realize that the energy flux ratios in Eq. (16)
apply to the time $t=t_{*}$. They change at later times toward 1/4 and 3/4 for
the pressure work and growth in kinetic energy; see Appendix A, which implies
$\alpha_{\rm vir}\approx 3/2$.

These ratios of the work terms implies that about one-third of the
gravitational energy goes into compressional heating and two-thirds go into
kinetic energy before eventually also being thermalized. At the reference time
$t_{*}=1.5/(c_{\rm s}k_{1})$, however, viscous dissipation contributes only
about 3%; see Runs A–E in Table LABEL:Tsummary.

In the kinematic regime of the dynamo, the work done by the Lorentz force is
very small; see Figure 4(b). Nevertheless, this term exceeds the Joule
dissipation when Re is large enough. Here, the work done against the Lorentz
force, $-W_{\rm L}$, leads to Joule heating and an increase in the magnetic
energy density, i.e.,

$-\langle\bm{u}\cdot(\bm{J}\times\bm{B})\rangle=\langle\mu_{0}\eta\bm{J}^{2}\rangle+\frac{{\rm
d}{}{\cal E}_{\rm M}}{{\rm d}{}t}.$ (17)

Therefore, based on the positivity of $(-W_{\rm L})-Q_{\rm M}$, i.e., the
positivity of $\gamma$, we argue that we can determine the threshold for
dynamo action to be around 25. This is the value for which the interpolated
line of $(-W_{\rm L}/Q_{\rm M})-1$ vs Re crosses zero, which corresponds to
marginal dynamo excitation; Figure 5. In fact, Run C with $\mbox{\rm
Re}_{k_{\rm f}}=20$ is close to the marginal point; see Table LABEL:Tsummary.

During the exponential growth phase of the Jeans instability, the velocity
grows at the rate $\sqrt{24}c_{\rm s}k_{1}$. During that period, $\mbox{\rm
Re}(t)$ increases rapidly with time, and so does also the difference $(-W_{\rm
L})-Q_{\rm M}$.

Figure 5:  Time dependence of $(-W_{\rm L}/Q_{\rm M})-1$ on $\mbox{\rm
Re}_{t}$ for Runs A (red), B (orange), C (green), D (blue), and E (black). The
values of $(-W_{\rm L}/Q_{\rm M})-1$ from Table LABEL:Tsummary are shown as
diamonds as a function of $\mbox{\rm Re}_{k_{\rm f}}$ and are seen to obey an
approximate fit given by $\ln(0.45\,\mbox{\rm Re}^{1/4})$; see the dashed-
dotted line.  Figure 6:  Characteristic wavenumbers $k_{\omega}$ (blue),
$k_{\bm{\omega}\cdot\bm{u}}$ (red), $k_{\,{\bm{\nabla}}\cdot\bm{u}}$ (solid
black), and $k_{p{\bm{\nabla}}\cdot\bm{u}}$ (dotted black).

Figure 7:  $B_{z}$ for Runs M2 (left) and M1 (right) at $tc_{\rm s}k_{1}=1.6$.
Note the large-scale concentrations (marked by white ellipses) and voids
(marked by black circles) for stronger fields, even though the stage of the
collapse is the same. For Run M1, one sees indications that the $2048^{3}$
resolution begins to become insufficient at $t=1.6$.

###  3.3 Dependence on the magnetic Reynolds number

As the magnetic Reynolds number increases, we expect a dynamo to become
stronger and thus, $\epsilon_{\rm M}^{\Delta}\equiv(-W_{\rm L}-Q_{\rm
M})/Q_{\rm M}$ to increase. However, the magnetic Reynolds number is time-
dependent, because $u_{\rm rms}$ increases. This raises the question whether
this dependence follows a similar trend that is seen by comparing different
runs.

In Figure 5 we plot $\epsilon_{\rm M}^{\Delta}(t)$ versus $\mbox{\rm Re}_{t}$
and compare with the values listed in Table LABEL:Tsummary for the nominal
time $t_{*}=1.5/(c_{\rm s}k_{1})$ versus $\mbox{\rm Re}_{k_{\rm f}}$. We see
that $\epsilon_{\rm M}^{\Delta}(t_{*})$ shows a rather shallow increase with
$\mbox{\rm Re}_{\rm M}(t_{*})$ of the form

$\epsilon_{\rm M}^{\Delta}\approx\ln(0.45\,\mbox{\rm Re}^{1/4}).$ (18)

The time-dependent tracks do approximately match this dependence at around
intermediate times, but all the lines are curved and shallower for early times
and steeper at late times, where the Jeans instability becomes dominant.

###  3.4 Compression and vorticity

The magnetic field growth seems to be strongly correlated with the rms
vorticity. Earlier studies showed that it is only the vortical part of the
flow that leads to dynamo action (Mee & Brandenburg, 2006). To quantify the
relative importance of vortical and irrotational or compressive contributions
to the velocity field, we show in Figure 6 the evolution of the characteristic
wavenumbers $k_{\omega}$, $k_{\bm{\omega}\cdot\bm{u}}$,
$k_{\,{\bm{\nabla}}\cdot\bm{u}}$, and $k_{p{\bm{\nabla}}\cdot\bm{u}}$.

We see that $k_{\omega}\approx 35$ and $k_{\bm{\omega}\cdot\bm{u}}\approx 20$
during the early phase when the collapse velocity is still subdominant. When
the collapse becomes dominant, $k_{\omega}$ and $k_{\bm{\omega}\cdot\bm{u}}$
rapidly decline and $k_{\,{\bm{\nabla}}\cdot\bm{u}}\approx 2.5$ prior to the
final collapse. Nevertheless, we always find
$k_{\omega}>k_{\,{\bm{\nabla}}\cdot\bm{u}}$, i.e., vorticity is still
important. This is partially because of the compression of vortices, as seen
in Figure 2, which enhances the vorticity. Furthermore, we find that
$k_{p{\bm{\nabla}}\cdot\bm{u}}\ll k_{\,{\bm{\nabla}}\cdot\bm{u}}$, except very
close to the collapse when the two are similar.

###  3.5 Dependence on resolution

For our reference models, we use a fairly high resolution of $2048^{3}$ mesh
points. Larger resolutions become easily impracticable and reach computational
memory limitations. However, the dependence on resolution is small and the
maximum run time before the collapse stops the calculation hardly changes at
all when increasing the resolution by a factor of two. By comparing Runs b and
B, we see that the values of $W_{\rm P}/W_{\rm J}$ and $\dot{\cal E}_{\rm
K}/W_{\rm J}$, as well as those of $\dot{\cal E}_{\rm M}/(-W_{\rm L})$ and
$Q_{\rm M}/(-W_{\rm L})$ agree with each other within 10%.

Figure 8:  Work terms $-W_{\rm L}$ (solid red) and $Q_{\rm M}$ (solid blue),
together with $W_{\rm J}$ (dashed red) and $-W_{\rm P}$ (dashed orange), as
well as $Q_{\rm K}$ (dashed blue) for (a) Run M3 with a relatively weak
magnetic field and (b) Run M1 with a strong magnetic field. Dotted lines
denote the slopes corresponding to the growth rates $1.6\,\sigma_{\rm J}$ and
$3.2\,\sigma_{\rm J}$ in panel (a) and $1.5\,\sigma_{\rm J}$ in panel (b). The
thick dashed-dotted line in panel (b) denotes $W_{\rm L}>0$. Note that the
work terms are in units of $\rho_{0}c_{\rm s}^{3}k_{1}$.

##  4 Results for strong magnetic fields

###  4.1 Earlier collapse with stronger fields

We now consider cases where the magnetic field is dynamically important. This
situation is of particular interest for dense pre-stellar cores, where the
measured (e.g., Crutcher, 2012) or inferred (e.g., Karoly et al., 2020; Pattle
et al., 2021) magnetic fields are particularly strong.

Figure 7 shows a comparison of the magnetic field patterns for Runs M3 and M1,
i.e., for weaker and stronger fields, respectively. For the stronger magnetic
fields in Run M1, the magnetic eddies appear to be organized in larger patches
that correspond to over- or underdense regions. For the stronger magnetic
fields in Run M1, the magnetic eddies appear earlier than in Run M3, which is
probably the result of the magnetically driven motion early on; see Figure 1.
This behavior is suggestive of an accelerated collapse process. This is an
important difference to the standard paradigm of magnetically controlled star
formation that employs a uniform magnetic field (Mestel & Spitzer, 1956;
Mouschovias & Spitzer, 1976; Shu, 1977). Instead, here the magnetic field is
turbulent and only has moderate large-scale coherence.

###  4.2 When $W_{\rm L}$ affects the collapse

In Figures 8(a) and (b) we show the work terms for Runs M3 and M1,
respectively. Again, we see that a turbulent magnetic field does not
systematically delay the collapse, and a strong field can even accelerate it.

Figure 9:  Comparison of kinetic, magnetic, and potential energy and density
spectra at $tc_{\rm s}k_{1}=1.7$ for (a) Run M1 with a strong turbulent
initial magnetic field and (b) Run I1 with a strong uniform initial magnetic
field. The dashed lines denoted with $E_{\bm{\omega}}(k)/k^{2}$ correspond to
kinetic energy spectra of the vortical part of the velocity. Here, the $E(k)$
are in units of $c_{\rm s}^{2}k_{1}^{-1}$. Note that both axis ranges in (a)
and (b) are the same.

We also see that the late-time exponential increase of $-W_{\rm L}$ and
$Q_{\rm M}$ changes with respect to the weak-field behavior from being twice
the rate of $W_{\rm J}$ to being equal to it; see Figure 8(b). Assuming the
change in $\bm{J}\times\bm{B}$ to be itself proportional to the change in
$\bm{u}$, we see that $-W_{\rm L}$ is quadratic in $\bm{u}$, which explains
the growth of $-W_{\rm L}$ at twice to Jeans rate. Once the field is strong,
$\bm{J}\times\bm{B}$ does not change much anymore, and so $-W_{\rm L}$ only
grows at the Jeans rate.

The gravitational collapse is primarily characterized by the increase in the
velocity $\bm{u}$ and not so much by a change in gravity. In reality, however,
the change in gravity does contribute about 50%, and so we find that $W_{\rm
J}$ increases at a rate of about $1.6\,\sigma_{\rm J}$, and $-W_{\rm L}$
increases at a rate of about $3.2\,\sigma_{\rm J}$.

Figure 10:  Comparison of relative helicity spectra with (a) turbulent and (b)
uniform initial magnetic fields. Blue (red) lines denote kinetic (magnetic)
relative helicity spectra. The solid lines are for Run B (I3), the dotted
lines are for Run D (I2), and the thick dashed lines are for Run E (I1) for
turbulent (uniform) initial magnetic fields.

###  4.3 Comparison with uniform initial magnetic fields

In Figure 9(a) we show kinetic and magnetic energy spectra together with
logarithmic density spectra and the normalized enstrophy spectra representing
the kinetic energy spectra of the vortical part, and compare with the case of
a uniform initial magnetic field in Figure 9(b). The two panels correspond to
Runs M1 and I1. Even though $\mbox{\rm Lu}_{t_{*}}$ is 25% higher in I1, the
spectral energies are lower. The density spectra show a rapid increase for
$k\leq k_{\rm J}$, which is associated with the Jeans instability, as was
already seen in the kinetic energy spectra; see Figure 3(a).

The magnetic energy spectra and the kinetic energy spectra of the vortical
part do not show the same increase, but there is a slight one, which is
different for the cases with a turbulent and an initially uniform field. To
understand this, it is useful to discuss now the kinetic and magnetic helicity
spectra for the two cases. We see that, for runs with an initially uniform
magnetic field at intermediate wavenumbers, the magnetic energy is either in
nearly perfect equipartition, or in super-equipartition with the kinetic
energy. However, it becomes subdominant at small wavenumbers, where the
behavior is affected by gravitational collapse. Therefore, the spectrum
resembles that of a small-scale dynamo, where the magnetic field is also in
super-equipartition at large $k$; see Haugen et al. (2003, 2004a). However,
this similarity should not be regarded in any way as evidence in favor of a
dynamo. It appears to be instead just a typical behavior of any type of
hydromagnetic turbulence.

Returning to the slight uprise of the magnetic field for $1\leq k/k_{1}<5$ in
the case of an initially uniform magnetic field, we argue that this is caused
by the tangling of the magnetic field by the collapsing gas motions. It is
therefore not due to an inverse cascade, which usually occurs only in the
absence of a mean magnetic flux through the domain. In the case of a turbulent
magnetic field, the build-up of vorticity at small wavenumbers could be caused
by the shear flows, which leads to what is known as a vorticity dynamo
(Elperin et al., 2003). In the case of an initially uniform magnetic field,
this vorticity dynamo is suppressed (Käpylä et al., 2009). However, the uprise
of the vortical part of the velocity field for $1\leq k/k_{1}<5$ appears to be
caused by the magnetic field and becomes weaker for Runs M2 and M3.

###  4.4 Helicity spectra

It is important to realize that, owing to the use of periodic boundary
conditions, an initially uniform magnetic field is equivalent to what is
sometimes described as an imposed magnetic field. This is simply because the
mean magnetic flux is preserved. A well-known difference between cases with
and without an imposed magnetic field is the fact that
$\langle\bm{A}\cdot\bm{B}\rangle$ is no longer conserved in the former case
(Berger, 1997). This is because now the magnetic helicity in the domain
interacts with the magnetic helicity on scales larger than the size of the
periodic domain, but this part is no longer included in the simulation; see
the discussion in Brandenburg & Matthaeus (2004). In simulations of decaying
turbulence, it has been found that the magnetic fluctuations decay more
rapidly when there is an imposed magnetic field (Brandenburg et al., 2020).

Figure 10 shows kinetic and magnetic helicity spectra for the cases with a
turbulent and an imposed field for runs with different magnetic field
strengths. The kinetic helicity spectra are similar in the two cases, but the
magnetic helicity spectra are not. In the case of an imposed magnetic field,
there is magnetic helicity of the same sign at all wavenumbers, although it is
less strong at small wavenumbers. By contrast, in the case with a turbulent
magnetic field with zero net flux, the magnetic helicity is predominantly of
opposite (negative) sign and relatively strong also at small wavenumbers,
except in the case with the strongest magnetic field (Run M1). This is caused
by magnetic helicity conservation, where a small-scale driving of magnetic
helicity of one sign causes automatically the appearance of magnetic helicity
of opposite sign at large scales; see also Brandenburg et al. (2019) for
similar results. The small-scale magnetic helicity does get slowly dissipated
at late times through finite microphysical magnetic diffusivity, leaving
predominantly the large-scale magnetic helicity of opposite sign behind.

###  4.5 Contributions to the Lorentz work

The work done against the Lorentz force serves as one of our main tools to
quantify dynamo action. As we have mentioned at the end of Sect. 2.2.1, the
work done against the Lorentz force can be subdivided into contributions from
the magnetic pressure gradient, the tension force, and the curvature force. In
Figure 11, we show that, for weak initial magnetic fields, the most important
contribution to the pseudo growth rate comes from the work done against the
curvature force, but later during the collapse, a more important contribution
comes from the compressional work done against the magnetic pressure gradient.

Figure 11:  Evolution of the pseudo growth rate $\gamma=\gamma_{\rm
c}+\gamma_{\|}+\gamma_{\perp}$ (solid lines), with contributions from the work
done against the curvature force ($\gamma_{\perp}$, red dashed lines), the
tension force ($\gamma_{\|}$, black dashed-triple-dotted lines), and the
magnetic pressure gradient ($\gamma_{\rm c}$, blue dotted lines) for (a) Run
B, (b) Run O2, (c) Run M1, and (d) Run I1. In panels (a) and (b),
$\gamma_{\|}=0$. In panels (c) and (d), $\gamma_{\|}\neq 0$, and the zero line
has been drawn as a straight black line. The vertical dashed-dotted lines
denote the critical time $t_{*}$ when the Mach number has recovered to the
original value of about 0.2.  Figure 12:  Evolution of the pseudo growth rate
$\gamma$ (black lines), with contributions from $\gamma_{\rm 2D}$ (blue lines)
and the residual $\gamma-\gamma_{\rm 2D}$ (red lines), for (a) Run B and (b)
Run M1.

In the runs with a strong magnetic field (turbulent or imposed), the value of
$\gamma$ is still negative at the critical time $t_{*}$, but the compressional
work done against the magnetic pressure gradient is positive, and it was
positive also during the earlier phase, especially in the case of a turbulent
magnetic field; see Figure 11(c). The contribution of $\gamma_{\rm c}$ is a
characteristic feature of amplification or at least sustenance of the magnetic
field in collapsing turbulence through compression.

In Figure 12 we plot the time dependences of $\gamma$, $\gamma_{\rm 2D}$, and
$\gamma_{\rm 3D}=\gamma-\gamma_{\rm 2D}$. We see that $\gamma_{\rm 2D}$ is
always close to zero, except during an early phase which can be associated
with 2-D tangling of the initial magnetic field. When $\gamma_{\rm 3D}$ is
included, the resulting pseudo growth rate is positive during much of the
early part of the evolution.

Based on the positive values of $\gamma_{\perp}$ and $\gamma_{\rm 3D}$ in the
cases of weak magnetic fields in Figures 11 and 12, we are led to suggest that
those runs do indeed host supercritical dynamos. When the magnetic field is
strong, however, $\gamma_{\perp}$ and $\gamma_{\rm 3D}$ are now negative,
suggesting that Run M1 cannot be classified as a dynamo.

For strong magnetic fields, only near the end of the collapse does
$\gamma_{\perp}$ become positive. On the other hand, $\gamma_{\rm 3D}$ becomes
then the dominant term during the collapse and $\gamma_{\rm 2D}$ becomes
negative; see Figure 12(b). This is probably caused by the strong alteration
of the flow by the magnetic field, making now $\gamma_{\rm 2D}$ strongly
negative. This increases the compression and tangling terms associated with
$\gamma_{\rm 2D}$, which then contribute to enhancing the kinetic energy
rather than the other way around (as in a dynamo). Nevertheless, since the
magnetic field is now increasing, $\gamma_{\rm 3D}$ becomes positive. From
Figures 11(c) and (d) we know, however, that this increase is manifestly
caused by compression, which must therefore be a 3-D compression, and not a
dynamo effect.

In summary, we are led to conclude that there is dynamo action in all cases
with weak magnetic fields prior to collapse, but probably no longer or not
very much during the actual collapse. When the magnetic field is already
strong, there is no longer dynamo action, but just 3-D compression. For large
Mach numbers, at late stages of the collapse, shocks form and cross each
other, which causes vorticity production (Porter et al., 2015), resulting then
also in dynamo action (Federrath et al., 2011a). However, this happens at such
small scales and such late times that this effect cannot be captured in direct
numerical simulations, even at a resolution of $2048^{3}$ mesh points. After
$tc_{\rm s}k_{1}=1.98$, the Jeans length is no longer resolved with 30 grid
cells – the minimum proposed by Federrath et al. (2011b); see Appendix D,
where we show the evolution of the maximum density during the collapse.

Figure 13:  Density spectra at $t=0.2$, 0.6, 1, and 1.4, shifted upward by an
$\exp(10t)$ factor to be able to see the shift of the oscillations in the
spectra in $k$. The dotted orange lines denote approximate fits of the form
given by Eq. (19).

###  4.6 Waves in the density spectra

Before concluding, let us comment on an interesting feature that we noticed in
spectra of the logarithmic density. At early and intermediate times, we see a
wavy structure in $E_{\ln\rho}(k)$; see Figure 13. In fact, this wavy
modulation is of the form $\cos k\xi(t)$, where $\xi(t)=c_{\rm s}t$ is the
distance a sound wave has propagated in the time $t$ since the initial
condition was applied. As time goes on, and as $\xi(t)$ therefore increases,
the waves appear to propagate toward smaller values of $k$ and are of
progressively shorter length in $k$ space. The changing phase of these waves
can be described by the formula

$E_{\ln\rho}(k,t)=E_{\rho}^{(0)}(k)\,\left[1+g(k,t)\,(1-\cos kc_{\rm
s}t)\right],$ (19)

where $E_{\rho}^{(0)}(k)\propto k^{-5/3}$ denotes the unmodulated spectrum,
and $g(t)=(140/k)\exp(-0.8\,t)$ is an empirically defined function (in code
units). Only at later times, the fit is going somewhat out of phase. Given the
agreement of our hypothetical modulation of the form $\cos k\xi(t)$ with the
actual spectrum, we can argue that the wavy structure is indeed caused by the
initial velocity perturbation launching sound waves from multiple locations in
the domain all at the same time, and that their characteristic scale increases
with time like $\xi(t)=c_{\rm s}t$. Similar waves have also been seen in
simulations of gravitational waves that are being initiated from an
instantaneous perturbation; see Figure 2 of Roper Pol et al. (2020). No
explanation for this phenomenon was offered there, but we have now confirmed
that it can be explained in a similar way, except that the relevant speed in
the expression for $\xi(t)$ is the speed of light.

##  5 Conclusions

In the present work, we have used the instantaneous excess of the work done by
the Lorentz force over the Joule dissipation as a quantity that characterizes
the dynamo. Under stationary conditions, the dynamo can easily be
characterized by the growth rate. However, a growth rate cannot be defined in
situations when the velocity itself decays or grows exponentially with time,
like we observe in our models.

The dynamo criterion based on the work terms results in a value of the
critical Reynolds number of about 25, which is smaller than the critical value
of about 35 for small-scale dynamo action (Haugen et al., 2004a), but larger
than the value for large-scale dynamo action in the presence of helicity of
below six (Brandenburg, 2009). Also, in the present case there is helicity, so
we do expect a critical value that is less than 35. However, there is a strong
contribution from irrotational motions that makes the dynamo harder to excite
and does itself not contribute to dynamo action (Mee & Brandenburg, 2006).
During the collapse, i.e., after $t=t_{*}$, and for weak magnetic fields (Runs
B and O2), $\gamma_{\perp}$ and $\gamma_{\rm 3D}$ show a slight decrease when
the magnetic field is weak, supporting the idea that this magnetic field
growth is not primarily caused by dynamo action, but just by compression.

Our investigation shows that the most important contribution to the growth of
a magnetic field comes from the work done against the curvature force,
although later during the collapse, there is an even more important
contribution from the compressional work done against the magnetic pressure
gradient. However, as we have argued above, this type of magnetic field
amplification happens also in one or two dimensions and should therefore not
be associated with dynamo action. By considering the decomposition into
$\gamma_{\rm 2D}$ and $\gamma_{\rm 3D}$ we have made an attempt of
distinguishing dynamo action from the type of non-dynamo amplification seen
also in two dimensions. Nevertheless, our dynamo criterion is not very
precise, as the pseudo growth rate changes behavior with different initial
conditions and a number of factors need to be considered in combination.

We stated in the introduction that the exact fraction of potential energy that
goes into turbulence is unknown. Our results now show that one-third of the
energy input from potential energy goes into compressional heating, and two-
third goes into the kinetic and magnetic energies of the turbulence. Thus, one
would expect that, at the end of the collapse, the sum of kinetic and magnetic
energy densities is twice the thermal energy density from compressive heating.
This is different from the virial theorem, which relates potential and kinetic
energies to each other. As explained in Sect. 3.2, however, since the two-
third contribution to the kinetic energy comes from potential energy, which
becomes more negative with time, it follows that the ratio of kinetic to
potential energy is 2/3 and thus, the virial parameter is 4/3. It would be
unity, if the contribution to the kinetic energy was half the Jeans work. At
later times, however, the fractional kinetic energy gain increases toward 3/4
of the Jeans work, which implies a virial parameter of about 3/2.

In all the simulations presented here, we have used an isothermal equation of
state. However, deviations from isothermality probably play an important role
during molecular cloud collapse. For example, Lee & Hennebelle (2018) caution
that, in simulations studying gravitational fragmentation of molecular clouds,
an isothermal equation of state cannot lead to converged results with
increasing numerical resolution. They propose that an adiabatic equation of
state at high densities, essentially accounting for the formation of the
Larson core is a more physically meaningful approach. In future work, it would
be interesting to perform simulations using an ideal gas equation of state
instead of an isothermal one. Such simulations could also allow for cooling,
which would further increase the density in regions of strong flow convergence
and counteract an otherwise singular collapse.

It would also be interesting to apply our analyses to earlier work that uses
Bonnor-Ebert spheres as initial conditions (Sur et al., 2010, 2012; Federrath
et al., 2011b). Bonnor-Ebert spheres are nonuniform equilibria in a non-
expanding finite space, so there is no $\rho_{0}$ term in Eq. (1). Our
simulations never go through such a state or reach any near-equilibrium state.
Starting with a Bonnor-Ebert sphere could lead to the collapse proceeding in a
different way from ours. If the collapse occurs sufficiently slowly, it could
be easier to achieve dynamo growth rates that remain faster than the collapse
rate for a longer time.

Another useful extension would be to compare with simulations that make use of
adaptive mesh refinement (see, for example, Federrath et al., 2010). Such
simulations would have varying accuracy in space, and it is currently unclear
how this affects the kinetic and magnetic energy spectra and other
diagnostics. Since the varying accuracy is not a concern in the present
simulations, they can be used as a benchmark. Another advantage of the present
simulations is the fact that the viscosity and magnetic diffusivity are fixed
and that therefore numerically converged and accurate results are possible.

## Acknowledgements

We thank the referee for useful comments and references. This work was
supported in part through the Swedish Research Council, grant 2019-04234. EN
acknowledges funding by the ERC Grant “Interstellar” (Grant agreement 740120)
and the Hellenic Foundation for Research and Innovation (Project number 224).
Nordita is supported in part by Nordforsk. We acknowledge the allocation of
computing resources provided by the Swedish National Infrastructure for
Computing (SNIC) at the PDC Center for High Performance Computing Stockholm
and the National Supercomputer Centre (NSC) at Linköping.

Code and data availability. The source code used for the simulations of this
study, the Pencil Code (Pencil Code Collaboration et al., 2021), is freely
available on https://github.com/pencil-code/. The DOI of the code is
https://doi.org/10.5281/zenodo.2315093. The simulation setup and the
corresponding data are freely available on
https://doi.org/10.5281/zenodo.5760126; see also
https://www.nordita.org/~brandenb/projects/SelfGrav for easier access to the
same material as on the Zenodo site.

##  Appendix A Virial parameter

In Sect. 2.2.1, we noted that $\alpha_{\rm vir}=2{\cal E}_{\rm K}/|{\cal
E}_{\rm P}|$ is expected to be around unity, but that its value can be
different at large Mach numbers and for strong magnetic fields. We have also
stated that a value of 4/3 is expected if 2/3 of the Jeans work goes into
building up kinetic energy. The purpose of this appendix is now to compare the
evolution of $\alpha_{\rm vir}$ for runs with strong magnetic field (Run M1),
larger Mach number (Run S), with a subsonic run with weak magnetic field (Run
B).

Initially, when the density is uniform, the potential energy density density
${\cal E}_{\rm P}=-\langle({\bm{\nabla}}\Phi)^{2}\rangle/8\pi G$ is
small,111Using the identity
${\bm{\nabla}}\cdot(\Phi{\bm{\nabla}}\Phi)=({\bm{\nabla}}\Phi)^{2}+\Phi\nabla^{2}\Phi$,
the potential energy density can also be written as ${\cal E}_{\rm
P}=-\langle\rho\Phi\rangle/2$, where we have made use of periodicity and the
fact that $\langle\Phi\rangle=0$. This yields $\dot{\cal E}_{\rm
P}=-\langle{\bm{\nabla}}\Phi\cdot{\bm{\nabla}}\dot{\Phi}\rangle/4\pi
G=-\langle\dot{\rho}\Phi\rangle$, which leads to Eq. (5) after using Eq. (3).
so $\alpha_{\rm vir}$ is large. However, when the collapse has started, deep
potential wells develop and $|{\cal E}_{\rm P}|$ increases with ${\cal E}_{\rm
P}<0$, so $\alpha_{\rm vir}$ drops and eventually settles at a value of around
$1.5$ for weak magnetic fields (Run B), and perhaps also for the supersonic
run before the collapse occurred (Run S); see Figure 14(a). This larger value
$\alpha_{\rm vir}\approx 1.5$ is caused primarily by the fact that the
fractional kinetic energy gain from the Jeans work in Eq. (16) increases with
time from 2/3 to 3/4; see Figure 14(b). At the same time, the fractional
pressure work decreases correspondingly from 1/3 to 1/4; see Figure 14(c). For
strong magnetic fields (Run M1), $\alpha_{\rm vir}$ continues to decrease
below unity; see Figure 14.

Figure 14:  Time evolution of (a) the virial parameter, (b) the fractional
kinetic energy gain, and (c) the fractional pressure work for Runs B, M1, and
S. In (a), the inset is a logarithmic representation of $\alpha_{\rm vir}$
over a larger range. At the end of Runs S and M1, the results are affected by
insufficient resolution and cannot be trusted.

##  Appendix B Growth of and relation between magnetic field and density

Figure 15:  Different scalings of the magnetic field strength as a function of
time (scaled up by $10^{16}$ for Runs B and S). Note that the units of
$\langle\bm{B}^{2}/\rho^{2n}\rangle^{1/2}$ are $c_{\rm
s}\mu_{0}^{1/2}\rho_{0}^{1/2-n}$.  Figure 16:  Logarithm of 2-D histograms
$P(\ln\rho,\ln|\bm{B}|)$ for Runs M1 and I1. The solid and dashed lines
correspond to $|\bm{B}|\sim\rho^{2/3}$ and $|\bm{B}|\sim\rho^{1/2}$,
respectively.  Figure 17:  Histograms
$P(\bm{B}\cdot{\bm{\nabla}}\ln\rho/|\bm{B}||{\bm{\nabla}}\ln\rho|)$ for Runs
M1 and I1, for different times.

The purpose of this appendix is to assess the role of the density in
determining a relevant measure of the magnetic field in collapsing turbulence.
During radial collapse, a uniform magnetic field is amplified such that the
ratio $|\bm{B}|/\rho^{2/3}$ is constant (Pudritz et al., 2006). Therefore, it
is customary to monitor the evolution of this quantity (Sur et al., 2010,
2012; Federrath et al., 2011b; Sharda et al., 2021). As a suitable volume
average, one can consider
$|\bm{B}|/\rho^{2/3}\sim\langle\bm{B}^{2}/\rho^{4/3}\rangle^{1/2}$. In Figure
15, we show that for Runs B and M1, the quantities
$\langle\bm{B}^{2}/\rho^{2n}\rangle^{1/2}/\rho_{0}^{n}$ are nearly the same
for different exponents $n$. For Run S, the differences for different $n$ are
somewhat larger, but the differences between the cases $n=1/2$ and $n=2/3$ are
still negligible. This justifies the use of $v_{\rm A}^{\rm rms}$ in Figure 1
of the main text, which corresponds to $n=1/2$ in Figure 15. In fact, also the
case $n=1/2$ has been discussed previously in the context of gravitational
collapse (Crutcher, 1999). In Figure 16, we present the logarithm of 2-D
histograms $P(\ln\rho,\ln|\bm{B}|)$ to show that most of the points in the
volume lie within elliptical islands stretched along the line
$|\bm{B}|\sim\rho^{2/3}$. They are normalized such that $\int
P(\ln\rho,\ln|\bm{B}|)\,{\rm d}{}\ln\rho\,{\rm d}{}\ln|\bm{B}|=1$. This
behavior has been observed in numerous simulations of self-gravitating
turbulence (e.g. Soler et al., 2013; Chen et al., 2016; Barreto-Mota et al.,
2021), as a result of the converging motions driven by gravity (Soler &
Hennebelle, 2017).

In Sect. 3.1, we mentioned that $\ln\rho$ lacks small-scale structure, and
that only ${\bm{\nabla}}\ln\rho$ displays noticeable small-scale variations.
In Figure 17 we present histograms of the cosine of the angle between $\bm{B}$
and the logarithmic density gradients,
$P(\bm{B}\cdot{\bm{\nabla}}\ln\rho/|\bm{B}||{\bm{\nabla}}\ln\rho|)$, for Runs
M1 and I1 at different times. They show that $\bm{B}$ is mostly perpendicular
to ${\bm{\nabla}}\ln\rho$, i.e., the magnetic field lines tend to be aligned
with the contours of $\ln\rho$.

##  Appendix C Comparison with the Taylor microscale Reynolds number

In Figure 18, we compare the evolution of the Taylor microscale Reynolds
number, $\mbox{\rm Re}_{\lambda}$, with other Reynolds numbers: $\mbox{\rm
Re}_{k_{\rm f}}$ and $\mbox{\rm Re}_{t}$, whose values at $t=t_{*}$ are given
in Table LABEL:Tsummary for $k_{\rm f}/k_{1}=10$, and the maximum of
$\mbox{\rm Re}_{k}$ over $k$, $\max_{k}(\mbox{\rm Re}_{k})$, which is at later
times dominated by the peak at small $k$; see Figure 9, where the peak is at
$k/k_{1}\approx 2$. Note that $\max_{k}(\mbox{\rm Re}_{k})$ begins to grow
exponentially shortly after $tc_{\rm s}k_{1}\approx 1.1$.

Figure 18:  Comparison of the evolution of the Taylor microscale Reynolds
number with other Reynolds numbers for Run B.

In Sect. 5, we compared our critical magnetic Reynolds number of 25 with the
earlier value of 35 by Haugen et al. (2004b), and we argued that the new value
is smaller because of helicity in the flow. We also stated that the results
depend on the Mach number (Haugen et al., 2004b). We can now compare with the
highly compressible, supersonic hydromagnetic turbulence simulations of
Federrath et al. (2014), who gave a critical magnetic Reynolds number of
$\approx 130$. This value is based on half the size of the computational
domain. To convert it to the normalization of Haugen et al. (2004a, b), it
should be divided by $2\pi$. This results in a value of $\approx 20$, which is
smaller than the values of 35 to 70 found by Haugen et al. (2004b) for Mach
numbers below and above unity, respectively. However, when they increased
their magnetic Prandtl number from unity to five, they found $\mbox{\rm
Re}_{\rm M}^{\rm crit}=25$ and $50$ for Mach numbers below and above unity,
respectively; compare Figs. 7 and 8 of Haugen et al. (2004b). Federrath et al.
(2014) used $\mbox{\rm Pr}_{\rm M}=10$, so their value of $\mbox{\rm Re}_{\rm
M}^{\rm crit}=20$ cannot directly be compared with those of Haugen et al.
(2004b) for $\mbox{\rm Pr}_{\rm M}=5$, but it seems a bit low. This question
cannot be fully clarified here and might depend on subtle numerical aspects.
The FLASH code (Fryxell et al., 2000) used by Federrath et al. (2014) is based
on a Riemann solver, which may affect the effective viscosity and magnetic
diffusivity.

##  Appendix D Maximum density and other density moments

Figure 19:  Evolution of the maximum density along with the minimum density
and intermediate moments of the density for $n=2$, $4$, and $12$. Note that
time increases toward the left. The time $tc_{\rm s}k_{1}=1.98$ for which the
Jeans length becomes resolved by less than 30 mesh points is marked by the
black symbol. This is when $(\rho_{\max}/\rho_{0})^{1/2}=410/30$.

As time goes on and the collapse proceeds, the density contrast increases.
This has implications for the nominal Jeans length, which decreases with
increasing maximum density. To get an idea of this, we plot in Figure 19 the
evolution of the maximum density. For completeness, we also plot the minimum
density and intermediate moments of the density,
$\langle\rho^{n}\rangle^{1/n}$ for $n=2$, $4$, and $12$. Note that

$\rho_{\max}=\lim_{n\to\infty}\langle\rho^{n}\rangle^{1/n},\quad\mbox{and}\quad\rho_{\min}=\lim_{n\to-\infty}\langle\rho^{n}\rangle^{1/n}.$
(20)

To show the values close to the collapse time $t_{\max}$ more clearly, we plot
the densities versus $(t_{\max}-t)\,c_{\rm s}k_{1}$ in a doubly-logarithmic
representation.

The nominal Jeans length is proportional to $\rho^{-1/2}$, although the
standard stability analysis breaks down if the density is no longer constant.
We see that the resolution criterion of Federrath et al. (2011b) is reached
when the square root of the density exceeds 1/30 of the initial Jeans length
of 410 grid cells for $\sigma_{\rm J}/c_{\rm s}k_{1}=5$; see Sect. 2.3.

## References

  * Alecian & Leorat (1988) Alecian G., Leorat J., 1988, A&A, 196, 1
  * André et al. (2007) André P., Belloche A., Motte F., Peretto N., 2007, A&A, 472, 519
  * Banerjee & Kritsuk (2018) Banerjee S., Kritsuk A. G., 2018, PhRvE, 97, 023107
  * Barreto-Mota et al. (2021) Barreto-Mota L., de Gouveia Dal Pino E. M., Burkhart B., Melioli C., Santos-Lima R., Kadowaki L. H. S., 2021, MNRAS, 503, 5425
  * Berger (1997) Berger M. A., 1997, J. Geophys. Res., 102, 2637
  * Binney & Tremaine (2008) Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition
  * Bonazzola et al. (1987) Bonazzola S., Heyvaerts J., Falgarone E., Perault M., Puget J. L., 1987, A&A, 172, 293
  * Bournaud et al. (2010) Bournaud F., Elmegreen B. G., Teyssier R., Block D. L., Puerari I., 2010, MNRAS, 409, 1088
  * Brandenburg (2009) Brandenburg A., 2009, ApJ, 697, 1206
  * Brandenburg & Matthaeus (2004) Brandenburg A., Matthaeus W. H., 2004, PhRvE, 69, 056407
  * Brandenburg et al. (2019) Brandenburg A., Kahniashvili T., Mandal S., Roper Pol A., Tevzadze A. G., Vachaspati T., 2019, PhRvF, 4, 024608
  * Brandenburg et al. (2020) Brandenburg A., Durrer R., Huang Y., Kahniashvili T., Mandal S., Mukohyama S., 2020, PhRvD, 102, 023536
  * Chen et al. (2016) Chen C.-Y., King P. K., Li Z.-Y., 2016, ApJ, 829, 84
  * Crutcher (1999) Crutcher R. M., 1999, ApJ, 520, 706
  * Crutcher (2012) Crutcher R. M., 2012, ARA&A, 50, 29
  * Dekel et al. (2009) Dekel A., et al., 2009, Natur, 457, 451
  * Elperin et al. (2003) Elperin T., Kleeorin N., Rogachevskii I., 2003, PhRvE, 68, 016311
  * Falco et al. (2013) Falco M., Hansen S. H., Wojtak R., Mamon G. A., 2013, MNRAS, 431, L6
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761, 156
  * Federrath et al. (2010) Federrath C., Banerjee R., Clark P. C., Klessen R. S., 2010, ApJ, 713, 269
  * Federrath et al. (2011a) Federrath C., Chabrier G., Schober J., Banerjee R., Klessen R. S., Schleicher D. R. G., 2011a, PhRvL, 107, 114504
  * Federrath et al. (2011b) Federrath C., Sur S., Schleicher D. R. G., Banerjee R., Klessen R. S., 2011b, ApJ, 731, 62
  * Federrath et al. (2014) Federrath C., Schober J., Bovino S., Schleicher D. R. G., 2014, ApJL, 797, L19
  * Federrath et al. (2021) Federrath C., Klessen R. S., Iapichino L., Beattie J. R., 2021, Nature Astronomy, 5, 365
  * Field et al. (2008) Field G. B., Blackman E. G., Keto E. R., 2008, MNRAS, 385, 181
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Haugen et al. (2003) Haugen N. E. L., Brandenburg A., Dobler W., 2003, ApJL, 597, L141
  * Haugen et al. (2004a) Haugen N. E., Brandenburg A., Dobler W., 2004a, PhRvE, 70, 016308
  * Haugen et al. (2004b) Haugen N. E. L., Brandenburg A., Mee A. J., 2004b, MNRAS, 353, 947
  * Hennebelle & Falgarone (2012) Hennebelle P., Falgarone E., 2012, A&AR, 20, 55
  * Hide & Palmer (1982) Hide R., Palmer T. N., 1982, Geophysical and Astrophysical Fluid Dynamics, 19, 301
  * Ibáñez-Mejía et al. (2017) Ibáñez-Mejía J. C., Mac Low M.-M., Klessen R. S., Baczynski C., 2017, ApJ, 850, 62
  * Jeans (1902) Jeans J. H., 1902, RSPTA, 199, 1
  * Johansen et al. (2007) Johansen A., Oishi J. S., Mac Low M.-M., Klahr H., Henning T., Youdin A., 2007, Natur, 448, 1022
  * Käpylä et al. (2009) Käpylä P. J., Mitra D., Brandenburg A., 2009, PhRvE, 79, 016302
  * Karoly et al. (2020) Karoly J., Soam A., Andersson B. G., Coudé S., Bastien P., Vaillancourt J. E., Lee C. W., 2020, ApJ, 900, 181
  * Klessen & Hennebelle (2010) Klessen R. S., Hennebelle P., 2010, A&A, 520, A17
  * Krumholz & Burkhart (2016) Krumholz M. R., Burkhart B., 2016, MNRAS, 458, 1671
  * Krumholz et al. (2018) Krumholz M. R., Burkhart B., Forbes J. C., Crocker R. M., 2018, MNRAS, 477, 2716
  * Lane et al. (2022) Lane H. B., Grudić M. Y., Guszejnov D., Offner S. S. R., Faucher-Giguère C.-A., Rosen A. L., 2022, MNRAS, 510, 4767
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Lee & Hennebelle (2018) Lee Y.-N., Hennebelle P., 2018, A&A, 611, A89
  * Mac Low & Klessen (2004) Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  * Mattsson & Hedvall (2022) Mattsson L., Hedvall R., 2022, MNRAS, 509, 3660
  * Mee & Brandenburg (2006) Mee A. J., Brandenburg A., 2006, MNRAS, 370, 415
  * Mestel & Spitzer (1956) Mestel L., Spitzer L. J., 1956, MNRAS, 116, 503
  * Moffatt (1978) Moffatt H. K., 1978, Magnetic Field Generation in Electrically Conducting Fluids. Cambridge University Press, Cambridge
  * Mouschovias & Spitzer (1976) Mouschovias T. C., Spitzer L. J., 1976, ApJ, 210, 326
  * Myers et al. (1991) Myers P. C., Ladd E. F., Fuller G. A., 1991, ApJL, 372, L95
  * Nelson et al. (2015) Nelson D., Genel S., Vogelsberger M., Springel V., Sijacki D., Torrey P., Hernquist L., 2015, MNRAS, 448, 59
  * Nordlund et al. (1992) Nordlund A., Brandenburg A., Jennings R. L., Rieutord M., Ruokolainen J., Stein R. F., Tuominen I., 1992, ApJ, 392, 647
  * Ntormousi & Hennebelle (2019) Ntormousi E., Hennebelle P., 2019, A&A, 625, A82
  * Passot et al. (1995) Passot T., Vazquez-Semadeni E., Pouquet A., 1995, ApJ, 455, 536
  * Pattle et al. (2021) Pattle K., et al., 2021, ApJ, 907, 88
  * Pencil Code Collaboration et al. (2021) Pencil Code Collaboration et al., 2021, JOSS, 6, 2807
  * Porter et al. (2015) Porter D. H., Jones T. W., Ryu D., 2015, ApJ, 810, 93
  * Pudritz et al. (2006) Pudritz R. E., Rogers C. S., Ouyed R., 2006, MNRAS, 365, 1131
  * Roper Pol et al. (2020) Roper Pol A., Mandal S., Brandenburg A., Kahniashvili T., Kosowsky A., 2020, PhRvD, 102, 083512
  * Schneider et al. (2013) Schneider N., et al., 2013, ApJL, 766, L17
  * Sharda et al. (2021) Sharda P., Federrath C., Krumholz M. R., Schleicher D. R. G., 2021, MNRAS, 503, 2014
  * Shu (1977) Shu F. H., 1977, ApJ, 214, 488
  * Shu (1992) Shu F. H., 1992, The physics of astrophysics. Volume II: Gas dynamics.
  * Soler & Hennebelle (2017) Soler J. D., Hennebelle P., 2017, A&A, 607, A2
  * Soler et al. (2013) Soler J. D., Hennebelle P., Martin P. G., Miville-Deschênes M. A., Netterfield C. B., Fissel L. M., 2013, ApJ, 774, 128
  * Sur (2019) Sur S., 2019, MNRAS, 488, 3439
  * Sur et al. (2010) Sur S., Schleicher D. R. G., Banerjee R., Federrath C., Klessen R. S., 2010, ApJL, 721, L134
  * Sur et al. (2012) Sur S., Federrath C., Schleicher D. R. G., Banerjee R., Klessen R. S., 2012, MNRAS, 423, 3148
  * Tennekes & Lumley (1972) Tennekes H., Lumley J. L., 1972, First Course in Turbulence
  * Truelove et al. (1997) Truelove J. K., Klein R. I., McKee C. F., Holliman John H. I., Howell L. H., Greenough J. A., 1997, ApJL, 489, L179
  * Xu & Lazarian (2020) Xu S., Lazarian A., 2020, ApJ, 899, 115