AppendixAppendix References

# Magnetic field fluctuations in anisotropic, supersonic turbulence

James R. Beattie1, Christoph Federrath1 and Amit Seta1
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia E-mail<EMAIL_ADDRESS>
(Accepted 2020 July 24. Received 2020 July 18; in original form 2020 May 28.)

###### Abstract

The rich structure that we observe in molecular clouds is due to the interplay
between strong magnetic fields and supersonic (turbulent) velocity
fluctuations. The velocity fluctuations interact with the magnetic field,
causing it too to fluctuate. Using numerical simulations, we explore the
nature of such magnetic field fluctuations,
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, over a wide range of turbulent
Mach numbers, $\operatorname{\mathcal{M}}=2-20$ (i.e., from weak to strong
compressibility), and Alfvén Mach numbers,
$\operatorname{\mathcal{M}_{\text{A0}}}=0.1-100$ (i.e., from strong to weak
magnetic mean fields, $B_{0}$). We derive a compressible quasi-static
fluctuation model from the magnetohydrodynamical (MHD) equations and show that
velocity gradients parallel to the mean magnetic field give rise to
compressible modes in sub-Alfvénic flows, which prevents the flow from
becoming two-dimensional, as is the case in incompressible MHD turbulence. We
then generalise an analytical model for the magnitude of the magnetic
fluctuations to include $\operatorname{\mathcal{M}}$, and find
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|=\delta
B=c_{s}\sqrt{\pi\rho_{0}}\operatorname{\mathcal{M}}\operatorname{\mathcal{M}_{\text{A0}}}$,
where $c_{s}$ is the sound speed and $\rho_{0}$ is the mean density of gas.
This new relation fits well in the strong $B$-field regime. We go on to study
the anisotropy between the perpendicular ($B_{\perp}$) and parallel
($B_{\parallel}$) fluctuations and the mean-normalised fluctuations, which we
find follow universal scaling relations, invariant of
$\operatorname{\mathcal{M}}$. We provide a detailed analysis of the morphology
for the $\delta B_{\perp}$ and $\delta B_{\parallel}$ probability density
functions and find that eddies aligned with $B_{0}$ cause parallel
fluctuations that reduce $B_{\parallel}$ in the most anisotropic simulations.
We discuss broadly the implications of our fluctuation models for magnetised
gases in the interstellar medium.

###### keywords:

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

††pubyear: 2020††pagerange: Magnetic field fluctuations in anisotropic,
supersonic turbulence–C

##  1 Introduction

Magnetised plasmas that undergo random velocity fluctuations are ubiquitous in
the Universe. For cool, molecular clouds, the birthplace of stars, the
velocity fluctuations are supersonic (Larson, 1981; Solomon et al., 1987;
Padoan et al., 1997; Ossenkopf & Mac Low, 2002; Elmegreen & Scalo, 2004; Heyer
& Brunt, 2004; Mac Low & Klessen, 2004; Krumholz & McKee, 2005; Ballesteros-
Paredes et al., 2007; Roman-Duval et al., 2011; Federrath, 2013; Schneider et
al., 2013; Chevance et al., 2020). But random velocity fluctuations alone
cannot explain all of the structure that we observe in these clouds. Indeed,
strong magnetic fields facilitate and help set the structure of the neutral
hydrogen clouds (Planck Collaboration et al., 2016a, b; Hennebelle & Inutsuka,
2019; Krumholz & Federrath, 2019; Clark & Hensley, 2019). For example,
magnetic fields create magnetosonic striations in star-forming clouds (Tritsis
& Tassis, 2016, 2018; Tritsis et al., 2018; Beattie & Federrath, 2020b) and in
general, contribute to changing the density dispersion of the clouds through
magnetic cushioning (Molina et al., 2012; Mandal et al., 2020). Magnetic
fields polarise the density structures in molecular clouds and hydrodynamical
shocks preferentially form across magnetic field lines (Clark et al., 2015;
Planck Collaboration et al., 2016a, b; Cox et al., 2016; Malinen et al., 2016;
Soler et al., 2017; Soler, 2019; Beattie & Federrath, 2020b; Clark & Hensley,
2019; Seifried et al., 2020). Magnetic fields inhibit cloud fragmentation, in
turn influencing the initial mass function for new-born stars in the modern
star-forming era (Price & Bate, 2008; Hennebelle et al., 2011; Federrath &
Klessen, 2012; Federrath & Banerjee, 2015; Krumholz & Federrath, 2019) and the
era of first stars (Sharda et al., 2020). They facilitate anisotropic cloud
collapse through flux-freezing of the initial mass-to-flux ratio, which causes
molecular clouds to preferentially collapse parallel to large scale magnetic
fields (Tritsis et al., 2015; Mocz et al., 2017; Mocz & Burkhart, 2018). At
smaller scales magnetic fields also facilitate the launching of protostellar
jets and outflows (Frank et al., 2014; Kuruwita et al., 2017; Gerrard et al.,
2019; Krumholz & Federrath, 2019; Kuruwita & Federrath, 2019), and influence
the accretion rate of newly forming stars (Kuruwita et al., 2020). This means,
if one wants to study the dynamics of star-forming molecular clouds (MCs), and
indeed calculate the potential that a cloud has to form stars, one must
understand the properties of magnetised, supersonic turbulence (Elmegreen &
Scalo, 2004; Mac Low & Klessen, 2004; McKee & Ostriker, 2007; Hennebelle &
Falgarone, 2012; Padoan et al., 2014; Federrath & Klessen, 2012, 2013;
Federrath, 2018).

###  1.1 Magnetised Turbulence

Magnetohydrodynamical (MHD) turbulence comes in a number of different
flavours, largely depending on the strength of the (kinetic) turbulence
compared to the strength of the magnetic field and whether or not the flow is
compressible (Beresnyak, 2019, for a modern review). One can encode this
information into the Alfvén Mach number,

$\operatorname{\mathcal{M}_{\text{A}}}=\frac{V}{V_{A}}=\sqrt{4\pi\rho}\frac{V}{B}$
(1)

where $V$, $V_{A}=B/\sqrt{4\pi\rho}$, $\rho$ and $B$ are the root-mean-squared
(rms) velocity, Alfvén velocity, density and magnetic field strength,
respectively. For $\operatorname{\mathcal{M}_{\text{A}}}<1$, $\rho
V^{2}\lesssim B^{2}$ the magnetic field contributes significantly to the
dynamics of the flow through the Lorentz force. This is the sub-Alfvénic
regime. For $\operatorname{\mathcal{M}_{\text{A}}}>1$, $\rho V^{2}\gtrsim
B^{2}$ the magnetic field plays a lesser role and the turbulent motions set,
for example the statistics of the flow (Padoan & Nordlund, 2011; Molina et
al., 2012; Beattie & Federrath, 2020b). This is called the super-Alfvénic
regime, and the trans-Alfvénic regime,
$\operatorname{\mathcal{M}_{\text{A}}}\sim 1$, separates the two. In our study
we explore a broad range of $\operatorname{\mathcal{M}_{\text{A}}}$, from
$0.1-100$, to include the two dynamically dissimilar regimes.

Turbulent magnetic fields can be separated into two components,

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)=\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}(t)+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t),$
(2)

where $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ is the total field,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$ is the ordered component of the
field, and $\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ is the turbulent (or
fluctuating) component of the field (sometimes called
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{\text{turb}}$ or $B_{\text{t}}$ in
some literature; Pillai et al. 2015; Federrath 2016a). The ordered component
of the magnetic field can extend over an entire molecular cloud, creating
coherent magnetic field structures at the parsec scale that evolve on much
longer dynamical times than the fluctuating component of the field (Bertrang
et al., 2014; Pillai et al., 2015; Federrath et al., 2016a; Hu et al., 2019).
We explore the case in our study where the ordered component does not change
in time at all, and the fluctuating component evolves self-consistently with
the MHD equations, in a statistically stationary state, hence, we can rewrite
Equation (2) as

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)=\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t),$
(3)

where
$\mathrm{\left\langle\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)\right\rangle}_{t}=|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|$
which means
$\mathrm{\left\langle\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)\right\rangle}_{t}=0$.
For this reason, we will refer to $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$
as the mean-field throughout this study.

###  1.2 The Alfvén Mach Number of MCs

Since the magnetic field has two distinct components one can talk about the
mean-field or fluctuating Alfvén Mach number. In this study, we find the
dependence of the magnetic field fluctuations,
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, on the Alfvén Mach number
defined with respect to the mean-field. We define the mean-field Alfvén Mach
number exactly as we defined the rms $\operatorname{\mathcal{M}_{\text{A}}}$
in Equation (1), but with mean-field components,
$\operatorname{\mathcal{M}_{\text{A0}}}=V/V_{\text{A}0}=(\sqrt{4\pi\rho_{0}}V)/B_{0}$,
where $B_{0}=|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|$ and $V_{\text{A}0}$
is the velocity of Alfvén waves along the mean-field. We explore a wide
parameter set that encompasses both the sub- and super-Alfvénic mean-field
regime, covering the parameter space of observed molecular clouds.

For example, for the central molecular zone cloud, G0.253+0.016 studied in
Federrath et al. (2016b), they measure a mean-field of $B_{0}=(2.07\pm
0.95)\,\text{mG}$ (Pillai et al., 2015), a mean volume density of
$\rho_{0}=(6.2\pm 3.3)\times 10^{-20}\,\text{g/cm${}^{3}$}$ and a velocity
dispersion of $V=(6.8\pm 0.2)\,\text{km/s}$ (see Federrath et al. 2016b for
details on assumptions and further references). Using these quantities and
Equation (1) one can calculate $V_{A0}=(24\pm 17)\,\text{km/s}$ which means
that $\operatorname{\mathcal{M}_{\text{A0}}}=0.3\pm 0.2$, placing it well into
the sub-Alfvénic mean-field regime. Furthermore, using the direction of the
velocity gradients to infer the magnetic field strength, Hu et al. (2019)
measure the Alfvén Mach number for five star-forming clouds in the Gould Belt:
Taurus ($\operatorname{\mathcal{M}_{\text{A0}}}=1.19\pm 0.02$), Perseus A
($1.22\pm 0.05$), L 1551 ($0.73\pm 0.13$), Serpens ($0.98\pm 0.08$) and NGC
1333 ($0.82\pm 0.24$). Hence, the trans to sub-Alfvénic magnetised turbulence
regime is of great importance for understanding the environment that stars
form in, which is the focus of this study.

Previous work by Federrath (2016a, herein called F16) explored the magnitude
of magnetic field fluctuations across a broad range of $B_{0}$ and $\delta B$
values, but for a single $\operatorname{\mathcal{M}}$ value. In our study we
will expand upon some of the key results in F16 and elucidate how the
$\operatorname{\mathcal{M}}$ number influences the fluctuating magnetic field.
Of particular importance for our study is the F16 analytical model for $\delta
B$ in the strong-field regime, where $B_{0}\gg\delta B$. By relating the
turbulent magnetic energy density, $e_{\text{m}}\approx B_{0}\delta B/(4\pi)$
and the turbulent kinetic energy $e_{\text{k}}=(\rho_{0}V^{2})/2$, F16 found
that

$\displaystyle\delta B=\mathcal{M}^{2}_{\text{A0}}B_{0}/2,$ (4)

by assuming that all of the turbulent magnetic energy was fed from the kinetic
energy in the plasma, i.e., $e_{\text{m}}=e_{\text{k}}$. We expand upon this
study significantly by exploring the fluctuations perpendicular and parallel
to $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$ separately in high-resolution
MHD simulations, by showing the mean normalised fluctuations are independent
of $\operatorname{\mathcal{M}}$, and by deriving a new model directly from the
induction equation in the compressible, ideal MHD model of plasmas. Hence our
study contributes significantly to better understanding the nature of magnetic
fluctuations in compressible, astrophysically relevant flows.

The study is organised into the following sections. First, in §2 we derive a
compressible quasi-static model for the fluctuations of the magnetic field
from the induction equation, which is relevant for understanding how the
magnetic anisotropy influences the dynamics of molecular clouds with a strong
mean magnetic field present. We use our model to show that velocity gradients
parallel to the magnetic field give rise to compressible modes in the
turbulence, and we provide an order of magnitude estimate for the fluctuations
as a function of $\operatorname{\mathcal{M}}$. In §3 we discuss the
simulations that we use to explore the magnetic field fluctuations and test
our models. In §4 we qualitatively explore how velocity gradients parallel to
the mean magnetic field create compressible modes in the velocity field and
discuss why this is relevant to the analysis of astrophysical observations. In
§5 we test our fluctuation model on the simulation data, and show it reduces
to the F16 analytical model for the strong-field fluctuations. We explore the
anisotropy and mean-normalised fluctuations, which show universal scaling laws
that do not depend upon $\operatorname{\mathcal{M}}$. We then explore the
morphology of the magnetic field probability density functions (PDFs) parallel
and perpendicular to the mean-field. Finally in §6 we summarise our key
findings.

##  2 Compressible Magnetic Fluctuation Model

###  2.1 Model derivation

In this section we derive a model for the amplitude of the fluctuations
directly from the induction equation in the ideal MHD framework. We use the
model derived in this section to explain the magnetic field fluctuations in
the remainder of the study. The ideal, isothermal MHD equations are

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}})$ $\displaystyle=0,$
(5) $\displaystyle\left(\frac{\partial}{\partial
t}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla\right)\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
$\displaystyle=\frac{(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{4\pi}-\nabla\left(c_{s}^{2}\rho+\frac{|\mathrm{{\mn@boldsymbol{\mathit{B}}}}|^{2}}{8\pi}\right)+\rho\mathrm{{\mn@boldsymbol{\mathit{F}}}},$
(6)
$\displaystyle\frac{\partial\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{\partial
t}$
$\displaystyle=\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}),$
(7) $\displaystyle\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
$\displaystyle=0,$ (8)

where $\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ is the fluid velocity, $\rho$
the density, $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ the magnetic field,
$c_{s}$ the sound speed and $\mathrm{{\mn@boldsymbol{\mathit{F}}}}$, a
stochastic driving function that gives rise to the turbulence (which could be
from, for example, supernova shocks permeating through the interstellar
medium, or internal to the MC, gravity, galactic-scale shocks 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; Colling et al. 2018; Schruba
et al. 2019; Lu et al. 2020). We consider a magnetic field of the form shown
in Equation (3), with constant mean component of the field in the
$z$-direction,

$\displaystyle\begin{pmatrix}B_{x}(t)\\\ B_{y}(t)\\\ B_{z}(t)\end{pmatrix}$
$\displaystyle=\begin{pmatrix}0\\\ 0\\\
B_{0}\end{pmatrix}+\begin{pmatrix}\delta B_{x}(t)\\\ \delta B_{y}(t)\\\ \delta
B_{z}(t)\end{pmatrix}.$ (9)

This gives rise to two natural length scales in the vector field: one
perpendicular to $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$, and one parallel
to $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$. This is because the mean
magnetic field defines an axis of symmetry in the turbulent flow (Cho et al.,
2002). We will denote the two scales with $\perp$ (perpendicular) and
$\parallel$ (parallel) subscripts throughout this study, respectively. To
develop some insight into the nature of the magnetic fluctuations we now
explore the time evolution of Equation (9) by propagating it through the
induction equation (Equation 7). By construction, one can show that the
induction equation only governs the time evolution of the fluctuating
component of the field, since the mean-field is time-independent and, since
the mean-field is only in the $z$ direction, the
$\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}})$
term can be significantly simplified (we show the full derivation of this in
Appendix A), to show that

$\displaystyle\left(\frac{\partial}{\partial
t}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla\right)\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}=$
$\displaystyle\overbrace{B_{0}\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}}^{\text{$v_{z}$-gradient}}+\overbrace{(\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{v}}}}}^{\text{advection}}$
$\displaystyle-\underbrace{(\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}})(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})}_{\text{compression}},$
(10)

which reveals that in the Lagrangian frame of the fluid, the magnetic field
fluctuations are influenced by (1) the velocity gradients along the mean
magnetic field (which we indicate with $\partial_{\parallel}$, i.e., they are
anisotropic), (2) the advection of the velocity by the fluctuations, and (3)
the compression of the velocity field, assuming
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$ is independent of time.
Disregarding the compression term,
$(\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}})(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})$,
and as $\partial_{t}\delta B\rightarrow 0$, this form of the induction
equation is known as the quasi-static approximation of the magnetic field
perturbation, which is only applicable in the case of incompressible plasma,
(Zikanov & Thess, 1998; Verma, 2017, and references therein), but here we have
generalised the equation for compressible MHD plasmas, relevant for
astrophysical phenomena.

For sub-Alfvénic molecular clouds, such as G0.253+0.016, $B_{0}\gg\delta B$
(in this case, $B_{0}>\delta B$ by an order of magnitude; Federrath et al.
2016b), and the time derivative of $\delta B$ is very small, i.e., much
smaller than the time scale of compression. Hence, by assuming $B_{0}\gg\delta
B$ and $(D/Dt)\delta B\approx
0$,111$\frac{D}{Dt}\equiv\frac{\partial}{\partial
t}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla$ is the material
derivative. Equation (2.1) simply becomes

$\displaystyle\frac{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}{B_{0}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})=\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}},$
(11)

and by taking the magnitudes of both sides,

$\displaystyle|\nabla\cdot\mathrm{{\mn<EMAIL_ADDRESS>(12)

Thus, it is primarily the velocity streams that run parallel with the mean
magnetic field that give rise to compressive modes in the clouds222Note that
these are in addition to the mixture of various modes caused by the turbulent
driving (Federrath et al., 2008, 2010).. We show the full derivation of
Equation (12) in Appendix B. This is consistent with other studies, where
hydrodynamical shocks appear perpendicular to the mean magnetic field, with
the compression happening in the parallel direction to the mean magnetic field
(Mocz & Burkhart, 2018; Beattie & Federrath, 2020b; Seifried et al., 2020). We
will discuss this further in §4.

Equation (12) shows the key difference between incompressible and compressible
MHD turbulence for the case when $B_{0}\gg\delta B$. If, for example,
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}=0$, the velocity gradient
along the field disappears, and hence the magnetic field fluctuations cannot
persist along the field. This is known as the two-dimensionalisation of the
three-dimensional flow (Alexakis, 2011; Verma, 2017), which for supersonic
turbulence, is forbidden to happen as shocks form and travel parallel to the
field from the parallel velocity gradient. This means magnetic field
fluctuations persist along the field, and the flow retains its three-
dimensional nature. We will see this is indeed the case for parallel
fluctuations in numerical simulations (demonstrated in §5).

###  2.2 An order of magnitude estimate for $\delta B$

We now consider the dimensionless form of Equation (2.1), where length scales,
$\hat{\mathrm{{\mn@boldsymbol{\mathit{\ell}}}}}=\mathrm{{\mn@boldsymbol{\mathit{\ell}}}}/L$,
velocities
$\hat{\mathrm{{\mn@boldsymbol{\mathit{v}}}}}=\mathrm{{\mn@boldsymbol{\mathit{v}}}}/V$,
and magnetic fields,
$\delta\hat{\mathrm{{\mn@boldsymbol{\mathit{B}}}}}=\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}/\delta
B$,
$\hat{\mathrm{{\mn@boldsymbol{\mathit{B}}}}}_{0}=\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}/B_{0}$,
are normalised by their magnitudes, and the time-scale in the derivative we
chose to be $\hat{t}=t/T$, where $T=L/V_{A0}$, is the time-scale of
fluctuations travelling along the mean magnetic field that spans the full
system scale $L$. By grouping $\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
and $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$ terms on the right hand side
of Equation (2.1) and by applying our non-dimensionalisation, we find,

$\displaystyle\frac{\delta
BV_{A0}}{L}\frac{\partial\delta\mathrm{{\mn@boldsymbol{\mathit{\hat{B}}}}}}{\partial\hat{t}}$
$\displaystyle+\frac{\delta
BV}{L}(\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}\cdot\hat{\nabla})\delta\mathrm{{\mn@boldsymbol{\mathit{\hat{B}}}}}=$
$\displaystyle\frac{B_{0}V}{L}\left[\hat{\partial}_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}-\mathrm{{\mn@boldsymbol{\mathit{\hat{B}}}}}_{0}(\hat{\nabla}\cdot\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}})\right]+$
$\displaystyle\frac{\delta
BV}{L}\left[(\delta\mathrm{{\mn@boldsymbol{\mathit{\hat{B}}}}}\cdot\hat{\nabla})\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}-\delta\mathrm{{\mn@boldsymbol{\mathit{\hat{B}}}}}(\hat{\nabla}\cdot\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}})\right].$
(13)

If we assume that
$\hat{\partial}_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}-\mathrm{{\mn@boldsymbol{\mathit{\hat{B}}}}}_{0}(\hat{\nabla}\cdot\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}})$
and
$(\delta\mathrm{{\mn@boldsymbol{\mathit{\hat{B}}}}}\cdot\hat{\nabla})\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}}-\delta\mathrm{{\mn@boldsymbol{\mathit{\hat{B}}}}}(\hat{\nabla}\cdot\mathrm{{\mn@boldsymbol{\mathit{\hat{v}}}}})$
are $\mathcal{O}(1)$, which is true by construction, then we can make an
order-of-magnitude estimate of the fluctuations,

$\displaystyle\frac{\delta BV_{A0}}{L}+\frac{\delta
BV}{L}\sim\frac{B_{0}V}{L}+\frac{\delta BV}{L},$ (14) $\displaystyle
V_{A0}+V\sim\frac{B_{0}}{\delta B}V+V,$ (15) $\displaystyle\delta B\sim
B_{0}\operatorname{\mathcal{M}_{\text{A0}}}.$ (16)

Since
$B_{0}=2c_{s}\sqrt{\pi\rho_{0}}(\operatorname{\mathcal{M}}/\operatorname{\mathcal{M}_{\text{A0}}})$
by definition, we ultimately find

$\displaystyle\delta B=2c_{s}\sqrt{\pi\rho_{0}}C\operatorname{\mathcal{M}},$
(17)

where $C$ is a proportionality factor, and is most likely a function of the
mean magnetic field since the fluctuations will change as the strength of the
mean-field changes (Federrath, 2016a, herein called F16). The key feature of
this analysis is that from the induction equation we derived that there is a
linear dependence between the magnetic field fluctuations ($\delta B$) and the
turbulent velocity fluctuations ($\propto\operatorname{\mathcal{M}}$). This is
consistent with the strong-field model for the fluctuations proposed by F16
(where $C=\operatorname{\mathcal{M}_{\text{A0}}}/2)$, which will be discussed
later in §5. We note also that from Equation (14), the ratio between $\delta
B$ and $B_{0}$ is independent of $\operatorname{\mathcal{M}}$, which we will
also explore further. Next we discuss the simulations that we perform to test
our models and explore the magnetic fluctuations.

Figure 1: A slice of $|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|$
compensated by $B_{0}$, at $t=5\,T$, where $B_{0}$ is the magnitude of the
mean of the magnetic field, which is orientated up the page,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}=B_{0}\hat{\mathrm{{\mn@boldsymbol{\mathit{z}}}}}$.
The turbulent Mach number increases from $\operatorname{\mathcal{M}}=2-20$,
from left to right, and the Alfvénic Mach number of the mean magnetic field
increases from $\operatorname{\mathcal{M}_{\text{A0}}}=0.1-100$ from top to
bottom, hence the simulation with the largest $B_{0}$ and
$\operatorname{\mathcal{M}}$ is in the top right corner, and the smallest
$B_{0}$ and $\operatorname{\mathcal{M}}$ in the bottom left corner. White
indicates that $|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|=B_{0}$, red
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|<B_{0}$ and black
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|>B_{0}$.

##  3 Numerical Simulations

Table 1: Main simulation parameters and derived quantities used throughout
this study.

| Turbulence | Mean Magnetic Field | Fluctuating Magnetic Field
---|---|---|---
Simulation | $\operatorname{\mathcal{M}}$ | $\operatorname{\mathcal{M}_{\text{A0}}}$ | $B_{0}$ | $\delta B_{\perp}\,[\mu\text{G}]$ | $\delta B_{\parallel}\,[\mu\text{G}]$ | $\delta B\,[\mu\text{G}]$
ID | $(\pm 1\sigma)$ | $(\pm 1\sigma)$ | $[\mu\text{G}]$ | $(\pm 1\sigma)$ | $(\pm 1\sigma)$ | $(\pm 1\sigma)$
(1) | (2) | (3) | (4) | (5) | (6) | (7)
M2Ma100 | 2.1 |  $\pm$ 0.1 | 103 |  $\pm$ 5 | 0.07 | 1.7 |  $\pm$ 0.1 | 1.2 |  $\pm$ 0.1 | 2.1 |  $\pm$ 0.2
M2Ma10 | 1.80 |  $\pm$ 0.08 | 9.0 |  $\pm$ 0.8 | 0.71 | 3.8 |  $\pm$ 0.2 | 2.6 |  $\pm$ 0.1 | 4.6 |  $\pm$ 0.2
M2Ma2 | 1.66 |  $\pm$ 0.05 | 1.7 |  $\pm$ 0.1 | 3.54 | 4.4 |  $\pm$ 0.4 | 2.7 |  $\pm$ 0.1 | 5.2 |  $\pm$ 0.4
M2Ma1 | 2.0 |  $\pm$ 0.1 | 0.98 |  $\pm$ 0.07 | 7.09 | 3.5 |  $\pm$ 0.3 | 2.0 |  $\pm$ 0.1 | 4.1 |  $\pm$ 0.3
M2Ma0.5 | 2.2 |  $\pm$ 0.2 | 0.54 |  $\pm$ 0.04 | 14.18 | 2.5 |  $\pm$ 0.3 | 1.5 |  $\pm$ 0.1 | 2.9 |  $\pm$ 0.3
M2Ma0.1 | 2.6 |  $\pm$ 0.2 | 0.131 |  $\pm$ 0.008 | 70.90 | 0.4 |  $\pm$ 0.04 | 0.8 |  $\pm$ 0.1 | 0.8 |  $\pm$ 0.1
M4Ma100 | 4.0 |  $\pm$ 0.2 | 101 |  $\pm$ 4 | 0.14 | 2.6 |  $\pm$ 0.3 | 1.8 |  $\pm$ 0.3 | 3.1 |  $\pm$ 0.4
M4Ma10 | 3.7 |  $\pm$ 0.1 | 9.2 |  $\pm$ 0.6 | 1.42 | 6.5 |  $\pm$ 0.4 | 4.4 |  $\pm$ 0.1 | 7.8 |  $\pm$ 0.4
M4Ma2 | 3.5 |  $\pm$ 0.1 | 1.73 |  $\pm$ 0.07 | 7.09 | 8.3 |  $\pm$ 0.4 | 4.9 |  $\pm$ 0.2 | 9.6 |  $\pm$ 0.4
M4Ma1 | 3.8 |  $\pm$ 0.3 | 0.95 |  $\pm$ 0.08 | 14.18 | 6.6 |  $\pm$ 0.6 | 3.6 |  $\pm$ 0.2 | 7.5 |  $\pm$ 0.7
M4Ma0.5 | 4.4 |  $\pm$ 0.2 | 0.54 |  $\pm$ 0.03 | 28.36 | 4.7 |  $\pm$ 0.6 | 2.7 |  $\pm$ 0.3 | 5.4 |  $\pm$ 0.7
M4Ma0.1 | 5.2 |  $\pm$ 0.4 | 0.13 |  $\pm$ 0.01 | 141.80 | 0.9 |  $\pm$ 0.2 | 1.0 |  $\pm$ 0.2 | 1.3 |  $\pm$ 0.3
M10Ma100 | 10.0 |  $\pm$ 0.3 | 100 |  $\pm$ 3 | 0.35 | 5.9 |  $\pm$ 0.8 | 4.0 |  $\pm$ 0.6 | 7 |  $\pm$ 1
M10Ma10 | 9.2 |  $\pm$ 0.4 | 9.2 |  $\pm$ 0.7 | 3.54 | 15.1 |  $\pm$ 0.6 | 9.9 |  $\pm$ 0.4 | 18 |  $\pm$ 1
M10Ma2 | 9.0 |  $\pm$ 0.4 | 1.8 |  $\pm$ 0.1 | 17.72 | 19 |  $\pm$ 1 | 11.7 |  $\pm$ 0.7 | 22 |  $\pm$ 1
M10Ma1 | 9.3 |  $\pm$ 0.5 | 0.93 |  $\pm$ 0.05 | 35.45 | 16 |  $\pm$ 1 | 8.2 |  $\pm$ 0.6 | 18 |  $\pm$ 2
M10Ma0.5 | 10.5 |  $\pm$ 0.4 | 0.52 |  $\pm$ 0.02 | 70.90 | 10 |  $\pm$ 1 | 5.8 |  $\pm$ 0.8 | 12 |  $\pm$ 1
M10Ma0.1 | 12 |  $\pm$ 1 | 0.125 |  $\pm$ 0.006 | 354.49 | 2.3 |  $\pm$ 0.3 | 2.5 |  $\pm$ 0.3 | 3.4 |  $\pm$ 0.5
M20Ma100 | 20 |  $\pm$ 1 | 101 |  $\pm$ 4 | 0.71 | 9.9 |  $\pm$ 0.8 | 6.8 |  $\pm$ 0.6 | 12 |  $\pm$ 1
M20Ma10 | 19 |  $\pm$ 1 | 9.3 |  $\pm$ 0.8 | 7.09 | 28 |  $\pm$ 2 | 18.5 |  $\pm$ 0.8 | 34 |  $\pm$ 2
M20Ma2 | 18 |  $\pm$ 1 | 1.8 |  $\pm$ 0.1 | 35.45 | 4.4 |  $\pm$ 0.3 | 2.7 |  $\pm$ 0.1 | 41 |  $\pm$ 2
M20Ma1 | 19 |  $\pm$ 1 | 0.93 |  $\pm$ 0.03 | 70.90 | 31 |  $\pm$ 3 | 16 |  $\pm$ 1 | 35 |  $\pm$ 3
M20Ma0.5 | 21 |  $\pm$ 1 | 0.53 |  $\pm$ 0.02 | 141.80 | 21 |  $\pm$ 3 | 11 |  $\pm$ 1 | 24 |  $\pm$ 3
M20Ma0.1 | 24 |  $\pm$ 1 | 0.119 |  $\pm$ 0.003 | 708.98 | 3.8 |  $\pm$ 0.4 | 4.1 |  $\pm$ 0.5 | 5.6 |  $\pm$ 0.6

  * •

Notes: For each simulation we extract 51
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ realisations at 0.1$\,T$ intervals,
where $T$ is the turbulent turnover time, between $5\,T$ and $10\,T$. All
$1\sigma$ fluctuations listed are from the time-averaging over the $5\,T$. We
show an example of the fluctuating component of the field in units of mean-
field in Figure 1. All simulations have $512^{3}$ grid cells. Column (1): the
simulation ID. Column (2): the rms turbulent Mach number,
$\operatorname{\mathcal{M}}=\sigma_{V}/c_{s}$. Column (3): the Alfvén Mach
number for the mean-$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ component,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$,
$\operatorname{\mathcal{M}_{\text{A}}}_{0}=(2c_{s}\operatorname{\mathcal{M}}\sqrt{\pi\rho_{0}})/|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|$,
where $\rho_{0}$ is the mean density, $c_{s}$ is the sound speed. Column (4):
the magnitude of the mean magnetic field,
$|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|=B_{0}$. Column (5): the
fluctuating component of the magnetic field that is perpendicular to the mean
magnetic field, where $B_{\perp}=B_{x}=B_{y}$, discussed in §5, and plotted in
Figure 4. Column (6): the same as column (5) but for the fluctuating component
parallel to the mean magnetic field, $B_{z}=B_{\parallel}$. Column (7): the
total magnitude of the fluctuations $\delta
B=|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|=\sqrt{\delta B_{x}^{2}+\delta
B_{y}^{2}+\delta B_{z}^{2}}$.

###  3.1 Turbulent MHD Model

In this study we analyse the magnetic field fluctuations in 24 high-
resolution, 3D turbulent, ideal magnetohydrodynamical (MHD) simulations with a
non-zero mean magnetic field and an isothermal equation of state,
$P=c_{s}^{2}\rho$, where $P$ is the pressure, with $c_{s}$ normalised to 1. We
use a modified version of the flash code, based on the public version 4.0.1
(Fryxell et al., 2000; Dubey et al., 2008) to solve the MHD equations (5–8) in
a periodic box with dimensions $L^{3}$, on a uniform grid with a uniform
resolution of $512^{3}$ grid cells, using the multi-wave, approximate Riemann
solver framework described in Bouchut et al. (2010), and implemented and
tested in Waagan et al. (2011). A comprehensive parameter set, including Mach
number, mean magnetic fields, and derived quantities for each of the
simulations is listed in Table 1.

###  3.2 Turbulent Driving, Density and Velocity Fields

The turbulent acceleration field $\mathrm{{\mn@boldsymbol{\mathit{F}}}}$ in
Equation (6) follows an Ornstein-Uhlenbeck (OU) process in time and is
constructed such that we can control the mixture of solenoidal and purely
compressive modes in $\mathrm{{\mn@boldsymbol{\mathit{F}}}}$ (see Federrath et
al. 2008; Federrath et al. 2009; Federrath et al. 2010 for a detailed
discussion of the turbulence driving). We choose to drive with a natural
mixture of the two modes (Federrath et al., 2010). We isotropically drive in
wavenumber space at $k\approx 2$, corresponding to real-space scales of
$\ell_{D}\approx L/2$. The driving amplitude is centred on $k=2$ and falls off
to zero with a parabolic spectrum towards $k=1$ on one side and $k=3$ on the
other. Thus, only large scales are driven and the turbulence on smaller scales
develops through the turbulent energy cascade driven from those large scales.
The auto-correlation timescale of $\mathrm{{\mn@boldsymbol{\mathit{F}}}}$ is
equal to $T=L/(2c_{s}\operatorname{\mathcal{M}})$, where $L$ is the system
scale, hence we use the auto-correlation timescale of
$\mathrm{{\mn@boldsymbol{\mathit{F}}}}$ to set the turbulent turnover
timescale for the desired turbulent Mach number on $L/2$ for each simulation.
We vary the sonic Mach number between $\operatorname{\mathcal{M}}=$ 2 and 20,
encompassing the range of observed $\operatorname{\mathcal{M}}$ values for the
turbulent interstellar medium (e.g., Schneider et al. 2013; Federrath et al.
2016b; Orkisz et al. 2017; Beattie et al. 2019). The initial velocity field is
set to $\mathrm{{\mn@boldsymbol{\mathit{v}}}}(x,y,z,t=0)=0$ and the density
field $\rho(x,y,z,t=0)=1$, i.e., the mean density, $\rho_{0}=1$. We run the
simulations for $10\,T$, where $T$ is the eddy turnover time. We note that the
physical evolution of these systems is fully described by the dimensionless
$\operatorname{\mathcal{M}}$ and $\operatorname{\mathcal{M}_{\text{A0}}}$
numbers, and all dimensional quantities, such as $L$, $\rho$, $v$, $B$, etc.,
can be scaled arbitrarily, as long as their chosen values leave
$\operatorname{\mathcal{M}}$ and $\operatorname{\mathcal{M}_{\text{A0}}}$
unchanged for a given simulation model.

Figure 2: The 3D geometry of magnetic field (blue) and velocity streamlines
(red) for the $\operatorname{\mathcal{M}}=20$,
$\operatorname{\mathcal{M}_{\text{A0}}}=10$ (top) and
$\operatorname{\mathcal{M}}=20$, $\operatorname{\mathcal{M}_{\text{A0}}}=0.1$
(bottom) simulations at $t=5\,T$. We see isotropic magnetic and velocity
structure for the super-Alfvénic turbulence, and highly-anisotropic magnetic
fields and vortex creation in the velocity streamlines for the sub-Alfvénic
turbulence. In incompressible flows the sub-Alfvénic turbulence becomes quasi
two-dimensional, but in our compressible supersonic simulations, shocks form
perpendicular to the vortices (Beattie & Federrath, 2020b), which sustain
parallel fluctuations in the magnetic field.

###  3.3 Magnetic Fields

The initial magnetic field, $\mathrm{{\mn@boldsymbol{\mathit{B}}}}(x,y,z)$ at
$t=0$, in Equations (6–8) is a uniform field with field lines threaded through
the $\hat{\mathrm{{\mn@boldsymbol{\mathit{z}}}}}$ direction of the
simulations. The total magnetic field is, as given in Equations 3 and 9 and
reiterated here for clarity,

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)=B_{0}\mathrm{{\mn@boldsymbol{\mathit{\hat{z}}}}}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t).$
(18)

Since we work in units of sound speed for velocity, and in units of $\rho_{0}$
for density, we chose to scale the units of the magnetic field for a typical
molecular cloud density (hydrogen number density of
$n_{\text{H}}=10^{3}\,\mathrm{cm}^{-3}$) and temperature ($10\,\mathrm{K}$)
corresponding to a sound speed of $c_{\mathrm{s}}=0.2\,\mathrm{km\,s}^{-1}$,
to

$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{B}}}}=\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{\text{sim}}\left(\frac{c_{s}}{0.2\,\text{km\,s}^{-1}}\right)\left(\frac{n_{\text{H}}}{10^{3}\,\text{cm}^{-3}}\right)^{1/2},$
(19)

where $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{\text{sim}}$ is the magnetic
field in simulation units, such that all magnetic fields are in units of
$\mu\text{G}$, as previously done in Mocz & Burkhart (2018), for example. We
list all of the magnetic field components in Table 1. Note that in our
simulations $0.01\lesssim B_{0}/(\mu\text{G})\lesssim 1000$, providing us with
$\sim$ five orders of magnitude to explore for the value of $B_{0}$ (column 4
in Table 1). This parameter space encompasses the wide variety of molecular
clouds that one might observe.

The fluctuating component of the field evolves self-consistently from the MHD
equations. However, we set $B_{0}$, and by using the definition of the Alfvén
velocity (defined previously in reference to Equation 1) and the turbulent
Mach number, we can set the target $\operatorname{\mathcal{M}_{\text{A0}}}$,

$\operatorname{\mathcal{M}_{\text{A0}}}=2c_{s}\sqrt{\pi\rho_{0}}\operatorname{\mathcal{M}}/B_{0}.$
(20)

We vary this value for each of the 24 simulations, spanning
$\operatorname{\mathcal{M}_{\text{A0}}}=0.1-100$, encompassing very weak and
strong mean magnetic fields. We extract 51 time realisations of the $B_{x}$,
$B_{y}$ and $B_{z}$ components of the total magnetic field between $5\leq
t/T\leq 10$ to ensure that the turbulence is in a statistically stationary
state (Federrath et al., 2009; Price & Federrath, 2010).

We show a slice of the magnitude of
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ in units of $B_{0}$, for each of
the 24 simulations at $t=5\,T$, in Figure 1. Red colours in the plot show
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|<B_{0}$, white
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|=B_{0}$ and black
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|>B_{0}$. In the sub-Alfvénic
regime we see $|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|<B_{0}$, in the
super-Alfvénic regime $|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|>B_{0}$
and $\operatorname{\mathcal{M}_{\text{A0}}}\approx 1$ marks the transition
between the two. Using the magnetic field components we construct the
variables $B_{\parallel}=B_{z}$ and $B_{\perp}=B_{x}=B_{y}$333Note that we are
defining these magnetic field quantities within the global mean magnetic field
frame, and not the local mean magnetic field frame. See Cho et al. (2002) for
the difference between the two frames.. We create PDFs for the two new
magnetic field variables, and average them over $5-10\,T$. In §5 we discuss
these PDFs in detail, but first we will discuss one of the main predictions
from our fluctuation model in Equation (2.1), namely that it is primarily the
velocity gradients parallel to the mean-field that create compressive modes in
the velocity field.

Figure 3: Top-left: The convergence of the velocity for a single $xz$ slice
from the M20MA0.1 simulation at $t=5\,T$. The direction of the mean magnetic
field is shown in the bottom left corner. We see highly
$-\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}>0$ (converging) structures
forming perpendicular to $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$. Fast
magnetosonic waves can also be seen travelling as compressive structures
perpendicular to $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$ (Tritsis &
Tassis, 2016). Top-right: The magnitude of the velocity gradient along the
magnetic field, $|\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}|$.
Equation (12) predicts that this plot and the magnitude of the convergence
should be the same as $\partial_{t}\delta B\rightarrow 0$. This is indeed
true, and the only significant difference between the two plots is the fast
magnetosonic waves, which will disappear as $B_{0}\gg\delta B$. Bottom-left:
The same as the top-left plot, but for a density slice, shown in units of mean
density, $\rho_{0}$. We see that the
$-\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}>0$ structures in the top
panel correspond to hydrodynamical shocks (Mocz & Burkhart, 2018; Beattie &
Federrath, 2020b), forming from mostly parallel velocity streams in the flow.
Bottom-right: Velocity streamlines integrated through the divergence of the
velocity, showing how parallel velocity channels feed into high-density
filaments, creating the velocity gradients parallel to the magnetic field.
Figure 4: Left: The perpendicular magnetic field fluctuations (Equation 23) as
a function of $\operatorname{\mathcal{M}}$, coloured by
$\operatorname{\mathcal{M}_{\text{A0}}}$. The model fits, $\delta
B_{\perp}=2c_{s}\sqrt{\pi\rho_{0}}C\operatorname{\mathcal{M}}$, (Equation 17)
are shown with dashed lines, also coloured by
$\operatorname{\mathcal{M}_{\text{A0}}}$. The $C$ parameter from the model is
shown in the inset plot, for each different
$\operatorname{\mathcal{M}_{\text{A0}}}$. We show a blue, dashed line for
$C=\operatorname{\mathcal{M}_{\text{A0}}}/2$, which is the strong magnetic
field model derived in F16 and written in terms of
$\operatorname{\mathcal{M}}$ and $\operatorname{\mathcal{M}_{\text{A0}}}$ in
Equation (25). Right: The same as the left panel but for the parallel magnetic
field fluctuations. Both plots reveal that there is a linear dependence of the
fluctuations on $\operatorname{\mathcal{M}}$, as predicted by our compressible
quasi-static model (Equation 2.1), and a close-to-linear dependence for
$\operatorname{\mathcal{M}_{\text{A0}}}<1$, i.e., for sub-Alfvénic mean-field
flows. Figure 5: The parallel (shown with triangles), perpendicular (shown
with circles) and total (shown with diamonds) fluctuations of the magnetic
field as a function of $\operatorname{\mathcal{M}}$, coloured by
$\operatorname{\mathcal{M}_{\text{A0}}}$. We plot the modified F16 magnetic
field fluctuation model, Equation (25), using dashed lines for
$\operatorname{\mathcal{M}_{\text{A0}}}=0.1$, $0.5$ and $1$. We show only the
sub-Alfvénic simulations where the modified F16 model is valid.

##  4 Compressible flows along the mean magnetic field

In §2 we showed how one could construct a compressive quasi-static model for
the magnetic field fluctuations (Equation 2.1) and that in the strong mean-
field regime the model reduces to a simple statement about compressibility and
velocity gradients along $B_{0}$, shown in Equation (12). This is a
significant result for astrophysical flows, since both strong, coherent mean
magnetic fields and velocity gradients perpendicular to filamentary structures
are measured in observations (Palmeirim et al., 2013; Beuther et al., 2015;
Federrath, 2016b; Shimajiri et al., 2019; Chen et al., 2020b). Here we show
some direct, qualitative evidence to support Equation (12).

We show a single time realisation of the full 3D geometry of the magnetic and
velocity (scaled by $c_{s}$) fields in Figure 2 for the M20MA10 (top) and
M20MA0.1 (bottom) simulations, where we refer to Table 1 for the respective
labels. M20MA10 shows an example of the tangled, isotropic magnetic fields
that we have in the high $\operatorname{\mathcal{M}_{\text{A0}}}$ simulations.
The magnetic field is dominated by
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, because the mean-field
component is too weak to impart any systematic, ordered structure on either
the velocity or total magnetic field. On the other hand, the M20MA0.1
simulation shows a highly-ordered, anisotropic magnetic field, dominated by
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$. Many of the velocity streamlines
form $\sim L/2$ eddies perpendicular to the magnetic field lines. If all of
the streamlines were organised into eddies the flow would be statistically the
same along $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$, hence becoming quasi
two-dimensional. However, we see some parallel velocity streams developing
towards the left corner of the plot. To explore this further we take an
$xz$-slice through this same time realisation.

We plot the convergence of the $xz$ slice of the velocity field for the
M20MA0.1 simulation in the top-left panel of Figure 3. Hydrodynamical shocks
that form perpendicular to the mean-field (see Mocz & Burkhart 2018 for the
shock-jump conditions, showing that they are hydrodynamical) can be clearly
identified as $-\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}>0$
structures, i.e., converging flows, by comparing the convergence with the
density structures in the bottom-left panel of Figure 3. Tracing the velocity
streamlines, shown in the bottom-right panel, we can see the most significant
shocks are coupled with a strong velocity channels along the mean magnetic
field. Furthermore, probing the magnitude of the velocity gradient along
$B_{0}$, shown in the top-left panel, shows a strong correlation between the
$-\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}>0$ structures and the
velocity gradient, consistent with Equation (12). The only difference between
the panels is the $\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
signatures of fast magnetosonic waves that form perpendicular to the field in
the convergence plot. However, these waves will cease as the mean-field
becomes strong and magnetic tension locks them in place and then
$|\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}|$ will indeed equal
$|\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}|$. To summarise,
what we show here is that as $B_{0}\gg\delta B$ parallel velocity gradients
cause compressive modes in the velocity field, which cause strong hydrodynamic
shocks to form across the field.

Indeed, this MHD phenomenon may contribute to the anisotropic cloud collapse
mechanism that is seen in sub-Alfvénic supersonic turbulence, simulations and
(most likely) observations (Tritsis et al., 2015; Mocz et al., 2017) since
velocity streams could be responsible for accreting material into dense, star-
forming filaments. Compressive modes that steepen velocity gradients could
even contribute to coherent velocity gradients and flows along and across the
filamentary structures in the clouds, like observed in star-forming and
infrared dark clouds (Palmeirim et al., 2013; Beuther et al., 2015; Shimajiri
et al., 2019) and investigated in great detail in Chen et al. (2020b, a). Chen
et al. (2020a) suggest that the velocity gradients are caused by gravitational
accretion onto a self-gravitating filament. Here we add that compressive modes
in highly-magnetised gases are coupled with velocity gradients along the
direction of the mean magnetic field, which are perpendicular to the densest
filaments, and parallel to striations (Beattie & Federrath, 2020b), and could
also contribute to the accretion process. Furthermore, Molina et al. (2012)
found, when measuring the density dispersion for sub-Alfvénic turbulence, that
more compressive modes appeared relative to the compressions induced by
turbulent forcing alone (i.e. the $b$ parameter that is an input in the
simulations). This may also be explained by the emergence of more compressive
modes when the turbulence is in this highly-magnetised state. Certainly this
result warrants detailed investigation in future studies and we will be
further exploring these velocity structures in an ensemble of higher-
resolution simulations in a forthcoming study on the two-point velocity
statistics of anisotropic, supersonic MHD turbulence.

To conclude this section we comment on the velocity gradient method, which is
summarised in Hu et al. (2019), and used to measure the
$\operatorname{\mathcal{M}_{\text{A}}}$ values that we report for the Gould
Belt in §1. This method has been significantly developed in the literature
(González-Casanova & Lazarian, 2017; Yuen & Lazarian, 2017b, a; Lazarian et
al., 2018, and references therein) and assumes that velocity gradients are
perpendicular to magnetic fields, which may be true for incompressible
$\operatorname{\mathcal{M}_{\text{A0}}}<1$ flows (i.e., through the two-
dimensionalisation due to the exponentially decaying parallel magnetic field
fluctuations; Verma 2017), but is not be strictly true for compressible flows,
as indicated in Equation (12) and the discussion above. Yuen & Lazarian
(2017a) attribute co-alignment of velocity and magnetic field vectors to
indicate gravitational collapse. We show that this need not be the case, and
for supersonic flows in a sub-Alfvénic mean-field regime, one can indeed find
strong velocity gradients parallel to the mean magnetic field, without the
presence of self-gravity. Next we discuss the magnetic field PDFs, and how the
magnetic field fluctuations depend upon $\operatorname{\mathcal{M}}$ and
$\operatorname{\mathcal{M}_{\text{A0}}}$.

##  5 Structure of the Magnetic Field

In this section we provide a detailed study of the magnetic field structure
for each of our 24 simulations, over a wide range of
$\operatorname{\mathcal{M}}$ and $\operatorname{\mathcal{M}_{\text{A0}}}$. We
first explore the magnitude and anisotropy of the fluctuations, followed by
mean-normalised fluctuations and finally the morphology and intermittency of
the magnetic field PDFs.

###  5.1 Magnetic field fluctuations

We construct the PDFs for both the $B_{\parallel}$ and $B_{\perp}$ components
of the magnetic field and average them over the time range $5-10\,T$. The
first moment and second moment of the PDFs have a physical interpretation. For
example, for the $B_{\parallel}$ PDF the first moment is,

$\displaystyle\mathrm{\left\langle
B_{\parallel}\right\rangle}=B_{0}=\int\operatorname{d}\\!{B_{\parallel}}\,B_{\parallel}p(B_{\parallel}),$
(21)

where $p(B_{\parallel})$ is the PDF for the parallel component of the magnetic
field, describes the mean-field value from Equation (20). Likewise, the
variance of the field is

$\displaystyle\mathrm{\left\langle
B_{\parallel}^{2}\right\rangle}-\mathrm{\left\langle
B_{\parallel}\right\rangle}^{2}=\delta
B_{\parallel}^{2}=\int\operatorname{d}\\!{B_{\parallel}}\,\left(B_{\parallel}-B_{0}\right)^{2}p(B_{\parallel}),$
(22)

which means the standard deviation of the field is exactly $\delta
B_{\parallel}$. Since the mean-field only has a $z$ component,
$\mathrm{\left\langle B_{\perp}\right\rangle}=0$, and hence the second moment
of the PDF is exactly equal to the fluctuations squared,

$\displaystyle\mathrm{\left\langle B_{\perp}^{2}\right\rangle}=\delta
B_{\perp}^{2}=\int\operatorname{d}\\!{B_{\perp}}\,B_{\perp}p(B_{\perp}).$ (23)

Less important in our study is the first moment of these PDFs, since we set
this as an input parameter in our simulations, as described in §3. Hence we
focus our analysis on the standard deviation of these distributions, which
tells us information about the magnitude of the total fluctuations across all
(perpendicular or parallel to the mean guide field) length scales in the
simulations. For a scale-dependent analysis of the fluctuations one would need
to calculate a two-point statistic, which we plan to explore in detail, in a
future study.

####  5.1.1 $\delta B_{\perp}$ and $\delta B_{\parallel}$ fluctuations

In Figure 4 we show the perpendicular and parallel fluctuation amplitudes of
the magnetic field as a function of $\operatorname{\mathcal{M}}$, coloured by
$\operatorname{\mathcal{M}_{\text{A0}}}$. Both types of fluctuations show a
linear dependence in $\operatorname{\mathcal{M}}$, as predicted by our order-
of-magnitude estimate of the fluctuations in Equation (17). Conceptually, the
$\operatorname{\mathcal{M}}$ dependence can be explained by the turbulent
motions contributing to tangling, twisting and perturbing the magnetic fields.
More quantitatively, in our quasi-static model for the fluctuations, shown in
Equation (2.1), each of the RHS terms depend linearly on the velocity field,
regardless of the relation between $B_{0}$ and $\delta B$, hence the linear
dependence upon $\operatorname{\mathcal{M}}$ is a somewhat universal feature
for the amplitude of the fluctuations. For a constant
$\operatorname{\mathcal{M}_{\text{A0}}}$, we fit Equation (17),

$\delta B=2c_{s}\sqrt{\rho_{0}\pi}C\operatorname{\mathcal{M}},$ (24)

where $C$ is the slope parameter. Note that we need not include an offset in
the fits, since when $\operatorname{\mathcal{M}}\rightarrow 0$ there is no
source for the magnetic field fluctuations, hence $\delta B\rightarrow 0$.
Clearly $C$ is a function of $\operatorname{\mathcal{M}_{\text{A0}}}$, which
is shown by how the slopes change for different
$\operatorname{\mathcal{M}_{\text{A0}}}$ in Figure 4, and plotted in the inset
figure for each panel. For $\operatorname{\mathcal{M}_{\text{A0}}}\lesssim 1$
we see an approximately linear dependence,
$C\approx\operatorname{\mathcal{M}_{\text{A0}}}/2$, for both types of
fluctuations, until reaching a turnover at
$\operatorname{\mathcal{M}_{\text{A0}}}\approx 1-2$, where reducing the mean-
field strength decreases the $C$ parameter. This is because as
$|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|$ shrinks the total magnetic field
energy strength also shrinks, and thus the magnitude of the fluctuations
decreases, at least until the turbulent small-scale dynamo mechanism can
exponentially grow the fluctuations (Batchelor, 1950; Kazantsev, 1968;
Schekochihin et al., 2002; Federrath et al., 2014; Seta et al., 2015; Seta et
al., 2020, and references therein), which is a regime we do not consider in
this study (but is studied in F16). Next we explain the linear dependence of
$\operatorname{\mathcal{M}_{\text{A0}}}$ in the strong-field regime.

####  5.1.2 Strong-field Model

Using the definition of the Alfvén Mach number,
$\operatorname{\mathcal{M}_{\text{A0}}}=V/V_{A0}$, the turbulent Mach number,
$\operatorname{\mathcal{M}}=V/c_{s}$, and the Alfvén velocity,
$V_{A0}=B_{0}/\sqrt{4\pi\rho_{0}}$ one can show that the strong-field F16
model (Equation 4) can be rewritten in terms of $\operatorname{\mathcal{M}}$,

$\displaystyle\delta
B=c_{s}\sqrt{\pi\rho_{0}}\operatorname{\mathcal{M}_{\text{A0}}}\operatorname{\mathcal{M}}.$
(25)

This modified version of the F16 model shows a linear dependence in
$\operatorname{\mathcal{M}}$, in agreement with our model derived from the
induction equation, and linear in $\operatorname{\mathcal{M}_{\text{A0}}}$.
Hence, the $C$ obtained is analytical in this regime,
$C=\operatorname{\mathcal{M}_{\text{A0}}}/2$. This means, as suggested in the
derivation of Equation (17), that $C$ can be thought of as the Alfvénic
control parameter, i.e., it changes based on the type of Alfvénic turbulence,
which is shown clearly in the inset plots of Figure 4.

We plot Equation (25) on the magnitudes of the perpendicular, parallel and
total fluctuations of the magnetic field in Figure 5, where we show only data
from the trans- to sub-Alfvénic regimes, since this is the regime where the
strong-field model is valid. In general, we find that the model captures the
nature of the fluctuations well, tracking either values between $B_{\perp}$
and $B_{\parallel}$ in the sub-Alfvénic simulations, or very close to the
total fluctuations in the trans-Alfvénic simulation,
$\operatorname{\mathcal{M}_{\text{A0}}}\sim 1$. Plotting both perpendicular
and parallel fluctuations in Figure 5 reveals the difference between their
magnitudes due to the anisotropy inherent to MHD turbulence with a mean-field
present (Goldreich & Sridhar, 1995; Lazarian & Vishniac, 1999; Cho & Vishniac,
2000; Cho & Lazarian, 2003; Boldyrev, 2006; Kowal et al., 2007; Burkhart et
al., 2014, 2015). In the following, we explore the anisotropy of the
fluctuations.

Figure 6: The ratio between $\delta B_{\parallel}$ and $\delta B_{\perp}$ as a
function of $\operatorname{\mathcal{M}_{\text{A0}}}$, coloured by
$\operatorname{\mathcal{M}}$. The anisotropy in the fluctuations is largest
for small $\operatorname{\mathcal{M}_{\text{A0}}}$, and disappears for large
$\operatorname{\mathcal{M}_{\text{A0}}}$, where $\delta B_{\parallel}/\delta
B_{\perp}=1/\sqrt{2}$ is completely isotropic, indicated by the blue, dashed
line.

###  5.2 Anisotropy of the magnetic fluctuations

In Figure 6 we show the ratio between the $B_{\parallel}$ and $B_{\perp}$
fluctuations. Since $\delta B_{\perp}=\sqrt{\delta B_{x}^{2}+\delta
B_{y}^{2}}$ when $\delta B_{x}=\delta B_{y}=\delta B_{\parallel}$ the ratio
$\delta B_{\parallel}/\delta B_{\perp}=1/\sqrt{2}$. For
$\operatorname{\mathcal{M}_{\text{A0}}}\gtrsim 2$ we see the ratio
monotonically approaching the $1/\sqrt{2}$ isotropic limit. We can understand
this from our model in Equation (2.1). When $\delta B\gg B_{0}$ we can
disregard the $\mathcal{O}(B_{0})$ terms in Equation (2.1),

$\displaystyle\frac{D}{Dt}\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}=$
$\displaystyle(\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}),$
(26)

where $\frac{D}{Dt}\equiv\left(\frac{\partial}{\partial
t}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla\right)$, which no longer
has the anisotropic term,
$B_{0}\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}$. This means
that each of the components of $\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
are affected by the advection and compression equally, causing isotropic
fluctuations, without any preferential direction (which we visualised in the
top panel of Figure 2). However, in the opposite case, where $B_{0}\gg\delta
B$, the $B_{0}\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
becomes large and we measure $\delta B_{\parallel}$ up to a factor 2 larger
than $\delta B_{\perp}$. This is because the
$B_{0}\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ from our
quasi-static model causes strong compressions across the field, which
increases the magnetic field fluctuations along the mean-field direction, as
shown in Figure 3, and discussed in detail in §2 and §4. The impact that these
shocks have on the parallel fluctuations is quite phenomenal, and we see a
significant growth in the absolute value of the anisotropy between
$\operatorname{\mathcal{M}_{\text{A0}}}=0.5$, which has stronger $B_{\perp}$
fluctuations, and $\operatorname{\mathcal{M}_{\text{A0}}}=0.1$.

An important feature seen from Figure 6 is that the ratio between magnetic
field fluctuations divides out most of the $\operatorname{\mathcal{M}}$
dependence, which we can see from each of the different
$\operatorname{\mathcal{M}}$ simulations falling approximately upon each
other. This suggests that the anisotropy in the magnetic field fluctuations
are controlled mostly by $\operatorname{\mathcal{M}_{\text{A0}}}$, although
some $\operatorname{\mathcal{M}}$ dependence is certainly present in the
$\operatorname{\mathcal{M}_{\text{A0}}}=0.1$ simulations. This motivates
another avenue of enquiry, the $B_{0}$ normalised fluctuations, which should
also be $\operatorname{\mathcal{M}}$ independent according to Equation (14).
We explore this in the next section.

Figure 7: The mean-field normalised fluctuations, $\delta B/B_{0}$, as a
function of $\operatorname{\mathcal{M}_{\text{A0}}}$, coloured by
$\operatorname{\mathcal{M}}$. The vertical blue line indicates when $(\delta
B/B_{0})=1$. The solid black line is the strong-field model, Equation ( 25),
and the dashed black line is our fit, $(\delta B/B_{0})=\mathcal{M}^{0.69\pm
0.05}_{\text{A0}}$, fit on $2\leq\operatorname{\mathcal{M}_{\text{A0}}}\leq
100$ data.

###  5.3 Mean-field normalised fluctuations

We show the mean-field normalised fluctuations in Figure 7 as a function of
$\operatorname{\mathcal{M}_{\text{A0}}}$, coloured by
$\operatorname{\mathcal{M}}$, for the total magnetic field fluctuations. We do
not find a large difference in the two fluctuating components of the field
(see Appendix C), so we focus our analysis on the total fluctuations. As we
saw in the Figure 6, there is only a very weak $\operatorname{\mathcal{M}}$
dependence left, which we show explicitly in Figure 10 in the Appendix, which
is the same as Figure 4, but for the mean normalised quantities). This plot
shows variations in $\operatorname{\mathcal{M}}$, regardless of the value for
$\operatorname{\mathcal{M}_{\text{A0}}}$, will not influence $\delta B/B_{0}$.
We can see why this is the case in our compressible quasi-static model, since
$\delta
B/B_{0}=C(\operatorname{\mathcal{M}_{\text{A0}}})\operatorname{\mathcal{M}_{\text{A0}}}$
as shown in Equations 14 – 17.

One of the key features of Figure 7 is the distinct kink at $\delta B\approx
B_{0}$ (indicated by the horizontal blue line in Figure 7), which happens at
$\operatorname{\mathcal{M}_{\text{A0}}}\approx 2$. The same transition can be
seen in Figure 1, where fluctuations are $\mathcal{O}(B_{0})$, and are
becoming more isotropic in nature, which is shown in Figure 6. Clearly this
marks a critical transition for the magnetic field evolution, where all of the
terms in our quasi-static model play a role, so very little can be deduced
using our model in this transition region. Instead, we propose a semi-
analytical model that encapsulates both the $\delta B/B_{0}<1$ and $\delta
B/B_{0}>1$ regimes,

$\displaystyle\frac{\delta B}{B_{0}}=$ $\displaystyle
C\operatorname{\mathcal{M}_{\text{A0}}},$ (27) $\displaystyle C=$
$\displaystyle\left\\{\begin{matrix}\operatorname{\mathcal{M}_{\text{A0}}}/2,\,\text{for}\;\delta
B/B_{0}<1,\\\\[5.0pt]
\alpha\mathcal{M}^{\beta}_{\text{A0}},\,\text{for}\;\delta B/B_{0}\geq
1,\end{matrix}\right.$ (28)

where $\alpha$ and $\beta$ are fit parameters. For $\delta B/B_{0}<1$ we use
the F16 model, but since no analytical model is known in the $\delta
B/B_{0}\geq 1$ regime, and our compressible quasi-static model does not give
much insight, we instead directly fit to the data. Using a linear least-
squares fit we determine $\alpha\approx 1$ and $1+\beta=0.69\pm 0.05$. This is
somewhat consistent with the qualitative model for the “intermediate regime"
in F16 ($\delta B\sim B_{0}^{1/3}\implies(\delta B/B_{0})\sim
B_{0}^{-2/3}\sim\mathcal{M}^{2/3}_{\text{A0}}$), but here we explicitly show
that the mean-field weighted fluctuations are independent of
$\operatorname{\mathcal{M}}$, regardless of the ratio between $\delta B$ and
$B_{0}$. Hence the relation is,

$\displaystyle\frac{\delta B}{B_{0}}=$
$\displaystyle\left\\{\begin{matrix}\mathcal{M}^{2}_{\text{A0}}/2,\,\text{for}\;\delta
B/B_{0}<1,\\\\[5.0pt] \mathcal{M}^{0.69\pm
0.05}_{\text{A0}},\,\text{for}\;\delta B/B_{0}\geq 1,\end{matrix}\right.$ (29)

which we show with the solid and dashed black line in Figure 7, respectively.
Indeed, we expect that our model will hold until the dynamo regime is reached
as $\delta B/B_{0}\gg 1$ and for all $\operatorname{\mathcal{M}}$, making it a
universal scaling relation for anisotropic, supersonic MHD turbulence. This is
the key result from this investigation.

The relative scatter around our relation is clearly set by
$\operatorname{\mathcal{M}}$, and becomes larger as the mean field weakens.
Calculating the standard deviation for each of the fixed
$\operatorname{\mathcal{M}_{\text{A0}}}$ values in Figure 7 we find that the
relative scatter changes systematically with
$\operatorname{\mathcal{M}_{\text{A0}}}$, ranging between $\approx 5-20\,\%$.
Now that we have discussed in great detail the absolute values, the anisotropy
and the mean-field weighted $\delta B$ we move on to give a deeper, and more
complete analysis of the magnetic field PDFs.

Figure 8: Variance ($=1$) and mean ($=0$) normalised PDFs of $B_{\perp}$ (left
column) and $B_{\parallel}$ (right column), for each of the 24 simulations
revealing the morphology of the distributions. Each panel shows simulations
with different $\operatorname{\mathcal{M}}$ values, and within each panel
shows the ensemble with different $\operatorname{\mathcal{M}_{\text{A0}}}$
values, with the legend shown at the top. We normalise by $\sigma_{B}=\delta
B$ values for each of the distributions because it makes it easier to compare
the shapes of the distributions without focusing on the amplitudes of the
fluctuations, which are discussed in detail in §5.1.1. We plot a standard
Gaussian, $\mathcal{N}(0,1)$, in red for comparative purposes. $1\sigma$
fluctuations are shown with the bands around the $5-10\,T$ time average.

###  5.4 Magnetic field distribution and intermittency

Magnetic field fluctuations are highly intermittent (Bruno et al., 2007; Seta
et al., 2020, specifically see §3 in Seta et al. 2020), i.e., the fluctuations
cannot be aptly described by the variance of the magnetic field distribution
that we have described in detail in the previous sections. In this section,
however, we systematically describe the shape and features of the
distributions for the $B_{\perp}$ and $B_{\parallel}$ components of the
turbulent magnetic field, for each of the simulations. We show the time-
averaged distributions in Figure 8. Note that we plot the distributions of the
normalised variable, $(B-\mathrm{\left\langle B\right\rangle})/\sigma_{B}$,
which enforces that all of the distributions have a mean of zero and variance
of one. This allows us to compare the shape of the distributions without being
obscured by the absolute magnitude of the magnetic field fluctuations, which
were studied in §5.1.1. We now consider the $B_{\perp}$ and $B_{\parallel}$
PDFs in the following two subsections.

####  5.4.1 $B_{\perp}$ PDF

In the left column of Figure 8 we show the perpendicular magnetic field
distributions. The distributions are reasonably symmetric, for all
$\operatorname{\mathcal{M}_{\text{A0}}}$ and $\operatorname{\mathcal{M}}$.
Indeed, as described earlier in §3.3, this is because the Lorentz force acts
symmetrically about the mean magnetic field, defining an axis of symmetry for
the fluctuations (Cho et al., 2002). We plot a standardised Normal
distribution, $\mathcal{N}(0,1)$, shown in dashed-red in each of the panels to
compare the distributions with purely Gaussian fluctuations. The most Gaussian
fluctuations are found in the $\operatorname{\mathcal{M}}=2$ simulations where
$\delta B\gtrsim B_{0}$, i.e., for
$\operatorname{\mathcal{M}_{\text{A0}}}\approx 2-10$, as shown in Figure 7.
This marks the transition between the
low-$\operatorname{\mathcal{M}_{\text{A0}}}$ and $\delta B\gg B_{0}$ flows.
The key difference between $\mathcal{N}(0,1)$ and the other distributions are
the long, extended tails and peaked mode, shown in both the low- and
high-$\operatorname{\mathcal{M}_{\text{A0}}}$ simulations for
$\operatorname{\mathcal{M}}=2$, and then present in all of the distributions
for $\operatorname{\mathcal{M}}\geq 4$. This is a signature of intermittency.

In the low-$\operatorname{\mathcal{M}_{\text{A0}}}$ simulations, where the
magnetic field is extremely strong and under large amounts of tension, which
scales as $1/\mathcal{M}^{2}_{\text{A0}}$, the intermittent tails of the
distributions can be associated with only large magnetic field perturbations.
These could be from strong hydrodynamical shocks, which are intermittent
events in the velocity field. Shock number densities (in a given volume, for
example) are a function of $\operatorname{\mathcal{M}}$, with more shocks and
with more power in each shock, being created for higher values of
$\operatorname{\mathcal{M}}$ (Gotoh, 1994; Girimaji & Zhou, 1995; Beattie &
Federrath, 2020a). This is why as the $\operatorname{\mathcal{M}}$ increases,
we see more intermittency in all of the distributions, regardless of
$\operatorname{\mathcal{M}_{\text{A0}}}$. In the
high-$\operatorname{\mathcal{M}_{\text{A0}}}$ simulations, like the
$\operatorname{\mathcal{M}_{\text{A0}}}=100$ ensemble, the mean-field is
extremely weak compared to the fluctuations. Hence these distributions are
only composed of the fluctuating magnetic field, which is purely intermittent
in nature. Next we discuss the $B_{\parallel}$ PDF.

Figure 9: Magnetic eddies shown in an $xy$ slice of the $B_{\parallel}$
component of the magnetic field for the M2MA0.1 simulation. Velocity
streamlines, with the thickness of the line proportional to the magnitude of
the velocity, are overlaid on the plot to reveal the eddy structure. We
subtract the mean field to reveal the regions of negative $\delta
B_{\parallel}$, corresponding to the extended log-linear tails in the
$B_{\parallel}$ PDFs from Figure 8, discussed in §5.4.2.

####  5.4.2 $B_{\parallel}$ PDF

In the right column of Figure 8 we show the parallel magnetic field PDFs. In
the high-$\operatorname{\mathcal{M}_{\text{A0}}}$ regime, the fluctuations are
mostly symmetrical about the mean, similar to the $B_{\perp}$ PDFs described
previously. However, for the sub-Alfvénic flows we see extended $\log$-linear
tails developing in the negative values of the PDFs, most obviously shown in
the top, right panel, for the M2MA0.1 simulation. This extended tail
corresponds to values of the total magnetic field, which are less than
$B_{0}$, i.e., $B_{\parallel}<B_{0}$, where $B_{\parallel}=B_{0}+\delta
B_{\parallel}$, are the parallel components of the magnetic field from
Equation (9). Clearly this implies that preferentially $\delta B_{\parallel}$
opposes the mean field, reducing the total magnetic field where the
fluctuation is present, i.e., in these local regions, the parallel components
of the vector that describes the fluctuation must have the form
$B_{\parallel}=B_{0}-\delta B_{\parallel}$. The tail is most prominent in the
morphology of the low-$\operatorname{\mathcal{M}}$ simulations.

By creating a mask that reveals just the magnetic field values in the tail of
the $B_{\parallel}$ PDF we find that eddies are responsible for creating
localised magnetic voids in the field. We show an example of the magnetised
eddies in the M2MA0.1 simulation, where we plot $B_{\parallel}-B_{0}$ for an
$xy$ slice through the magnetic field, overlaid with velocity streamlines in
Figure 9. The eddies form low-density and magnetic pressure regions in the
fluid, and with a fixed $B_{0}$ this means that the fluctuating component of
the field must oppose the mean-field to reduce the local magnetic pressure,
giving rise to the extended $B_{\parallel}<B_{0}$ tails in the PDFs. However,
as the magnetic fields become more isotropic with increasing
$\operatorname{\mathcal{M}}$ (Figure 6) the tails become less prominent as the
low pressure regions are mixed through the flow.

The remaining structure in the distribution is then most likely due to the
turbulent fluctuations, which we showed in §5.1.1 are linear in
$\operatorname{\mathcal{M}}$. Hence for large $\operatorname{\mathcal{M}}$
values we mostly see the velocity interactions, and the interaction between
$\delta B_{\perp}$ and $\delta B_{\parallel}$ no longer dominates the
morphology of the PDFs.

##  6 Summary and Key Findings

In this study, we explore the amplitude of magnetic field fluctuations that
are perpendicular, $\delta B_{\perp}$, and parallel, $\delta B_{\parallel}$,
to the mean magnetic field, $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$, in an
ensemble of 24 simulations, across a large range of Alfvénic mean-field Mach
numbers, $\operatorname{\mathcal{M}_{\text{A0}}}=0.1-100$ and root-mean-
squared turbulent Mach numbers, $\operatorname{\mathcal{M}}=2-20$,
encompassing realistic values for molecular clouds, the birthplaces of stars.
First, we derive a new compressible quasi-static fluctuation model, which we
use to explain much of the phenomena discussed in this study and show how,
even in the strong mean magnetic field regime the flow remains three-
dimensional. Next we explore how the nature of the fluctuations perpendicular
and parallel to the mean magnetic field change as a function of
$\operatorname{\mathcal{M}}$ and $\operatorname{\mathcal{M}_{\text{A0}}}$, and
generalise an analytical model for the fluctuations in the strong mean-field
regime. This is followed by an investigation of the anisotropy between $\delta
B_{\perp}$ and $\delta B_{\parallel}$, and the mean-field normalised
fluctuations. Finally, we explore the shape of the probability density
functions of magnetic fluctuations. We summarise the key findings below:

  * •

We derive a compressible quasi-static model for the magnetic field
fluctuations, shown in Equation (2.1) with full derivation in Appendix A,
which predicts that for sub-Alfvénic flows, compressions (e.g., shocks) in the
velocity field are associated with gradients along the mean magnetic field,

$|\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}|=|\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}|,$

which could be a contributing factor for anisotropic cloud collapse (Tritsis
et al., 2015; Mocz et al., 2017, etc.) and coherent velocity structures
perpendicular to the principle axis of filaments, observed in real molecular
clouds (Chen et al., 2020b). We use the non-dimensional form of the equation
to predict that the magnitude of the magnetic field fluctuations are linear in
$\operatorname{\mathcal{M}}$,

$\delta
B=2c_{s}\sqrt{\pi\rho_{0}}C(\operatorname{\mathcal{M}_{\text{A0}}})\operatorname{\mathcal{M}},$

where $c_{s}$ is the sound speed, $\rho_{0}$ is the mean density and $C$ is
the proportionality factor, which changes between the $B_{0}\gg\delta B$ and
$B_{0}\ll\delta B$ regimes. We show in Figure 4 that the linear dependence
predicted for all $\operatorname{\mathcal{M}_{\text{A0}}}$ holds true.

  * •

We plot the velocity divergence and velocity streamlines for one of the
$\operatorname{\mathcal{M}_{\text{A0}}}=0.1$ simulations in Figure 3,
revealing the compressive modes coupled with parallel velocity gradients,
$|\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}|=|\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}|$.
We show that, even in the absence of self-gravity, shocks that form
perpendicular to the mean magnetic field turn into converging flows that are
accreted upon by coherent parallel velocity channels, consistent with our
compressible quasi-static model.

  * •

Our new fluctuation model reduces to the Federrath (2016a) analytical model
for strong-field magnetic field fluctuations by setting
$C=\operatorname{\mathcal{M}_{\text{A0}}}/2$. We further rewrite the Federrath
(2016a) model to incorporate the explicit dependence upon
$\operatorname{\mathcal{M}}$,

$\delta
B=c_{s}\sqrt{\pi\rho_{0}}\operatorname{\mathcal{M}_{\text{A0}}}\operatorname{\mathcal{M}},$

and we show in Figure 5 that it is in good agreement with the measured
fluctuations for the sub-Alfvénic simulations, across all
$\operatorname{\mathcal{M}}$.

  * •

We find the ratio between the $\delta B_{\perp}$ and $\delta B_{\parallel}$
components of the magnetic field are similar for different
$\operatorname{\mathcal{M}}$, but do depend upon
$\operatorname{\mathcal{M}_{\text{A0}}}$ (shown in Figure 6). As
$\operatorname{\mathcal{M}_{\text{A0}}}$ gets larger the fluctuations become
isotropic, but for $\operatorname{\mathcal{M}_{\text{A0}}}\sim 0.1$ $\delta
B_{\parallel}\sim 2\delta B_{\perp}$. This is because the strong velocity
gradients along the field, when $B_{0}\gg\delta B$, leads to large values of
$\delta B_{\parallel}$. This is a key distinguishing feature from
incompressible MHD turbulence with a strong mean magnetic field, where
parallel fluctuations exponentially decay and the flow becomes quasi two-
dimensional (Alexakis, 2011; Verma, 2017). We show that the flow remains very
much three-dimensional for compressible MHD turbulence.

  * •

We show that the mean-field normalised fluctuations, $\delta B/B_{0}$, are
independent of $\operatorname{\mathcal{M}}$ in Figure 7. We propose a semi-
analytic model using the Federrath (2016a) model in the strong-field regime
and a least-squares fit in the weak-field regime,

$\frac{\delta
B}{B_{0}}=\left\\{\begin{matrix}\operatorname{\mathcal{M}_{\text{A0}}}^{2}/2,\,\delta
B/B_{0}<1,\\\\[5.0pt] \mathcal{M}^{0.69\pm 0.05}_{\text{A0}},\,\delta
B/B_{0}\geq 1,\end{matrix}\right.$

which is valid until the turbulence reaches the dynamo growth regime, $(\delta
B/B_{0})\gg 1$, and independent of $\operatorname{\mathcal{M}}$, making it a
universal feature of anisotropic, supersonic MHD turbulence.

  * •

We calculate the time-averaged ($5-10\,T$, where $T$ is the large-scale
turbulent eddy turnover time) $B_{\perp}$ and $B_{\parallel}$ distributions,
shown in Figures 8, and discuss their morphology in §5.4. We find distinct
signatures of intermittency in most of the distributions and an extended tail
into the negative values for the $B_{\parallel}$ distribution. This
corresponds to $B_{\parallel}$ fluctuations opposing the mean-field, which we
show is caused by eddies that cause local, low-magnetised pressure regions in
the fluid.

Data analysis and visualisation software used in this study: numpy (Oliphant,
2006), matplotlib (Hunter, 2007), cython (Behnel et al., 2011), visit (Childs
et al., 2012), scipy (Virtanen et al., 2020).

## Acknowledgements

We thank the reviewer, Pierre Lesaffre, for providing a thorough and
constructive review of our study. We thank the Australian National University
and Research School of Astronomy and Astrophysics for the support in place for
staff and students during the COVID19 isolation, which is when most of this
study was written. J. R. B. acknowledges funding from the Australian National
University, specifically the Deakin PhD and Dean’s Higher Degree Research
(theoretical physics) Scholarships. C. F. acknowledges funding provided by the
Australian Research Council (Discovery Project DP170100603, and 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, pr48pi and GCS Large-scale project 10391), the Partnership for
Advanced Computing in Europe (PRACE grant pr89mu), the Australian National
Computational Infrastructure (grant ek9), and the Pawsey Supercomputing Centre
with funding from the Australian Government and the Government of Western
Australia, in the framework of the National Computational Merit Allocation
Scheme and the ANU Allocation Scheme. The simulation software FLASH was in
part developed by the DOE-supported Flash Centre for Computational Science at
the University of Chicago.

## Data Availability

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

## References

  * Alexakis (2011) Alexakis A., 2011, Phys. Rev. E, 84, 056330
  * Ballesteros-Paredes et al. (2007) Ballesteros-Paredes J., Klessen R. S., Mac Low M.-M., Vazquez-Semadeni E., 2007, Protostars and Planets V, pp 63–80
  * Batchelor (1950) Batchelor G. K., 1950, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 201, 405
  * Beattie & Federrath (2020a) Beattie J. R., Federrath C., 2020a, MNRAS, submitted
  * Beattie & Federrath (2020b) Beattie J. R., Federrath C., 2020b, MNRAS, 492, 668
  * Beattie et al. (2019) Beattie J. R., Federrath C., Klessen R. S., Schneider N., 2019, MNRAS, 488, 2493
  * Behnel et al. (2011) Behnel S., Bradshaw R., Citro C., Dalcin L., Seljebotn D. S., Smith K., 2011, Computing in Science & Engineering, 13, 31
  * Beresnyak (2019) Beresnyak A., 2019, Living Reviews in Computational Astrophysics, 5, 2
  * Bertrang et al. (2014) Bertrang G., Wolf S., Das H. S., 2014, A&A, 565, A94
  * Beuther et al. (2015) Beuther H., Ragan S. E., Johnston K., Henning T., Hacar A., Kainulainen J. T., 2015, A&A, 584, A67
  * 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
  * Bruno et al. (2007) Bruno R., Carbone V., Chapman S., Hnat B., Noullez A., Sorriso-Valvo L., 2007, Physics of Plasmas, 14, 032901
  * Brunt et al. (2009) Brunt C. M., Heyer M. H., Mac Low M. M., 2009, A&A, 504, 883
  * Burkhart et al. (2014) Burkhart B., Lazarian A., Leão I. C., de Medeiros J. R., Esquivel A., 2014, ApJ, 790, 130
  * Burkhart et al. (2015) Burkhart B., Lee M.-Y., Murray C. E., Stanimirović S., 2015, ApJ, 811, L28
  * Chen et al. (2020a) Chen C.-Y., Mundy L. G., Ostriker E. C., Storm S., Dhabal A., 2020a, MNRAS,
  * Chen et al. (2020b) Chen M. C.-Y., et al., 2020b, The Astrophysical Journal, 891, 84
  * Chevance et al. (2020) Chevance M., et al., 2020, arXiv e-prints, p. arXiv:2004.06113
  * Childs et al. (2012) Childs H., et al., 2012, in , High Performance Visualization–Enabling Extreme-Scale Scientific Insight. pp 357–372
  * Cho & Lazarian (2003) Cho J., Lazarian A., 2003, MNRAS, 345, 325
  * Cho & Vishniac (2000) Cho J., Vishniac E. T., 2000, ApJ, 539, 273
  * Cho et al. (2002) Cho J., Lazarian A., Vishniac E. T., 2002, ApJ, 564, 291
  * Clark & Hensley (2019) Clark S. E., Hensley B. S., 2019, ApJ, 887, 136
  * Clark et al. (2015) Clark S. E., Hill J. C., Peek J. E. G., Putman M. E., Babler B. L., 2015, Phys. Rev. Lett., 115, 241302
  * Colling et al. (2018) Colling C., Hennebelle P., Geen S., Iffrig O., Bournaud F., 2018, A&A, 620, A21
  * Cox et al. (2016) Cox N. L. J., et al., 2016, A&A, 590, A110
  * 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
  * Elmegreen & Scalo (2004) Elmegreen B. G., Scalo J., 2004, ARA&A, 42, 211
  * Federrath (2013) Federrath C., 2013, MNRAS, 436, 1245
  * Federrath (2015) Federrath C., 2015, MNRAS, 450, 4035
  * Federrath (2016a) Federrath C., 2016a, Journal of Plasma Physics, 82, 535820601
  * Federrath (2016b) Federrath C., 2016b, MNRAS, 457, 375
  * Federrath (2018) Federrath C., 2018, Physics Today, 71, 38
  * Federrath & Banerjee (2015) Federrath C., Banerjee S., 2015, MNRAS, 448, 3297
  * Federrath & Klessen (2012) Federrath C., Klessen R. S., 2012, ApJ, 761
  * Federrath & Klessen (2013) Federrath C., Klessen R. S., 2013, ApJ, 763, 51
  * 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. (2014) Federrath C., Schober J., Bovino S., Schleicher D. R. G., 2014, ApJ, 797, L19
  * Federrath et al. (2016a) Federrath C., Klessen R. S., Iapichino L., Hammer N. J., 2016a, preprint,  (arXiv:1607.00630)
  * Federrath et al. (2016b) Federrath C., et al., 2016b, ApJ, 832, 143
  * 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
  * Frank et al. (2014) Frank A., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 451 (arXiv:1402.3553), doi:10.2458/azu_uapress_9780816531240-ch020
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Gerrard et al. (2019) Gerrard I. A., Federrath C., Kuruwita R., 2019, MNRAS, 485, 5532
  * Girimaji & Zhou (1995) Girimaji S. S., Zhou Y., 1995, Technical report, Spectrum and energy transfer in steady burgers turbulence
  * Goldreich & Sridhar (1995) Goldreich P., Sridhar S., 1995, ApJ, 438, 763
  * González-Casanova & Lazarian (2017) González-Casanova D. F., Lazarian A., 2017, ApJ, 835, 41
  * Gotoh (1994) Gotoh T., 1994, Physics of Fluids, 6, 3985
  * Grisdale et al. (2017) Grisdale K., Agertz O., Romeo A. B., Renaud F., Read J. I., 2017, MNRAS, 466, 1093
  * Hennebelle & Falgarone (2012) Hennebelle P., Falgarone E., 2012, Astronomy and Astrophysics Review, 20
  * Hennebelle & Inutsuka (2019) Hennebelle P., Inutsuka S.-i., 2019, Frontiers in Astronomy and Space Sciences, 6, 5
  * Hennebelle et al. (2011) Hennebelle P., Commerçon B., Joos M., Klessen R. S., Krumholz M., Tan J. C., Teyssier R., 2011, A&A, 528, A72
  * Heyer & Brunt (2004) Heyer M. H., Brunt C. M., 2004, ApJ, 615, L45
  * Hu et al. (2019) Hu Y., et al., 2019, Nature Astronomy, 3, 776
  * 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
  * Kazantsev (1968) Kazantsev A. P., 1968, Soviet Journal of Experimental and Theoretical Physics, 26, 1031
  * Körtgen et al. (2017) Körtgen B., Federrath C., Banerjee R., 2017, MNRAS, 472, 2496
  * Kowal et al. (2007) Kowal G., Lazarian A., Beresnyak A., 2007, ApJ, 658, 423
  * 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 & McKee (2005) Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250
  * Kuruwita & Federrath (2019) Kuruwita R. L., Federrath C., 2019, Monthly Notices of the Royal Astronomical Society, 486, 3647
  * Kuruwita et al. (2017) Kuruwita R. L., Federrath C., Ireland M., 2017, MNRAS, 470, 1626
  * Kuruwita et al. (2020) Kuruwita R., Federrath C., Haugbølle T., 2020, arXiv e-prints, p. arXiv:2004.07523
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Lazarian & Vishniac (1999) Lazarian A., Vishniac E. T., 1999, ApJ, 517, 700
  * Lazarian et al. (2018) Lazarian A., Yuen K. H., Ho K. W., Chen J., Lazarian V., Lu Z., Yang B., Hu Y., 2018, ApJ, 865, 46
  * 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
  * Mac Low & Klessen (2004) Mac Low M. M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  * Malinen et al. (2016) Malinen J., et al., 2016, MNRAS, 460, 1934
  * Mandal et al. (2020) Mandal A., Federrath C., Körtgen B., 2020, Monthly Notices of the Royal Astronomical Society, 493, 3098
  * McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
  * Mocz & Burkhart (2018) Mocz P., Burkhart B., 2018, MNRAS, 480, 3916
  * Mocz et al. (2017) Mocz P., Burkhart B., Hernquist L., McKee C. F., Springel V., 2017, ApJ, 838, 40
  * 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/
  * Orkisz et al. (2017) Orkisz J. H., et al., 2017, A&A, 599, A99
  * Ossenkopf & Mac Low (2002) Ossenkopf V., Mac Low M.-M., 2002, A&A, 390, 307
  * Padoan & Nordlund (2011) Padoan P., Nordlund Å., 2011, ApJ, 730, 40
  * Padoan et al. (1997) Padoan P., Nordlund P., Jones B. J. T., 1997, Commmunications of the Konkoly Observatory Hungary, 100, 341
  * Padoan et al. (2014) Padoan P., Federrath C., Chabrier G., Evans II N. J., Johnstone D., Jørgensen J. K., McKee C. F., Nordlund Å., 2014, Protostars and Planets VI, pp 77–100
  * Palmeirim et al. (2013) Palmeirim P., et al., 2013, A&A, 550, A38
  * Pillai et al. (2015) Pillai T., Kauffmann J., Tan J. C., Goldsmith P. F., Carey S. J., Menten K. M., 2015, The Astrophysical Journal, 799, 74
  * Planck Collaboration et al. (2016a) Planck Collaboration et al., 2016a, A&A, 586, A137
  * Planck Collaboration et al. (2016b) Planck Collaboration et al., 2016b, A&A, 586, A138
  * Price & Bate (2008) Price D. J., Bate M. R., 2008, MNRAS, 385, 1820
  * Price & Federrath (2010) Price D. J., Federrath C., 2010, MNRAS, 406, 1659
  * Roman-Duval et al. (2011) Roman-Duval J., Federrath C., Brunt C., Heyer M., Jackson J., Klessen R. S., 2011, ApJ, 740, 120
  * Schekochihin et al. (2002) Schekochihin A. A., Boldyrev S. A., Kulsrud R. M., 2002, ApJ, 567, 828
  * Schneider et al. (2013) Schneider N., et al., 2013, ApJ, 766, L17
  * Schruba et al. (2019) Schruba A., Kruijssen J. M. D., Leroy A. K., 2019, The Astrophysical Journal, 883, 2
  * Seifried et al. (2020) Seifried D., Walch S., Weis M., Reissl S., Soler J. D., Klessen R. S., Joshi P. R., 2020, arXiv e-prints, p. arXiv:2003.00017
  * Seta et al. (2015) Seta A., Bhat P., Subramanian K., 2015, Journal of Plasma Physics, 81, 395810503
  * Seta et al. (2020) Seta A., Bushby P. J., Shukurov A., Wood T. S., 2020, Physical Review Fluids, 5, 043702
  * Sharda et al. (2020) Sharda P., Federrath C., Krumholz M. R., 2020, arXiv e-prints, p. arXiv:2002.11502
  * Shimajiri et al. (2019) Shimajiri Y., André P., Palmeirim P., Arzoumanian D., Bracco A., Könyves V., Ntormousi E., Ladjelate B., 2019, A&A, 623, A16
  * Soler (2019) Soler J. D., 2019, A&A, 629, A96
  * Soler et al. (2017) Soler J. D., et al., 2017, A&A, 603, A64
  * Solomon et al. (1987) Solomon P. M., Rivolo A. R., Barrett J., Yahil A., 1987, ApJ, 319, 730
  * Tritsis & Tassis (2016) Tritsis A., Tassis K., 2016, MNRAS, 462, 3602
  * Tritsis & Tassis (2018) Tritsis A., Tassis K., 2018, Science, 360, 635
  * Tritsis et al. (2015) Tritsis A., Panopoulou G. V., Mouschovias T. C., Tassis K., Pavlidou V., 2015, MNRAS, 451, 4384
  * Tritsis et al. (2018) Tritsis A., Federrath C., Schneider N., Tassis K., 2018, MNRAS, 481, 5275
  * Verma (2017) Verma M. K., 2017, Reports on Progress in Physics, 80, 087001
  * 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
  * Yuen & Lazarian (2017a) Yuen K. H., Lazarian A., 2017a, arXiv e-prints, p. arXiv:1703.03026
  * Yuen & Lazarian (2017b) Yuen K. H., Lazarian A., 2017b, ApJ, 837, L24
  * Zikanov & Thess (1998) Zikanov O., Thess A., 1998, Journal of Fluid Mechanics, 358, 299–333

##  Appendix A The compressible quasi-static fluctuations model

The general MHD induction equation is

$\frac{\partial\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{\partial
t}=\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}),$
(30)

but with a magnetic field decomposition that we show in Equation (9) it is
possible to simplify the equation, and gain some physical intuition for the
fluctuating component of $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$. First we
write
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)=\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)$,

$\frac{\partial\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{\partial
t}=\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times[\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}]),$
(31)

since $\partial_{t}\mathrm{{\mn@boldsymbol{\mathit{B_{0}}}}}=0$. Expanding the
cross product twice we find,

$\displaystyle\frac{\partial\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{\partial
t}=\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0})+\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}).$
(32)

Next we use the identity
$\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{A}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}})=\mathrm{{\mn@boldsymbol{\mathit{A}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{B}}}})-\mathrm{{\mn@boldsymbol{\mathit{B}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{A}}}})+(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{A}}}}-(\mathrm{{\mn@boldsymbol{\mathit{A}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
to expand each of the two terms. The first term is,

$\displaystyle\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0})=$
$\displaystyle\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0})-\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})$
$\displaystyle+(\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{v}}}}-(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0},$
(33)

which we can simplify because the magnetic field is divergence free, and
because $\partial_{x_{i}}\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}=0$, hence,

$\displaystyle\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0})=-\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})+(\mathrm{{\mn<EMAIL_ADDRESS>(34)

Since
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}=B_{0}\hat{\mathrm{{\mn@boldsymbol{\mathit{z}}}}}$
the second term can be written in terms of just the velocity gradient in the
$z$ direction,

$\displaystyle\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0})=B_{0}\partial_{z}\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}).$
(35)

Now we turn our attention to the fluctuating term in Equation (32). Using the
same identity as we used for the mean-field term and immediately dropping the
$\nabla\cdot\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ term we find,

$\displaystyle\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}})=-\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})+(\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{v}}}}-(\mathrm{{\mn<EMAIL_ADDRESS>(36)

Now combing the two Equation (35) and (36) to construct Equation (32),

$\displaystyle\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0})+$
$\displaystyle\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}})=$
$\displaystyle
B_{0}\partial_{z}\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})+$
$\displaystyle(\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{v}}}}-(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla)\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}-\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}),$
(37)

which can be simplified to reveal,

$\displaystyle\frac{\partial\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{\partial
t}+(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla)\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}=$
$\displaystyle
B_{0}\partial_{z}\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\mathrm{{\mn@boldsymbol{\mathit{B}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})+(\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{v}}}},$
(38)

or more succinctly, using
$\frac{D}{Dt}=\partial_{t}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla$
as the Lagrangian derivative, i.e. the derivative in the frame co-moving with
the fluid,

$\displaystyle\frac{D\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{Dt}=$
$\displaystyle
B_{0}\partial_{z}\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\mathrm{{\mn@boldsymbol{\mathit{B}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})+(\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{v}}}},$
(39)

The first term, $B_{0}\partial_{z}\mathrm{{\mn@boldsymbol{\mathit{v}}}}$, (we
use the notation
$B_{0}\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ in the main
text) tells us that the fluctuations change when the velocity gradient along
the mean magnetic field changes. This term encodes the anisotropy of the
fluctuations into the induction equation. Thus, when $B_{0}$ increases (or
equivalently $\operatorname{\mathcal{M}_{\text{A0}}}$ decreases) the
anisotropy in fluctuations increases and the fluctuations become more
isotropic. This is also demonstrated via results from our simulations in
Figure 6. The second term,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}})$,
tells us that the change in the fluctuations scale with the compression of the
velocity field weighted by the magnetic field. The last term,
$(\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
is the advection of the velocity from the magnetic field fluctuations.

##  Appendix B The compressible quasi-static fluctuations model in the limit
of large $B_{0}$

For $B_{0}\gg\delta B$ Equation (2.1) simplifies significantly. The time
derivative of $\delta B$, on the LHS, is approximately zero, and so are the
terms that have strict $\delta B$ dependence on the RHS. Hence

$\displaystyle 0$
$\displaystyle=B_{0}\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}-\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}),$
(40) $\displaystyle\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}$
$\displaystyle=\frac{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}{B_{0}}(\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}).$
(41)

The $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}/B_{0}$ term is just the unit
vector of $\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$,
$\hat{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, which contributes nothing
to the magnitude, but sets direction of the compressive modes. Hence the
magnitude of the compression is equal to the magnitude of the velocity
gradient along the magnetic field,

$\displaystyle|\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}|=|\partial_{\parallel}\mathrm{{\mn@boldsymbol{\mathit{v}}}}|,$
(42)

as shown in Equation (12), in the main text.

##  Appendix C Weak $\operatorname{\mathcal{M}}$ dependence of mean-field
weighted fluctuations

Figure 10 shows how the perpendicular (left) and parallel (right) magnetic
field fluctuations only very weakly depend upon $\operatorname{\mathcal{M}}$.
We discuss this result in §5.3.

Figure 11 shows the $\delta B_{\parallel}/B_{0}$ (triangles) and $\delta
B_{\perp}/B_{0}$ fluctuations. We show that they follow a similar trend as the
magnitude of the full 3D fluctuations,
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|/B_{0}$, as discussed in §5.3.

Figure 10: The same as Figure 4, but for the mean-field weighted fluctuations
as discussed in §5.3. Figure 11: The same as Figure 7, but for $\delta
B_{\parallel}/B_{0}$ (triangles) and $\delta B_{\perp}/B_{0}$ (circles). We
see both of the magnetic field components follow a similar trend as
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|/B_{0}$, as discussed in §5.3.