AppendixAppendix References

# Filaments and striations: anisotropies in observed, supersonic, highly-
magnetised turbulent clouds

James R. Beattie1 and Christoph Federrath1
1Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Stars form in highly-magnetised, supersonic turbulent molecular clouds. Many
of the tools and models that we use to carry out star formation studies rely
upon the assumption of cloud isotropy. However, structures like high-density
filaments in the presence of magnetic fields, and magnetosonic striations
introduce anisotropies into the cloud. In this study we use the two-
dimensional (2D) power spectrum to perform a systematic analysis of the
anisotropies in the column density for a range of Alfvén Mach numbers
($\operatorname{\mathcal{M}_{\text{A}}}=0.1$–$10$) and turbulent Mach numbers
($\operatorname{\mathcal{M}}=2$–$20$), with 20 high-resolution, three-
dimensional (3D) turbulent magnetohydrodynamic simulations. We find that for
cases with a strong magnetic guide field, corresponding to
$\operatorname{\mathcal{M}_{\text{A}}}<1$, and
$\operatorname{\mathcal{M}}\lesssim 4$, the anisotropy in the column density
is dominated by thin striations aligned with the magnetic field, while for
$\operatorname{\mathcal{M}}\gtrsim 4$ the anisotropy is significantly changed
by high-density filaments that form perpendicular to the magnetic guide field.
Indeed, the strength of the magnetic field controls the degree of anisotropy
and whether or not any anisotropy is present, but it is the turbulent motions
controlled by $\operatorname{\mathcal{M}}$ that determine which kind of
anisotropy dominates the morphology of a cloud.

###### keywords:

turbulence – magnetohydrodynamics – ISM: clouds – ISM: kinematics and dynamics
– ISM: magnetic fields

††pubyear: 2019††pagerange: Filaments and striations: anisotropies in
observed, supersonic, highly-magnetised turbulent clouds–11

##  1 Introduction

Stars are born inside high-density regions of molecular H2 clouds (MCs) that
have undergone collapse. But cloud collapse is not a straightforward process.
MCs are anisotropic, fragmented and undergo supersonic, turbulent motions with
dynamically relevant magnetic fields (Elmegreen & Falgarone, 1996; Elmegreen &
Scalo, 2004; Mac Low & Klessen, 2004; Lazarian & Cho, 2004; McKee & Ostriker,
2007; Crutcher, 2012; Federrath & Klessen, 2012; Padoan et al., 2014;
Hennebelle & Inutsuka, 2019). The supersonic turbulent motions play a dual
role in forming stars by both mixing the H2 gas, which increases the critical
density required for collapse and by creating high-density structures where
the cloud can cool and form dense protostellar cores (Krumholz & McKee, 2005;
Federrath & Klessen, 2012; Federrath & Banerjee, 2015). The role of magnetic
fields, however, is less understood. Certainly the presence of magnetic fields
suppresses the density fluctuations caused by the turbulent motions (e.g.
Molina et al. 2012) and reduces the star-formation rate by a factor of a few
(Padoan & Nordlund, 2011; Federrath & Klessen, 2012; Federrath, 2015; Krumholz
& Federrath, 2019). However, magnetic fields may also be responsible for
overcoming resistance to star formation, for example, by magnetic
reconnection, ambipolar diffusion or magnetic breaking. The complex interplay
of these physical processes is still not fully understood (Mouschovias et al.,
2006; Lazarian et al., 2012; Hennebelle & Inutsuka, 2019). Magnetic fields
also may play a role in feedback, for example, by driving protostellar jets
that inject mass, momentum and energy into the medium, possibly stirring
turbulent motions (Frank et al., 2014; Gerrard et al., 2019; Krumholz &
Federrath, 2019; Kuruwita & Federrath, 2019). They also seem to play a role in
the orientation of filamentary structures and striations in the MCs which may
not directly influence the star-forming potential of the cloud, but certainly
influences the global structure of the cloud through the density PDF (Planck
Collaboration et al., 2016a, b; Cox et al., 2016; Malinen et al., 2016; Soler
et al., 2017; Tritsis & Tassis, 2016; Tritsis et al., 2018; Soler, 2019).
High-density filamentary structures house the local conditions necessary to
overcome the turbulent motions, and observations find protostars located at
junctions (or hubs) between intersecting filaments (André et al., 2010;
Men’shchikov et al., 2010; Schneider et al., 2013; Arzoumanian et al., 2018;
Tokuda et al., 2018; Trevino-Morales et al., 2019; Roy et al., 2019; Xu et
al., 2019). Indeed, high-density filaments may even play a role in setting the
stellar initial mass function (André et al., 2019). Striations, in contrast to
high-density filaments, most likely form through fast magnetosonic waves,
largely independent of the turbulent properties of the cloud (Tritsis &
Tassis, 2016; Tritsis et al., 2018), consisting of thin channels of compressed
gas that run coincident with the magnetic field direction (Cox et al., 2016;
Malinen et al., 2016; Soler et al., 2017). In this study, we explore how both
striations and high-density filaments form and contribute to the anisotropy in
the column density.

###  1.1 Anisotropy in MHD Turbulence

Magnetic fields locally111locally in this context means at length scales where
the magnetic tension,
$([\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{B}}}}]\mathrm{{\mn@boldsymbol{\mathit{B}}}})/(4\pi)$,
can be neglected, and that are comparable to turbulent eddy scales (Burkhart
et al. 2015 and references therein). impart anisotropic behaviour on the
velocity and density structures that are present in the turbulent clouds
(Goldreich & Sridhar, 1995; Cho & Vishniac, 2000; Cho & Lazarian, 2003; Kowal
et al., 2007; Burkhart et al., 2014). This is because the magnetic field does
not disappear at any length scale and flows along the field lines do not feel
the pressures and tension from the Lorentz force. What this means is that,
locally, length scales perpendicular and parallel to the magnetic field lines
become dynamically dissimilar in the cloud. Pioneering work from Goldreich &
Sridhar (1995) showed that the anisotropy in the energy spectrum for an
incompressible, turbulent magnetohydrodynamical (MHD) fluid with equal energy
contributions from the magnetic and turbulent forces is,

$\frac{\operatorname{k_{\parallel}}}{\operatorname{k_{\perp}}}\sim\left(L\operatorname{k_{\perp}}\right)^{-1/3},$
(1)

where $\operatorname{k_{\parallel}}$ and $\operatorname{k_{\perp}}$ are
$k$-space scales parallel and perpendicular, respectively, to the mean, local
magnetic field, and $L$ is the cloud scale. The anisotropy is introduced from
a critical balance between the turnover time of turbulent eddies perpendicular
to the magnetic field, $T_{\perp}\sim\ell_{\perp}/V_{\ell}$, where
$V_{\ell}\sim\ell_{\perp}^{1/3}$ is the velocity scaling for eddies
perpendicular to the field-lines, and the period of the Alfvén waves
travelling along the field-lines,
$T_{\parallel}\sim\ell_{\parallel}/V_{\text{A}}$, where
$V_{\text{A}}\sim|\mathrm{{\mn@boldsymbol{\mathit{B}}}}|/\sqrt{\rho}$ is the
Alfvén velocity, where $|\mathrm{{\mn@boldsymbol{\mathit{B}}}}|$ is the
strength of the magnetic field and $\rho$ is the density. Putting this
together,

$\ell_{\parallel}/V_{\text{A}}\sim\ell_{\perp}/V_{\ell}\sim\ell_{\perp}/\ell_{\perp}^{1/3}\implies\ell_{\parallel}\sim\ell_{\perp}^{2/3}.$
(2)

Furthermore, this means that the anisotropy in the energy spectrum,
$\operatorname{k_{\parallel}}\sim\operatorname{k_{\perp}}^{2/3}$, is a
function of scale. The anisotropy grows slowly with
$\operatorname{k_{\perp}}$, i.e. for small length scales in the turbulence we
expect stronger anisotropies. This results in the turbulent eddies becoming
stretched along the magnetic field and small-scale eddies are stretched more
than large-scale eddies. The Goldreich & Sridhar (1995) model has been
numerically validated in the incompressible regime (Cho & Vishniac, 2000;
Maron & Goldreich, 2001; Cho et al., 2002) and extended for compressible
turbulence in the sub-Alfvénic regime
$(\operatorname{\mathcal{M}_{\text{A}}}<1$, where
$\operatorname{\mathcal{M}_{\text{A}}}=(c_{s}\operatorname{\mathcal{M}})/V_{A}$
is the Alfvén Mach number, $\operatorname{\mathcal{M}}=V_{L}/c_{s}$ is the
sonic Mach number, where $c_{s}$ is the sound speed in the cloud and $V_{L}$
is the velocity at cloud scale $L$), and in the super-Alfvénic regime
$(\operatorname{\mathcal{M}_{\text{A}}}>1$; Lazarian & Vishniac 1999; Cho &
Lazarian 2003). However, the scaling laws predicted in the Goldreich & Sridhar
(1995) model are still debated (see discussion in Perez & Boldyrev 2007, and
spectral analysis in Boldyrev 2006, for example) and will only be resolved
when high-resolution simulations can separate the relevant scaling ranges
without ambiguity.

The anisotropic behaviour of MHD turbulence is of great interest (e.g. Bigot
et al. 2008; Esquivel & Lazarian 2011; Burkhart et al. 2014; Verdini et al.
2015; Tritsis et al. 2018; Hennebelle & Inutsuka 2019; Xu et al. 2019). It is
not only one of the key differences between MHD turbulence and the Kolmogorov
(1941) hydrodynamical (HD) turbulence, and leads to many interesting
astrophysical phenomena, but also because the ansiotropy can be used as a tool
for measuring cloud kinematics and may end up playing an important role in the
creation of filaments (Hennebelle & Inutsuka, 2019; Xu et al., 2019). For
example, properties of anisotropic density structures, striations, have been
used to measure magnetic field strength in the cloud (Tritsis & Tassis, 2016;
Tritsis et al., 2018), the velocity anisotropies (centroids in PPV space) have
been used to determine the Alfvénic regime that the turbulence is in, first by
Esquivel & Lazarian (2011) and then extended by Burkhart et al. (2014). Xu et
al. (2019) provides a model for low-density filament orientation caused by
anisotropies in the density field of subsonic turbulence. In our study we
further explore the anisotropic density structures induced by strong magnetic
fields and strong turbulence, systematically across a large parameter range by
comparing the three-dimensional (3D) simulation data with line-of-sight
projections.

###  1.2 Line-of-sight Projections

Line-of-sight (LOS) position-position (PP) projection data, i.e. the column
density, $\Sigma$, can be obtained through molecular lines, dust emission, or
dust extinction observations of the MCs (Schneider et al., 2015). The
intrinsically 3D density structures appear significantly altered in LOS
projections. For example both high and low densities are truncated in the LOS
projection (i.e. the variance of the density distribution decreases; Brunt et
al. 2010b; Federrath & Klessen 2013) and density variations on large scales
are smoothed less than on small scales (Beattie et al., 2019a). Star formation
models, and in general models of turbulence are however constructed and tested
upon the 3D density and velocity data. Hence reconstructing 3D kinematic or
dynamical information from 2D projections is a highly-nontrivial but important
task (Larson, 1981; Federrath et al., 2010; Brunt et al., 2010b, a; Beaumont
et al., 2013; Ginsburg et al., 2013; Brunt & Federrath, 2014; Kainulainen et
al., 2014; Tritsis & Tassis, 2018; Beattie et al., 2019a). There are a number
of analytic models using relations in $k$-space to extract meaningful PPP
cloud observables, such as the density dispersion and type of turbulent
driving (Brunt et al. 2010a, b; Brunt & Federrath 2014; who assumed cloud
isotropy) empirical methods used to determine the PPP fractal dimension and
turbulent Mach number of the clouds from the statistics of the column density
(Sanchez et al., 2005; Beattie et al., 2019a, b), and correlations in the
(bi)spectrum between PP and PPP data (Burkhart et al., 2009).

In this study we explore the anisotropic behaviour using the 2D power spectrum
of the column density in supersonic, sub-Alfvénic and super-Alfvénic MHD
turbulent regimes with many time-realisations in 20, high-resolution, 3D MHD
simulations. Our study is organised into the following sections. In §2 we
discuss the 20 supersonic, turbulent MHD cloud models that we use. In §3 we
briefly discuss the structures that are present in the column density maps for
each of the simulations. Next, in §4 we define the 2D power spectrum that we
use to explore the anisotropic structures in each of the column density maps
and show how the different regimes give rise to different types of anisotropic
behaviours. In §5 we calculate the $\operatorname{k_{\perp}}$ and
$\operatorname{k_{\parallel}}$ cascades from the 2D power spectra. In §6 we
quantify the average anisotropy in the power spectra and reveal the
$\operatorname{\mathcal{M}}$ and $\operatorname{\mathcal{M}_{\text{A}}}$
dependency of the anisotropy. Finally, in §7 we summarise the key results of
the study.

Table 1: Main simulation parameters and derived quantities.

|  | Mean $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field Components | rms Quantities | Derived Quantities
---|---|---|---|---
Simulation | Grid | $\operatorname{\mathcal{M}_{\text{A}}}_{0}$ |  $|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$| | $\operatorname{\mathcal{M}}$ | $\operatorname{\mathcal{M}_{\text{A}}}$ | $|\mathrm{{\mn@boldsymbol{\mathit{B}}}}|\,\text{[$\mu$G]}$ | $|\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}|\,\text{[$\mu$G]}$ | $\mathrm{\left\langle a\right\rangle}_{\operatorname{k_{\perp}}}$
ID | Res. | $(\pm 1\sigma)$ | [$\mu$G] | $(\pm 1\sigma)$ | $(\pm 1\sigma)$ | $(\pm 1\sigma)$ | $(\pm 1\sigma)$ | $(\pm 1\sigma)$
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9)
M2Ma10 | $512^{3}$ | 9.0 $\pm$ 0.8 | 0.71 | 1.80 |  $\pm$ 0.08 | 2.4 |  $\pm$ 0.2 | 2.6 |  $\pm$ 0.1 | 2.0 |  $\pm$ 0.1 | 1.03 |  $\pm$ 0.11
M2Ma2 | $512^{3}$ | 1.7 $\pm$ 0.1 | 3.54 | 1.66 |  $\pm$ 0.05 | 1.32 |  $\pm$ 0.07 | 4.5 |  $\pm$ 0.1 | 1.0 |  $\pm$ 0.1 | 0.88 |  $\pm$ 0.02
M2Ma1 | $512^{3}$ | 0.98 $\pm$ 0.07 | 7.09 | 2.0 |  $\pm$ 0.1 | 0.94 |  $\pm$ 0.07 | 7.36 |  $\pm$ 0.02 | 0.27 |  $\pm$ 0.02 | 0.61 |  $\pm$ 0.06
M2Ma0.5 | $512^{3}$ | 0.54 $\pm$ 0.04 | 14.18 | 2.2 |  $\pm$ 0.2 | 0.54 |  $\pm$ 0.04 | 14.26 |  $\pm$ 0.01 | 0.08 |  $\pm$ 0.01 | 0.5 |  $\pm$ 0.1
M2Ma0.1 | $512^{3}$ | 0.133 $\pm$ 0.008 | 70.90 | 2.6 |  $\pm$ 0.2 | 0.131 |  $\pm$ 0.01 | 70.902 |  $\pm$ 0.001 | 0.002 |  $\pm$ 0.001 | 0.57 |  $\pm$ 0.07
M4Ma10 | $512^{3}$ | 9.2 $\pm$ 0.6 | 1.42 | 3.7 |  $\pm$ 0.1 | 2.8 |  $\pm$ 0.2 | 4.6 |  $\pm$ 0.1 | 3.2 |  $\pm$ 0.1 | 1.00 |  $\pm$ 0.03
M4Ma2 | $512^{3}$ | 1.73 $\pm$ 0.07 | 7.09 | 3.5 |  $\pm$ 0.1 | 1.43 |  $\pm$ 0.06 | 8.6 |  $\pm$ 0.1 | 1.5 |  $\pm$ 0.1 | 0.94 |  $\pm$ 0.02
M4Ma1 | $512^{3}$ | 0.95 $\pm$ 0.08 | 14.18 | 3.8 |  $\pm$ 0.3 | 0.98 |  $\pm$ 0.07 | 14.61 |  $\pm$ 0.06 | 0.43 |  $\pm$ 0.06 | 0.80 |  $\pm$ 0.08
M4Ma0.5 | $512^{3}$ | 0.54 $\pm$ 0.03 | 28.36 | 4.4 |  $\pm$ 0.2 | 0.54 |  $\pm$ 0.03 | 28.49 |  $\pm$ 0.03 | 0.13 |  $\pm$ 0.03 | 0.7 |  $\pm$ 0.1
M4Ma0.1 | $512^{3}$ | 0.13 $\pm$ 0.01 | 141.80 | 5.2 |  $\pm$ 0.4 | 0.13 |  $\pm$ 0.01 | 141.800 |  $\pm$ 0.002 | 0.000 |  $\pm$ 0.002 | 0.8 |  $\pm$ 0.1
M10Ma10 | $512^{3}$ | 9.2 $\pm$ 0.7 | 3.54 | 9.2 |  $\pm$ 0.4 | 3.1 |  $\pm$ 0.2 | 10.6 |  $\pm$ 0.4 | 7.1 |  $\pm$ 0.4 | 1.05 |  $\pm$ 0.01
M10Ma2 | $512^{3}$ | 1.8 $\pm$ 0.1 | 17.72 | 9.0 |  $\pm$ 0.4 | 1.5 |  $\pm$ 0.1 | 21.2 |  $\pm$ 0.4 | 3.5 |  $\pm$ 0.4 | 1.03 |  $\pm$ 0.03
M10Ma1 | $512^{3}$ | 0.93 $\pm$ 0.05 | 35.45 | 9.3 |  $\pm$ 0.5 | 0.91 |  $\pm$ 0.04 | 36.3 |  $\pm$ 0.1 | 0.9 |  $\pm$ 0.1 | 1.13 |  $\pm$ 0.09
M10Ma0.5 | $512^{3}$ | 0.52 $\pm$ 0.02 | 70.90 | 10.5 |  $\pm$ 0.4 | 0.52 |  $\pm$ 0.02 | 71.13 |  $\pm$ 0.06 | 0.23 |  $\pm$ 0.06 | 1.1 |  $\pm$ 0.2
M10Ma0.1 | $512^{3}$ | 0.125 $\pm$ 0.006 | 354.49 | 12 |  $\pm$ 1 | 0.126 |  $\pm$ 0.006 | 354.500 |  $\pm$ 0.002 | 0.010 |  $\pm$ 0.002 | 1.1 |  $\pm$ 0.3
M20Ma10 | $512^{3}$ | 9.3 $\pm$ 0.8 | 7.09 | 19 |  $\pm$ 1 | 3.4 |  $\pm$ 0.3 | 19.7 |  $\pm$ 0.7 | 12.6 |  $\pm$ 0.7 | 1.03 |  $\pm$ 0.01
M20Ma2 | $512^{3}$ | 1.8 $\pm$ 0.1 | 35.45 | 18 |  $\pm$ 1 | 1.58 |  $\pm$ 0.09 | 41.2 |  $\pm$ 0.4 | 5.8 |  $\pm$ 0.4 | 1.08 |  $\pm$ 0.03
M20Ma1 | $512^{3}$ | 0.93 $\pm$ 0.03 | 70.90 | 19 |  $\pm$ 1 | 0.92 |  $\pm$ 0.03 | 72.5 |  $\pm$ 0.3 | 1.6 |  $\pm$ 0.3 | 1.2 |  $\pm$ 0.12
M20Ma0.5 | $512^{3}$ | 0.53 $\pm$ 0.02 | 141.80 | 21 |  $\pm$ 1 | 0.52 |  $\pm$ 0.02 | 142.2 |  $\pm$ 0.1 | 0.4 |  $\pm$ 0.1 | 1.2 |  $\pm$ 0.2
M20Ma0.1 | $512^{3}$ | 0.119 $\pm$ 0.003 | 708.98 | 24 |  $\pm$ 1 | 0.119 |  $\pm$ 0.003 | 708.993 |  $\pm$ 0.003 | 0.013 |  $\pm$ 0.003 | 1.3 |  $\pm$ 0.3

  * •

Notes: For each simulation we extract 51 realisations at 0.1$\,T$ intervals,
where $T$ is the turbulent turnover time, between $5\,T$ and $10\,T$. All
$1\sigma$ fluctuations listed are from the time-averaging over $5\,T$. Column
(1): the simulation ID. Column (2): the native resolution of the 3D simulation
grid. Column (3): the Alfvén Mach number for the
mean-$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ component,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$,
$\operatorname{\mathcal{M}_{\text{A}}}_{0}=(2c_{s}\operatorname{\mathcal{M}}\sqrt{\pi\rho_{0}})/|\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}|$,
where $\rho_{0}$ is the mean density, $c_{s}$ is the sound speed and
$\operatorname{\mathcal{M}}$ is the turbulent Mach number. Column (4): the
strength of the mean magnetic field in direction $\hat{z}$,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}=2c_{s}\sqrt{\pi\rho_{0}}(\operatorname{\mathcal{M}}/\operatorname{\mathcal{M}_{\text{A}}}_{0})\hat{\mathrm{{\mn@boldsymbol{\mathit{z}}}}}$.
Column (5): the rms turbulent Mach number. Column (6): the rms Alfvén Mach
number,
$\operatorname{\mathcal{M}_{\text{A}}}=(2c_{s}\operatorname{\mathcal{M}}\sqrt{\pi\rho_{0}})/|\mathrm{{\mn@boldsymbol{\mathit{B}}}}|$.
Column (7): the rms magnetic field strength. Column (8): the turbulent
magnetic field component $\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ from
the relation
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}=\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$.
Column (9): the average anisotropy over all $\operatorname{k_{\perp}}$length
scales in the zero-mean transformed column density, $\Sigma/\Sigma_{0}-1$,
where $\Sigma_{0}$ is the mean column density. $\operatorname{k_{\perp}}$
corresponds to wave vectors perpendicular to the magnetic field direction.

##  2 Magnetised, Turbulent Molecular Cloud Models

###  2.1 MHD Model

In this study we analyse the column densities of 20 high-resolution, 3D
turbulent, ideal magnetohydrodynamical (MHD) simulations of quiescent (non-
star-forming), supersonic molecular clouds, with no self-gravitation and
isothermal equation of state (EOS). A comprehensive parameter set, including
mean-field, root-mean-squared (RMS) and derived components for each of the
molecular cloud models is listed in Table 1. We use a modified version of
flash based on version 4.0.1 (Fryxell et al., 2000; Dubey et al., 2008) to
solve the compressible, ideal MHD equations,

$\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mathrm{{\mn@boldsymbol{\mathit{v}}}})$ $\displaystyle=0,$
(3) $\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
P_{*}+\rho\mathrm{{\mn@boldsymbol{\mathit{F}}}},$ (4)
$\displaystyle\frac{\partial\mathrm{{\mn@boldsymbol{\mathit{B}}}}}{\partial
t}$
$\displaystyle=\nabla\times(\mathrm{{\mn@boldsymbol{\mathit{v}}}}\times\mathrm{{\mn@boldsymbol{\mathit{B}}}}),$
(5) $\displaystyle\nabla\cdot\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
$\displaystyle=0,$ (6)

where
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}=\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}+\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$
is the magnetic field, made from mean-field,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$, and fluctuating,
$\delta\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, components, $\rho$ is the
density field, $\mathrm{{\mn@boldsymbol{\mathit{v}}}}$ is the velocity field
and $P_{*}=c_{s}^{2}\rho+|\mathrm{{\mn@boldsymbol{\mathit{B}}}}|^{2}/(8\pi)$
is the pressure, which is the sum of the thermal pressure $c_{s}^{2}\rho$,
where $c_{s}$ is the sound speed (i.e. the isothermal EOS), which is
normalised to $c_{s}=1$ in our simulations so all velocities are in units of
$\operatorname{\mathcal{M}}$, and magnetic pressure,
$|\mathrm{{\mn@boldsymbol{\mathit{B}}}}|^{2}/(8\pi)$. The reader should note
that since we have an isothermal EOS we need not include the energy equation
in our MHD model, since the isothermal system is completely closed, even in
the absence of the energy equation. We solve the MHD equations in a 3D box
with dimensions $L\times L\times L$, on a uniform grid with resolution
$512^{3}$, and with periodic boundary conditions.

###  2.2 Turbulent Driving

The $\mathrm{{\mn@boldsymbol{\mathit{F}}}}$ in Equation 4 is an Ornstein-
Uhlenbeck (OU) process that satisfies the stochastic differential equation,

$\operatorname{d}\\!{\hat{\mathrm{{\mn@boldsymbol{\mathit{F}}}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}},t)=F_{0}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})\mathscr{P}_{ij}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})\operatorname{d}\\!{\mathrm{{\mn@boldsymbol{\mathit{W}}}}}(t)-\hat{\mathrm{{\mn@boldsymbol{\mathit{F}}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}},t)\frac{\operatorname{d}\\!{t}}{T},$
(7)

where $\mathrm{{\mn@boldsymbol{\mathit{W}}}}$(t) is a Wiener process which
creates a random Gaussian increment for the forcing field,
$\mathrm{{\mn@boldsymbol{\mathit{F}}}}$. The increment is then projected onto
the forcing field in $k$-space at $k\sim 2$, i.e. $\sim 1/2$ the box-size,
using a projection tensor, $\mathscr{P}_{ij}$,

$\mathscr{P}_{ij}=\zeta\left(\delta_{ij}+\frac{k_{i}k_{j}}{|k|^{2}}\right)+(1-\zeta)\frac{k_{i}k_{j}}{|k|^{2}},$
(8)

where $\delta_{ij}$ is the Kronecker delta tensor. We control the contribution
from each of the driving modes through the $\zeta$ parameter, e.g. $\zeta=1$
is for purely solenoidal driving in $\mathrm{{\mn@boldsymbol{\mathit{F}}}}$,
and $\zeta=0$ produces purely compressive driving (see Federrath et al. 2008;
Federrath et al. 2009; Federrath et al. 2010 for a detailed discussion of the
driving). Here we choose a natural mixture of the two modes, $\zeta=0.5$
(Federrath et al., 2010). The $T$ in Equation 7 is the autocorrelation
timescale of the OU process. For example, if we only consider the second term
in Equation 7,

$\hat{\mathrm{{\mn@boldsymbol{\mathit{F}}}}}(\mathrm{{\mn@boldsymbol{\mathit{k}}}},t)\sim
F_{0}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})\exp\left\\{-t/T\right\\},$ (9)

the autocorrelation timescale defines the characteristic time that the field
has lost $1/e$ of its previous structure. The timescale is equal to
$T=L/(2c_{s}\operatorname{\mathcal{M}})$, hence we use $T$ to set the desired
turbulent Mach number of the cloud model. For the 20 simulations we vary
$\operatorname{\mathcal{M}}$ between 2 and 20, encompassing the range of
observed $\operatorname{\mathcal{M}}$ values for molecular clouds (e.g.
Schneider et al. 2013; Federrath et al. 2016; Orkisz et al. 2017).

###  2.3 Initial Conditions and Magnetic Fields

The initial magnetic field $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$ in
Equations 4–6 is set to a uniform value with field lines threaded through the
$\hat{\mathrm{{\mn@boldsymbol{\mathit{z}}}}}$ direction of the cloud. We set
the mean-field component of $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$, i.e.,
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}=2c_{s}\sqrt{\pi\rho_{0}}(\operatorname{\mathcal{M}}/\operatorname{\mathcal{M}_{\text{A}}}_{0})=B_{z}$,
where $\operatorname{\mathcal{M}_{\text{A}}}_{0}$ is the mean-field Alfvén
Mach number. Hence by setting $\operatorname{\mathcal{M}}$ and
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}$ we can set the target
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$. We explore the supersonic gas
dynamics of the clouds in the highly-magnetised, sub-Alfvénic regime,
$\operatorname{\mathcal{M}_{\text{A}}}_{0}\approx 0.1$ and 0.5, and in the
trans- and super-Alfvénic regimes,
$\operatorname{\mathcal{M}_{\text{A}}}_{0}\approx 1$, 2 and 10. These values
allow us to explore the transition between the two regimes, which is important
to see how the anisotropic behaviour of the column densities changes as the
magnetic and thermal energies contribute in different proportions to the cloud
dynamics.

We run the simulations to 10$\,T$, where $T$ is eddy the turnover time,
equivalent to the autocorrelation timescale discussed in §2.2. We extract
column density maps, $\Sigma$, from the volume density along the
$\hat{\mathrm{{\mn@boldsymbol{\mathit{y}}}}}$ direction every $0.1\,T$,
between $5\,T$ and $10\,T$, ensuring that the turbulent cloud is in a fully-
developed, stationary statistical regime (Federrath et al., 2009; Price &
Federrath, 2010). We however perform the same analysis outlined in the
forthcoming sections, §4 and §6, on different column density projection angles
in Appendix C. It also ensures that the initial conditions,
$\rho(x,y,z,t=0)=\rho_{0}=1$, where $\rho_{0}$ is the mean density of the
cloud, $\mathrm{{\mn@boldsymbol{\mathit{v}}}}(t=0)=0$ and
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}_{0}(t=0)=B_{z}\hat{\mathrm{{\mn@boldsymbol{\mathit{z}}}}}$,
do not influence the statistics of the flow that we analyse. Our simulations
are dimensionalised so that $\rho(x,y,z)$ is in units of mean density,
$\rho_{0}$. The reason for this is so that our results do not depend upon
choice of $\rho_{0}$ and are solely determined by the
$\operatorname{\mathcal{M}}$ and $\operatorname{\mathcal{M}_{\text{A}}}_{0}$,
i.e. we could scale $\rho_{0}$ to any value as long as we also scale
$|\mathrm{{\mn@boldsymbol{\mathit{B}}}}|$, $c_{s}$, and $T$, the
autocorrelation time, such that $\operatorname{\mathcal{M}}$ and
$\operatorname{\mathcal{M}_{\text{A}}}$ stay the same. We extract column
density maps, which are in the same dimensionlised form as, $\rho(x,y,z)$,
every $0.1\,T$, which results in a total of 51 realisations of the column
density for each of the simulations.

##  3 Column Density Maps

Figure 1 shows a single realisation of the rich structure of the column
density at $5\,T$. From left-to-right in the plot,
$\operatorname{\mathcal{M}}$ increases from 2 to 20, and from top-to-bottom,
$\operatorname{\mathcal{M}_{\text{A}}}$ increases from 0.1 to 10. Hence the
strongest $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field and turbulence is in
the top-right corner, and the weakest in the bottom-left corner. In the strong
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field,
low-$\operatorname{\mathcal{M}}$ models we find significant striations
_parallel_ to the magnetic field. By contrast, in the strong
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field,
high-$\operatorname{\mathcal{M}}$ models, high-density filamentary structures
and low-density voids develop _perpendicular_ to the magnetic field. The
parallel striations dominate the morphology of the density field for the
$\operatorname{\mathcal{M}}=2$ simulation, while the high-density filamentary
structures dominate the column density for the $\operatorname{\mathcal{M}}\geq
4$ simulations. This suggests that there is an interplay between
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$ and $\operatorname{\mathcal{M}}$
on the mode of alignment for filaments with the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field. The strong-field simulations
are the most interesting for the purpose of exploring the anisotropic
behaviour, whereas the weak-field simulations, i.e.
high-$\operatorname{\mathcal{M}_{\text{A}}}_{0}$, qualitatively appear to be
in an isotropic regime, where anisotropic structures are present, but they are
orientated randomly, contributing to no systematic, preferred direction. To
better understand and quantify how the anisotropic structures in the sub-
Alfvénic regime contribute to the observed properties of the column density we
now explore the 2D power spectrum.

Figure 1: A single realisation at $t/T=5$ of the $xz$-projected column
densities from the 20 simulations analysed in this study. The column density
maps are shown on a logarithmic scale and in units of the mean column density,
$\log_{10}\left(\Sigma/\Sigma_{0}\right)$, with increasing
$\operatorname{\mathcal{M}}$ (shown in the top-left corner) from left to right
and decreasing $\operatorname{\mathcal{M}_{\text{A}}}_{0}$ (shown in the top-
right corner) from top to bottom. The magnetic field is oriented along the
$z$-axis (up the page). An animation of this figure, revealing the dynamics,
is available in the online version. Figure 2: 2D power spectrum of the zero-
mean column density defined in Equation (10), averaged over 5$\,T$. Each panel
shows a different simulation, with the sonic and Alfvén Mach number labelled
at the top of each panel as in Figure 1. The colour range is logarithmically
scaled from red to yellow, with the range shown at the bottom of each panel.
Smoothed iso-surfaces of power are plotted every half order of magnitude,
shown with red contours. The
$\mathrm{{\mn@boldsymbol{\mathit{k}}}}_{\perp},\mathrm{{\mn@boldsymbol{\mathit{k}}}}_{\parallel}$
origin is shown with the black coordinate system at the centre of each of the
power spectra. Ellipses are fit to the iso-surfaces, shown in blue, to
approximate the anisotropy parameter discussed in the main text.

##  4 The 2D Power Spectra

The power spectrum allows us to reveal order out of complicated and stochastic
structures that are mostly intangible in real space. There are, however, also
techniques for directly analysing real-space structures, such as structure
functions, topological (Appleby et al., 2018; Henderson et al., 2019), and
fractal methods (Scalo, 1990; Elmegreen & Falgarone, 1996; Stutzki et al.,
1998; Kowal et al., 2007; Federrath et al., 2009; Roman-Duval et al., 2010;
Donovan Meyer et al., 2013; Konstandin et al., 2016; Beattie et al., 2019b).
The power spectrum has a special role in the study of turbulence, since
turbulence models rely heavily upon understanding flow characteristics on
different length scales in real space, e.g., on the driving scale, in the
inertial (for incompressible flows) scaling range (cascade) of turbulence, and
on the dissipation scale (Kolmogorov, 1941; Burgers, 1948). The benefit of
using the power spectrum to describe these scales is that we can organise them
into $k$-modes, symmetrical (in the isotropic case) about an origin, rather
than distributed through real space. However, before constructing the power
spectrum we first must consider the quantity that we are taking the power
spectra of, i.e. in this study, the column density.

We transform the column density into a zero-mean field using the transform,

$\xi(x,y)=\Sigma(x,y)/\Sigma_{0}-1,$ (10)

where $\Sigma_{0}$ is the mean of the column density. Zero-mean fields are
useful for constructing a power spectrum of the field. This is because one can
use the power spectrum to measure moments of the field (Brunt et al., 2010a,
b; Konstandin et al., 2016), e.g., Parsevel’s theorem applied to a stochastic
field $f$ is,

$\int_{\forall k}\operatorname{d}\\!{k}\,\mathscr{P}_{f}=\mathrm{\left\langle
f^{2}\right\rangle},$ (11)

where $\mathrm{\left\langle\ldots\right\rangle}$ is the ensemble average
operator. We can relate Parsevel’s theorem to the 1${}^{\text{st}}$ and
2${}^{\text{nd}}$ moments of the field by adding and subtracting the square of
the 1${}^{\text{st}}$ moment, $\mathrm{\left\langle f\right\rangle}^{2}$,

$\int_{\forall
k}\operatorname{d}\\!{k}\,\mathscr{P}_{f}(k)=\mathrm{\left\langle
f\right\rangle}^{2}+\mathrm{\left\langle
f^{2}\right\rangle}-\mathrm{\left\langle
f\right\rangle}^{2}=\mu^{2}+\sigma^{2},$ (12)

where $\mu$ is the mean (1${}^{\text{st}}$ moment) and $\sigma^{2}$ is the
variance. For a zero-mean field, $g$, where $\mu=0$

$\int_{\forall k}\operatorname{d}\\!{k}\,\mathscr{P}_{g}(k)=\sigma^{2}.$ (13)

Zero-mean fields are therefore of special interest, since the integral of the
power spectrum for a zero-mean field is exactly the variance of the field.
Brunt et al. (2010a) and Federrath et al. (2016) use this property to
reconstruct 3D volume-density distributions from the 2D column density
distribution (note however that Brunt et al. (2010a) removes the mean in
$k$-space rather than in real space). All the power spectra that we construct
in this study will be of the zero-mean transformed column density, as defined
in Equation (10).

Figure 3: The anisotropic structures that dominate the morphology of the 2D
power spectrum for the column density revealed in real space. We show the two
sub-Alfvénic cases: simulations M2Ma0.1 (top panels) and M20Ma0.1 (bottom
panels). In the first column we show the column densities for each simulation,
discussed in §3. In the second column we show Fourier transforms (Equation 14)
of the zero-mean transformed column densities (Equation 10) convolved with the
annulus, $\mathcal{E}_{\sigma}$ (discussed in more detail in Appendix A),
derived from the elliptic fits shown in Figure 2. In the third column we show
the column density structure corresponding to the modes in each of the
annulus, i.e., the real-space structures responsible for the anisotropies in
the 2D power spectra discussed in §4. For M2Ma0.1 we find highly orientated,
thin striations that follow the direction of the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field. By contrast, for M20Ma0.1, we
additionally see high-density filaments aligned perpendicular to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field. An animation of this figure
over the full 10$\,T$ is available in the online version.

###  4.1 Construction of the Power Spectra

The zero-mean transformed column density, $\xi(x,y)$, is a discrete field with
a finite number of resolution elements. Hence we must use the discrete Fourier
transform, $\mathscr{F}[\ldots]$, to map the density into $k$-space,

$\hat{\xi}(\mathrm{{\mn@boldsymbol{\mathit{k}}}}_{xy})=\mathscr{F}\left[\xi(\mathrm{{\mn@boldsymbol{\mathit{\ell}}}}_{xy})\right]=\frac{1}{N^{2}}\displaystyle\sum^{N-1}_{x=0}\sum^{N-1}_{y=0}\xi(\mathrm{{\mn@boldsymbol{\mathit{\ell}}}}_{xy})\exp\left\\{-i\frac{\mathrm{{\mn@boldsymbol{\mathit{\ell}}}}_{xy}\cdot\mathrm{{\mn@boldsymbol{\mathit{k}}}}_{xy}}{N}\right\\},$
(14)

where $\mathrm{{\mn@boldsymbol{\mathit{\ell}}}}_{xy}$ is the real space vector
with coordinates $(x,y)$, $\mathrm{{\mn@boldsymbol{\mathit{k}}}}_{xy}$ is the
corresponding $k$-space vector,
$\mathrm{{\mn@boldsymbol{\mathit{k}}}}_{xy}=(k_{x},k_{y})$, $N$ is the grid
size in one dimension (i.e., $N=512$ for all simulations) and
$i\equiv\sqrt{-1}$. The 2D column-density power spectrum is the square of the
Fourier transform, hence,

$\mathscr{P}_{\xi}(\mathrm{{\mn@boldsymbol{\mathit{k}}}})=\hat{\xi}\cdot\hat{\xi}^{*}=\mathscr{P}_{\Sigma/\Sigma_{0}-1},$
(15)

where $\hat{\xi}^{*}$ is the complex conjugate of $\hat{\xi}$. We show the
time-averaged (across all available turnover time realisations) 2D power
spectra of the zero-mean column density for each simulation in Figure 2. Shown
in red are contours of equal power, plotted every half order of magnitude in
power.

For each of the closed contours we fit an ellipse through five degrees of
freedom following the generalised eigensystem method outlined in Fitzgibbon et
al. (1996). The five degrees of freedom are translations in
$\operatorname{k_{\parallel}}$ and $\operatorname{k_{\perp}}$, rotations of
the ellipse in the
$\operatorname{k_{\parallel}}$-$\operatorname{k_{\perp}}$-plane, and the
modulus of the major and minor principle axes of the ellipse. The purpose of
the fit is to approximate the power spectrum with a morphology that captures
the main anisotropic features at each $\operatorname{k_{\perp}}$ scale. In
particular, we want to know how the power spectrum varies along
$\operatorname{k_{\parallel}}$, compared with $\operatorname{k_{\perp}}$, at
each of the equi-power contours. The elliptic fits lead to a canonical
anisotropy parameter,

$a\equiv\frac{e_{\parallel}}{e_{\perp}},$ (16)

where $e_{\parallel}$ is the principle axis of the ellipse along the
$\operatorname{k_{\parallel}}$ direction, and $e_{\perp}$ is the principle
axis along the $\operatorname{k_{\perp}}$ direction. This means that for
$a<1$, $e_{\parallel}<e_{\perp}$, hence the power spectrum is elongated along
the $\operatorname{k_{\perp}}$ axis. Likewise, $a>1$ implies
$e_{\parallel}>e_{\perp}$, which corresponds to elongation along the
$\operatorname{k_{\parallel}}$ axis. For $a=1$ (circle), the power spectrum is
isotropic. We note that $a$ may depend upon the length scale, i.e., what is
isotropic on large length scales need not be isotropic on small length scales
or vice versa. This means that $a$ may change for different ellipses, at
different $\operatorname{k_{\perp}}$ scales in the same power spectra. We show
this by choosing to write $a$ as a function of $\operatorname{k_{\perp}}$,
$a(\operatorname{k_{\perp}})$, consistent with the MHD theory outlined in
§1.1.

In the following sections we discuss the results from the elliptic fits to the
power spectra. In §4.2 we focus on the super-Alfvénic simulations, and in §4.3
we discuss the sub-Alfvénic simulations.

###  4.2 Super-Alfvénic Regime

We start our analysis in the bottom eight panels of Figure 2. In these
simulations, the magnetic field is weak, placing the turbulence in the super-
Alfvénic regime ($\operatorname{\mathcal{M}_{\text{A}}}_{0}>1$). Our elliptic
fits to the contours reveal that on all $k$-scales the power spectrum remains
almost entirely circular, i.e., the ratio between $e_{\parallel}\sim
e_{\perp}$, and thus, $a\sim 1$. Perfect rotational symmetry of the power
spectrum indicates isotropic statistics, reminiscent of hydrodynamical (non-
magnetized) turbulence. We can state this mathematically as

$\displaystyle\lim_{\operatorname{\mathcal{M}_{\text{A}}}_{0}>1}a(\operatorname{k_{\perp}})\sim
1,$ (17)

which means that in the super-Alfvénic regime, the anisotropy factor $a\sim 1$
for all length scales $\operatorname{k_{\perp}}$. This situation changes with
increasing magnetic guide field, introducing scale-dependent anisotropies, as
discussed in the next subsection.

###  4.3 Trans and Sub-Alfvénic Regime

Now we move onto the trans ($\operatorname{\mathcal{M}_{\text{A}}}_{0}\sim 1$)
and sub-Alfvénic ($\operatorname{\mathcal{M}_{\text{A}}}_{0}<1$) turbulent
regimes, corresponding to the top 12 panels in Figure 2. On quick inspection
one can note two distinct shapes emerging in these regimes. Ellipses stretched
along the $\operatorname{k_{\perp}}$ axis ($a<1$) are present in the fits for
the $\operatorname{\mathcal{M}}\lesssim 4$ simulations, and ellipses stretched
along the $\operatorname{k_{\parallel}}$ axis ($a>1$) for the
$\operatorname{\mathcal{M}}\gtrsim 4$ simulations. These are two distinct
anisotropies that we will discuss in the next two subsections in detail. We
begin with the low-$\operatorname{\mathcal{M}}$ simulations.

####  4.3.1 $\operatorname{\mathcal{M}}\lesssim 4$

For $\operatorname{\mathcal{M}}\lesssim 4$ (top-left corner of Figure 2) the
power spectrum is elongated along the $\operatorname{k_{\perp}}$ axis. This
means that length scales along the magnetic field carry more power on the same
length scales perpendicular to the magnetic field, corresponding to an
anisotropic factor $a<1$, as described in §4.2. We observe the same kind of
anisotropy at both large $k$, $\sim k_{\text{max}}$ (small spatial scales) and
small $k$, $\sim k_{\text{min}}$ (box-size spatial scales). This is especially
obvious in the M2Ma0.1 and M2Ma0.5 simulations, where the anisotropy grows on
smaller $k$ scales (the 2D power spectrum becomes more elongated on large
spatial scales, corresponding to small $k$). We state this concisely,

$a(\operatorname{k_{\perp}}\sim k_{\text{min}})<1,$ (18)

at small $k$, but also at large $k$, hence

$a(\operatorname{k_{\perp}}>k_{\text{min}})<1,$ (19)

but the anisotropy is strongest at small $k$, hence,

$a(\operatorname{k_{\perp}}\sim
k_{\text{min}})<a(\operatorname{k_{\perp}}>k_{\text{min}})<1.$ (20)

Returning to Figure 1 briefly, we can hypothesise that in real space this
corresponds to the density striations parallel to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field dominating the morphology of the
column density on $L$ scales (the box-size) and throughout all of the smaller
scales local to $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$. Striations are most
likely formed through fast magnetosonic waves moving perpendicular to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field (Tritsis & Tassis, 2016; Tritsis
et al., 2018), which is why we see them dominating some of the sub-Alfvénic
density structures. The reason why these structures lead to stronger
anisotropies at small $k$ is because they span coherently across the entire
box-length. Figure 11 in the Appendix shows the 3D density structures for the
M1Ma0.1 and M20Ma0.1 simulations. In 3D, the striations appear to be sheets of
density that are sheared along the field lines and then twisted tightly around
them by the turbulent eddies. We will explore these 3D structures in future
studies but for now we keep our focus on the column-density structures.

To test our hypothesis we use the elliptic fits (discussed in §4.1) in
$k$-space to extract the real-space structures. The M2Ma0.1 simulation is
shown in the top row of Figure 3. The first column is the column density map,
the same as in Figure 1. The second column is the real component of the 2D
Fourier transform of $\xi$ (Equation 14), the zero-mean transformed column
density (Equation 10) with a selected annular region, $\mathcal{E}_{\sigma}$
defined by our elliptic fits222Here $\sigma$ corresponds to the size of a
Gaussian smoothing kernel that we apply to the annular region before doing the
real-space inversion. We do this to avoid the Gibbs phenomenon at the edges of
the annulus, i.e., $\mathcal{E}_{\sigma}=\mathcal{E}*\mathcal{N}(\sigma)$,
where $\mathcal{E}$ is the unsmoothed annulus and $\mathcal{N}(\sigma)$ is a
Gaussian smoothing function. We pick the smoothing kernel,
$\sigma=2\operatorname{d}\\!{x}$, where $\operatorname{d}\\!{x}$ is the size
of a grid cell, for our problem. For more details on this method we refer to
Appendix A.. We perform the inverse Fourier transform on only the $k$-modes
inside of $\mathcal{E}_{\sigma}$, which are shown in the third column of
Figure 3. This figure reveals that the $a<1$ anisotropy in the power spectrum
is indeed caused by magnetosonic striations running parallel to
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$.

The striations revealed in Figure 3 are strongly aligned with the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field, consistent with previous
observational analyses (Cox et al., 2016; Malinen et al., 2016; Tritsis et
al., 2018), which indicates a trans- to sub-Alfvénic flow regime. Indeed, the
striations are such a highly anisotropic feature of the column density and so
strongly aligned with the $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field that
the 2D power spectrum becomes significantly pinched along the
$\operatorname{k_{\perp}}$ direction. This is evident in the red contours in
the top-left four panels in Figure 2. As $\operatorname{k_{\perp}}$ grows, so
do the pinched ends of the contours. What this means is that at smaller real
space length scales the densities in the cloud are becoming completely
dominated by the magnetic field, i.e., the small-scale structures (Equation
19) along the $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field essentially become
uncorrelated with the hydrodynamics. This is because when a power spectrum
lacks diagonal modes, the dynamics in each of the directions is not
correlated, i.e., the dynamics along the magnetic field is independent of the
dynamics perpendicular to the field (Tennekes & Lumley, 1972).

As the Mach number increases (i.e., as we move right across the panels in
Figure 2), we quickly see the transition and development of a different kind
of anisotropy. The transition occurs at $\operatorname{\mathcal{M}}\sim 4$
(the second column of Figure 2), where the power spectrum is full of
uncorrelated parallel and perpendicular anisotropic structures resulting in
nearly square-shaped contours. The transition occurs at large $k$ first, with
small $k$ (large structures), still being significantly aligned with the
magnetic field. This is because there are still strongly-aligned, parallel
striations on the box-scale, however, gas is becoming more compressed along
the field lines, which generates turbulent shocks, perpendicular to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field. As $\operatorname{\mathcal{M}}$
continues to increase we see this new, $a>1$ anisotropy beginning to dominate
the morphology. We now turn our attention to this type of anisotropy.

####  4.3.2 $\operatorname{\mathcal{M}}\gtrsim 4$

For $\operatorname{\mathcal{M}}\gtrsim 4$ the power contours in Figure 2 form
an $a>1$ anisotropy, where more power is in the $\operatorname{k_{\perp}}$
direction than on the same length scale in the $\operatorname{k_{\parallel}}$
direction. This means that our elliptic fits on the column density power
spectrum become elongated along the $\operatorname{k_{\parallel}}$ direction
at large $k$, i.e.,

$a(\operatorname{k_{\perp}}>k_{\text{min}})>1.$ (21)

However, at small $k$, our elliptic fits are near circular, i.e., on large
scales close to $L$, the power spectra return to being isotropic,

$a(\operatorname{k_{\perp}}\sim k_{\text{min}})\sim 1,$ (22)

in contrast with the $\operatorname{\mathcal{M}}\lesssim 4$, sub-Alfvénic
simulations discussed in §4.3.1, where the anisotropic effects of the magnetic
field do not disappear on any scale, large or small. Indeed, Equation 22 is
true for all high-$\operatorname{\mathcal{M}}$ simulations, where the large-
scale turbulent energies overcome the preferential orientation imparted by the
magnetic fields,

$\displaystyle\lim_{\operatorname{\mathcal{M}}\gg
1}a(\operatorname{k_{\perp}}\sim k_{\text{min}})\sim 1.$ (23)

However, the most important feature of the $\operatorname{\mathcal{M}}\gtrsim
4$ power spectra is the $a>1$ anisotropy. The elongation along the
$\operatorname{k_{\parallel}}$ axis, at large $k$ indicates that there are
structures perpendicular to the magnetic field that do not span the entire box
length. We reveal what these structures are in the right bottom row of Figure
3, using an annulus derived from our elliptic fits, as discussed in §4.3.1,
but for the M20Ma0.1 simulation. The bottom right panel reveals the density
structures that are causing the $a>1$ anisotropy. We see many, high-density,
filamentary structures perpendicular to the magnetic field direction,
consistent with orientation studies of the Gould Belt molecular clouds (Planck
Collaboration et al., 2016a, b). Some parallel striations are still present,
since the magnetic field is still dynamically important, but the perpendicular
filaments dominate the morphology of the power spectrum.

Federrath (2016) and Xu et al. (2019) independently suggested that the width
of the high-density filaments,

$\ell\propto L/\operatorname{\mathcal{M}}^{2},$ (24)

using two different models. Federrath (2016) uses the turbulent sonic scale
and Xu et al. (2019) use the Rankine-Hugoniot jump conditions to motivate
Equation (24). In $k$-space, this means that
$\operatorname{k_{\parallel}}\sim\operatorname{\mathcal{M}}^{2}$, and so once
the filamentary structures become locked perpendicular to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field, increasing
$\operatorname{\mathcal{M}}$ grows the anisotropy along the
$\operatorname{k_{\parallel}}$ direction. This is because as the perpendicular
filaments become thinner they contribute to more stretching along
$\operatorname{k_{\parallel}}$. Furthermore, since the number of filaments
present in isothermal turbulence is in someway controlled by
$\operatorname{\mathcal{M}}$ (Beattie & Federrath, 2019), we expect that as
$\operatorname{\mathcal{M}}$ increases, more filaments develop, stretching the
power spectrum further. Hence for observations of MCs with highly-orientated,
high-density filaments, the magnetic field must be dynamically important, and
at the same time, the turbulent motions must be significantly supersonic,
i.e., $\operatorname{\mathcal{M}}\gtrsim 4$.

###  4.4 Summary

We find that the anisotropy parameter $a$, defined in Equation (16), has two
dominant modes in the column density. The first mode, $a<1$, is elongation
along the $\operatorname{k_{\perp}}$ direction, which is due to striations
that run parallel to the $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field,
examined in §4.3.1. The $a<1$ anisotropy persists across all $k$-scales,
similar to the magnetic field, and indeed, seems to be largely controlled by
the magnetic field. The second mode, $a>1$, is elongation along the
$\operatorname{k_{\parallel}}$ direction in the highly-turbulent, strong
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field regime, discussed in §4.3.2.
This anisotropy coincides with the development of high-density filaments that
form perpendicular to the magnetic field, and that are largely controlled by
the turbulent Mach number (e.g., in their width). Thus, our results strongly
suggest that MC observations revealing highly-orientated, high-density
filaments are in a sub-Alfvénic, highly-supersonic regime
($\operatorname{\mathcal{M}}>4$). On large scales, the $a>1$ anisotropy
disappears, coinciding with large-scale isotropic, turbulent motions. More
generally, for all simulations in the super-Alfvénic limit
($\operatorname{\mathcal{M}_{\text{A}}}_{0}>1$), the density and flow
structures are largely isotropic, regardless of $\operatorname{\mathcal{M}}$,
as discussed in §4.2.

It should be noted that our ensemble of simulations do not include any with
gravitational dynamics, and the discussed anisotropic structures emerge purely
in the presence of strongly-orientated magnetic fields and turbulent motions.
Indeed, this suggests that these structures can be a product of magnetic
fields and turbulence alone, however one might expect to observe more extreme
and time-dependant anisotropic structures developing in clouds governed by
gravito-magnetohydrodynamics. This is because gas that is accreted onto the
high-density filamentary structures will cause gravitationally unstable
filaments to collapse and create thin, orientated structures, further
contributing to the stretching of the 2D power spectrum. We leave this
investigation for future studies of the time-dependent anisotropy in
collapsing molecular cloud systems.

##  5 The 1D Power Spectra

Typical in spectral analysis is the reduction of the 2D power spectrum (or 3D
power spectrum) into the 1D azimuthally-averaged power spectrum. However, the
process of azimuthally averaging the power spectrum assumes that the
turbulence is isotropic, i.e., $\mathrm{{\mn@boldsymbol{\mathit{k}}}}$ is
rotationally symmetric in $k$-space, hence the power spectrum only depends
upon the modulus of $\mathrm{{\mn@boldsymbol{\mathit{k}}}}$. For the
anisotropic 2D power spectra discussed in §4, however, no such rotational
symmetry exists in the sub-Alfvénic regime (at least in the LOS perpendicular
to the magnetic field). For this reason, we decompose the 2D power spectra
into $\operatorname{k_{\perp}}$ and $\operatorname{k_{\parallel}}$ components,
with respect to the orientation of the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field. This has been done before in
Bigot et al. (2008), for example.

Table 2: Power-law scaling exponents for the perpendicular and parallel
cascades in the zero-mean transformed column density.

Simulation | $\alpha_{\perp}$ | $\alpha_{\parallel}$ | $\alpha_{\parallel}/\alpha_{\perp}$
---|---|---|---
ID | $(\pm 1\sigma)$ | $(\pm 1\sigma)$ | $(\pm 1\sigma)$
M2Ma10 | $-3.1\pm 0.3$ | $-3.3\pm 0.2$ | $1.1\pm 0.1$
M2Ma2 | $-3.1\pm 0.2$ | $-3.4\pm 0.2$ | $1.1\pm 0.1$
M2Ma1 | $-2.7\pm 0.2$ | $-3.8\pm 0.2$ | $1.4\pm 0.1$
M2Ma0.5 | $-1.7\pm 0.2$ | $-4.2\pm 0.2$ | $2.5\pm 0.3$
M2Ma0.1 | $-1.4\pm 0.2$ | $-4.2\pm 0.3$ | $3.0\pm 0.5$
M4Ma10 | $-2.8\pm 0.1$ | $-3.0\pm 0.2$ | $1.1\pm 0.1$
M4Ma2 | $-2.9\pm 0.2$ | $-3.3\pm 0.2$ | $1.1\pm 0.1$
M4Ma1 | $-2.7\pm 0.2$ | $-3.1\pm 0.1$ | $1.2\pm 0.1$
M4Ma0.5 | $-1.5\pm 0.2$ | $-3.7\pm 0.3$ | $2.5\pm 0.3$
M4Ma0.1 | $-1.7\pm 0.2$ | $-3.6\pm 0.3$ | $2.1\pm 0.3$
M10Ma10 | $-2.8\pm 0.2$ | $-2.5\pm 0.2$ | $0.9\pm 0.1$
M10Ma2 | $-2.3\pm 0.2$ | $-2.6\pm 0.3$ | $1.1\pm 0.2$
M10Ma1 | $-2.4\pm 0.2$ | $-2.5\pm 0.2$ | $1.0\pm 0.1$
M10Ma0.5 | $-2.0\pm 0.2$ | $-2.5\pm 0.3$ | $1.3\pm 0.2$
M10Ma0.1 | $-1.8\pm 0.2$ | $-2.3\pm 0.2$ | $1.3\pm 0.2$
M20Ma10 | $-2.3\pm 0.2$ | $-2.4\pm 0.2$ | $1.0\pm 0.1$
M20Ma2 | $-1.9\pm 0.1$ | $-2.4\pm 0.2$ | $1.3\pm 0.1$
M20Ma1 | $-2.2\pm 0.2$ | $-2.3\pm 0.2$ | $1.0\pm 0.1$
M20Ma0.5 | $-2.0\pm 0.1$ | $-2.3\pm 0.2$ | $1.2\pm 0.1$
M20Ma0.1 | $-1.7\pm 0.2$ | $-2.2\pm 0.2$ | $1.3\pm 0.2$

  * •

Notes: Column (1): the simulation ID. A full list of the simulation parameters
for each ID is shown in Table 1. Column (2): the perpendicular (to the
magnetic field) power-law scaling exponent,
$\mathscr{P}_{\xi}\propto\operatorname{k_{\perp}}^{\alpha_{\perp}}$ for the
zero-mean transformed column density (Equation 10). Column (3): the same as
column (2) but the parallel power-law scaling exponent,
$\mathscr{P}_{\xi}\propto\operatorname{k_{\parallel}}^{\alpha_{\parallel}}$.
Column (4): the ratio between the two scaling exponents from columns (2) and
(3).

Figure 4: Slices of the 2D power spectra across the $\operatorname{k_{\perp}}$
(left panels) and $\operatorname{k_{\parallel}}$ (right panels) directions. We
fit slopes on power between $5\leq k\leq 20$, illustrated in black, and
tabulated in Table 2. The plots are coloured by $\operatorname{\mathcal{M}}$,
indicated in the legend of the bottom-right panel. Figure 5: Left: power-law
scaling exponents, $\alpha$, from the 1D sliced power spectra of the zero-mean
transformed column density plotted as a function of turbulent Mach number. The
crosses indicate the scaling exponent for the perpendicular component of the
1D power spectrum (Equation 26), and the triangles the parallel scaling
exponent (Equation 25). The markers are coloured by the Alfvén Mach number,
where yellow is weak $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field, and dark-
blue is strong $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field. Right: the same
as the left panel but for the ratio between the parallel and perpendicular
scaling exponents. We see the scaling exponents converging to the same value
in the high-$\operatorname{\mathcal{M}}$ limit, indicating the development of
an isotropic 2D power spectrum.

###  5.1 Construction

The simplest way to decompose the 2D power spectrum into 1D
$\operatorname{k_{\perp}}$ and $\operatorname{k_{\parallel}}$ components is to
take slices through the 2D power spectrum. We slice at
$\mathscr{P}_{\xi}(\operatorname{k_{\perp}}=0,\operatorname{k_{\parallel}})$,
for the $\operatorname{k_{\parallel}}$ spectrum and
$\mathscr{P}_{\xi}(\operatorname{k_{\perp}},\operatorname{k_{\parallel}}=0)$,
for the $\operatorname{k_{\perp}}$ spectrum. Next we fit power laws between
wavenumbers $5\leq k\leq 20$, which is a conservative estimate of where the
scaling range might be found in the velocity structure for our simulations. We
denote the power law parallel to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field

$\mathscr{P}_{\xi}\propto\operatorname{k_{\parallel}}^{\alpha_{\parallel}},$
(25)

and perpendicular to the $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field,

$\mathscr{P}_{\xi}\propto\operatorname{k_{\perp}}^{\alpha_{\perp}},$ (26)

where $\alpha_{\parallel}$ and $\alpha_{\perp}$ are the scaling exponents in
the parallel and perpendicular cascades, respectively. We tabulate the
exponents from the power-law fits, along with the ratio between them in Table
2. Unlike the single-point statistics for the anisotropy parameters calculated
in §6, the power-law exponents computed here encode all the anisotropy
information from length scales within the (approximate) cascades.

###  5.2 1D Power Spectrum Results

Figure 4 shows that the slope of the parallel component of the power spectrum
becomes shallower with increasing $\operatorname{\mathcal{M}}$. This is
because parallel to the $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field the
densities do not feel the Lorentz force and are able to create higher-$k$
modes in the column density, corresponding to the development of density
fluctuations and strong shocks (Kim & Ryu, 2005; Kowal et al., 2007). Indeed,
the strong shocks that develop from parallel compression along the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field lines are what creates the high-
density filaments perpendicular to the field, consistent with our results from
§4. We plot the scaling exponents in Figure 5 to better visualise the
$\operatorname{\mathcal{M}}$ and $\operatorname{\mathcal{M}_{\text{A}}}_{0}$
dependencies. The triangle markers, corresponding to the parallel scaling
exponents, have values of $\sim-4$ at low-$\operatorname{\mathcal{M}}$, and
are completely distinct from the perpendicular scaling exponents, shown with
cross markers. However, with increasing $\operatorname{\mathcal{M}}$, both of
the scaling exponents become centralised around $\sim-2$, consistent with the
reoccurring theme of isotropy in the high-$\operatorname{\mathcal{M}}$ limit,
discussed in the context of the 2D power spectrum in §4.2.

Unlike the component of the column density parallel to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field, the perpendicular component
feels very little of the turbulent motions. This is shown by tracing the
$\alpha_{\perp}$ (cross markers) in the left panel of Figure 5 with the same
magnetic field strength (i.e., same colour) and observing how they stay
relatively constant with changing $\operatorname{\mathcal{M}}$. The spread of
the scaling exponents seems to be determined by
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$, where steeper values correspond
to weaker fields, similar to the average anisotropic parameters discussed in
§6. This is consistent with the “magnetic cushioning" principle described in
Molina et al. (2012), whereby the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field suppresses fluctuations in the
density field, resulting in a steeper slope in the power spectrum.

In the right-hand panel of Figure 5 we show the ratio between the parallel and
perpendicular scaling exponents as a function of $\operatorname{\mathcal{M}}$.
In the low-$\operatorname{\mathcal{M}_{\text{A}}}_{0}$,
low-$\operatorname{\mathcal{M}}$ power spectra, to the left of the panel,
$\alpha_{\parallel}$ is $\sim 2.5$, steeper than $\alpha_{\perp}$. This means
that most of the high-$k$ density fluctuations are transverse to the field
lines, which could indeed be associated with the striations discussed in §4.3.
However, as $\operatorname{\mathcal{M}}$ increases, unlike the 2D power
spectrum, there is no distinct signature of the high-density filaments.

Figure 6: Left Panel: the average anisotropy, $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}$ (Equation 27), in the zero-mean
transformed column density (Equation 10) as a function of
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$, coloured by
$\operatorname{\mathcal{M}}$. The red, dashed line indicates perfect isotropy
at $\mathrm{\left\langle a\right\rangle}_{\operatorname{k_{\perp}}}=1$. Above
the isotropic boundary are the clouds with high-density filaments dominating
the morphology of the power spectrum, which we examine in §4.3.2, and clouds
below are dominated by magnetosonic striations, which we discuss in §4.3.1.
Right Panel: the same as the left panel, but as a function of
$\operatorname{\mathcal{M}}$ and coloured by
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$. We fit power laws to data with
the same $\operatorname{\mathcal{M}_{\text{A}}}_{0}$. The fits are indicated
with the dotted lines, and labelled by the power law, $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}\sim\operatorname{\mathcal{M}}^{\beta}$.
Decreasing $\operatorname{\mathcal{M}_{\text{A}}}_{0}$ steepens the power law,
increasing the overall anisotropy, but whether or not the average anisotropy
is above or below unity is controlled by $\operatorname{\mathcal{M}}$.

##  6 The Average Anisotropy in the Column Density

In §4 we qualitatively described the anisotropy present in the column
densities as a function of scale in the 2D power spectra. Next in §5 we
explored how the 1D power spectra perpendicular and parallel to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field encode some (but not all)
information about the anisotropic structures found in §4. Now we summarise the
anisotropy for each model using a single, scale-averaged anisotropy parameter.
We define our scale-averaged anisotropy parameter as

$\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}=\frac{1}{\operatorname{k_{\perp}}_{\text{max}}}\displaystyle\int\displaylimits^{\operatorname{k_{\perp}}_{\text{max}}}_{0}\operatorname{d}\\!{\operatorname{k_{\perp}}}\,a(\operatorname{k_{\perp}}),$
(27)

where we integrate across all anisotropy parameters for all available
$\operatorname{k_{\perp}}$-scales, and divide by the maximum
$\operatorname{k_{\perp}}$, $\operatorname{k_{\perp}}_{\text{max}}$, which is
$\operatorname{k_{\perp}}$ at the last closed contour in the 2D power
spectrum. We show examples of the $a(\operatorname{k_{\perp}})$ curves for the
two extreme anisotropic simulations ($\operatorname{\mathcal{M}}=2$ and
$\operatorname{\mathcal{M}}=20$) in Appendix B. We do this for each power
spectrum, creating a single number that represents the average anisotropy of
the 2D power spectrum, tabulated in column (9) of Table 1, shown with
$1\sigma$ uncertainties defined as

$\sigma_{\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}}=\sqrt{\mathrm{\left\langle
a^{2}\right\rangle}_{\operatorname{k_{\perp}}}-\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}^{2}}.$ (28)

This is similar to what has been done in Appleby et al. (2018), where they use
Minkowski tensors to approximate the anisotropy of dark matter fields. Our
anisotropy parameter is analogous to taking the average of the ratio between
the two eigenvalues of the $(W_{2}^{1,1})_{ij}$ tensor (equation 15 in Appleby
et al. 2018) over all excursion sets on the 2D power spectrum.

###  6.1 $\mathrm{\left\langle a\right\rangle}_{\operatorname{k_{\perp}}}$
Results

Figure 7: The power-law exponent, $\beta$, for the $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}$ – $\operatorname{\mathcal{M}}$
relation, shown in Equation 29 and illustrated in Figure 6, as a function of
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$. We see a smooth transition across
the values of $\operatorname{\mathcal{M}_{\text{A}}}_{0}$ that is bounded
above by $\beta\sim 0.4$ and below by $\beta\sim 0$.

We show the average anisotropy parameter $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}$, defined in Equation (28), as a
function of $\operatorname{\mathcal{M}_{\text{A}}}_{0}$ and
$\operatorname{\mathcal{M}}$ in Figure 6. In the left plot we see the average
anisotropy as a function of $\operatorname{\mathcal{M}_{\text{A}}}_{0}$,
coloured by $\operatorname{\mathcal{M}}$. There are two key features in this
plot. First, as $\operatorname{\mathcal{M}_{\text{A}}}_{0}$ increases, the
average anisotropy disappears (i.e., $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}$ approaches unity). Second, for a
fixed $\operatorname{\mathcal{M}_{\text{A}}}_{0}$, the spread is controlled by
$\operatorname{\mathcal{M}}$. To understand this relationship we show the
average anisotropy as a function of $\operatorname{\mathcal{M}}$ in the right-
hand panel of Figure 6. For a fixed
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$ we find that there exits a power
law in $\operatorname{\mathcal{M}}$,

$\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}\propto\operatorname{\mathcal{M}}^{\beta}.$
(29)

The power laws may not strictly hold for the very high- and
low-$\operatorname{\mathcal{M}}$ limits (e.g., there may be an anisotropic
saturation), but they provide a reasonable approximation for the intermediate
regime covered by the simulations. We show the fits with the dotted lines,
coloured by $\operatorname{\mathcal{M}_{\text{A}}}_{0}$ in Figure 6. At
$\operatorname{\mathcal{M}_{\text{A}}}_{0}\approx 10$ (shown in yellow) we
find that the average anisotropy is approximately constant in
$\operatorname{\mathcal{M}}$, $\beta\sim 0$, shown to the right of Figure 7,
hence in the weak $\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field limit, no
anisotropy exists for all $\operatorname{\mathcal{M}}$, consistent with our
previous qualitative investigations in §4. However, in the strong-
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field regime we see the slope of the
power laws increase until reaching a saturation at $\beta\sim 0.4$ when
$\operatorname{\mathcal{M}_{\text{A}}}_{0}\sim 0.1-0.5$, shown to the left of
Figure 7. As we move along the $\operatorname{\mathcal{M}}$-axis, we capture
the transition from the striation-dominated regime, $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}<1$, to the high-density, filament-
dominated regime, $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}>1$. This shows that the type of
anisotropy, $a>1$ or $a<1$, is controlled by $\operatorname{\mathcal{M}}$, but
when (at what $\operatorname{\mathcal{M}}$) the transition occurs is set by
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$, which is a key result from this
study. Indeed, this means that one may be able to use the anisotropy in the 2D
power spectrum to determine important cloud parameters, such as the mean
component of the magnetic field or the sonic and Alfvénic Mach numbers, which
will be explored in a future study. Next we move onto another measure of the
anisotropy, the perpendicular and parallel 1D power spectra of the column
density.

##  7 Summary and Key Findings

In this study we use 51 time-realisations from 20 high-resolution, turbulent
magnetohydrodynamical (MHD) simulations with mean-field Alfvén Mach numbers,
$0.1\leq\operatorname{\mathcal{M}_{\text{A}}}_{0}\leq 10$, and turbulent Mach
numbers, $2\leq\operatorname{\mathcal{M}}\leq 20$ (simulation parameters are
shown in Table 1). We create two-dimensional (2D) projections (column
densities) perpendicular to the threaded magnetic fields, shown in Figure 1.
We construct 2D power spectra for each of the zero-mean transformed column
densities (Equation 10), and average them over five turbulent crossing times
to show how the average equi-power structures are organised in $k$-space,
shown in Figure 2. We calculate the anisotropy for each simulation by fitting
ellipses to the contours of the 2D power spectra and calculate the ratio, $a$,
between the minor and major principle axis of the fitted ellipses. We reveal
two distinct anisotropies, $a<1$, corresponding to magnetosonic striations,
and $a>1$, corresponding to high-density filaments, distinguished in Figure 3,
and discussed in detail in §4. Next we create 1D slices of the 2D power
spectra, shown in Figure 4, and discuss how the perpendicular and parallel
scaling exponents in the slices reveal some similarities to our 2D spectral
analysis in §4 and depend upon the magnetic field and turbulence parameters.
In §6 we average the anisotropy parameter, $a$, across all spatial scales in
each simulation and show how it depends upon $\operatorname{\mathcal{M}}$ and
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$ in Figure 6. In the following, we
list the keys results from this study.

  * •

The 2D power spectrum of the column density (Figure 2) reveals two distinct
types of anisotropy, $a$, in the sub-Alfvénic MHD turbulence regime, neither
of which are present in the super-Alfvénic regime. Both types of anisotropies
exhibit scale-dependency, discussed in detail in §4.3.1 and §4.3.2.

  * •

The column density of sub-Alfvénic, $\operatorname{\mathcal{M}}\lesssim 4$
turbulence is dominated by magnetosonic striations running parallel to the
magnetic field (top panels in Figure 3). This results in significant
stretching of the 2D power spectrum along the $\operatorname{k_{\perp}}$ axis,
which is defined as an $a<1$ anisotropy using our anisotropy definition
(§4.3.1), and shallows the slope of the perpendicular cascade in the 1D
perpendicular power spectrum (§5). The anisotropy from the striations becomes
stronger on large (box) scales, since the long structures, parallel to the
$\mathrm{{\mn@boldsymbol{\mathit{B}}}}$-field, span across the full box-
length. The presence of striations in the observed column densities of
molecular clouds may therefore mean that the cloud has a dynamically important
magnetic field.

  * •

For the sub-Alfvénic, $\operatorname{\mathcal{M}}\gtrsim 4$ simulations, high-
density filaments, perpendicular to the magnetic field stretch the 2D power
spectrum of the column density along the $\operatorname{k_{\parallel}}$ axis
(bottom panels of Figure 3), resulting in an $a>1$ anisotropy, distinct from
the $a<1$ anisotropy in $\operatorname{\mathcal{M}}\lesssim 4$ turbulence
(§4.3.2). For this regime, the anisotropy is only present on scales smaller
than the box-length. At the box-length the turbulence returns to the isotropic
regime. Indeed, molecular cloud observations that reveal highly-orientated,
high-density filaments are therefore likely to be in the sub-Alfvénic, highly-
supersonic ($\operatorname{\mathcal{M}}\gtrsim 4$) dynamical regime.

  * •

The anisotropy parameter averaged over all spatial scales,
$\mathrm{\left\langle a\right\rangle}_{\operatorname{k_{\perp}}}$,
approximately follows a power law for fixed magnetic field strength in
$\operatorname{\mathcal{M}}$, $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}\propto\operatorname{\mathcal{M}}^{\beta}$
(see §6 and Figure 6). The value of $\beta$ changes from $\sim 0-0.4$ (see
Figure 7) for $\operatorname{\mathcal{M}_{\text{A}}}_{0}\approx 10-0.1$,
respectively. This means that $\operatorname{\mathcal{M}_{\text{A}}}_{0}$
largely controls the degree of anisotropy, but it is
$\operatorname{\mathcal{M}}$ that controls the transition between $a<1$ and
$a>1$, i.e., the stretching of the 2D power spectrum either along the
$\operatorname{k_{\perp}}$ or $\operatorname{k_{\parallel}}$ axis,
respectively.

The quantification of anisotropies in the column density may be useful for
constraining cloud parameters such as
$\operatorname{\mathcal{M}_{\text{A}}}_{0}$ and $\operatorname{\mathcal{M}}$
in follow-up studies.

## Acknowledgements

We thank the anonymous reviewer for the useful suggestions that improved this
study. J. R. B. thanks Jonah Hansen for helpful editorial comments. C. F.
acknowledges funding provided by the Australian Research Council (Discovery
Project DP170100603, and Future Fellowship FT180100495), and the Australia-
Germany Joint Research Cooperation Scheme (UA-DAAD). We further acknowledge
high-performance computing resources provided by the Leibniz Rechenzentrum and
the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale
project 10391), the Partnership for Advanced Computing in Europe (PRACE grant
pr89mu), the Australian National Computational Infrastructure (grant ek9), and
the Pawsey Supercomputing Centre with funding from the Australian Government
and the Government of Western Australia, in the framework of the National
Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme. The
simulation software flash was in part developed by the DOE-supported Flash
Centre for Computational Science at the University of Chicago.

## References

  * André et al. (2010) André P., et al., 2010, A&A, 518, L102
  * André et al. (2019) André P., Arzoumanian D., Könyves V., Shimajiri Y., Palmeirim P., 2019, arXiv e-prints, p. arXiv:1907.13448
  * Appleby et al. (2018) Appleby S., Chingangbam P., Park C., Hong S. E., Kim J., Ganesan V., 2018, ApJ, 858, 87
  * Arzoumanian et al. (2018) Arzoumanian D., Shimajiri Y., Inutsuka S.-i., Inoue T., Tachihara K., 2018, Publications of the Astronomical Society of Japan, 70
  * Beattie & Federrath (2019) Beattie J. R., Federrath C., 2019, MNRAS, submitted
  * Beattie et al. (2019a) Beattie J. R., Federrath C., Klessen R. S., 2019a, MNRAS, 487, 2070
  * Beattie et al. (2019b) Beattie J. R., Federrath C., Klessen R. S., Schneider N., 2019b, MNRAS, 488, 2493
  * Beaumont et al. (2013) Beaumont C. N., Offner S. S. R., Shetty R., Glover S. C. O., Goodman A. A., 2013, ApJ, 777, 173
  * Bigot et al. (2008) Bigot B., Galtier S., Politano H., 2008, Phys. Rev. E, 78, 066301
  * Boldyrev (2006) Boldyrev S., 2006, Phys. Rev. Lett., 96, 115002
  * Brunt & Federrath (2014) Brunt C. M., Federrath C., 2014, MNRAS, 442, 1451
  * 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 et al. (2009) Burkhart B., Falceta-Gonçalves D., Kowal G., Lazarian A., 2009, ApJ, 693, 250
  * Burkhart et al. (2014) Burkhart B., Lazarian A., Leão I. C., de Medeiros J. R., Esquivel A., 2014, ApJ, 790, 130
  * Burkhart et al. (2015) Burkhart B., Lee M.-Y., Murray C. E., Stanimirović S., 2015, ApJ, 811, L28
  * Cho & Lazarian (2003) Cho J., Lazarian A., 2003, MNRAS, 345, 325
  * Cho & Vishniac (2000) Cho J., Vishniac E. T., 2000, ApJ, 539, 273
  * Cho et al. (2002) Cho J., Lazarian A., Vishniac E. T., 2002, ApJ, 564, 291
  * Cox et al. (2016) Cox N. L. J., et al., 2016, A&A, 590, A110
  * Crutcher (2012) Crutcher R. M., 2012, Annual Review of Astronomy and Astrophysics, 50, 29
  * Donovan Meyer et al. (2013) Donovan Meyer J., et al., 2013, ApJ, 772, 107
  * Dubey et al. (2008) Dubey A., et al., 2008, in Pogorelov N. V., Audit E., Zank G. P., eds, Astronomical Society of the Pacific Conference Series Vol. 385, Numerical Modeling of Space Plasma Flows. p. 145
  * Elmegreen & Falgarone (1996) Elmegreen B. G., Falgarone E., 1996, ApJ, 471, 816
  * Elmegreen & Scalo (2004) Elmegreen B. G., Scalo J., 2004, ARA&A, 42, 211
  * Esquivel & Lazarian (2011) Esquivel A., Lazarian A., 2011, ApJ, 740, 117
  * Federrath (2015) Federrath C., 2015, MNRAS, 450, 4035
  * Federrath (2016) Federrath C., 2016, 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. (2008) Federrath C., Klessen R. S., Schmidt W., 2008, ApJ, 688, L79
  * Federrath et al. (2009) Federrath C., Klessen R. S., Schmidt W., 2009, ApJ, 692, 364
  * Federrath et al. (2010) Federrath C., Roman-Duval J., Klessen R., Schmidt W., Mac Low M. M., 2010, A&A, 512
  * Federrath et al. (2016) Federrath C., et al., 2016, ApJ, 832, 143
  * Fitzgibbon et al. (1996) Fitzgibbon A. W., Pilu M., Fisher R. B., 1996, in Proceedings of 13th International Conference on Pattern Recognition. pp 253–257 vol.1, doi:10.1109/ICPR.1996.546029
  * Frank et al. (2014) Frank A., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 451 (arXiv:1402.3553), doi:10.2458/azu_uapress_9780816531240-ch020
  * Fryxell et al. (2000) Fryxell B., et al., 2000, ApJS, 131, 273
  * Gerrard et al. (2019) Gerrard I. A., Federrath C., Kuruwita R., 2019, MNRAS, 485, 5532
  * Ginsburg et al. (2013) Ginsburg A., Federrath C., Darling J., 2013, ApJ, 779
  * Goldreich & Sridhar (1995) Goldreich P., Sridhar S., 1995, ApJ, 438, 763
  * Henderson et al. (2019) Henderson R., Makarenko I., Bushby P., Fletcher A., Shukurov A., 2019, Journal of the American Statistical Association, 0, 1
  * Hennebelle & Inutsuka (2019) Hennebelle P., Inutsuka S.-i., 2019, Frontiers in Astronomy and Space Sciences, 6, 5
  * Kainulainen et al. (2014) Kainulainen J., Federrath C., Henning T., 2014, Science, 344, 183
  * Kim & Ryu (2005) Kim J., Ryu D., 2005, ApJ, 630, L45
  * Kolmogorov (1941) Kolmogorov A. N., 1941, Doklady Akademii Nauk Sssr, 30, 301
  * Konstandin et al. (2016) Konstandin L., Schmidt W., Girichidis P., Peters T., Shetty R., Klessen R. S., 2016, MNRAS, 460, 4483
  * Kowal et al. (2007) Kowal G., Lazarian A., Beresnyak A., 2007, ApJ, 658, 423
  * Krumholz & Federrath (2019) Krumholz M. R., Federrath C., 2019, Frontiers in Astronomy and Space Sciences, 6, 7
  * Krumholz & McKee (2005) Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250
  * Kuruwita & Federrath (2019) Kuruwita R. L., Federrath C., 2019, Monthly Notices of the Royal Astronomical Society, 486, 3647
  * Larson (1981) Larson R. B., 1981, MNRAS, 194, 809
  * Lazarian & Cho (2004) Lazarian A., Cho J., 2004, Astrophysics and Space Science, 292, 29
  * Lazarian & Vishniac (1999) Lazarian A., Vishniac E. T., 1999, ApJ, 517, 700
  * Lazarian et al. (2012) Lazarian A., Esquivel A., Crutcher R., 2012, ApJ, 757, 154
  * Mac Low & Klessen (2004) Mac Low M. M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125
  * Malinen et al. (2016) Malinen J., et al., 2016, MNRAS, 460, 1934
  * Maron & Goldreich (2001) Maron J., Goldreich P., 2001, ApJ, 554, 1175
  * McKee & Ostriker (2007) McKee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
  * Men’shchikov et al. (2010) Men’shchikov A., et al., 2010, A&A, 518, L103
  * Molina et al. (2012) Molina F. Z., Glover S. C. O., Federrath C., Klessen R. S., 2012, MNRAS, 423, 2680
  * Mouschovias et al. (2006) Mouschovias T. C., Tassis K., Kunz M. W., 2006, ApJ, 646, 1043
  * Orkisz et al. (2017) Orkisz J. H., et al., 2017, A&A, 599, A99
  * Padoan & Nordlund (2011) Padoan P., Nordlund Å., 2011, ApJ, 730, 40
  * Padoan et al. (2014) Padoan P., Federrath C., Chabrier G., Evans II N. J., Johnstone D., Jørgensen J. K., McKee C. F., Nordlund Å., 2014, Protostars and Planets VI, pp 77–100
  * Perez & Boldyrev (2007) Perez J. C., Boldyrev S., 2007, ApJ, 672, L61
  * Planck Collaboration et al. (2016a) Planck Collaboration et al., 2016a, A&A, 586, A137
  * Planck Collaboration et al. (2016b) Planck Collaboration et al., 2016b, A&A, 586, A138
  * Price & Federrath (2010) Price D. J., Federrath C., 2010, MNRAS, 406, 1659
  * Roman-Duval et al. (2010) Roman-Duval J., Jackson J. M., Heyer M., Rathborne J., Simon R., 2010, ApJ, 723, 492
  * Roy et al. (2019) Roy A., et al., 2019, arXiv e-prints, p. arXiv:1903.12608
  * Sanchez et al. (2005) Sanchez N., Alfaro E. J., Perez E., 2005, ApJ, 625, 849
  * Scalo (1990) Scalo J., 1990, in Capuzzo-Dolcetta R., Chiosi C., di Fazio A., eds, Astrophysics and Space Science Library Vol. 162, Physical Processes in Fragmentation and Star Formation. pp 151–176, doi:10.1007/978-94-009-0605-1_12
  * Schneider et al. (2013) Schneider N., et al., 2013, ApJ, 766, L17
  * Schneider et al. (2015) Schneider N., et al., 2015, A&A, 575, A79
  * Soler (2019) Soler J. D., 2019, A&A, 629, A96
  * Soler et al. (2017) Soler J. D., et al., 2017, A&A, 603, A64
  * Stutzki et al. (1998) Stutzki J., Bensch F., Heithausen A., Ossenkopf V., Zielinsky M., 1998, A&A, 336, 697
  * Tennekes & Lumley (1972) Tennekes H., Lumley J., 1972, A first course in turbulence. Cambridge, Massachusetts, MIT Press
  * Tokuda et al. (2018) Tokuda K., et al., 2018, arXiv e-prints, p. arXiv:1811.04400
  * Trevino-Morales et al. (2019) Trevino-Morales S. P., et al., 2019, arXiv e-prints, p. arXiv:1907.03524
  * Tritsis & Tassis (2016) Tritsis A., Tassis K., 2016, MNRAS, 462, 3602
  * Tritsis & Tassis (2018) Tritsis A., Tassis K., 2018, Science, 360, 635
  * Tritsis et al. (2018) Tritsis A., Federrath C., Schneider N., Tassis K., 2018, MNRAS, 481, 5275
  * Verdini et al. (2015) Verdini A., Grappin R., Hellinger P., Landi S., MÃŒller W. C., 2015, ApJ, 804, 119
  * Xu et al. (2019) Xu S., Ji S., Lazarian A., 2019, arXiv e-prints, p. arXiv:1905.06341

##  Appendix A Inversion into Real-space

We can access anisotropic features from the 2D power spectrum and put them
back into real-space by convolving a top hat mask with the same elliptic
properties of the contours to the 2D Fourier transform of the density field,
and then inverting it. First we define the mask,

$\mathcal{E}_{\sigma}=\mathcal{N}(\sigma)*\mathcal{E}(e_{\parallel,1},e_{\perp,1},e_{\parallel,2},e_{\perp,2}),$
(30)

where, $\mathcal{E}$, is the annular region defined between two ellipses with
principle axes ($e_{\parallel,1},e_{\perp,1}$) and
($e_{\parallel,2},e_{\perp,2}$), $\mathcal{N}$ is a Gaussian smoothing
function, which smooths the annulus with kernel size $\sigma$ and $*$ is the
convolution operator. We smooth the annulus to avoid the Gibbs phenomenon at
the discontinuous edge of the annulus, which will otherwise be present after
taking the inverse Fourier transform. We choose
$\sigma=2\,\operatorname{d}\\!{x}$, where $\operatorname{d}\\!{x}$ is the size
of a grid-cell of the $\mathcal{E}$ mask. We choose this because it
sufficiently smooths the discontinuous edge of the mask, whilst minimising the
amount of overflow into neighbouring grid-cells (which will correspond to
neighbouring $k$-modes). Next we take the convolution with the 2D Fourier
transform, shown in Equation 14, of the column density, $\hat{\xi}$, where we
remind the reader, $\xi=\Sigma/\Sigma_{0}-1$ is the zero-mean transformed
column density, and invert the convolved Fourier image of the column density
into real-space, i.e.,

$\displaystyle\xi_{\mathcal{E}}$
$\displaystyle=\mathscr{F}^{-1}\left[\hat{\xi}*\mathcal{E}_{\sigma}\right]$
(31)
$\displaystyle=\displaystyle\sum^{N-1}_{k_{x}=0}\sum^{N-1}_{k_{y}=0}\left[\hat{\xi}*\mathcal{E}_{\sigma}\right]\exp\left\\{i\frac{\mathrm{{\mn<EMAIL_ADDRESS>(32)

Since the inverse Fourier transform is a complex-valued function we expect
there to be a real and imaginary component to $\xi_{\mathcal{E}}$. We take
just the real structures and plot them in the third row of Figure 3. In Figure
3 we use the annulus $18\leq\operatorname{k_{\perp}}\leq 58\cup
10\leq\operatorname{k_{\parallel}}\leq 27$ and
$18\leq\operatorname{k_{\perp}}\leq 62\cup
17\leq\operatorname{k_{\parallel}}\leq 84$, for the M2Ma0.1 and M20Ma0.1
simulation, respectively. These regions are extracted directly from the
elliptic fits, shown in Figure 2. This method reveals what kind of structures
in the column density maps are responsible for the anisotropic features in the
2D power spectra.

Figure 8: The anisotropy parameter, $a$, shown in Equation 16, as a function
of wave vector perpendicular to the guide field, $\operatorname{k_{\perp}}$
for the two extreme anisotropic simulations examined in §4, coloured by mean-
field Alvén Mach number, $\operatorname{\mathcal{M}_{\text{A}}}_{0}$. The
Goldreich & Sridhar (1995)(GS1995) scaling, Equation 1, is shown in black,
which is the anisotropy predicted in the energy spectrum for incompressible
MHD turbulence. We find that the spectra for the column density especially in
the $\operatorname{\mathcal{M}}=20$ simulation do not follow the GS1995
scaling.

##  Appendix B Anisotropy in the column density as a function of
$\operatorname{k_{\perp}}$

We show the anisotropy, $a$, as a function of scale,
$\operatorname{k_{\perp}}$ in Figure 8 for the $\operatorname{\mathcal{M}}=2$
(top plot) and $\operatorname{\mathcal{M}}=20$ (bottom plot) simulations,
varying from $\operatorname{\mathcal{M}_{\text{A}}}_{0}=0.1-10$. These are
illustrative of the curves that we use to average over in §6 to create the
single spatially-averaged statistic, $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}$, for each simulation. The
Goldreich & Sridhar (1995) (GS1995) scaling, shown in Equation 1, is
illustrated in black, at the left of each of the plots.

The GS1995 scaling is steeper than the observed decline in the
$\operatorname{\mathcal{M}}=2$ simulations and in the wrong direction for the
$\operatorname{\mathcal{M}}=20$ simulations. This result is not entirely
surprising, since the critical balance in GS1995 is concerned with the energy
spectrum of incompressible MHD turbulence, and in this study we are concerned
with the power spectrum of the column density in highly-compressible MHD
turbulence – two very different quantities.

##  Appendix C Anisotropy in the Column Density as a function of LOS
projection Angles

We explore the anisotropy as a function of line-of-sight (LOS) viewing angle,
$\phi$, with respect to the perpendicular viewing angle to the magnetic guild
field, $\mathrm{{\mn@boldsymbol{\mathit{B}}}}=B_{z}\hat{z}$, as explored in
the main text. We call this $\phi=0$. We focus our analysis on the most
anisotropic simulations, M2Ma0.1 and M20Ma0.1. We show the projections and 2D
power spectra for $\phi=0,\,\pi/6,\,\pi/3$ and $\pi/2$, in Figure 9, following
the methods from §4, where $\phi=\pi/2$ is the projection that runs along the
guide field.
Figure 9 reveals that the power spectrum becomes more isotropic as our viewing
angle moves between the perpendicular and parallel viewing angles with respect
to the guide field. We calculate the average anisotropy factor, as we did in
§6 and plot it as a function of $\phi$ in Figure 10. We find that, for the
highly-supersonic simulation, M20Ma0.1 there is a smooth transition towards
isotropy as the viewing angle coincides with the direction of the magnetic
field. However, for the M2Ma0.1 there is an abrupt change into the isotropic
regime only when the viewing angle is directly aligned with the magnetic
field. This means that the anisotropy caused by the striations is more robust
to changes in th viewing angle than the anisotropy caused by the perpendicular
filaments.

Figure 9: The column density maps and 2D power spectra for the M2Ma0.1 (top
two rows) and M20Ma0.1 (bottom two rows) simulations shown at line-of-sight
(LOS) projection angles $\phi=0,\,\pi/6,\,\pi/3$ and $\pi/2$, from left to
right, respectively. $\phi=0$ corresponds to the LOS perpendicular to the
magnetic guide field and $\phi=\pi/2$ the LOS parallel to the guide field. As
the LOS moves towards the parallel direction we find that the anisotropies
discussed in the main text disappear, revealing approximately isotropic
behaviour around the magnetic guide field in the column density. Figure 10:
The average anisotropy, $\mathrm{\left\langle
a\right\rangle}_{\operatorname{k_{\perp}}}$, as described in §6, as a function
of line-of-sight (LOS) viewing angle, $\phi$, where $\phi=0$ and $\phi=\pi/2$
correspond to the LOS viewing angle perpendicular and parallel to the magnetic
guide field, respectively. The red line indicates perfect isotropy. The top
panel is for the M2Ma0.1 simulation and the bottom for M20Ma0.1. We show
examples of the column density maps for each $\phi$ in Figure 9.

##  Appendix D 3D density visualisations

Figure 11: The gas density shown in 3D, for the two simulations, M2Ma0.1 and
M20Ma0.1, exhibiting extreme anisotropy along the direction of the magnetic
field lines. Both 3D densities reveal that striations are most likely long,
thin, sheets that twist around the magnetic field lines due to the turbulent
motions. Furthermore, in the M20Ma0.1 simulation we see high-density filaments
forming parallel to the magnetic fields.