AppendixAppendix References

# The density distribution and the physical origins of density intermittency
in sub- to trans-Alfvénic supersonic turbulence

James R. Beattie${}^{\href https://orcid.org/0000-0001-9199-7771\,1}$, Philip
Mocz${}^{\href https://orcid.org/0000-0001-6631-2566\,2\divideontimes,3}$,
Christoph Federrath${}^{\href https://orcid.org/0000-0002-0706-2306\,1,4}$ and
Ralf S. Klessen${}^{\href https://orcid.org/0000-0002-0560-3172\,5,6}$
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
2Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane,
Princeton, NJ 08544, USA
3Lawrence Livermore National Laboratory, 4 Ivy Lane, 7000 East Ave, Livermore,
CA 94550, USA
4Australian Research Council Centre of Excellence in All Sky Astrophysics
(ASTRO3D), Canberra, ACT 2611, Australia
5Universität Heidelberg, Zentrum für Astronomie, Institut für Theoretische
Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany
6Universität Heidelberg, Interdisziplinäres Zentrum für Wissenschaftliches
Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany
⋇Einstein Fellow E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The probability density function (PDF) of the logarithmic density contrast,
$s=\ln(\rho/\rho_{0})$, with gas density $\rho$ and mean density $\rho_{0}$,
for hydrodynamical supersonic turbulence is well-known to have significant
signatures of intermittency that monotonically increase with the turbulent
Mach number, $\operatorname{\mathcal{M}}$. By studying the mass- and volume-
weighted $s$-PDF for an ensemble of 16 sub- to trans-Alfv́enic mean-field,
supersonic, isothermal turbulence simulations, relevant to molecular gas in
the cool interstellar medium, we show that a more intricate picture emerges
for the intermittency of $s$. Using four independent measures of the
intermittency we find hydrodynamical-like intermittency in the highly
magnetised plasma for $\mathcal{M}\lesssim 4$. However, for
$\mathcal{M}\gtrsim 4$, the signatures of intermittency disappear, leaving
approximately lognormal $s$-statistics – exactly the opposite of
hydrodynamical turbulence in the high-$\mathcal{M}$ limit. To understand the
$\mathcal{M}\lesssim 4$ intermittency we use one-dimensional (1D) pencil beams
to explore the dynamics along and across the mean magnetic field,
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. We discuss
kinetic, density and magnetic field fluctuations from the pencil beams, and
identify physical sources of intermittency as single, strong shocks coupled to
fast magnetosonic compressions that form along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ and create large,
volume-poor under-densities. These under-densities contribute significantly to
the skewness of the $s$-PDF. We confirm this result independently using 1D
fluid shock simulations. We discuss the Gaussianisation of the
$\mathcal{M}\gtrsim 4$ $s$-fields through the lens of two phenomenologies: the
self-similarity of the $s$-field and homogenisation of the dynamical
timescales between the over- and under-dense regions in the compressible gas.

###### keywords:

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

††pubyear: 2021††pagerange: The density distribution and the physical origins
of density intermittency in sub- to trans-Alfvénic supersonic turbulence–A

##  1 Introduction

The PDF of the gas density is a valuable tool for understanding the nature of
compressible turbulence and star formation in the interstellar medium of
galaxies (Vazquez-Semadeni, 1994; Padoan et al., 1997; Passot & Vázquez-
Semadeni, 1998; Federrath et al., 2008b; Federrath et al., 2009; Brunt et al.,
2010a, b; Krumholz & McKee, 2005; Hennebelle et al., 2011; Federrath &
Klessen, 2012; Burkhart & Lazarian, 2012; Konstandin et al., 2012b; Molina et
al., 2012; Hopkins, 2013; Nolan et al., 2015; Federrath & Banerjee, 2015;
Squire & Hopkins, 2017; Pan et al., 2019; Mocz & Burkhart, 2019; Menon et al.,
2021; Menon et al., 2020; Khullar et al., 2021). Simple models of molecular
clouds (MCs) in the interstellar medium (ISM), which are supersonic and
magnetised, and have not yet started to collapse under their own self-gravity
have an approximately normal volume-weighted $s$-PDF111Or equivalently, a
lognormal $\rho$-PDF., where $s\equiv\ln(\rho/\rho_{0})$, and $\rho$ is the
cloud density, with $\rho_{0}$ the value of the volume-weighted mean. It
follows

$\displaystyle p_{\text{N},V}(s;\sigma_{s,V}^{2})$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma_{s,V}^{2}}}\exp\left\\{-\frac{(s-s_{0,V})^{2}}{2\sigma_{s,V}^{2}}\right\\},$
(1) $\displaystyle s_{0,V}$ $\displaystyle=-\frac{\sigma_{s,V}^{2}}{2},$ (2)
$\displaystyle\sigma^{2}_{s,V}$
$\displaystyle=f(\operatorname{\mathcal{M}},\operatorname{\mathcal{M}_{\text{A}}},b,\gamma,\Gamma).$
(3)

The log-density variance, $\sigma_{s}^{2}$, captures the density variations
induced by different physical processes in a MC. It is a function of (i) the
turbulent Mach number,

$\displaystyle\operatorname{\mathcal{M}}=\sigma_{V}/c_{s},$ (4)

where $\sigma_{V}$ is the velocity dispersion on system scale $L$, and $c_{s}$
is the sound speed (Vazquez-Semadeni, 1994; Padoan et al., 1997; Passot &
Vázquez-Semadeni, 1998; Price et al., 2011; Konstandin et al., 2012b), (ii)
the Alfvén Mach number,

$\displaystyle\operatorname{\mathcal{M}_{\text{A}}}=\sigma_{V}/V_{\rm
A}=c_{s}\operatorname{\mathcal{M}}/V_{\rm A},$ (5)

where $V_{\text{A}}=B/\sqrt{4\pi\rho}$ is the root-mean-squared Alfvén wave
velocity and $B$ is the magnetic field (Padoan & Nordlund, 2011; Molina et
al., 2012; Beattie et al., 2021), (iii) the turbulent driving parameter, $b$,
which captures the influence of compressive
($\nabla\times\mathrm{{\mn@boldsymbol{\mathit{f}}}}=0$) or solenoidal
($\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{f}}}}=0$) forcing modes on the
density fluctuations, where $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ is a
source of turbulent forcing (Federrath et al., 2008b, 2010; Menon et al.,
2021), (iv) the adiabatic index $\gamma$ (Nolan et al., 2015), and (v) the
polytropic index $\Gamma$ (Federrath & Banerjee, 2015). Likewise, the mass-
weighted distribution is,

$\displaystyle p_{\text{N},M}(s;\sigma_{s,M}^{2})$
$\displaystyle=\frac{\exp\left\\{s\right\\}}{\sqrt{2\pi\sigma_{s,M}^{2}}}\exp\left\\{-\frac{(s-s_{0,M})^{2}}{2\sigma_{s,M}^{2}}\right\\},$
(6) $\displaystyle s_{0,M}$ $\displaystyle=\frac{\sigma_{s,V}^{2}}{2},$ (7)
$\displaystyle\sigma^{2}_{s,M}$ $\displaystyle=\sigma^{2}_{s,V},$ (8)

which shows how, in a lognormal density-fluctuation theory, the mass-weighted
and volume-weighted distributions are intrinsically linked through the mean
(Equation 7) and variance (Equation 8) of the volume-weighted $s$-PDF (Li et
al., 2003). In hydrodynamical, supersonic turbulence it is well-known that
there is a significant difference between $\sigma^{2}_{s,M}$ and
$\sigma^{2}_{s,V}$ (an excess of $\sigma^{2}_{s,V}$), which in essence is from
the emergence of non-Gaussian, intermittent events in the fluid (Passot &
Vázquez-Semadeni, 1998; Kritsuk et al., 2007; Federrath et al., 2010;
Federrath, 2013; Hopkins, 2013; Squire & Hopkins, 2017; Mocz & Burkhart,
2019).

The $s$-statistics and intermittency of strongly magnetised turbulence, where
the flow is globally anisotropic about a coherent, strong mean-field has not
been studied in great detail, even though it is indeed a relevant flow regime
for either some regions inside of MCs (Li et al., 2013; Federrath et al.,
2016; Hu et al., 2019; Heyer et al., 2020; Skalidis & Tassis, 2020; Hwang et
al., 2021; Hoang et al., 2021) or even perhaps most MCs (Li, 2021). In this
study we probe the physics of the $s$-PDFs in simulations with
$2\lesssim\operatorname{\mathcal{M}}\lesssim 20$ and
$0.1\lesssim\operatorname{\mathcal{M}_{\text{A0}}}\lesssim 2$, where
$\operatorname{\mathcal{M}_{\text{A0}}}=c_{s}\operatorname{\mathcal{M}}\sqrt{4\pi\rho_{0}}/|\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}|$,
revealing the nature of voids, rarefactions and over-dense structures in the
trans- to sub-Alfvénic mean-field turbulence, i.e., turbulence where the
energy in the mean magnetic field is larger than that of the kinetic turbulent
energy. We show that the $s$-intermittency for compressible, hydrodynamical
turbulence is different from turbulence in the highly-magnetised mean-field
regime. Specifically, we find that Equation 8 approximately holds for moderate
to high-$\operatorname{\mathcal{M}}$ magnetised flows, which is not the case
for purely hydrodynamical turbulence due to strong intermittency. We discuss a
range of phenomenological models for why this may be the case and, for the
first time, identify physical sources of $s$-intermittency in real-space,
confirming that an important source of intermittency in these flows are strong
shocks that produce intermittent, deep, volume-poor rarefactions in the gas
density along the coherent magnetic field.

This study is organised as follows. First, in §2 we revisit the lognormal
model for the density fluctuations. In §3 we extend the discussion to
$s$-fluctuation models that include intermittency effects. In §4 we outline
the supersonic MHD turbulence simulations that we use to explore the
logarithmic density fluctuations. In §5 we analyse the morphology of the full
three-dimensional (3D) volume-weighted and mass-weighted $s$-PDFs from the
simulation data, including fitting lognormal and non-lognormal models for the
volume-weighted PDFs. Next, in §6 we compute four independent measures of
$s$-intermittency, and show how the highly-magnetised regime has a more
complicated intermittency structure than the hydrodynamical regime. In §7 we
use one-dimensional (1D) probes of the low-$\operatorname{\mathcal{M}}$
turbulence to reveal the parallel and perpendicular to mean magnetic field
$s$, velocity, and magnetic field dynamics, and use this 1D analysis to
identify the rarest (intermittent) events in the turbulence, in real-space. In
§8 we discuss the lack of intermittency in high-$\operatorname{\mathcal{M}}$,
magnetised turbulence. Finally in §9 we summarise the key results of this
study.

##  2 The lognormal density-PDF

###  2.1 The lognormal model

Lognormal models for the PDF of turbulent $\delta=\rho/\rho_{0}$ fluctuations
originate from Vazquez-Semadeni (1994). They consider $\delta$ in the self-
similar, hierarchical structure of the cool, isothermal ISM, where global
pressure is negligible $(\operatorname{\mathcal{M}}\gg 1)$ and self-gravity is
yet to dominate the dynamics of the region. The lognormal PDF is motivated by
assuming that for time, $t_{n}$, the density can be expressed as a
multiplicative interaction through independent jumps $\delta$ in the Eulerian
frame of the gas,

$\displaystyle\frac{\rho(t_{n})}{\rho_{0}}=\delta_{n}\delta_{n-1}\ldots\delta_{1}\delta_{0}\frac{\rho(t_{0})}{\rho_{0}}=\left(\prod^{n}_{i=0}\delta_{i}\right)\frac{\rho(t_{0})}{\rho_{0}},$
(9)

where $\rho(t_{0})/\rho_{0}=1$ is the initial density in units of the mean.
This means that under the log-transformation, fluctuations become additive,

$\displaystyle
s(t_{n})=\ln\frac{\rho(t_{n})}{\rho_{0}}=\sum^{n}_{i=0}\ln\delta_{i}=\sum^{n}_{i=0}s_{i},$
(10)

If the density fluctuations are

  1. 1.

generated by the same underlying probability distribution, and

  2. 2.

are independent from one another (i.e., not temporally correlated),

then the central limit theorem states that the distribution of the logarithmic
density fluctuations should tend towards a Gaussian distribution as the number
of fluctuations increase. This framework lets us understand the nature of a
single (Eulerian) density fluctuation changing in time. However, Vazquez-
Semadeni (1994) further argued that since the hydrodynamical equations are
self-similar in space (i.e., invariant to arbitrary length scaling) the
fluctuations should be lognormal on all scales. With these assumptions, the
lognormal distribution should aptly describe the density fluctuations on any
length scale of the turbulence.

###  2.2 Key issues with the lognormal model

####  2.2.1 Mass conservation

Hopkins (2013) articulates that having a lognormal distribution on all scales
in the flow violates mass conservation. We show this by considering the
$\rho/\rho_{0}$-PDF, $p_{L}(\delta)$, on the system scale, $L$. This
distribution is the convolution of all PDFs from scales below the system
scale, $L/\Gamma$, $L/\Gamma^{2}$, $\ldots$, $L/\Gamma^{n}$, $\ldots$, where
$\Gamma\gtrsim 1$, in the turbulence, which is simply an application of the
Law of Total Probability222This is simply because
$(\rho/\rho_{0})_{L/\Gamma^{n}}$ on each scale is a random variable composed
of the (volume-weighted) linear some of densities on smaller scales,
$(\rho/\rho_{0})_{L/\Gamma^{n}}\propto(\rho/\rho_{0})_{L/\Gamma^{n+1}}+(\rho/\rho_{0})_{L/\Gamma^{n+2}}+\ldots$
This, by definition, leads to the PDF on each scale being an infinite
convolution of PDFs from scales below it. (Castaing, 1996). Writing this
formerly,

$\displaystyle p_{L}(\delta)$
$\displaystyle=p_{L/\Gamma}(\delta)\otimes\ldots\otimes
p_{L/\Gamma^{n}}(\delta)\otimes\ldots=\bigotimes_{n=1}^{\infty}p_{L/\Gamma^{n}}(\delta),$
(11)

where $\otimes$ is the convolution operator, $f(x)\otimes
g(x)=\int^{\infty}_{-\infty}\operatorname{d}\\!{t}\,f(x-t)g(t)$. Now we assume
that all of the PDFs, $p_{L/\Gamma}(\delta)\ldots
p_{L/\Gamma^{n}}(\delta)\ldots$, on the RHS of Equation 11 are lognormal
distributions and inquire about the LHS of the equation. This is a well-known
and historical problem in statistics, motivated by problems from a broad range
of disciplines, from telecommunications to the biosciences (see for example,
Wu et al., 2005; Lo, 2012; Ben Hcine & Bouallegue, 2015). The key conclusion
one makes is that the LHS of Equation 11 cannot possibly be a lognormal
distribution. In fact, because no scale is special in the turbulence, and
instead of picking $p_{L}(\delta)$, we pick any arbitrary
$p_{L/\Gamma^{n}}(\delta)$, then it is clear that only one of the infinite
$p_{L/\Gamma^{n}}(\delta)$ can possibly be lognormal. This is directly related
to mass conservation. It follows by assuming that all
$p_{L/\Gamma^{n}}(\delta)$ are lognormal, then the total mass is,

$\displaystyle M_{1}=\int_{0}^{\infty}\operatorname{d}\\!{\delta}\,\delta
p_{L}(\delta),$ (12)

which ought to equal

$\displaystyle
M_{2}=\int_{0}^{\infty}\operatorname{d}\\!{\delta}\,\delta\bigotimes_{n=1}^{\infty}p_{L/\Gamma^{n}}(\delta).$
(13)

However, because $\bigotimes_{n=1}^{\infty}p_{L/\Gamma^{n}}(\delta)$ does not
converge to a lognormal distribution, and because we also have assumed that
_all_ PDFs are, including $p_{L}(\delta)$, then $M_{1}\neq M_{2}$, and hence a
lognormal model on all scales violates mass conservation.

####  2.2.2 The sonic scale

This is not the only reason why the lognormal distribution, which describes a
scale-free random process, cannot be strictly “the right model”333Note that it
can, however, be very close to a lognormal model, and for practical purposes
this might suffice, but theoretically the $s$-PDF can not be _exactly_ a
normal model. for the density fluctuations in supersonic turbulence. Recent
numerical experiments by Federrath et al. (2021) reveal that supersonic
turbulence is not scale-free, and has a characteristic length, the sonic
scale, $\ell_{s}$, for which

$\displaystyle\left(\left\langle|\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{r}}}})-\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{r}}}}+\ell_{s})|^{2}\right\rangle_{\mathrm{{\mn@boldsymbol{\mathit{r}}}}}/2\right)^{1/2}=c_{s},$
(14)

where
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}(\mathrm{{\mn@boldsymbol{\mathit{r}}}})$
is the velocity field of the turbulence. Furthermore, $c_{s}$ is the sound
speed and
$\left\langle\ldots\right\rangle_{\mathrm{{\mn@boldsymbol{\mathit{r}}}}}$ is
the ensemble average over positions, $\mathrm{{\mn@boldsymbol{\mathit{r}}}}$
where the turbulence above $\ell_{s}$ behaves like Burgers (1948) turbulence
(the turbulence of interacting sawtooth waves), and the scales below
Kolmogorov (1941) turbulence (the turbulence of inertially interacting
eddies). This means that one can not model the small-scale turbulence and
scale the results to the system-scale, as in Vazquez-Semadeni (1994). In fact,
through direct measurements of the $s$-PDF below and above $\ell_{s}$
(Federrath et al., 2021) show that on scales below $\ell_{s}$ the $s$-PDF
becomes peaked and kurtotic, and on scales above, negatively skewed. It is
clear that the $s$-PDF is non-Gaussian, even on the subsonic scales of the
turbulence.

##  3 Beyond lognormality

In both ISM observations and in high-resolution numerical studies of
supersonic turbulence, we find that the $s$-PDF deviates from perfect
lognormality (Kritsuk et al., 2007; Federrath et al., 2008b, 2010; Price et
al., 2011; Konstandin et al., 2012a; Federrath, 2013; Hopkins, 2013; Pan et
al., 2019; Menon et al., 2021). The reason for these deviations (ignoring
gravity and power-law tails due to gravitationally bound, collapsing, dense
structures; see Klessen 2000; Federrath et al. 2008a; Kritsuk et al. 2011;
Federrath 2013; Girichidis et al. 2014; Mocz et al. 2017; Burkhart 2018;
Khullar et al. 2021, or due to non-isothermal effects; see Passot & Vázquez-
Semadeni 1998; Nolan et al. 2015; Federrath & Banerjee 2015) are due to
intermittency, in reference to intermittent, i.e., rare events that are non-
Gaussian in nature. More quantitatively, Federrath et al. (2010) described
three different realisations that intermittency manifests itself in turbulent
flows:

  1. 1.

non-Gaussian wings / non-Gaussian higher-order moments of PDFs for turbulent
quantities, e.g., $\rho$, $\mathrm{{\mn@boldsymbol{\mathit{v}}}}$,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, their derivatives, and combinations
of the quantities (Burkhart et al., 2009; Hopkins, 2013; Mocz & Burkhart,
2019; Seta et al., 2020; Mohapatra et al., 2020; Beattie et al., 2020);

  2. 2.

anomalous scaling of the higher-order structure functions of the velocity
field (She & Leveque, 1994; Castaing, 1996; Kowal et al., 2007; Schmidt et
al., 2008; Konstandin et al., 2012a; Hopkins, 2013);

  3. 3.

structures with intense vorticity,
$\nabla\times\mathrm{{\mn@boldsymbol{\mathit{v}}}}$, and energy dissipation,
such as strong shocks and filamentary structures,
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}<0$, and rarefied regions
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}>0$ (Kritsuk et al., 2007;
Federrath, 2013; Squire & Hopkins, 2017; Park & Ryu, 2019a; Mocz & Burkhart,
2019; Beattie & Federrath, 2020; Yoffe & McComb, 2021).

In supersonic turbulence, let alone magnetised supersonic turbulence, the
intermittent structures described above dominate the flow, for example through
populations of interacting shocks. This has led to large efforts to understand
intermittency, which is still one of the most fundamental problems in
turbulence, and to characterise it in the ISM (Falgarone et al., 1995, 2009;
Hily-Blant & Falgarone, 2009; Falgarone et al., 2015). It means to properly
model the density field one must include intermittency in the recipe. We now
discuss two models that are based upon the life and happenings of shocked
regions in the turbulence, and give some insight into the processes
underpinning intermittency in the turbulent density field.

Figure 1: Left: The H2013 $s$-PDF model shown for different intermittency
parameters: $0\,(p_{\text{N},V})\leq T\leq 0.25$. This parameter range spans
from no intermittency (lognormal, shown with black, dashed line style, $T=0$)
to strong intermittency ($T=0.25$). We plot in the variable $s/\sigma_{s}$ to
draw attention to just the non-Gaussian morphology of the distribution. Right:
The same as left plot but for the MB2019 $s$-PDF model, plotted for different
$f$ values, corresponding to changing the dynamical timescale of over-dense
and under-dense regions in the turbulence. $f=0\iff p_{\text{N},V}$ is the
lognormal model, and $f=10.0$ is the extreme case where over-dense regions
exist on 16$\times$ shorter timescales than under-dense regions.

###  3.1 Continuous log-Poisson model

Hopkins (2013) (hereafter called H2013) considers the Castaing (1996) model to
describe a general steady-state result from a class of multiplicative random
relaxation processes that produce random multiplicative changes to the
density. The essence of the H2013 model is that an exponential PDF describes
the continuous jump in $s$ between two nested scales in the turbulence, the
$(s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}})$-PDF. Then, using an infinite
convolution series, as in Equation 11, of
$(s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}})$-PDFs from each $L/\Gamma^{n}$, the
$s$-PDF is derived to be:

$\displaystyle p_{\rm{H}}(s;\sigma_{s,V}^{2},T)$
$\displaystyle=I_{1}(2\sqrt{\lambda
u})\sqrt{\frac{\lambda}{T^{2}u}}\exp\left\\{-(\lambda+u)\right\\},$ (15)
$\displaystyle u$ $\displaystyle\equiv\frac{\lambda}{1+T}-\frac{s}{T},\,u\geq
0,$ (16) $\displaystyle s_{0,V}$
$\displaystyle=-\frac{\sigma_{s,V}^{2}}{2}(1+T)^{-1},$ (17)
$\displaystyle\sigma_{s,V}^{2}$ $\displaystyle=2\lambda T^{2},$ (18)

where $I_{1}(x)$ is the modified Bessel function of the first kind, $\lambda$
is associated with the probability of having a density jump (in a Poisson
manner) between two arbitrary, nested scales in the turbulence,
$L/\Gamma^{n}$, $L/\Gamma^{n+1}$, and $T$ is the magnitude of the mean
logarithmic density jump between these scales,

$\displaystyle T=\mathrm{\left\langle
s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}}\right\rangle}_{L/\Gamma^{n}}=\mathrm{\left\langle\ln\frac{(\rho/\rho_{0})_{L/\Gamma^{n}}}{(\rho/\rho_{0})_{L/\Gamma^{n+1}}}\right\rangle}_{L/\Gamma^{n}},$
(19)

where $\mathrm{\left\langle\ldots\right\rangle}$ is the ensemble average.
Hence $T$ can be thought of as an mean $s$-fluctuation between nested scales,
averaged over all scales in the turbulence (Squire & Hopkins, 2017). Large
values (large in this context is not much greater than 0, noting that $T$ has
scales similar to that of $s$, as shown in left panel of Figure 1) of $T$
implies that the $\mathrm{\left\langle
s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}}\right\rangle}_{L/\Gamma^{n}}$ is
significant, due to, for example, low-volume filling, high-mass over-dense
regions. We show illustrations of the H2013 PDF in the left panel of Figure 1
for different $T$ values, showing how larger values of $T$ increase the
negative skewness of the PDF. For small values of $T$, there are only
negligible density jumps (on average) between any two nested scales. This
means as $T\rightarrow 0$, the $s$-field becomes perfectly self-similar,
$\mathrm{\left\langle
s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}}\right\rangle}_{L/\Gamma^{n}}=0$, and can
be approximated by Gaussian-like statistics. Hence $T$ can also be thought of
as the average deviation away from perfect self-similarity in the $s$-field,
across all scales. Yet another way of interpreting $T$ is by rewriting
Equation 18 in terms of the mass-weighted variance,

$\displaystyle\frac{\sigma_{s,V}^{2}}{\sigma_{s,M}^{2}}$
$\displaystyle=\left(1+T\right)^{3},$ (20) $\displaystyle T$
$\displaystyle=\left(\frac{\sigma_{s,V}}{\sigma_{s,M}}\right)^{2/3}-1.$ (21)

As we mentioned in §1, for a lognormal density-fluctuation theory
$\sigma_{s,V}^{2}=\sigma_{s,M}^{2}$, corresponding to $T=0$ in the above
expression. Therefore non-zero values of $T$ lead a modified relationship
between the $\sigma_{s,V}^{2}$ and $\sigma_{s,M}^{2}$, which defines a non-
lognormal map between the mass-weighted and volume weighted distributions. We
will discuss this in more detail in §3.

The H2013 and Castaing (1996) models were motivated for isotropic, homogeneous
turbulence, and hence may not be phenomenologically appropriate for MHD. Still
we find that the PDF, at least empirically, captures the morphology of the $s$
fluctuations well over a wide range of $\operatorname{\mathcal{M}}$ and for
$\operatorname{\mathcal{M}_{\text{A0}}}\lesssim 2$. In this study we explore
how well this PDF can be extrapolated to the highly-magnetised, highly-
supersonic regime and using the $T$ statistic, we explore how the
intermittency of $s$ behaves in the presence of strong magnetic fields and
high-$\operatorname{\mathcal{M}}$, which was a parameter range not explored in
the original H2013 analysis (large $s$ variance, strong
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$-field).

###  3.2 Langevin model

Recent work by Mocz & Burkhart (2019), herein MB2019, model the density-PDF of
hydrodynamical density fluctuations using a Markov process framework. MB2019
construct a Langevin model,

$\displaystyle s(t+\operatorname{d}\\!{t})$
$\displaystyle=s(t)+A(s)\operatorname{d}\\!{t}+\mathcal{N}(0,1)\sqrt{D(s)\operatorname{d}\\!{t}},$
(22) $\displaystyle A(s)$
$\displaystyle=-\frac{s-s_{0,V}}{\tau_{A}}\left[1+H(s-s_{0,V})\frac{3f}{2}\right],$
(23) $\displaystyle D(s)$ $\displaystyle=\frac{2\sigma_{s,V}^{2}}{\tau_{D}},$
(24)

where $A(s)$ is the deterministic or advective term in the model, $D(s)$ is
the stochastic or diffusive term, and $\mathcal{N}(0,1)$ is a standard normal
distribution. The stochastic term, $D(s)$, contains the turbulent fluctuations
$\sigma_{s}^{2}$, which change on a dynamical timescale

$\displaystyle\tau_{D}=\ell_{0}/(c_{s}\operatorname{\mathcal{M}}),$ (25)

where $\ell_{0}$ is the turbulent driving scale. The deterministic term,
$A(s)$, encodes how logarithmic density dynamically fluctuates about $s_{0}$
(as a mean-reverting random walk) on timescales

$\displaystyle\tau_{A}=\frac{\tau_{A,0}}{[1+H(s-s_{0,V})\frac{3f}{2}]},$ (26)

where $H(s-s_{0,V})$ is the Heaviside function, $\tau_{A,0}$ is the
characteristic timescale for the fluctuations and $f$ is a constant that
encodes how high-density structures, such as the density contrast caused by a
shock, live on shorter timescales than the rest of the density fluctuations in
the fluid (Robertson & Goldreich, 2018; Scannapieco & Safarzadeh, 2018). For
example, for $s>s_{0,V}$, $\tau_{A}$ is reduced by $1+3f/2$, and hence $f$
becomes the fitting parameter for how much shorter the dynamical timescales
are for the shocked density structures compared to rarefied and mean-density
regions in the turbulence. MB2019 found $f\approx 0.15$ using 1D simulations
over a large range of $\operatorname{\mathcal{M}}$, which means shocked
regions operate on dynamical timescales $\approx 20\%$ shorter than the low-
density regions. In this study, we test if $f$ is constant or rather varies
over a large range of $\operatorname{\mathcal{M}}$ and
$\operatorname{\mathcal{M}_{\text{A0}}}$, extending the analysis of MB2019.

Figure 2: A 2D slice through the 3D logarithmic density field in the $x=L/2$
plane parallel to the mean magnetic field (the direction of
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ is indicated in the
top-left panel). Red indicates the over-densities $(\rho/\rho_{0}>1)$ and blue
the under-densities $(\rho/\rho_{0}<1)$; the density range is indicated in
each panel. The plots are organised such that the weakest turbulence (smallest
root-mean-squared velocities) is on the left,
$\operatorname{\mathcal{M}}\approx 2$, and strongest on the right,
$\operatorname{\mathcal{M}}\approx 20$, and the strongest mean magnetic field
is on the top, $\operatorname{\mathcal{M}_{\text{A0}}}\approx 0.1$ and weakest
on the bottom, $\operatorname{\mathcal{M}_{\text{A0}}}\approx 2.0$. Overlaid
on each plot is an iterative line integral convolution through the
corresponding 2D slice of the magnetic field, revealing the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field topology. Qualitatively, the top
panels show density structures that are highly-anisotropic and stretched along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ whilst in the
bottom panels the structures become more isotropic as the total
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ becomes dominated by the turbulent
component of the field which is, among other things, advected with the
supersonic, turbulent velocities.

The PDF of MB2019’s Langevin model, $p_{\rm{MB}}$, defines a solution to the
steady-state (time-independent) Fokker-Planck equation,

$\displaystyle\frac{\partial}{\partial
s}\left[A(s)p_{\rm{MB}}(s)\right]=\frac{1}{2}\frac{\partial^{2}}{\partial
s^{2}}\left[D(s)p_{\rm{MB}}(s)\right],$ (27)

which has a solution of the form,

$\displaystyle p_{\rm{MB}}(s)$
$\displaystyle\propto\exp\left\\{-\frac{(s-s_{0,V})^{2}[1+f(s-s_{0,V})H(s-s_{0,V})]}{2\sigma_{s,V}^{2}(\tau_{A}/\tau_{D})}\right\\},$
(28)

which we show for a number of different $f$ parameters in the right panel of
Figure 1, where $f=0$ corresponds to no difference between the high- and low-
density fluctuation timescales, i.e., a lognormal PDF, up to $f=10$, which
corresponds to high-density structures that operate on timescales 10 times
shorter than the low-density regions. It is obvious that the value of $f$
encodes a $3^{\rm rd}$ moment, skewness (discussed more in §6), into the PDF
through the $\mathcal{O}([s-s_{0}]^{3})$ term in the exponential. This
provides a similar morphology to the H2013 model but, with a softer high-
density tail (exponential rather than Bessel function).

We will use the MB2019 and H2013 non-lognormal models as building blocks to
interpret results in §5. To summarise, MB2019’s Langevin model describes
skewness in terms of the difference in dynamical timescales between over-dense
and under-dense regions in the turbulence. In contrast, the H2013 model
describes a continuous log-Poisson model process, which quantifies the
intermittency with the parameter $T$, which is phenomenologically associated
with the deviation from perfect self-similarity of the $s$-field, and can
alternatively be seen as a modification to the mass-weighted and volume-
weighted variance relation. Finally, we note that both of the models we
discuss were constructed to address the hydrodynamical phenomenology of
intermittency. We believe both models, which in essence rely upon the
statistics and timescales of shocked density regions, are able to provide some
insight into general supersonic turbulent flows, because large populations of
shocked regions form in the density field, with or without a magnetic field
present (Lehmann et al., 2016; Park & Ryu, 2019b; Beattie & Federrath, 2020;
Beattie et al., 2020; Beattie et al., 2021). Before discussing sub-Alfvénic
$s$-PDF data, we turn to the details of the numerical experiments that we use
in our study.

##  4 Supersonic turbulence Simulations

Figure 3: An example of the M2MA01 simulation (Table 1) shown at $t/\tau=6$
with a section of the magnetic field (blue) and velocity streamlines (magma)
surrounded by a box of $\ln(\rho/\rho_{0})$ slices in each direction, and with
volume-rendered over-dense regions $(\rho/\rho_{0}\gg 1)$, shown in black. The
velocity structure reveals a deformed vortex twisting around the magnetic
field. Typical of sub-Alfvénic mean-field turbulence, the magnetic field is
extremely ordered, showing no turbulent component at all via the streamline
visualisation. This is because all of the energy is contained in the
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ (direction
indicated at the bottom-right of the plot) and in the sub-Alfvénic mean-field
regime the turbulent velocity fluctuations are too weak to significantly bend
the field. However, because $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ is not
fixed at the boundary, the turbulent motions can advect the magnetic fields in
the
$\mathrm{{\mn@boldsymbol{\mathit{r}}}}\perp\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$
plane, establishing supersonic vortical motions with dynamical timescales
$\tau=\ell_{0}/(c_{s}\operatorname{\mathcal{M}})$, where $\ell_{0}$ is the
energy injection or driving scale. Figure 4: The volume-weighted (top) and
mass-weighted (bottom) variances of the logarithmic density field for each of
the simulations in Table 1, as a function of eddy turnover time, $t/\tau$,
coloured by $\operatorname{\mathcal{M}}$ and with different line styles for
$\operatorname{\mathcal{M}_{\text{A0}}}$. We find that almost all simulations
have become approximately stationary between $5\leq t/\tau\leq 10$, which is
the interval where we average over to develop all of the statistics in our
study. $\operatorname{\mathcal{M}}$ (changing colour) plays the largest role
in setting the value of the variance, and the magnetic field (changing line
style) only plays a significant role for the high-$\operatorname{\mathcal{M}}$
simulations, because it preferentially suppresses small-scale fluctuations
that are introduced in high-$\operatorname{\mathcal{M}}$ flows (Beattie et
al., 2021), coincident with when the turbulence is becoming more globally
isotropic. Table 1: Main simulation parameters.

Simulation ID | $\operatorname{\mathcal{M}}\,(\pm 1\sigma)$ | $\operatorname{\mathcal{M}_{\text{A0}}}\,(\pm 1\sigma)$ | $N^{3}$
---|---|---|---
M2Ma01 | 2.6 |  $\pm$ 0.2 | 0.131 |  $\pm$ 0.008 | $512^{3}$
M4Ma01 | 5.2 |  $\pm$ 0.4 | 0.13 |  $\pm$ 0.01 | $512^{3}$
M10Ma01 | 12 |  $\pm$ 1 | 0.125 |  $\pm$ 0.006 | $512^{3}$
M20Ma01 | 24 |  $\pm$ 1 | 0.119 |  $\pm$ 0.003 | $512^{3}$
M2Ma05 | 2.2 |  $\pm$ 0.2 | 0.54 |  $\pm$ 0.04 | $512^{3}$
M4Ma05 | 4.4 |  $\pm$ 0.2 | 0.54 |  $\pm$ 0.03 | $512^{3}$
M10Ma05 | 10.5 |  $\pm$ 0.5 | 0.52 |  $\pm$ 0.02 | $512^{3}$
M20Ma05 | 21 |  $\pm$ 1 | 0.53 |  $\pm$ 0.02 | $512^{3}$
M2Ma1 | 2.0 |  $\pm$ 0.1 | 0.98 |  $\pm$ 0.07 | $512^{3}$
M4Ma1 | 3.8 |  $\pm$ 0.3 | 0.95 |  $\pm$ 0.08 | $512^{3}$
M10Ma1 | 9.3 |  $\pm$ 0.5 | 0.93 |  $\pm$ 0.05 | $512^{3}$
M20Ma1 | 18.8 |  $\pm$ 0.7 | 0.93 |  $\pm$ 0.03 | $512^{3}$
M2Ma2 | 1.7 |  $\pm$ 0.1 | 1.7 |  $\pm$ 0.1 | $512^{3}$
M4Ma2 | 3.5 |  $\pm$ 0.1 | 1.73 |  $\pm$ 0.07 | $512^{3}$
M10Ma2 | 9.0 |  $\pm$ 0.4 | 1.8 |  $\pm$ 0.1 | $512^{3}$
M20Ma2 | 18 |  $\pm$ 1 | 1.8 |  $\pm$ 0.1 | $512^{3}$
M2Ma$\infty$ | 2.15 |  $\pm$ 0.07 |  | $\infty$ | $576^{3}$
M4Ma$\infty$ | 4.1 |  $\pm$ 0.1 |  | $\infty$ | $576^{3}$
M10Ma$\infty$ | 10.1 |  $\pm$ 0.5 |  | $\infty$ | $576^{3}$
M20Ma$\infty$ | 20.2 |  $\pm$ 0.7 |  | $\infty$ | $576^{3}$

  * •

Notes: For each simulation we extract 51 realisations at 0.1$\,\tau$
intervals, where $\tau=L/(2c_{s}\operatorname{\mathcal{M}})$ is the turbulent
turnover time, between $5$-$10\,\tau$. All $1\sigma$ fluctuations listed in
the table and study are from time-averaging quantities over the $5\,\tau$.
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{A0}}}=(2c_{s}\operatorname{\mathcal{M}}\sqrt{\pi\rho_{0}})/|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|$,
where $\rho_{0}$ is the mean density and $c_{s}$ is the sound speed. Column
(4): the number of grid cells in the discretisation of the spatial domain,
$L^{3}$.

###  4.1 Ideal (M)HD model and turbulent driving

In this study we analyse the $s$-PDFs using high-resolution, 3D turbulent,
ideal magnetohydrodynamical (MHD) simulations. The corresponding equations are

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}})$ $\displaystyle=0,$
(29) $\displaystyle\rho\left(\frac{\partial}{\partial
t}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla\right)\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}}}}\cdot\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{8\pi}\right)+\rho\mathrm{{\mn@boldsymbol{\mathit{f}}}},$
(30)
$\displaystyle\frac{\partial\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{\partial
t}$
$\displaystyle=\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}),\;\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{B}}}}=0,$
(31)

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}}}}$ the
turbulent acceleration field. The latter follows an Ornstein-Uhlenbeck process
with finite auto-correlation time and is constructed such that we control the
mixture of solenoidal and purely compressive modes in
$\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ (see Federrath et al. 2010 for a
detailed discussion of the turbulence driving). We choose to drive equal
energy in both modes, i.e., a natural mixture of solenoidal and compressive
driving modes. The driving is isotropic in wavenumber space, centred on
$|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|=2$, corresponding to real-space
scales of $\ell_{0}=L/2$, and falling off to zero with a parabolic spectrum
within $1\leq|\mathrm{{\mn@boldsymbol{\mathit{k}}}}L/2\pi|\leq 3$. Thus,
energy is injected only on large scales and the turbulence on smaller scales
develops self-consistently through the turbulent energy cascade. The auto-
correlation timescale of $\mathrm{{\mn@boldsymbol{\mathit{f}}}}$ is equal to
$\tau=\ell_{0}/\sigma_{V}=L/(2c_{s}\operatorname{\mathcal{M}})$, hence, we use
this timescale to set the desired turbulent Mach number on $\ell_{0}=L/2$ for
each simulation. We vary the sonic Mach number between
$\operatorname{\mathcal{M}}=2$ and 20, encapsulating the range of observed
$\operatorname{\mathcal{M}}$ values for turbulent MCs (e.g., Schneider et al.
2013; Federrath et al. 2016; Orkisz et al. 2017; Beattie et al. 2019).

We use a modified version of flash based on version 4.0.1 (Fryxell et al.,
2000; Dubey et al., 2008) to solve the MHD equations in a periodic box with
dimensions $L^{3}$, on a uniform grid with resolution $512^{3}-576^{3}$, using
the multi-wave, approximate Riemann solver framework described in Bouchut et
al. (2010) and implemented in flash in Waagan et al. (2011). We utilise up to
approximately 6000 compute cores in parallel, for roughly one million compute
hours for the most supersonic and magnetised simulations, running at almost
$99\%$ efficiency on the Gadi supercomputer hosted by the National Computing
Infrastructure Australia. For details on the performance of the modified flash
code we refer the readers to Federrath et al. (2021) and for more details
about the current simulations Beattie & Federrath (2020), Beattie et al.
(2020), and Beattie et al. (2021). With a grid resolution of $\gtrsim 512^{3}$
many previous studies have found that the 1-point density statistics are well
converged (e.g. Kowal et al., 2007; Kitsionas et al., 2009; Federrath et al.,
2010; Price & Federrath, 2010; Kritsuk et al., 2011; Federrath, 2013;
Mohapatra et al., 2021), hence using a $\gtrsim 512^{3}$ computational grid is
an appropriate resolution for discussing converged $s$ statistics.

###  4.2 Initial conditions and processing

The initial velocity field is set to
$\mathrm{{\mn@boldsymbol{\mathit{v}}}}(x,y,z,t=0)=(0,0,0)$, with units
$c_{s}=1$, and the density field $\rho(x,y,z,t=0)=\rho_{0}$, with units
$\rho_{0}=1$. The magnetic field is composed out of a fluctuating,
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)$, and a mean-field,
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, as
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)=B_{0}\mathrm{{\mn@boldsymbol{\mathit{\hat{z}}}}}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)$.
The $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field is initialised with
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ threaded through
$\hat{\mathrm{{\mn@boldsymbol{\mathit{z}}}}}$ of the simulations. We consider
systems where $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$
evolves on much larger timescales than
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}(t)$. This could be, for example,
a cool $H_{2}$ region that is pierced by a galactic scale mean-field. The
mean-field would evolve on galactic timescales, $\sim\mathcal{O}(\rm Gyr)$,
and the fluctuating field, akin to the turbulence in the region, on much
shorter timescales, $\sim\mathcal{O}(\rm Myr)$. Hence, the mean, coherent
field is approximately frozen in magnitude in the turbulence. We capture this
by fixing $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ to be
constant throughout the entire simulation runtime, i.e.
$\partial_{t}\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}=\partial_{x_{i}}\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}=0$,
$x_{i}\in\left\\{x,y,z\right\\}$. Even though the mean-field gradients are
zero, field lines, which are dominated by
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, are able to walk
stochastically (magnetic field line wander, Howes & Bourouaine, 2017) in the
plane $\perp\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. The
turbulent field, $\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, evolves self-
consistently from the MHD equations, with
$\mathrm{\left\langle\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}\right\rangle}_{t}=0$,
and
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|\propto\operatorname{\mathcal{M}}\operatorname{\mathcal{M}_{\text{A0}}}$
(Federrath, 2016a; Beattie et al., 2020; Skalidis & Tassis, 2020). $B_{0}$ is
set by using the definition of the Alfvén velocity and
$\operatorname{\mathcal{M}}$,
$B_{0}=2c_{s}\sqrt{\pi\rho_{0}}\operatorname{\mathcal{M}}/\operatorname{\mathcal{M}_{\text{A0}}}$,
where $\operatorname{\mathcal{M}_{\text{A0}}}$ is the desired Alfvénic Mach
number of $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. For the
magnetised simulations we fix this value between
$0.1\lesssim\operatorname{\mathcal{M}_{\text{A0}}}\lesssim 2.0$, ensuring that
the mean-field is sufficiently strong compared to the turbulence (Beattie et
al., 2020; Beattie & Federrath, 2020). From previous experiments, in the
highly-magnetised, supersonic regime,
$|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|/|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|=\mathcal{M}_{A0}^{2}/2$
(Beattie et al., 2020). For the hydrodynamical simulations we set
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}=0$, which implies
$\operatorname{\mathcal{M}_{\text{A0}}}=\infty$. We use these simulations to
compare with the MHD case throughout the study, but they are not the focus of
the study.

We run the simulations from $t/\tau=0-10$. We construct the $s$-field and bin
the data into volume-weighted and mass-weighted distributions. We also extract
the volume- and mass-weighted variance and volume-weighted skewness of the
$s$-field data. We show the time evolution of the volume- and mass-weighted
variance in Figure 4 across the entire run time of the simulations, $0\leq
t/\tau\leq 10$. All results in this study, including $s$-PDF fits, will be
based on time-averages across 51 realisations, within $5\leq t/\tau\leq 10$ to
gather data only when the turbulence is in a statistically stationary state,
unless explicitly indicated otherwise. We show this region with the grey band
in Figure 4. Compared to hydrodynamical turbulence, which takes roughly
$2\tau$ to become stationary (Federrath et al., 2009; Price & Federrath, 2010)
strong mean-field MHD turbulence takes longer to reach a statistically
stationary state, $\sim 5\tau$. Next we discuss the results from the volume-
weighted and mass-weighted $s$-PDFs.

Figure 5: The volume-weighted $s$-PDFs for each of the MHD simulations, shown
with $1\sigma$ fluctuations from averaging over 51 time realisations between
$5\leq t/\tau\leq 10$. The H2013 fit, $p_{\text{H}}(s;\sigma_{s,V}^{2},T)$,
Equations 15-18, is shown with dashed-dots, the MB2019 fit,
$p_{\text{MB}}(s;[\tau_{\rm A}/\tau_{\rm D}]\sigma_{s,V}^{2},f)$, Equation 28,
with dashes, and the lognormal fit, $p_{\rm N}(s;\sigma_{s,V}^{2})$, Equations
1-2, is shown with dots. Overall, the $p_{\text{H}}(s;\sigma_{s}^{2},T)$ and
$p_{\text{MB}}(s;[\tau_{\rm A}/\tau_{\rm D}]\sigma_{s,V}^{2},f)$ fits work
best for low-$\operatorname{\mathcal{M}}$ simulations, where there is an
extended low-$s$ tail, and all fits work equally well at
high-$\operatorname{\mathcal{M}}$. There are significant temporal fluctuations
in the low-$s$ tail, which indicates the volume of the rarefaction waves and
voids is a highly-volatile quantity, compared to the volume of the highest
density structures. We show the $\rho/\rho_{0}=1$ line in purple, which marks
the transition between over- and under-densities. Figure 6: The same as Figure
5 but for the mass-weighted $s$-PDFs. The lognormal fit,
$p_{\text{N},M}(s;\sigma_{s,M}^{2})$, Equations 6-8, is shown with dots. Most
of the lognormal PDFs overestimate the high-$s$ tail, and underestimate the
low-$s$ tail. The high-$s$ tail has significant temporal fluctuations,
revealing that the high-mass, high-density structures: filaments, over-dense
sheets and shocked regions in the turbulence, fluctuate significantly in mass.
Conversely, the lowest density structures, relative to the high-density tail,
do not fluctuate significantly.

##  5 Morphology of the 3D density-PDFs

In §3 we discussed two non-lognormal volume-weighted $s$-PDF models that have
been used to describe logarithmic density fluctuations in supersonic
turbulence. In Figure 5 we plot the volume-weighted $s$-PDFs for each of the
magnetised simulations from Table 1. We fit a H2013 model,
$p_{\text{H}}(s;\sigma_{s}^{2},T)$, Equations 15-18, MB2019 model,
$p_{\text{MB}}(s;\sigma_{s}^{2},T)$, Equations 28 and a lognormal model,
$p_{\text{N},V}(s;\sigma_{s}^{2})$, Equations 1-3 to the $s$-PDFs. The H2013
and MB2019 fits are qualitatively much better than the lognormal model for the
low-$\operatorname{\mathcal{M}}$ simulations, which show strong intermittent
behaviour, for all $\operatorname{\mathcal{M}_{\text{A0}}}$. This demonstrates
the robustness of both models for describing the $s$-fluctuations, with or
without a magnetic field. We draw purple dashed lines where
$\ln\rho/\rho_{0}=0$. This shows that the low-$\operatorname{\mathcal{M}}$
intermittency is not necessarily from an excess of volume-weighted high-
density structures, but from an excess of low-density rarefactions and density
voids. As we discussed in §5, for $\operatorname{\mathcal{M}}\gtrsim 4$ the
fits become, qualitatively, equally as good as the lognormal model. For the
mass-weighted $s$-PDFs, shown in Figure 6, we fit a lognormal model,
$p_{\text{N},M}(s;\sigma_{s}^{2})$, Equations 6-8. In the next subsections we
split the discussion of the volume- and mass-weighted PDFs and their fits into
the low- and high-density tail.

###  5.1 High-density tail

The high-density tail of the PDFs in Figure 5 and 6 trace the shocks,
filaments, sheets and other over-dense structures in the supersonic turbulence
(Robertson & Goldreich, 2018; Federrath et al., 2021). These structures make-
up a small fraction of the total volume, and most of which is filled with
rarefactions and voids that are found in the low-density tail between
$1/\operatorname{\mathcal{M}}\leq\rho/\rho_{0}\leq 1$. However, most of the
mass is found between $1\leq\rho/\rho_{0}\leq\operatorname{\mathcal{M}}$, so
the high-density tail contains most of the mass, but little of the volume
(Robertson & Goldreich, 2018).

Because the high-density tail is set by the shocked gas, it is no surprise
that Pan et al. (2019) found the tail is amplified by the
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ term in the continuity
equation and suppressed by the $\nabla P=c_{s}^{2}\nabla\rho$ term in the
momentum equation for hydrodynamical turbulence. Since
$\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}\sim
v/L=c_{s}\operatorname{\mathcal{M}}/L$ for a fixed $L$ the high-density tail
ought to increases systematically with $\operatorname{\mathcal{M}}$. However,
we find that it is asymptotic at high-$\operatorname{\mathcal{M}}$, which
Beattie et al. (2021) attributes to the total volume of the shocks forming
along $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ reducing
until they no longer have a significant contribution to the spread of the PDF.
Beattie et al. (2021) suggests that at this point the largest contribution is
from weakly compressible MHD shocks that form from field line compressions
(fast magnetosonic shocks), which set the maximum limit of the spread of the
PDF, along with the type of turbulent driving.

We find that the lognormal model (shown with dotted markers) systematically
overestimates the PDF in the high-density tail for the
$\operatorname{\mathcal{M}}=2-4$ simulations, regardless of
$\operatorname{\mathcal{M}_{\text{A0}}}$. This is likely from a mixture of the
strong magnetic pressure, $\nabla B^{2}/(8\pi)$, and the magnetic tension,
$(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}$.
Both the pressure and tension terms scale with $\mathcal{M}_{\rm A0}^{-2}$,
and when $\operatorname{\mathcal{M}_{\text{A0}}}$ is sufficiently small, the
pressure acts to homogenise the density, and the tension acts to smooth out
any curvature in the magnetic field, which in turn reduces any large contrasts
in the density through $B\propto\rho$ flux-freezing (Landau & Lifshitz, 1959;
Mocz & Burkhart, 2018; Yuen & Lazarian, 2020). The overall result is that the
magnetic field limits both high- and low-density fluctuations (Nordlund &
Padoan, 1999; Molina et al., 2012; Hennebelle, 2013; Mocz & Burkhart, 2018;
Beattie et al., 2021).

In this
$|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|\gg|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|$
regime, the strongest shocks are only able to form along the mean magnetic
field (Beattie & Federrath, 2020; Beattie et al., 2020; Beattie et al., 2021).
This results in less available volume where strong shocks are able to form and
therefore fewer over-densities. Because there are fewer shocked regions there
is less dense material that is able to contribute to the high-density tail of
the PDFs. We explore this in more detail in §7 using one-dimensional pencil
beams that probe the flow. However, as $\operatorname{\mathcal{M}}$ increases
to $\operatorname{\mathcal{M}}\gtrsim 10$, the $s$-field on large scales
becomes more isotropic, allowing for more space where shocks can form, filling
up again the high-density tail of the PDF which is then fit well by the
lognormal model. Shock frequencies have previously been found to increases
with $\operatorname{\mathcal{M}}$ and decrease with
$|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|$ (Park & Ryu, 2019b). Beattie &
Federrath (2020) observed large-scale mixing and isotropisation in the density
field by studying the 2D power spectra, which were isotropic on large scales,
and then anisotropic on the scales of the individual shocks for
$\operatorname{\mathcal{M}}\gtrsim 10$.

In the mass-weighted PDFs, shown in Figure 6, we find that the lognormal model
performs best for the high-density tails at low-$\operatorname{\mathcal{M}}$,
low-$\operatorname{\mathcal{M}_{\text{A0}}}$ and worst in the
high-$\operatorname{\mathcal{M}}$, trans-Alfvénic flows. Similar to the
volume-weighted PDF, the lognormal model over-predicts the amount of high-
density structures in the turbulence. Unlike the volume-weighted PDF, the
$1\sigma$ temporal fluctuations (indicated by using transparent bands) from
the mass-weighted PDF are strongest in the high-density tail. By tracking and
analysing individual over-dense regions Robertson & Goldreich (2018) showed
that the temporal fluctuations come from the transient nature of the shock
life-cycle, where shocks form, accumulate, lose mass and are torn apart within
a fraction of $t_{\rm cross}=L/c_{s}$. These are found to a much lesser degree
in the volume-weighted PDFs, which demonstrates how the highest density
structures in the turbulence must dynamically fluctuate through time, in mass,
but be relatively stable in volume. The temporal fluctuations in the high-
density bins increase with $\operatorname{\mathcal{M}}$, supporting the ideas
of Padoan et al. (1997), Robertson & Goldreich (2018) and Mocz & Burkhart
(2018), which state that the internal structure and geometrical properties of
in situ shocks is set by the properties of the turbulence.

###  5.2 Low-density tail

The low-density tail captures the density rarefaction waves and voids in the
turbulence (Kritsuk et al., 2007; Federrath et al., 2010). These are coupled
to the shocks because as the shocks compress the fluid they also evacuate
large regions in the turbulence. These regions occupy the largest volumes in
the flow, but have the smallest mass (Robertson & Goldreich, 2018). For the
volume-weighted PDF the low-density temporal fluctuations are significant,
shown by the large $1\sigma$ values in the low-density tails in Figure 5. This
can be interpreted to mean that the volumes of the low-density regions
fluctuate throughout the life-cycle of the under-densities and hence the
rarefactions are themselves strong sources of intermittency (Passot & Vázquez-
Semadeni, 1998; Kritsuk et al., 2007; Federrath et al., 2010).

For the mass-weighted PDFs in Figure 6, the temporal fluctuations are small in
the low-density tail, hence, we find an interesting and understandable
symmetry between the low-density and high-density regions. The voids and
rarefactions in the turbulence fluctuate in volume, and much less so in mass,
and the high-density regions fluctuate in mass, and much less so in volume.
This is one of the key results from our analysis of the time-averaged PDFs,
unrelated to the Gaussianity.

Similarly to what we found previously for the high-density tail in §5.1, as
$\operatorname{\mathcal{M}}$ increases the lognormal model fits the $s$-PDF
better. At low-$\operatorname{\mathcal{M}}$, and for most
$\operatorname{\mathcal{M}_{\text{A0}}}$, the low-density tail under-predicts
the amount of rarefactions and voids in the turbulence. We have until now
restricted most of our attention to the lognormal fit. In the next section we
consider the more general H2013 and MB2019 fits, which have associated
intermittency parameters, $T$, that encodes the deviation away from perfect
self-similarity through different scales in the $s$-field, and $f$, that
captures the difference between the dynamical times of shocked regions and
rarefactions, respectively.

Figure 7: Top: The volume-weighted skewness of the logarithmic density,
$\mathcal{S}_{s}$ as a function of turbulent eddy turnover time between $0\leq
t/\tau\leq 10$. The arrow describes the general trend we find in the skewness,
which is as $\operatorname{\mathcal{M}}$ increases the skewness tends towards
zero, corresponding to a more lognormal $s$ field. Bottom: $\mathcal{S}_{s}$
as a function of $\operatorname{\mathcal{M}}$, averaged over $5\leq t/\tau\leq
10$, with different line styles for $\operatorname{\mathcal{M}_{\text{A0}}}$,
as indicated in the legend from the top panel. Figure 8: Left: The H2013
intermittency parameter, $T$, as a function of $\operatorname{\mathcal{M}}$,
coloured by $\operatorname{\mathcal{M}_{\text{A0}}}$, split into three
categories: sub-Alfvénic ($\operatorname{\mathcal{M}_{\text{A0}}}<1.0$,
orange), trans-Alfvénic
($1\lesssim\operatorname{\mathcal{M}_{\text{A0}}}\lesssim 2$, aqua) and
hydrodynamical ($\operatorname{\mathcal{M}_{\text{A0}}}=\infty$, black) with
the H2013 hydrodynamical fit shown with the grey band for different driving
parameters, $b=0.3-1$. The Squire & Hopkins (2017) fit is shown with the black
dot-dashed line, which is based on the shape and density jump of
hydrodynamical shocks. A fit to the sub-Alfvénic models,
$T(\operatorname{\mathcal{M}})=0.27^{0.39}_{0.19}\operatorname{\mathcal{M}}^{-0.81^{-0.60}_{-1.03}}$,
is shown with the orange band and trans-Alfvénic,
$T(\operatorname{\mathcal{M}})=1.66^{2.51}_{1.10}\operatorname{\mathcal{M}}^{-3.24^{-2.70}_{-3.78}}$,
with the aqua band. The bands encapsulate the parameter uncertainties for each
of the fits. Right: The mass-weighted $s$ dispersion as a function of volume-
weighted $s$ dispersion, coloured the same as the left plot. The lognormal
model defines the 1:1 line shown in black. The $T=0.05$ and $T=0.5$ lines show
the cases for weak intermittency and strong intermittency, respectively, from
the H2013 model. We show a band that encapsulated the H2013 fit with the grey
band. In both plots, the $\operatorname{\mathcal{M}}=2-4$ simulations show
intermittency comparable to the H2013 fits. However, for
$\operatorname{\mathcal{M}}\gtrsim 4$ the logarithmic density fields become
significantly more Gaussian, with shrinking $T$ and dispersion relation closer
to the 1:1 lognormal line, illustrating a key difference between
magnetohydrodynamical and hydrodynamical $s$ statistics.

##  6 Intermittency in the s-PDF

We qualitatively observed in the volume- and mass-weighted logarithmic
density-PDFs that the moderate- to high-$\operatorname{\mathcal{M}}$
simulations are fit well by a lognormal distribution, and the
low-$\operatorname{\mathcal{M}}$ simulations show significant deviation from
the lognormal fit in both the low- and high-density tail. Now we aim to
quantify these non-Gaussian features, focusing primarily on the volume-
weighted PDF where the intermittency parameters computed from the H2013 and
MB2019 models are valid.

###  6.1 Skewness

We quantify the intermittency of the volume-weighted $s$-PDFs in four ways.
First, we calculate the $3^{\rm rd}$ central moment, the skewness of $s$. It
is defined as

$\displaystyle\mathcal{S}_{s}=\frac{\mathrm{\left\langle\left(s-\mathrm{\left\langle
s\right\rangle}\right)^{3}\right\rangle}}{\sigma_{s,V}^{3}}.$ (32)

This gives $\mathcal{S}_{s}=0$ for the normal distribution,
$\mathcal{S}_{s}>0$ for a distribution with an elongated tail towards the
high-density tail and $\mathcal{S}_{s}<0$ for an elongated low-density tail.
Because we compute the central moment the magnitude of $\mathcal{S}_{s}$ is
weighted with respect to the volume-weighted variance, $\sigma_{s}^{2}$. Hence
for $|\mathcal{S}_{s}|\geq 1$ the skewness is just as, or more important for
determining the morphology of the $s$-PDF than $\sigma_{s}$.

We plot the temporal evolution of the skewness for the volume-weighted
distribution in the top panel of Figure 7, coloured by different
$\operatorname{\mathcal{M}}$ and illustrated with different line styles for
different $\operatorname{\mathcal{M}_{\text{A0}}}$, using the same style as
Figure 4. $\mathcal{S}_{s}$ becomes stationary for $t\gtrsim 5\tau$, which we
demonstrate in the top panel. Before $t=5\tau$, $\mathcal{S}_{s}$ the skewness
increases in magnitude with $\operatorname{\mathcal{M}}$, which is the
opposite trend we found in averaged PDFs. We show the $\mathcal{S}_{s}$
averaged over $5\leq t/\tau\leq 10$ in the bottom panel of Figure 7. Tracing
the different line styles through the plot, which correspond to different
$\operatorname{\mathcal{M}_{\text{A0}}}$, we can immediately see that the
absolute value of the skewness decreases and tends towards
$\mathcal{S}_{s}\approx 0$ as $\operatorname{\mathcal{M}}\gg 1$ for the sub-
Alfvénic simulations, and $\mathcal{S}_{s}\approx 0.2$ for the trans-Alfvénic
simulations. For low-$\operatorname{\mathcal{M}}$ $\mathcal{S}_{s}$ varies
between $-0.3$ and $-0.5$, regardless of
$\operatorname{\mathcal{M}_{\text{A0}}}$, meaning that the peak of the PDF is
shifted towards the higher densities, with an elongated tail into the low-
densities, consistent with our qualitative findings in §5, and the type of
intermittency described in the H2013 and MB2019 models. The trans-Alfvénic
simulations show some deviation to higher values than $\mathcal{S}_{s}=0$ but
because $|\mathcal{S}_{s}|$ decreases with $\operatorname{\mathcal{M}}$ all of
the simulations share the same trend towards more Gaussian statistics. Since
we find $|\mathcal{S}_{s}|<1$, $\mathcal{S}_{s}$ is never more important than
$\sigma_{s}^{2}$ for describing the morphology of the distribution in these
strong mean-field, supersonic flows, i.e. a lognormal model is not strictly a
bad empirical approximation for the $s$-PDF. However, because any non-zero
$\mathcal{S}_{s}$ is a result of asymmetry in the PDF, models that compute
astrophysical properties, such as star-formation rate or efficiency from the
high-density tail (Krumholz & McKee, 2005; Hennebelle et al., 2011; Padoan &
Nordlund, 2011; Federrath & Klessen, 2012, 2013; Burkhart, 2018), which is
significantly overestimated by the lognormal model (Figures 5 and 6) may
overestimate the amount of shocked and over-dense gas, especially for
lower-$\operatorname{\mathcal{M}}$, magnetised MCs.

Our $\mathcal{S}_{s}$ results may seem to contradict previous findings of
Kowal et al. (2007) and Burkhart et al. (2009), but they analysed the higher-
order statistics of not the logarithmic densities, but the linear densities,
$\rho/\rho_{0}$. The higher-order statistics of $\rho/\rho_{0}$ and $s$ need
not be alike; and in fact we have shown here they are significantly different.
Our findings are consistent with those of Molina et al. (2012), who found that
magnetised gas densities are much more lognormal than in hydrodynamical
turbulence, but did not quantify it in detail.

Figure 9: The MB2019 intermittency parameter, $f$, written in terms of the
reduction factor between the dynamical timescales for the rarefactions and the
shocked gas density, $1+3f/2$, as shown in Equation 26, as a function of
$\operatorname{\mathcal{M}}$. At low-$\operatorname{\mathcal{M}}$ we find that
the timescales of shocked gas evolve approximately twice as fast as low-
density gas, but as $\operatorname{\mathcal{M}}$ increases the timescales
homogenise, resulting in more Gaussian $s$-PDFs.

###  6.2 H2013 intermittency parameter

To understand the transition from the low-$\operatorname{\mathcal{M}}$
intermittency to the high-$\operatorname{\mathcal{M}}$ Gaussianity we see both
in the skewness statistics and directly from the PDFs we plot $T$ as a
function of $\operatorname{\mathcal{M}}$, coloured by
$\operatorname{\mathcal{M}_{\text{A0}}}$, in the right panel of Figure 8. We
overlay an empirical model for the hydrodynamical turbulence that H2013 fit to
the simulation data in grey, $T(\operatorname{\mathcal{M}})\ \approx
0.05\,b\operatorname{\mathcal{M}}$. The width of the band demonstrates the
impact of the turbulent driving parameter on the fit (where $b$ falls in the
range $b=0.3-1$, corresponding to solenoidal to compressive turbulence,
respectively; see Federrath et al. 2010). We also show the Squire & Hopkins
(2017) model for $T$ with the black dot-dashed line, which is based upon the
volume and mass conservation, the size and density contrast of hydrodynamical
shocks,
$T(\operatorname{\mathcal{M}})=\kappa(1-\operatorname{\mathcal{M}}^{-2})$,
where $\kappa=0.2$ is associated with the shock width being a fraction of
$\ell_{\rm s}$ in the turbulence (Federrath, 2016b; Xu et al., 2019; Federrath
et al., 2021). We also fit our own empirical power-law models of the form
$T(\operatorname{\mathcal{M}})\propto\operatorname{\mathcal{M}}^{\alpha}$ to
the sub-Alfvénic and trans-Alfvénic data. The bands are associated with the
uncertainty in the power-law parameter estimates.

Consistent with our qualitative observations of the PDFs, we find that
regardless of $\operatorname{\mathcal{M}_{\text{A0}}}$, the
$\operatorname{\mathcal{M}}\lesssim 4$ simulations have H2013-like
intermittency, shown by how the low-$\operatorname{\mathcal{M}}$ $T$ values
cluster around the grey model band. We will explore the origin of these
intermittent structures in the next section, but for now we just state that
the logarithmic density must have some non-Gaussian structures, which lead to
a similar value of $T$ as hydrodynamical simulations (shown in black) of the
same $\operatorname{\mathcal{M}}$. For all of the MHD simulations $T$ peaks in
the low-$\operatorname{\mathcal{M}}$ and then decreases towards $T=0$, a
Gaussian field, which is the opposite trend found in the hydrodynamical
simulations. Note that this does not imply that either H2013 or Squire &
Hopkins (2017) are incorrect, rather that the conclusions made in those
studies should not be applied to highly-supersonic, magnetised turbulent
flows, assuming that the H2013 PDF can be used to interpret our MHD data. The
sub-Alfvénic $T$ values (orange) are fit by a power law,
$T(\operatorname{\mathcal{M}})=0.27^{0.39}_{0.19}\operatorname{\mathcal{M}}^{-0.81^{-0.60}_{-1.03}}$
and the trans-Alfvénic (aqua) by
$T(\operatorname{\mathcal{M}})=1.66^{2.51}_{1.10}\operatorname{\mathcal{M}}^{-3.24^{-2.70}_{-3.78}}$.
Hence, the sub-Alfvénic simulations seem to be more intermittent than the
trans-sonic $\operatorname{\mathcal{M}_{\text{A0}}}$ simulations, consistent
with what we found in the skewness. We interpret this with the H2013
phenomenology as follows: the average jump in density between nested scales in
the turbulence, $T$ (Equation 19) scales with
$\sim\operatorname{\mathcal{M}}^{-1}$ in the sub-Alfvénic simulations, and
$\sim\operatorname{\mathcal{M}}^{-3}$ in the trans-Alfvénic simulations. We
discuss this in much more detail in §8.

One of the key motivations for H2013 was to address the discrepancy between
the mass-weighted and volume-weighted variances (Equation 20), which, as we
highlighted in §1, should be identically equal,
$\sigma^{2}_{s,V}=\sigma^{2}_{s,M}$, if a lognormal theory describes the
density fluctuations. As another independent measure of the intermittency we
therefore plot $\sigma_{s,M}$ as a function of $\sigma_{s,V}$ in the left
panel of Figure 8. We plot a grey band that encapsulates the empirical
relation that H2013 fit to the hydrodynamical data, the weak ($T=0.05$, blue,
dotted) and strong ($T=0.5$, red, dotted) intermittency isocontours and the
$\sigma_{s,V}=\sigma_{s,M}$ line in black. We find that the
low-$\operatorname{\mathcal{M}}$ density structures follow more closely the
H2013 intermittency relation, but as $\operatorname{\mathcal{M}}$ increases
the density field becomes lognormal, tending towards the one-to-one lognormal
line with a small amount of scatter, consistent with our previous
measurements.

Figure 10: An example of profiles from two typical pencil beams along (left)
and across (right) $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$
in the M2MA01 simulation. The plot is organised such that the logarithmic
density is the blue solid line, corresponding to the left, blue axis, the
kinetic velocities are either solid or dashed red lines, corresponding to the
right red axis, in units of $c_{s}$, and Alfvén velocity $\delta
v_{A}=|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|/\sqrt{4\pi\rho}$, is a
dot-dashed red line, also sharing the red axis with the kinetic velocity.
Integer values of the global $\sigma_{s}$ for the $s$ field are shown as
horizontal lines across the plots, indicating how extreme the $s$-fluctuations
are. The $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ direction
is indicated in the top-left of each panel. Left: The parallel pencil beam
reveals large-scale, over-densities with an exponential tail, caused by
converging flows (in the grid frame) along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, consistent with
type C filament formation (Abe et al., 2020). The structure of the over-
density is qualitatively similar to atmospheric shocks studied in Robertson &
Goldreich (2018) and Mocz & Burkhart (2018). A sawtooth shockwave travels down
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ and is coupled to a
volume-filling rarefaction wave, travelling in the opposite direction. At the
shockwave interface fast magnetosonic (travelling at $\sim c_{s}$ with the
mean-field subtracted) compression waves are excited that also travel down
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. Right: The
perpendicular pencil beam reveals high-frequency logarithmic density features
that are correlated with velocity streams along the field, $\delta
v_{\parallel}/c_{s}$. Large-scale vortical structures are seen in the $\delta
v_{\perp}/c_{s}$ profile, with two vortices developing on
$[-L_{\perp}/2,L_{\perp}/2]$. The shear Alfvén wave fluctuations, $\delta
v_{\rm A,\perp}/c_{s}$, are extremely small-scale compared to the kinetic
fluctuations. Figure 11: Positions of the 1D pencil beams shown in a 2D slice
of the logarithmic density, which are randomly chosen and represent a typical,
average pencil beam through the turbulence. The along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$-field pencil beam
(left panel of Figure 10) is shown with the vertical line and the across
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$-field pencil beam
(right panel of Figure 10) is shown with the horizontal line. The along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$-field pencil beam
has been shifted into the centre of the 2D map for ease of viewing.

###  6.3 MB2019 intermittency parameter

Using the MB2019 fits we extract the intermittency parameter, $f$ (Equation
26), and plot $1+3f/2$ as a function of $\operatorname{\mathcal{M}}$, with
different line styles for $\operatorname{\mathcal{M}_{\text{A0}}}$ in Figure
9. We use $1+3f/2$ because this is the reduction factor between the dynamical
timescales for $s>s_{0}$, $\tau_{\rm A}$ and $s\leq s_{0}$, $\tau_{\rm A,0}$
structures in the MHD turbulence. The horizontal blue line shows the case
where both the under-density and over-density dynamical timescale are the
same, i.e. when $f=0$. At low-$\operatorname{\mathcal{M}}$, $\tau_{\rm
A}\approx\tau_{\rm A,0}/2$, hence the shocked structures are operating on
timescales twice as fast as the voids and rarefactions. This results in the
most non-Gaussian $s$-PDFs at low-$\operatorname{\mathcal{M}}$, similar to
what we found using the skewness and $T$ parameter. Consistent with our
previous measurements, as $\operatorname{\mathcal{M}}$ increases the field
becomes more Gaussian, which can be interpreted as the two dynamical
timescales becoming equal. The sub-Alfvénic simulations take the longest to
reach Gaussian-like statistics, and the trans-Alfvénic, the fastest, the same
as what was found in the $T(\operatorname{\mathcal{M}})$ plots in Figure 8.
Again, we discuss and compare these results with the other statistics in more
detail in §8.

###  6.4 Summary of §6

To summarise this section, we compute the logarithmic density intermittency
with four independent measures: (i) the skewness, $\mathcal{S}_{s}$ (Equation
32), (ii) the H2013 $T$ parameter (Equation 19), (iii) the
$\sigma_{s,M}-\sigma_{s,V}$ relation (Equation 20) and (iv) the MB2019 $f$
parameter (Equation 26). All statistics point to the conclusion that
$s$-intermittency in strong mean-field MHD turbulence is significantly
different, and in fact opposite (decreases with $\operatorname{\mathcal{M}}$)
compared to purely hydrodynamical turbulence (increases with
$\operatorname{\mathcal{M}}$), which is a key result of this study. However,
at low-$\operatorname{\mathcal{M}}$ ($\operatorname{\mathcal{M}}\lesssim 4$)
there considerable H2013-like intermittency present in the $s$-field, i.e.,
$T$ values that are comparable to $T$ values computed for hydrodynamical
turbulence. We now turn our attention to identifying what physics leads to
this low-$\operatorname{\mathcal{M}}$ intermittency by identifying the
intermittent structures in real-space, following the philosophy from Imara et
al. (2021) – even theoretical astrophysics is, at heart, an observational
endeavour that relies upon us being able to visualise the structures we are
interested in.

##  7 The physics of low-$\operatorname{\mathcal{M}}$ intermittency

Beattie & Federrath (2020) showed that the density structures in
$\operatorname{\mathcal{M}}\approx 2-4$, sub- to trans-Alfvénic turbulence are
anisotropic on all length scales. On average, equi-power surfaces in the power
spectrum reveal that they are stretched along the mean magnetic field444Note
that on average for $\operatorname{\mathcal{M}}\gtrsim 4$ turbulence the
density anisotropy is across and not along the mean magnetic field due to
formation of perpendicular, high-density filaments; see Figures 2 and 3 in
Beattie & Federrath (2020).. For this reason, this regime can naturally be
decomposed into two domains: across ($\ell_{\perp}$) and along
($\ell_{\parallel}$)
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, which on average
defines a symmetrical axis ($\delta v_{x}\sim\delta v_{y}$, $\delta
B_{x}\sim\delta B_{y}$) in the turbulence555From time-averaging experiments we
find the symmetry is an average property of the turbulence and events, for
example, strongly interacting vortices, can lead to velocity streams in
preferential directions across the field, “breaking” the symmetry around
<EMAIL_ADDRESS>Since the M2MA01
simulation is the most anisotropic (Beattie & Federrath, 2020; Beattie et al.,
2020), this is where such a simple domain decomposition will work the best.
Also, since the M2MA01 is significantly intermittent (of order the equivalent
hydrodynamical simulation), we focus our analysis in this section primarily on
the local dynamics and intermittency of this simulation.

To analyse the real-space intermittent structure we take pencil beams (1D
slices) in the $\ell_{\perp}$ and $\ell_{\parallel}$ directions to examine
examples of the frozen-in-time local dynamics along each direction. In Figure
10 we show pencil beams for the logarithmic density (blue), kinetic turbulent
velocity, $\delta v$, (red, solid, dashed) and Alfvén turbulent velocity,
$\delta v_{\rm
A}=|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|/\sqrt{4\pi\rho}$ (red, dot-
dashed) profiles. To make sure the beams we plot are representative of the
average, typical dynamics we check the other $1024$ possible pencil beams and
find that these represent that average dynamics in the M2MA01 simulation
faithfully (i.e., they are not outliers or peculiar in any way). The left
panel illustrates a beam along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ (magnetic field
direction is indicated in the top, left) and the right panel, across. To help
the reader contextualise the pencil beams, in Figure 11 we indicate the
positions of each of the pencil beams in the frame of a 2D logarithmic density
slice. The vertical dashed line shows the position in the left panel of Figure
10 and horizontal in the right panel.

###  7.1 Along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$-field pencil beam

In the left panel we see an in situ sawtooth shock coupled to an over-density
forming perpendicular to
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ (solid, blue) from
converging kinetic velocities (solid, red line) along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, in the simulation
frame. Converging flows are a popular hypothesis for filament formation in
molecular clouds (e.g. Chen et al., 2020; Bonne et al., 2020), and in
particular for sub-Alfvénic, compressible gases (Padoan & Nordlund, 1999; Chen
& Ostriker, 2014; Abe et al., 2020). However, the interplay between the
magnetic field and gravity orienting the filaments with respect to the
magnetic field is still debated in the literature (Planck Collaboration et
al., 2016; Soler et al., 2017; Tritsis et al., 2018; Mocz & Burkhart, 2018;
Heyer et al., 2020; Pillai et al., 2020; Barreto-Mota et al., 2021;
Girichidis, 2021). We show here that perpendicular filaments caused by
converging flows along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ can be seeded
without the presence of gravity, such that dense filaments are primarily
created (or at least seeded) by turbulent compression (Federrath, 2016b). The
formation scenario is as follows: the strong
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$-field acts by
constraining the turbulence along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ (suppressing
diagonal velocity modes), compressible modes, which are injected isotropically
from our turbulent driving source, are self-organised along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, and the gas
density along the magnetic field becomes shocked, forming over-dense regions
across the $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. Beattie
et al. (2020) found that the self-organisation of the compressible modes in
the velocity was a natural repercussion of the strong
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, where in the limit
of
$|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|\gg|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|$
the total velocity divergence reduces to just the along-the-field component,

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

where $\partial_{\parallel}$ is the derivative along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, and there is only
negligible velocity divergence across the field,
$\nabla_{\perp}\cdot\mathrm{{\mn@boldsymbol{\mathit{v}}}}_{\perp}\approx 0$,
and the compressive motions are along the $\ell_{\parallel}$ beam. It is also
worth noting how similar this phenomenology is to the model for weakly-
compressible MHD turbulence proposed by Bhattacharjee et al. (1998).

We highlight and emphasise here that in situ turbulent over-densities are not
morphologically the same as simple shock tube models, which is shown both by
the over-density (solid, blue line) and the shock (solid, red line) in the
left panel of Figure 10, which have exponential tails, and sawtooth structure,
respectively. This demonstrates that in situ over-densities have internal
structures (e.g., scale heights), i.e., are not just discontinuities, and that
in the strong $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$
regime Burgers (1948)-like turbulence develops along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. The over-density
we highlight is qualitatively similar to the hydrodynamical atmospheric shocks
studied in Robertson & Goldreich (2018) and extended to MHD in Mocz & Burkhart
(2018). Robertson & Goldreich (2018) explains the exponential tail by
balancing the ram pressure, which is $\propto\rho v^{2}$, with the ambient
pressure gradient $\nabla P$. We also note how the over-density is
approximately consistent with the simple hydrodynamical model proposed by
Padoan & Nordlund (2011) for the thickness of the post-shock layers. Like the
Robertson & Goldreich (2018) model, Padoan & Nordlund (2011) uses a pressure
balance to derive a characteristic width of a shock. For the sawtooth shock we
find, $\ln\rho_{\parallel}/\rho_{0}=\ln\operatorname{\mathcal{M}}^{2}\sim
1.4$, and shock width
$L_{\parallel}/\operatorname{\mathcal{M}}^{2}=L_{\parallel}/4$
($\operatorname{\mathcal{M}}=2$ for the M2MA01 simulation), which defines a
rectangle that captures the basic geometry of the over-density captured in the
pencil beam. Many models for the logarithmic density variance (e.g. Molina et
al., 2012; Federrath, 2015; Nolan et al., 2015; Beattie et al., 2021) rely
upon this model for relating the shock jump relations to the variance. Beattie
et al. (2021) in particular modelled the variance along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ with hydrodynamical
shock jump and width conditions, which is supported here.

At the sawtooth shock front (at roughly $-L_{\parallel}/8$ in the left panel
of Figure 10) some fast magnetosonic666Note here that the MHD waves speeds in
the parallel direction are $v=\delta v_{A}+v_{A0}$, but we show the mean-
subtracted speeds, $\delta v_{A}$ in Figure 10. This means that the waves we
find in the parallel direction travel at (group or phase) velocities
$v\approx\sqrt{c_{s}^{2}+v_{A0}^{2}}$, consistent with fast magnetosonic
waves, which travel at these speeds isotropically in the fluid, at least for
trans-sonic turbulence (Makwana & Yan, 2020). compression ($\delta
B_{\parallel})$ waves are excited (shown with dot-dashed, red line), which
propagate down $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$.
These magnetic compression waves become amplified in the sub-Alfvénic regime
when shock formation becomes highly oriented along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ (Beattie et al.,
2020), favouring the excitation of high-amplitude, fast magnetosonic
compression waves compared to perpendicular shear Alfvén waves. Beattie et al.
(2020) found that the compression waves preferentially oppose the direction of
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, reducing the total
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, and allowing pressure-supported
vortices to be maintained in the turbulence through $\sim\nabla B^{2}$. Away
from the shock front the $\delta B_{\parallel}$ gradient becomes very small,
$|\partial B/\partial\ell_{\parallel}|\sim 0$. This has repercussions for the
magnetic tension. The magnetic tension written in normal
$\mathrm{{\mn@boldsymbol{\mathit{\hat{e}}}}}_{n}$ and tangential
$\mathrm{{\mn@boldsymbol{\mathit{\hat{e}}}}}_{t}$
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field coordinates is,

$\displaystyle(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}=B\frac{\partial
B}{\partial\lambda}\mathrm{{\mn@boldsymbol{\mathit{\hat{e}}}}}_{t}-\kappa
B^{2}\mathrm{{\mn@boldsymbol{\mathit{\hat{e}}}}}_{n},$ (34)

where $\lambda$ is the arc length parameterisation variable along
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, and $\kappa$ is the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field curvature. If gradients along
the field are small, as they are everywhere except at the shock interface in
the left panel of Figure 10, then clearly

$\displaystyle(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}\approx-\kappa
B^{2}\mathrm{{\mn@boldsymbol{\mathit{\hat{e}}}}}_{n}\implies(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}\perp\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}},$
(35)

and the tension becomes purely a restoring force for the Alfvénic
fluctuations. At the shock interface the parallel
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field gradient is very large compared
to the amplitudes of the shear Alfvén waves as we show in the right panel of
Figure 10, and

$\displaystyle(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}\approx
B\frac{\partial
B}{\partial\lambda}\mathrm{{\mn@boldsymbol{\mathit{\hat{e}}}}}_{t}\implies(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn<EMAIL_ADDRESS>(36)

This demonstrates that in the sub-Alfvénic mean-field turbulence regime the
$(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
term becomes anisotropic. This has repercussions for both the energy cascade
in MHD turbulence, where
$(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
is found to play a role in suppressing the kinetic energy transfer (Grete et
al., 2020), and the nature of solenoidal modes in the velocity field, where
$(\mathrm{{\mn@boldsymbol{\mathit{B}}}}\cdot\nabla)\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
is found to be a strong source of vorticity generation when the
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ is strong (Lim et
al., 2020).

Figure 12: The same as the left panel of Figure 10 but for the events in the
turbulence that give rise to the largest under-densities across $5\leq
t/\tau\leq 10$, i.e. a proxy for the events in the turbulence that contribute
the most to the H2013-like intermittency in the $s$-PDF. The panels are
organised such that top-left panel corresponds to the 1D pencil beam at
$t/\tau=5$, and the bottom-right panel is $t/\tau=9.8$. We find that most of
the extreme under-densities ($s_{\parallel}=\ln(\rho_{\parallel}/\rho_{0})$,
blue, solid lines) correspond to single, strong sawtooth shocks (red, solid
lines) that evacuate mass from large-in-volume regions along the direction of
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. Sometimes these
intermittent events are also coupled with large parallel magnetic field
fluctuations caused by fast magnetosonic waves (red, dashed lines). Figure 13:
The space-time diagrams of $s$-fluctuations for the four 1D shock experiments:
(i) single shock, (ii) two shocks, (iii) two shocks, constant driving, (iv)
many shocks, constant driving. We use these experiments to explore what kinds
of 1D shock dynamics lead to the largest non-Gaussian components in the
$s$-field, analogous to the dynamics that are constrained to along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ in the 3D sub-
Alfvénic mean-field turbulence. $t/t_{\rm cross}=0$ is at the bottom of the
space-times, and $t/t_{\rm cross}=2$ is at the top. Figure 14: Skewness of $s$
as a function of $t/t_{\rm cross}$ for each of our four 1D shock experiments.
The generation of single shocks results in the extreme, negatively skewed
$s$-PDFs at the time that the shock is generated. Many-shock simulations give
rise to stationary, slightly negatively-skewed distributions that fluctuate
through $\mathcal{S}_{s}=0$.

###  7.2 Across
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$-field pencil beam

Turning our attention to the pencil beam perpendicular to
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ in the right panel
of Figure 10, we find significant amounts of high-frequency density structure,
which do not correlate with the perpendicular kinetic velocity (dashed, red
line). The perpendicular velocities trace the structure of two large,
supersonic, rigid body ($v_{\perp}/c_{s}\propto\ell/L_{\perp}$) counter-
rotating vortices in the plane $\ell_{\perp}$ around
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. These vortices are
visualised in 3D, through the velocity streamlines in Figure 3, and previously
in the 3D sub-Alfvénic turbulence renderings in Figure 2 in Beattie et al.
(2020). The large density contrasts, up to roughly
$\operatorname{\mathcal{M}}^{2}$, are more strongly correlated with
$v_{\parallel}/c_{s}$ (red, solid line), which we observed were caused from
converging flows along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. We conclude that
the largest over-densities in the strong-$B_{0}$ MHD regime are caused by
convergent flows along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, consistent with
the anisotropic $\rho/\rho_{0}$-variance model proposed by Beattie et al.
(2021). The shear Alfvén wave fluctuations (dot-dashed, red line) are
extremely weak and on small scales compared to both the kinetic velocities and
the magnetic compression waves in the left panel, consistent with the ratios
between $\delta B_{\parallel}/\delta B_{\perp}>2$ found for the sub-Alfvénic
regime in Beattie et al. (2020). This suggests that models like Goldreich &
Sridhar (1995), which rely on a critical balance between the turbulent eddies
and shear Alfvén wave packets may only be appropriate on very small scales in
supersonic turbulence, where the flow becomes strongly Alfvénic (with respect
to the Aflvénic fluctuations and not the mean field). This illustrates what is
a commonly not considered by some of the astrophysical turbulence community –
sub-Alfvénic mean-field turbulence does not mean that the shear Alfvén waves
are dominating the fluid dynamics. In fact, the perpendicular fluctuating
field is clearly super-Alfvénic with respect to the kinetic turbulence, and
all of the magnetic energy is stored in
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$.

###  7.3 Intermittent events and 1D shock experiments

Now that we understand how to interpret the pencil beams, we turn our
attention to using them to seek the events that contribute to the $s$-PDF
intermittency. The H2013-like intermittency manifests itself in the volume-
weighted $s$-PDF by elongating the low-$s$ tail, and truncating the high-$s$
tail. As a proxy for understanding what kind of flow events give rise to the
elongated low-$s$ tail we create pencil beams that intersect through the grid
element in the simulation that has the minimum value of $s$, $s_{\rm min}$. We
construct parallel pencil beams, because as we noted in the previous section,
the largest density contrasts are formed by compressive motions along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. We show these
pencil beams in Figure 12 for the density, and kinetic and fast magnetosonic
field fluctuations (with the same colouring and line style scheme as Figure
10), for time realisations between $t/\tau=5-9.8$, where $t/\tau=5$ is shown
in the top-left corner, and $t/\tau=9.8$ in the bottom-right corner of the
plot. We annotate the skewness, $\mathcal{S}_{s}$, in the top-left corner of
the plot.

Our first conclusion is that many of the 49 visualised extreme low-$s$ events
are formed by single, strong shocks, which is especially apparent in the time
series from $t/\tau=6.2-7.5$. In this time-correlated777Note that there need
not be any time correlation, by construction, but clearly some of the most
under-dense events are so extreme that they define $s_{\rm min}$ for a full
$\sim\tau=L/(2c_{s}\operatorname{\mathcal{M}})$. set of pencil beams we see a
$\operatorname{\mathcal{M}}\approx 3$ shock (a $\sim 1.5\sigma$ event in the
velocity field for the M2MA01 simulation) moving down the field, evacuating a
large void in the turbulence, where $s_{\rm min}$ resides. The strong shock
also facilitates an equally large fast magnetosonic modes. Clearly this hints
at density, magnetic and velocity intermittency being linked through these
strong-shock events, and explains the physical origins of the fast
magnetosonic wave intermittency studied in Ho & Lazarian (2021), i.e., as a
side effect of shock generation in the supersonic plasma. The median, 16th and
84th percentiles for the skewness of all of the intermittent beams is
$\mathcal{S}_{s}=-0.70^{-1.35}_{-0.31}$, almost a factor of 2 more skewed than
the average $\mathcal{S}_{s}$ calculated for the global $s$-field in M2MA01 in
§6. This demonstrates that these events contribute significantly to the non-
Gaussian features of the $s$-PDF. The most skewed beam is at $t/\tau=5.5$,
with $\mathcal{S}_{s}=-1.9$, where a strong, almost
$\operatorname{\mathcal{M}}\approx 4$ ($2\sigma$ event in the velocity) shock
compresses an under-dense region into a thin, low-volume rarefaction.

A few of the panels show multi-shock interactions along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ give rise to
$s_{\rm min}$, like, for example, at $t/\tau=5.0$ and $t/\tau=8.4$, where two
sawtooth shocks travel along the field. To confirm that strong, single shocks,
give rise to the most non-Gaussian features we perform four simple numerical
experiments. We give details about the numerical experiments in §A. To
summarise the experiments, we follow the 1D setup in Mocz & Burkhart (2018),
solving the compressible Euler equations on a periodic domain $x\in[0,2L]$ and
creating shocks by perturbing the velocity field with a Gaussian pulse that
has an amplitude which sets the root-mean-squared
$\operatorname{\mathcal{M}}$. We choose $\operatorname{\mathcal{M}}\approx 2$,
mimicking the characteristic shock velocities from the M2MA01 simulation. The
four experiments are:

  1. 1.

a single shock travelling across the domain,

  2. 2.

two shocks colliding at the centre of the domain,

  3. 3.

two shocks repeatedly injected,

  4. 4.

multiple shocks repeatedly injected.

We show the space-time diagrams for $s$, in each of the experiments in Figure
13, from left-to-right for experiments (i)-(iv), and the $s$-profiles as a
function of $t/t_{\rm cross}$ in Figure 16. In Figure 13, the over-densities
are coloured in red, which trace the shock fronts and the under-densities in
blue, which are almost stationary in the first three experiments. We evolve
the experiments for $2t_{\rm cross}$, which is shown by going up the vertical
axis as indicated with the $t/t_{\rm cross}$ annotation. In the experiments
(iii) and (iv) one can see the that repeated pulses in the $s$-field give rise
to intricate, interacting shock networks, analogous to Burgers (1948)
turbulence (Mocz & Burkhart, 2018).

We compute the $\mathcal{S}_{s}$ for each of the experiments as a function of
$t/t_{\rm cross}$, shown in Figure 14. We only take statistics for $0\leq
t/t_{\rm cross}\leq 0.8$, the time for which the periodic boundaries in the
first three experiments do not play a role in the flow. We find that the
largest $|\mathcal{S}_{s}|$ values, $\mathcal{S}_{s}\approx-5$, are associated
with the single shock experiment, consistent with our previous intermittent
pencil beam analysis. For both the one-shock (i) and two-shock (ii), (iii)
experiments, the initial generation of the over-dense regions, where both the
under-density and over-density occupy the smallest volumes, have the most
negative $\mathcal{S}_{s}$, similar to the features we saw in the most non-
Gaussian event in Figure 12 at $t/\tau=5.5$. The interaction between the two
shocks in experiment (ii), which happens at $t/t_{\rm corr}\approx 0.4$,
results in $\mathcal{S}_{s}>0$ as the shocked gas fills the volume of the
domain after the interaction. The two-shock experiment with repeated
injections almost matches experiment (ii), but the $\mathcal{S}_{s}$ reduces
after each pulse in the driving, before slowly relaxing back to a more
Gaussian field. The many-shock experiment settles into a stationary state
quickly, fluctuating around $\mathcal{S}_{s}=0$, slightly tending towards a
negative value of $\mathcal{S}_{s}$, similar to our global $\mathcal{S}_{s}$
statistics we computed in §6. Even in this experiment it is probably the rare,
single, strong shocks that we can see traced by the red, over-dense gas, that
contribute to the non-Gaussian components by creating the strongest under-
densities in the $s$-field.

The key conclusion from this analysis is that we find strong, single shocks in
the 3D, low-$\operatorname{\mathcal{M}}$, sub-Alfvénic mean-field turbulence
that give rise to some of the lowest values in $s_{\rm min}$. We show that it
is these strong, single shocks that give rise to the greatest non-Gaussian
contributions to the $s$-PDF. We have now solidified our understanding of the
$s$-intermittency for the $\operatorname{\mathcal{M}}\lesssim 4$ sub-Alfvénic
mean-field turbulence, and now turn focus on the
$\operatorname{\mathcal{M}}\gtrsim 4$ regime.

##  8 Phenomenology of high-$\operatorname{\mathcal{M}}$ Gaussianisation

In this section we outline a number of reasonable explanations for the
high-$\operatorname{\mathcal{M}}$ lognormality that we observed in the trans-
and sub-Alvénic mean-field regime in Figures 7, 8 and 9. The arguments for
this are based upon the self-similarity of $s$ and timescale phenomenologies
developed by H2013 and MB2019, respectively.

###  8.1 Self-similarity interpretation

Because the H2013 $s$-PDF model is characterised by the $T$ parameter
(Equations 19 & 20), discussed in §3, which is an average measure of how much
the $s$-field in the turbulence deviates from perfect self-similarity at each
neighbouring spatial scale, $\mathrm{\left\langle
s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}}\right\rangle}_{L/\Gamma^{n}}$, we call
this the self-similarity (SS) phenomenology. The H2013 SS phenomenology may be
the most robust of the phenomenologies we discuss, because the formulation
relies upon the Castaing (1996) method of probing probabilities for very
general steady-state classes of multiplicative-random-relaxation processes
within the turbulence cascade. Under this interpretation, the decreasing $T$
with high-$\operatorname{\mathcal{M}}$ in the magnetised flows is a simple
statement about the $s$-field becoming closer to perfect self-similarity. As
Squire & Hopkins (2017) points out, this does not mean that the underlying
statistics are Gaussian (we have shown that this cannot be the case if one
considers mass conservation or supersonic turbulence, which is not scale-free
in §2.2), but rather, it means that the flow is full of small-in-volume,
small-amplitude density contrasts, making $\mathrm{\left\langle
s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}}\right\rangle}_{L/\Gamma^{n}}$ approach
zero, and hence decreasing $T$.

We believe the $T$ trends we find in Figure 8 can be at least conceptually
understood through (1) the typical shock-jump conditions in MHD turbulence and
(2) the isotropic turbulent mixing of the flow as
$\operatorname{\mathcal{M}_{\text{A0}}}$ increases. Beattie et al. (2021)
found that the magnetic field preferentially suppresses the
$\rho/\rho_{0}$-variance at high-$\operatorname{\mathcal{M}}$ and attributed
this to the magnetic field being able to more easily smooth out small-scale
fluctuations, which are flux-frozen to the field (Landau & Lifshitz, 1959),
and are introduced into the flow as $\operatorname{\mathcal{M}}$ increases
(Kim & Ryu, 2005). One can see this immediately by considering a slab of
isothermal gas, threaded by magnetic fields. Perpendicular to the magnetic
field $B\propto\rho$ (Hennebelle & Inutsuka, 2019), and by using the regular
Rankine-Hugoniot jump conditions (Landau & Lifshitz, 1959) one can show that
the density jump for these magnetised fluctuations is

$\displaystyle\frac{\rho}{\rho_{0}}=\dfrac{1}{2}\left[\sqrt{\left(1+2\frac{\mathcal{M}^{2}_{\text{A0}}}{(b\operatorname{\mathcal{M}})^{2}}\right)^{2}+8\mathcal{M}^{2}_{\text{A0}}}-\left(1+2\frac{\mathcal{M}^{2}_{\text{A0}}}{(b\operatorname{\mathcal{M}})^{2}}\right)\right]$
(37)

(Molina et al., 2012; Mocz & Burkhart, 2019; Beattie et al., 2021). We plot
$\rho/\rho_{0}$ as a function of $\operatorname{\mathcal{M}}$ in Figure 15
(black), as well as the regular hydrodynamical shock jump conditions
$\rho/\rho_{0}=(b\operatorname{\mathcal{M}})^{2}$ (red) for comparison. The
magnetised gas density jump diverges from the hydrodynamical behaviour at
high-$\operatorname{\mathcal{M}}$ because the small-scale density fluctuations
are suppressed by the magnetic field. This significantly reduces the amplitude
of the fluctuations (both for over- and under-densities; see e.g., Figure 3 in
Mocz & Burkhart 2018), which means that $\mathrm{\left\langle
s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}}\right\rangle}_{L/\Gamma^{n}}$ is also
reduced, giving rise to a more self-similar (small $T$), Gaussian-like set of
global $s$-statistics for the high-$\operatorname{\mathcal{M}}$ simulations
(Squire & Hopkins, 2017), regardless of
$\operatorname{\mathcal{M}_{\text{A0}}}$.

Figure 15: The theoretical density jump for a hydrodynamical (black) and
magnetised (red) shocked region of gas as a function of
$\operatorname{\mathcal{M}}$ derived using the Rankine-Hugoniot jump
conditions. The magnetised gas density jump diverges from the hydrodynamical
value at high-$\operatorname{\mathcal{M}}$ as the small-scale density
fluctuations are suppressed by the magnetic field. The
high-$\operatorname{\mathcal{M}}$ limit, shown in blue, shows that stronger
fields suppress these fluctuations more efficiently.

In Figure 8 we found that the $T$ parameter depends upon
$\operatorname{\mathcal{M}_{\text{A0}}}$:
high-$\operatorname{\mathcal{M}_{\text{A0}}}$ simulations become lognormal,
$T\sim\operatorname{\mathcal{M}}^{-3}$ faster than
low-$\operatorname{\mathcal{M}_{\text{A0}}}$,
$T\sim\operatorname{\mathcal{M}}^{-2}$. We believe this could be due to the
global transition from an weakly-isotropic to strongly-anisotropic $s$-field.
For the sub-Alfvénic mean-field turbulence, the high-energy
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ provides a domain,
$\ell_{\parallel}$, in the turbulence along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, where
hydrodynamical-like shocks, as we showed in §7, increasing $\mathcal{S}_{s}$
and hence $T$. But for the trans-Alfvénic simulations,
$\operatorname{\mathcal{M}_{\text{A0}}}=1-2$,
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ and
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ are at almost energy
equipartition (Beattie et al., 2020, see Figure 7). The transition is
qualitatively illustrated in Figure 2 through the difference between the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field lines in the top row, where the
field is highly-coherent across the whole box, dominated by the
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, and bottom row,
where the field is becoming tangled as the
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ grows. Once there is significant
magnetic energy in $\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ there is no
longer a region in the turbulence, $\ell_{\parallel}$, where hydrodynamical-
like shocks can form (left panel of Figure 10), since the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-fields are tangled into random
orientations all through space. This means magnetised $s$-fluctuations are
more uniformly distributed through the fluid, and hence, there are many more
small amplitude fluctuations in the fluid than in the highly-sub-Alfvénic
flows. As we discussed above, for the high-$\operatorname{\mathcal{M}}$
simulations, having many small amplitude fluctuations leads to lower values of
$T$, and more, self-similar, $\mathrm{\left\langle
s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}}\right\rangle}_{L/\Gamma^{n}}\approx 0$,
Gaussian-like $s$-statistics.

Finally, we note that our interpretation relies upon measuring $T$ via the
$s$-PDF fits, and interpreting this parameter in the context of SS. One way to
directly measure the SS would be to compute the scale-dependent
$(s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}})$-PDF, which we do not do in this study,
but we think is feasible to do with high-resolution turbulence data that was
recently produced in Federrath et al. (2021). One can compare this with the
$T$ computed from the $s$-PDF with the direct measurement,
$\mathrm{\left\langle
s_{L/\Gamma^{n}}-s_{L/\Gamma^{n+1}}\right\rangle}_{L/\Gamma^{n}}$. This would
allow us to understand better the exact physics of $T$, and how we can use it
to interpret the $s$-PDF intermittency.

###  8.2 Timescale interpretation

The central tenant of the MB2019 $s$-PDF phenomenology is timescales that
determine the lifetimes of high-density and low-density structures in the
turbulence are different (see §3). High-density objects, $s>s_{0}$, live on
shorter timescales compared to lower-density objects, $s\leq s_{0}$, reduced
by a factor of $1+3f/2$ (see Equation 26) which amounts to skewness in the
$s$-PDF. It is true that the magnetic field can act to shield density
fluctuations, prolonging the life of the fluctuations. Hennebelle (2013) found
this was the case for $\approx 55-65\%$ of the dense, magnetised filamentary
structures in a suite of MHD box simulations. The radial Lorentz force was
found to mostly point towards the filaments, opposing thermal pressure and
preventing the fluctuation from expanding, prolonging the life of the over-
density. This could be in part responsible for the Gaussianisation of the
statistics, coupled with $\operatorname{\mathcal{M}}$ increasing, increasing
the number of shocks and getting larger and more complete samples of shock
lifetimes. The $\operatorname{\mathcal{M}_{\text{A0}}}$ dependency could be
explained in the same fashion as the previous SS phenomenology, i.e., sub-
Alfvénic flows allow for hydrodynamical-like shocks to form along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$, which live on
shorter timescales than magnetised shocks, whereas shocks in more trans-
Alfvénic, isotropic flows, always feel the magnetic field, prolonging the
lifetimes for all of the density fluctuations. This is consistent with what we
find in Figure 9, where we directly compute $1+3f/2$ for our simulations. It
shows a spread in $1+3f/2$ for different
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field strengths, where
$\operatorname{\mathcal{M}_{\text{A0}}}=2$ is about 25% slower than
$\operatorname{\mathcal{M}_{\text{A0}}}=0.1$, at
$\operatorname{\mathcal{M}}\approx 4$, before becoming completely homogenised
at $\operatorname{\mathcal{M}}\gtrsim 4$.

Robertson & Goldreich (2018) has probably done the most robust measurements of
lifetimes for over-dense regions using clustering and filtering techniques on
passive tracer particles. The lifetimes of dense shocked regions in the
supersonic isothermal turbulence simulation are short,
$\mathcal{O}(10^{-3}L/c_{s})$. In the frame of turbulent turnover times on the
driving scale that is
$\sim\mathcal{O}(10^{-2}L/[2\operatorname{\mathcal{M}}c_{s}])$ for
$\operatorname{\mathcal{M}}\gtrsim 4$. Based on our observations of the shocks
in the magnetised flows it is very reasonable to suggest that the shock
lifetimes are well beyond that of
$\sim\mathcal{O}(10^{-2}L/[2\operatorname{\mathcal{M}}c_{s}])$. For example,
we find some of the extreme shock events that we identified in Figure 12
persist for $\mathcal{O}(L/[2\operatorname{\mathcal{M}}c_{s}])$. Regardless, a
more detailed study of the lifetimes of over-dense objects, such as the one
performed for hydrodynamical simulations in Robertson & Goldreich (2018) would
be needed to test this directly, not relying on the indirect $f$ parameter to
determine how large of an affect the magnetic field has on the lifetimes of
the dense regions with different $\operatorname{\mathcal{M}}$ and
$\operatorname{\mathcal{M}_{\text{A0}}}$.

##  9 Summary and key findings

We provide a detailed analysis of the 1-point logarithmic density
($s\equiv\ln\rho/\rho_{0}$) statistics for magnetised supersonic turbulence in
the sub- to trans-Alfvénic mean-field
($\operatorname{\mathcal{M}_{\text{A0}}}\lesssim 2$) and supersonic
($\operatorname{\mathcal{M}}>1$) regime, relevant to the cool interstellar
medium. In §§1-2 we motivate and discuss lognormal models for the $s$
probability density function (PDF) and the fundamental assumptions and
problems associated with lognormal fluctuation models. Next in §3 we focus on
two models that go beyond the lognormal fluctuation picture, motivated by
relaxing the self-similarity assumption (H2013) and formulating the
$s$-fluctuations as a Langevin model with non-symmetric timescales for high-
and low-$s$ fluctuations (MB2019). In §4 we describe the numerical methods and
setup for the 16 MHD simulations we use to compute the density statistics.
Next in §5 we discuss in detail the morphology of the volume-weighted and
mass-weighted $s$-PDFs, focusing on the high- and low-$s$ tails and the
temporal fluctuations from the time-averaging. In §6 we provide a detailed
analysis of the $s$-intermittency based on four independent statistical
measures and identify, unlike hydrodynamical turbulence, that the sub- to
trans-Alfvénic mean-field turbulence becomes more lognormal with increasing
$\operatorname{\mathcal{M}}$. In §7 we use 1D pencil beams to probe the 3D
turbulence in the low-$\operatorname{\mathcal{M}}$ regime and analyse the
across and along $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$
dynamics, and show that the low-$\operatorname{\mathcal{M}}$ intermittency is
caused by strong, single shocks that form along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. We run 1D shock
simulations that support our conclusions. Finally in §8 we discuss two
different phenomenological explanations for the lack of $s$ intermittency in
the high-$\operatorname{\mathcal{M}}$ regime: self-similarity and inhomogenous
timescales between high- and low-density structures. We list the key results
below:

  * •

We show that the mass and volume-weighted statistics (Figures 4 & 7) only
become stationary in trans-to-sub-Alfvénic mean-field turbulence after about
five turbulent turnover times, $t=\ell_{0}/(c_{s}\operatorname{\mathcal{M}})$,
where $\ell_{0}$ is the driving scale, more than double the time required for
hydrodynamical turbulence to reach a stationary state.

  * •

The temporal fluctuations in the volume-weighted (Figure 5) and mass-weighted
(Figure 6) $s$-PDFs show how the low-density, high volume-filling structures,
such as rarefactions and voids, fluctuate significantly in volume, and much
less in mass. In contrast, the highest-density and lowest volume-filling
structures, such as shocked regions, filaments, and over-dense sheets,
fluctuate significantly in mass, but not in volume. This paints a picture of
magnetised turbulence where the densest regions are gaining and losing mass
over short timescales, using neighbouring voids as reservoirs, which are
expanding and contracting in volume as mass is exchanged between them.

  * •

We compare the non-lognormal models, H2013 (Equation 15-18) and MB2019
(Equation 28), and the lognormal model (Equation 1) for the volume-weighted
$s$-PDF. We find that unlike hydrodynamical turbulence, which monotonically
increases the non-Gaussian features of the $s$-PDF with
$\operatorname{\mathcal{M}}$, the H2013 or MB2019 are required to capture the
PDF morphology at $\operatorname{\mathcal{M}}\lesssim 4$, but for
$\operatorname{\mathcal{M}}\gtrsim 4$ the lognormal model is empirically (but
not theoretically) sufficient.

  * •

Motivated by the complex non-Gaussian behaviour of the $s$-PDF we calculate
four independent measures of the global, volume-weighted intermittency (1) the
skewness (Equation 32, Figure 7), (2) the mass-weighted versus volume-weighted
$s$ standard deviation relation (right panel of Figure 8), (3) the H2013 $T$
parameter (Equation 19, left panel of Figure 8) and (4) the MB2019 $f$
parameter (Equation 26). All four independent statistics confirm that
hydrodynamical-like intermittency is found for
low-$\operatorname{\mathcal{M}}$ highly-magnetised simulations, but the
intermittency monotonically reduces as a function of
$\operatorname{\mathcal{M}}$, with some dependence upon
$\operatorname{\mathcal{M}_{\text{A0}}}$, completely opposite to the
hydrodynamical turbulence studied in H2013 and Squire & Hopkins (2017).

  * •

Using 1D pencil beams (Figure 10, §7) we explore the
low-$\operatorname{\mathcal{M}}$ local parallel and perpendicular dynamics
along $\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$. Along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ large-scale
exponential shocks form through convergent flows and sawtooth shocks. We
stress, that unlike simple shock-tube experiments, the over-densities have an
internal structure, similar to over-densities described in Robertson &
Goldreich (2018). The shocks excite magnetic compression waves at the shock-
front and give rise to the largest over-densities in the flow. This also leads
to anisotropies in the magnetic tension. Across
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ we find large-
scale, rigid body vortices and very small-amplitude shear Alfvén waves
(compared to the kinetic and fast magnetosonic compression waves). We comment
that in sub-Alfvénic mean-field turbulence shear Alfvén waves can only
possibly play a dominant role in the dynamics on very small scales and the
large scales are dominated by sawtooth shocks in the gas velocity (i.e., akin
to large-scale Burgers (1948) turbulence).

  * •

Using our pencil beams we search for the low-$s$ intermittent events that are
captured by the H2013-type intermittency. We look for the most extreme under-
dense events through 49 time realisations (Figure 12), and find that the
largest contributors to the non-Gaussian $s$-statistics are strong, single
shocks ($1.5-2\sigma$ events in the velocity) that form along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ and create volume-
poor $s$ rarefactions, which are seemingly also coupled to rare velocity and
magnetosonic field fluctuation events. We show explicitly that it is single
shocks and not multi-shock interactions along
$\operatorname{\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}}$ that are the
largest contributors to the non-Gaussian components of $s$ using 1D numerical
shock experiments that we show in Figures 13,14 and 16.

  * •

For the high-$\operatorname{\mathcal{M}}$ flows, we discuss the two
phenomenologies that we use to understand why the $s$-statistics become more
Gaussian wiht increasing $\operatorname{\mathcal{M}}$ (§8). These are the
self-similarity (H2013, Squire & Hopkins 2017) and inhomogeneous timescale
(Mocz & Burkhart, 2019) interpretations of the $s$-statistics. We conclude it
is most likely that the self-similarity of the flow is maintained through
suppression of over- and under-dense $s$-fluctuations by
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$. This effect is stronger when the
flows are trans-Alfvénic and are becoming dominated by isotropic magnetised
$s$-fluctuations. We also conjecture that the magnetic field plays are role in
facilitating more uniform dynamical timescales between the low- and high-
density structures in the gas.

## Acknowledgements

J. R. B. thanks Christoph Federrath’s and Mark Krumholz’s research groups for
the many productive discussions, in particular, Shyam H. Menon, and
acknowledges financial support from the Australian National University, via
the Deakin PhD and Dean’s Higher Degree Research (theoretical physics)
Scholarships, the Research School of Astronomy and Astrophysics, via the Joan
Duffield Research Scholarship and the Australian Government via the Australian
Government Research Training Program Fee-Offset Scholarship. P. M.
acknowledges support for this work provided by NASA through Einstein
Postdoctoral Fellowship grant number PF7-180164 awarded by the Chandra X-ray
Centre, which is operated by the Smithsonian Astrophysical Observatory for
NASA under contract NAS8-03060. C. F. acknowledges funding provided by the
Australian Research Council (Future Fellowship FT180100495), and the
Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). R. S. K.
acknowledges financial support from the German Research Foundation (DFG) via
the Collaborative Research Center (SFB 881, Project-ID 138713538) ’The Milky
Way System’ (subprojects A1, B1, B2, and B8). He also thanks for funding from
the Heidelberg Cluster of Excellence STRUCTURES in the framework of Germany’s
Excellence Strategy (grant EXC-2181/1 - 390900948) and for funding from the
European Research Council via the ERC Synergy Grant ECOGAL (grant 855130). We
further acknowledge high-performance computing resources provided by the
Australian National Computational Infrastructure (grant ek9) in the framework
of the National Computational Merit Allocation Scheme and the ANU Merit
Allocation Scheme, and by the Leibniz Rechenzentrum and the Gauss Centre for
Supercomputing (grants pr32lo, pr48pi, pr74nu, pn73fi).

The simulation software, flash, was in part developed by the DOE-supported
Flash Centre for Computational Science at the University of Chicago. Data
analysis and visualisation software used in this study: C++ (Stroustrup,
2013), cython (Behnel et al., 2011), visit (Childs et al., 2012), numpy
(Oliphant, 2006; Harris et al., 2020), matplotlib (Hunter, 2007), scipy
(Virtanen et al., 2020), scikit-image (van der Walt et al., 2014), pandas
(pandas development team, 2020).

## Data Availability

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

## References

  * Abe et al. (2020) Abe D., Inoue T., Inutsuka S.-i., Matsumoto T., 2020, arXiv e-prints, p. arXiv:2012.02205
  * Barreto-Mota et al. (2021) Barreto-Mota L., de Gouveia Dal Pino E. M., Burkhart B., Melioli C., Santos-Lima R., Kadowaki L. H. S., 2021, MNRAS, 503, 5425
  * Beattie & Federrath (2020) Beattie J. R., Federrath C., 2020, MNRAS, 492, 668
  * Beattie et al. (2019) Beattie J. R., Federrath C., Klessen R. S., Schneider N., 2019, MNRAS, 488, 2493
  * Beattie et al. (2020) Beattie J. R., Federrath C., Seta A., 2020, MNRAS, 498, 1593
  * Beattie et al. (2021) Beattie J. R., Mocz P., Federrath C., Klessen R. S., 2021, MNRAS, 504, 4354
  * Behnel et al. (2011) Behnel S., Bradshaw R., Citro C., Dalcin L., Seljebotn D. S., Smith K., 2011, Computing in Science & Engineering, 13, 31
  * Ben Hcine & Bouallegue (2015) Ben Hcine M., Bouallegue R., 2015, arXiv e-prints, p. arXiv:1502.03619
  * Bhattacharjee et al. (1998) Bhattacharjee A., Ng C. S., Spangler S. R., 1998, ApJ, 494, 409
  * Bonne et al. (2020) Bonne L., et al., 2020, A&A, 641, A17
  * Bouchut et al. (2010) Bouchut F., Klingenberg C., Waagan K., 2010, Numerische Mathematik, 115, 647
  * Brunt et al. (2010a) Brunt C. M., Federrath C., Price D. J., 2010a, MNRAS, 403, 1507
  * Brunt et al. (2010b) Brunt C. M., Federrath C., Price D. J., 2010b, MNRAS, 405, L56
  * Burgers (1948) Burgers J., 1948, Advances in Applied Mechanics, 1, 171
  * Burkhart (2018) Burkhart B., 2018, ApJ, 863, 118
  * Burkhart & Lazarian (2012) Burkhart B., Lazarian A., 2012, ApJ, 755, L19
  * Burkhart et al. (2009) Burkhart B., Falceta-Gonçalves D., Kowal G., Lazarian A., 2009, ApJ, 693, 250
  * Castaing (1996) Castaing B., 1996, J. Phys. II France, 6, 105
  * Chen & Ostriker (2014) Chen C.-Y., Ostriker E. C., 2014, ApJ, 785, 69
  * Chen et al. (2020) Chen M. C.-Y., et al., 2020, The Astrophysical Journal, 891, 84
  * Childs et al. (2012) Childs H., et al., 2012, in , High Performance Visualization–Enabling Extreme-Scale Scientific Insight. Taylor & Francis, pp 357–372
  * 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
  * Falgarone et al. (1995) Falgarone E., Pineau des Forets G., Roueff E., 1995, A&A, 300, 870
  * Falgarone et al. (2009) Falgarone E., Pety J., Hily-Blant P., 2009, A&A, 507, 355
  * Falgarone et al. (2015) Falgarone E., Momferratos G., Lesaffre P., 2015, The Intermittency of ISM Turbulence: What Do the Observations Tell Us?. Springer Berlin Heidelberg, Berlin, Heidelberg, pp 227–252, doi:10.1007/978-3-662-44625-6_9, https://doi.org/10.1007/978-3-662-44625-6_9
  * 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 & 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. (2008a) Federrath C., Glover S. C. O., Klessen R. S., Schmidt W., 2008a, Physica Scripta Volume T, 132, 014025
  * Federrath et al. (2008b) Federrath C., Klessen R. S., Schmidt W., 2008b, 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. (2016) Federrath C., et al., 2016, ApJ, 832, 143
  * Federrath et al. (2021) Federrath C., Klessen R. S., Iapichino L., Beattie J. R., 2021, Nature Astronomy
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Girichidis (2021) Girichidis P., 2021, arXiv e-prints, p. arXiv:2106.12596
  * Girichidis et al. (2014) Girichidis P., Konstandin L., Whitworth A. P., Klessen R. S., 2014, ApJ, 781, 91
  * Goldreich & Sridhar (1995) Goldreich P., Sridhar S., 1995, ApJ, 438, 763
  * Grete et al. (2020) Grete P., O’Shea B. W., Beckwith K., 2020, ApJ, 889, 19
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hennebelle (2013) Hennebelle P., 2013, A&A, 556, A153
  * 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 et al. (2020) Heyer M., Soler J. D., Burkhart B., 2020, MNRAS, 496, 4546
  * Hily-Blant & Falgarone (2009) Hily-Blant P., Falgarone E., 2009, A&A, 500, L29
  * Ho & Lazarian (2021) Ho K. W., Lazarian A., 2021, ApJ, 911, 53
  * Hoang et al. (2021) Hoang T. D., et al., 2021, arXiv e-prints, p. arXiv:2108.10045
  * Hopkins (2013) Hopkins P. F., 2013, MNRAS, 430, 1880
  * Howes & Bourouaine (2017) Howes G. G., Bourouaine S., 2017, Journal of Plasma Physics, 83, 905830408
  * 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
  * Hwang et al. (2021) Hwang J., et al., 2021, ApJ, 913, 85
  * Imara et al. (2021) Imara N., Forbes J. C., Weaver J. C., 2021, arXiv e-prints, p. arXiv:2108.00014
  * Khullar et al. (2021) Khullar S., Federrath C., Krumholz M. R., Matzner C. D., 2021, arXiv e-prints, p. arXiv:2107.00725
  * Kim & Ryu (2005) Kim J., Ryu D., 2005, ApJ, 630, L45
  * Kitsionas et al. (2009) Kitsionas S., et al., 2009, A&A, 508, 541
  * Klessen (2000) Klessen R. S., 2000, ApJ, 535, 869
  * Kolmogorov (1941) Kolmogorov A. N., 1941, Doklady Akademii Nauk Sssr, 30, 301
  * Konstandin et al. (2012a) Konstandin L., Federrath C., Klessen R. S., Schmidt W., 2012a, Journal of Fluid Mechanics, 692, 183
  * Konstandin et al. (2012b) Konstandin L., Girichidis P., Federrath C., Klessen R. S., 2012b, ApJ, 761, 149
  * Kowal et al. (2007) Kowal G., Lazarian A., Beresnyak A., 2007, ApJ, 658, 423
  * Kritsuk et al. (2007) Kritsuk A. G., Norman M. L., Padoan P., Wagner R., 2007, ApJ, 665, 416
  * Kritsuk et al. (2011) Kritsuk A. G., et al., 2011, ApJ, 737, 13
  * Krumholz & McKee (2005) Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250
  * Landau & Lifshitz (1959) Landau L., Lifshitz E., 1959, Fluid Mechanics: Landau and Lifshitz: Course of Theoretical Physics. Butterworth-Heinemann
  * Lehmann et al. (2016) Lehmann A., Federrath C., Wardle M., 2016, MNRAS, 463, 1026
  * Li (2021) Li H.-B., 2021, Galaxies, 9
  * Li et al. (2003) Li Y., Klessen R. S., Mac Low M.-M., 2003, ApJ, 592, 975
  * Li et al. (2013) Li H.-b., Fang M., Henning T., Kainulainen J., 2013, Monthly Notices of the Royal Astronomical Society, 436, 3707
  * Lim et al. (2020) Lim J., Cho J., Yoon H., 2020, ApJ, 893, 75
  * Lo (2012) Lo C. F., 2012, Journal of Applied Mathematics, 2012, 838397
  * Makwana & Yan (2020) Makwana K. D., Yan H., 2020, Phys. Rev. X, 10, 031021
  * Menon et al. (2020) Menon S. H., Federrath C., Kuiper R., 2020, MNRAS, 493, 4643
  * Menon et al. (2021) Menon S. H., Federrath C., Klaassen P., Kuiper R., Reiter M., 2021, MNRAS, 500, 1721
  * Miyoshi & Kusano (2005) Miyoshi T., Kusano K., 2005, Journal of Computational Physics, 208, 315
  * Mocz & Burkhart (2018) Mocz P., Burkhart B., 2018, MNRAS, 480, 3916
  * Mocz & Burkhart (2019) Mocz P., Burkhart B., 2019, ApJ, 884, L35
  * Mocz et al. (2017) Mocz P., Burkhart B., Hernquist L., McKee C. F., Springel V., 2017, ApJ, 838, 40
  * Mohapatra et al. (2020) Mohapatra R., Federrath C., Sharma P., 2020, MNRAS, 493, 5838
  * Mohapatra et al. (2021) Mohapatra R., Federrath C., Sharma P., 2021, MNRAS, 500, 5072
  * Molina et al. (2012) Molina F. Z., Glover S. C. O., Federrath C., Klessen R. S., 2012, MNRAS, 423, 2680
  * Nolan et al. (2015) Nolan C. A., Federrath C., Sutherland R. S., 2015, MNRAS, 451, 1380
  * Nordlund & Padoan (1999) Nordlund Å. K., Padoan P., 1999, in Franco J., Carraminana A., eds, Interstellar Turbulence. p. 218 (arXiv:astro-ph/9810074)
  * 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
  * Padoan & Nordlund (1999) Padoan P., Nordlund Å., 1999, ApJ, 526, 279
  * 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
  * Pan et al. (2019) Pan L., Padoan P., Nordlund Å., 2019, ApJ, 881, 155
  * Park & Ryu (2019a) Park J., Ryu D., 2019a, ApJ, 875, 2
  * Park & Ryu (2019b) Park J., Ryu D., 2019b, ApJ, 875, 2
  * Passot & Vázquez-Semadeni (1998) Passot T., Vázquez-Semadeni E., 1998, Phys. Rev. E, 58, 4501
  * Pillai et al. (2020) Pillai T. G. S., et al., 2020, Nature Astronomy
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 586, A137
  * Price & Federrath (2010) Price D. J., Federrath C., 2010, MNRAS, 406, 1659
  * Price et al. (2011) Price D. J., Federrath C., Brunt C. M., 2011, ApJ, 727, 1380
  * Robertson & Goldreich (2018) Robertson B., Goldreich P., 2018, ApJ, 854, 88
  * Scannapieco & Safarzadeh (2018) Scannapieco E., Safarzadeh M., 2018, ApJ, 865, L14
  * Schmidt et al. (2008) Schmidt W., Federrath C., Klessen R., 2008, Phys. Rev. Lett., 101, 194505
  * Schneider et al. (2013) Schneider N., et al., 2013, ApJ, 766, L17
  * Seta et al. (2020) Seta A., Bushby P. J., Shukurov A., Wood T. S., 2020, arXiv e-prints, p. arXiv:2003.07997
  * She & Leveque (1994) She Z.-S., Leveque E., 1994, Phys. Rev. Lett., 72, 336
  * Skalidis & Tassis (2020) Skalidis R., Tassis K., 2020, arXiv e-prints, p. arXiv:2010.15141
  * Soler et al. (2017) Soler J. D., et al., 2017, A&A, 603, A64
  * Squire & Hopkins (2017) Squire J., Hopkins P. F., 2017, MNRAS, 471, 3753
  * Stroustrup (2013) Stroustrup B., 2013, The C++ Programming Language, 4th edn. Addison-Wesley Professional
  * Tritsis et al. (2018) Tritsis A., Federrath C., Schneider N., Tassis K., 2018, MNRAS, 481, 5275
  * Vazquez-Semadeni (1994) Vazquez-Semadeni E., 1994, ApJ, 423, 681
  * 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
  * Wu et al. (2005) Wu J., Mehta N., Zhang J., 2005, in GLOBECOM ’05. IEEE Global Telecommunications Conference, 2005.. pp 3413–3417, doi:10.1109/GLOCOM.2005.1578407
  * Xu et al. (2019) Xu S., Ji S., Lazarian A., 2019, ApJ, 878, 157
  * Yoffe & McComb (2021) Yoffe S. R., McComb W. D., 2021, arXiv e-prints, p. arXiv:2107.09112
  * Yuen & Lazarian (2020) Yuen K. H., Lazarian A., 2020, ApJ, 898, 66
  * pandas development team (2020) pandas development team T., 2020, pandas-dev/pandas: Pandas, doi:10.5281/zenodo.3509134, https://doi.org/10.5281/zenodo.3509134
  * van der Walt et al. (2014) van der Walt S., et al., 2014, PeerJ, 2, e453

##  Appendix A 1D shock experiments

Figure 16: The time-evolution of logarithmic density profiles from the four 1D
shock experiments, (1) single shock, (2) two shocks, (3) two shocks, constant
driving, (4) many shocks, constant driving, that we use to explore the
low-$\operatorname{\mathcal{M}}$ intermittency. The profiles are coloured such
that black corresponds to $t/t_{\rm cross}=0$, and pink $t/t_{\rm cross}\gg
0$.

We perform 1D shock experiments with a similar setup to the one described in
Mocz & Burkhart (2018). We solve the compressible Euler equations in 1D,

$\displaystyle\frac{\partial\rho}{\partial t}+\frac{\partial}{\partial x}\rho
v_{x}$ $\displaystyle=0,$ (38)
$\displaystyle\rho\left(\frac{\partial}{\partial
t}+v_{x}\frac{\partial}{\partial x}\right)v_{x}$
$\displaystyle=c_{s}^{2}\frac{\partial}{\partial x}\rho,$ (39)

based on the multi-state Harten–Lax–van Leer approximate Riemann method
described in Miyoshi & Kusano (2005). We solve the equations on a periodic