key: cord-0323423-ds76rjkr
authors: Beattie, James R.; Federrath, Christoph; Seta, Amit
title: Magnetic field fluctuations in anisotropic, supersonic turbulence
date: 2020-07-28
journal: nan
DOI: 10.1093/mnras/staa2257
sha: 55420a1b96ffdb7243abf8afef20fae065d353dd
doc_id: 323423
cord_uid: ds76rjkr

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

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

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

where V , VA = B/ √ 4πρ, ρ and B are the root-mean-squared (rms) velocity, Alfvén velocity, density and magnetic field strength, respectively. For MA < 1, ρV 2 B 2 the magnetic field contributes significantly to the dynamics of the flow through the Lorentz force. This is the sub-Alfvénic regime. For MA > 1, ρV 2 B 2 the magnetic field plays a lesser role and the turbulent motions set, for example the statistics of the flow (Padoan & Nordlund 2011; Molina et al. 2012; Beattie & Federrath 2020b) . This is called the super-Alfvénic regime, and the trans-Alfvénic regime, MA ∼ 1, separates the two. In our study we explore a broad range of MA, from 0.1 − 100, to include the two dynamically dissimilar regimes.

Turbulent magnetic fields can be separated into two components,

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

where B(t) t = |B0| which means δB(t) t = 0. For this reason, we will refer to B0 as the mean-field throughout this study.

Since the magnetic field has two distinct components one can talk about the mean-field or fluctuating Alfvén Mach number. In this study, we find the dependence of the magnetic field fluctuations, δB, on the Alfvén Mach number defined with respect to the mean-field. We define the meanfield Alfvén Mach number exactly as we defined the rms MA in Equation (1), but with mean-field components, MA0 = V /VA0 = ( √ 4πρ0V )/B0, where B0 = |B0| and VA0 is the velocity of Alfvén waves along the mean-field. We explore a wide parameter set that encompasses both the sub-and super-Alfvénic mean-field regime, covering the parameter space of observed molecular clouds. For example, for the central molecular zone cloud, G0.253+0.016 studied in , they measure a mean-field of B0 = (2.07 ± 0.95) mG (Pillai et al. 2015) , a mean volume density of ρ0 = (6.2 ± 3.3) × 10 −20 g/cm 3 and a velocity dispersion of V = (6.8 ± 0.2) km/s (see for details on assumptions and further references). Using these quantities and Equation (1) one can calculate VA0 = (24 ± 17) km/s which means that MA0 = 0.3 ± 0.2, placing it well into the sub-Alfvénic mean-field regime. Furthermore, using the direction of the velocity gradients to infer the magnetic field strength, Hu et al. (2019) measure the Alfvén Mach number for five star-forming clouds in the Gould Belt: Taurus (MA0 = 1.19 ± 0.02), Perseus A (1.22 ± 0.05), L 1551 (0.73 ± 0.13), Serpens (0.98 ± 0.08) and NGC 1333 (0.82 ± 0.24). Hence, the trans to sub-Alfvénic magnetised turbulence regime is of great importance for understanding the environment that stars form in, which is the focus of this study.

Previous work by Federrath (2016a, herein called F16) explored the magnitude of magnetic field fluctuations across a broad range of B0 and δB values, but for a single M value. In our study we will expand upon some of the key results in F16 and elucidate how the M number influences the fluctuating magnetic field. Of particular importance for our study is the F16 analytical model for δB in the strongfield regime, where B0 δB. By relating the turbulent magnetic energy density, em ≈ B0δB/(4π) and the turbulent kinetic energy e k = (ρ0V 2 )/2, F16 found that

by assuming that all of the turbulent magnetic energy was fed from the kinetic energy in the plasma, i.e., em = e k . We expand upon this study significantly by exploring the fluctuations perpendicular and parallel to B0 separately in high-resolution MHD simulations, by showing the mean normalised fluctuations are independent of M, and by deriving a new model directly from the induction equation in the compressible, ideal MHD model of plasmas. Hence our study contributes significantly to better understanding the nature of magnetic fluctuations in compressible, astrophysically relevant flows.

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

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

where v is the fluid velocity, ρ the density, B the magnetic field, cs the sound speed and F , a stochastic driving function that gives rise to the turbulence (which could be from, for example, supernova shocks permeating through the interstellar medium, or internal to the MC, gravity, galactic-scale shocks or ambient pressure from the galactic environment; Brunt et al. 2009; Elmegreen 2009; Federrath 2015; Krumholz & Burkhart 2016; Grisdale et al. 2017; Jin et al. 2017; Körtgen et al. 2017; Federrath et al. 2017; Colling et al. 2018; Schruba et al. 2019; Lu et al. 2020 ). We consider a magnetic field of the form shown in Equation (3), with constant mean component of the field in the z-direction,

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

which reveals that in the Lagrangian frame of the fluid, the magnetic field fluctuations are influenced by (1) the velocity gradients along the mean magnetic field (which we indicate with ∂ , i.e., they are anisotropic), (2) the advection of the velocity by the fluctuations, and (3) the compression of the velocity field, assuming B0 is independent of time. Disregarding the compression term, (B0 + δB)(∇ · v ), and as ∂tδB → 0, this form of the induction equation is known as the quasi-static approximation of the magnetic field perturbation, which is only applicable in the case of incompressible plasma, (Zikanov & Thess 1998; Verma 2017 , and references therein), but here we have generalised the equation for compressible MHD plasmas, relevant for astrophysical phenomena.

For sub-Alfvénic molecular clouds, such as G0.253+0.016, B0 δB (in this case, B0 > δB by an order of magnitude; , and the time derivative of δB is very small, i.e., much smaller than the time scale of compression. Hence, by assuming B0 δB and (D/Dt)δB ≈ 0, 1 Equation (10) simply becomes

and by taking the magnitudes of both sides,

Thus, it is primarily the velocity streams that run parallel with the mean magnetic field that give rise to compressive modes in the clouds 2 . We show the full derivation of Equation (12) in Appendix B. This is consistent with other studies, where hydrodynamical shocks appear perpendicular to the mean magnetic field, with the compression happening in the parallel direction to the mean magnetic field (Mocz & Burkhart 2018; Beattie & Federrath 2020b; Seifried et al. 2020 ). We will discuss this further in §4. Equation (12) shows the key difference between incompressible and compressible MHD turbulence for the case when B0

δB. If, for example, ∇·v = 0, the velocity gradient along the field disappears, and hence the magnetic field fluctuations cannot persist along the field. This is known as the two-dimensionalisation of the three-dimensional flow (Alexakis 2011; Verma 2017) , which for supersonic turbulence, is forbidden to happen as shocks form and travel parallel to the field from the parallel velocity gradient. This means magnetic field fluctuations persist along the field, and the flow retains its three-dimensional nature. We will see this is indeed the case for parallel fluctuations in numerical simulations (demonstrated in §5).

1 D Dt ≡ ∂ ∂t + v · ∇ is the material derivative. 2 Note that these are in addition to the mixture of various modes caused by the turbulent driving (Federrath et al. 2008; Federrath et al. 2010 ).

We now consider the dimensionless form of Equation (10), where length scales,ˆ = /L, velocitiesv = v /V , and magnetic fields, δB = δB/δB,B0 = B0/B0, are normalised by their magnitudes, and the time-scale in the derivative we chose to bet = t/T , where T = L/VA0, is the time-scale of fluctuations travelling along the mean magnetic field that spans the full system scale L. By grouping δB and B0 terms on the right hand side of Equation (10) and by applying our non-dimensionalisation, we find,

If we assume that∂ v −B0(∇ ·v ) and (δB ·∇)v − δB(∇ ·v ) are O(1), which is true by construction, then we can make an order-of-magnitude estimate of the fluctuations,

VA0

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

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

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

The initial magnetic field, B(x, y, z) at t = 0, in Equations (6-8) is a uniform field with field lines threaded through theẑ direction of the simulations. The total magnetic field is, as given in Equations 3 and 9 and reiterated here for clarity,

Since we work in units of sound speed for velocity, and in units of ρ0 for density, we chose to scale the units of the magnetic field for a typical molecular cloud density (hydrogen number density of nH = 10 3 cm −3 ) and temperature (10 K) corresponding to a sound speed of cs = 0.2 km s −1 , to

where Bsim is the magnetic field in simulation units, such that all magnetic fields are in units of µG, as previously done in Mocz & Burkhart (2018) , for example. We list all of the magnetic field components in Table 1 . Note that in our simulations 0.01 B0/(µG) 1000, providing us with ∼ five orders of magnitude to explore for the value of B0 (column 4 in Table 1 ). This parameter space encompasses the wide variety of molecular clouds that one might observe. The fluctuating component of the field evolves selfconsistently from the MHD equations. However, we set B0, and by using the definition of the Alfvén velocity (defined previously in reference to Equation 1) and the turbulent Mach number, we can set the target MA0,

We vary this value for each of the 24 simulations, spanning MA0 = 0.1 − 100, encompassing very weak and strong mean magnetic fields. We extract 51 time realisations of the Bx, By and Bz components of the total magnetic field between 5 ≤ t/T ≤ 10 to ensure that the turbulence is in a statistically stationary state (Federrath et al. 2009; Price & Federrath 2010 ).

We show a slice of the magnitude of δB in units of B0, for each of the 24 simulations at t = 5 T , in Figure 1 . Red colours in the plot show |δB| < B0, white |δB| = B0 and black |δB| > B0. In the sub-Alfvénic regime we see |δB| < B0, in the super-Alfvénic regime |δB| > B0 and MA0 ≈ 1 marks the transition between the two. Using the magnetic field components we construct the variables B = Bz and B ⊥ = Bx = By 3 . We create PDFs for the two new magnetic field variables, and average them over 5−10 T . In §5 we discuss these PDFs in detail, but first we will discuss one of the main predictions from our fluctuation model in Equation (10), namely that it is primarily the velocity gradients parallel to the mean-field that create compressive modes in the velocity field.

In §2 we showed how one could construct a compressive quasi-static model for the magnetic field fluctuations (Equation 10) and that in the strong mean-field regime the model reduces to a simple statement about compressibility and velocity gradients along B0, shown in Equation (12). This is a significant result for astrophysical flows, since both strong, coherent mean magnetic fields and velocity gradients perpendicular to filamentary structures are measured in observations (Palmeirim et al. 2013; Beuther et al. 2015; Shimajiri et al. 2019; Chen et al. 2020b ). Here we show some direct, qualitative evidence to support Equation (12). We show a single time realisation of the full 3D geometry of the magnetic and velocity (scaled by cs) fields in Figure 2 for the M20MA10 (top) and M20MA0.1 (bottom) simulations, where we refer to Table 1 for the respective labels. M20MA10 shows an example of the tangled, isotropic magnetic fields that we have in the high MA0 simulations. The magnetic field is dominated by δB, because the mean-field component is too weak to impart any systematic, ordered structure on either the velocity or total magnetic field. On the other hand, the M20MA0.1 simulation shows a highlyordered, anisotropic magnetic field, dominated by B0. Many of the velocity streamlines form ∼ L/2 eddies perpendicular to the magnetic field lines. If all of the streamlines were organised into eddies the flow would be statistically the same along B0, hence becoming quasi two-dimensional. However, we see some parallel velocity streams developing towards the left corner of the plot. To explore this further we take an xz-slice through this same time realisation.

We plot the convergence of the xz slice of the velocity field for the M20MA0.1 simulation in the top-left panel of Figure 3 . Hydrodynamical shocks that form perpendicular to the mean-field (see Mocz & Burkhart 2018 for the shock-jump conditions, showing that they are hydrodynamical) can be clearly identified as −∇ · v > 0 structures, i.e., converging flows, by comparing the convergence with the density structures in the bottom-left panel of Figure 3 . Tracing the velocity streamlines, shown in the bottom-right panel, we can see the most significant shocks are coupled with a strong velocity channels along the mean magnetic field. Furthermore, probing the magnitude of the velocity gradient along B0, shown in the top-left panel, shows a strong correlation between the −∇ · v > 0 structures and the velocity gradient, consistent with Equation (12). The only difference between the panels is the ∇ · v signatures of fast magnetosonic waves that form perpendicular to the field in the convergence plot. However, these waves will cease as the mean-field becomes strong and magnetic tension locks them in place and then |∇ · v | will indeed equal |∂ v |.

To summarise, what we show here is that as B0 δB parallel velocity gradients cause compressive modes in the velocity field, which cause strong hydrodynamic shocks to form across the field.

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

To conclude this section we comment on the velocity The magnitude of the velocity gradient along the magnetic field, |∂ v |. Equation (12) predicts that this plot and the magnitude of the convergence should be the same as ∂tδB → 0. This is indeed true, and the only significant difference between the two plots is the fast magnetosonic waves, which will disappear as B 0 δB. Bottom-left: The same as the top-left plot, but for a density slice, shown in units of mean density, ρ 0 . We see that the −∇·v > 0 structures in the top panel correspond to hydrodynamical shocks (Mocz & Burkhart 2018; Beattie & Federrath 2020b) , forming from mostly parallel velocity streams in the flow. Bottom-right: Velocity streamlines integrated through the divergence of the velocity, showing how parallel velocity channels feed into high-density filaments, creating the velocity gradients parallel to the magnetic field. gradient method, which is summarised in Hu et al. (2019) , and used to measure the MA values that we report for the Gould Belt in §1. This method has been significantly developed in the literature (González-Casanova & Lazarian 2017; Yuen & Lazarian 2017b,a; Lazarian et al. 2018 , and references therein) and assumes that velocity gradients are perpendicular to magnetic fields, which may be true for incompressible MA0 < 1 flows (i.e., through the twodimensionalisation due to the exponentially decaying parallel magnetic field fluctuations; Verma 2017), but is not be strictly true for compressible flows, as indicated in Equa-tion (12) and the discussion above. Yuen & Lazarian (2017a) attribute co-alignment of velocity and magnetic field vectors to indicate gravitational collapse. We show that this need not be the case, and for supersonic flows in a sub-Alfvénic mean-field regime, one can indeed find strong velocity gradients parallel to the mean magnetic field, without the presence of self-gravity. Next we discuss the magnetic field PDFs, and how the magnetic field fluctuations depend upon M and MA0. 

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

We construct the PDFs for both the B and B ⊥ components of the magnetic field and average them over the time range 5 − 10 T . The first moment and second moment of the PDFs have a physical interpretation. For example, for the B PDF the first moment is,

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

which means the standard deviation of the field is exactly δB . Since the mean-field only has a z component, B ⊥ = 0, and hence the second moment of the PDF is exactly equal to the fluctuations squared,

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

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

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

Using the definition of the Alfvén Mach number, MA0 = V /VA0, the turbulent Mach number, M = V /cs, and the Alfvén velocity, VA0 = B0/ √ 4πρ0 one can show that the strong-field F16 model (Equation 4) can be rewritten in terms of M,

This modified version of the F16 model shows a linear dependence in M, in agreement with our model derived from the induction equation, and linear in MA0. Hence, the C obtained is analytical in this regime, C = MA0 /2. This means, as suggested in the derivation of Equation (17), that C can be thought of as the Alfvénic control parameter, i.e., it changes based on the type of Alfvénic turbulence, which is shown clearly in the inset plots of Figure 4 . We plot Equation (25) on the magnitudes of the perpendicular, parallel and total fluctuations of the magnetic field in Figure 5 , where we show only data from the transto sub-Alfvénic regimes, since this is the regime where the strong-field model is valid. In general, we find that the model captures the nature of the fluctuations well, tracking either values between B ⊥ and B in the sub-Alfvénic simulations, or very close to the total fluctuations in the trans-Alfvénic simulation, MA0 ∼ 1. Plotting both perpendicular and parallel fluctuations in Figure 5 reveals the difference between their magnitudes due to the anisotropy inherent to MHD turbulence with a mean-field present (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999; Cho & Vishniac 2000; Cho & Lazarian 2003; Boldyrev 2006; Kowal et al. 2007; Burkhart et al. 2014 Burkhart et al. , 2015 . In the following, we explore the anisotropy of the fluctuations.

In Figure 6 we show the ratio between the B and B ⊥ fluctuations. Since δB ⊥ = δB 2 x + δB 2 y when δBx = δBy = δB the ratio δB /δB ⊥ = 1/ √ 2. For MA0 2 we see the ratio monotonically approaching the 1/ √ 2 isotropic limit. We can understand this from our model in Equation (10). When δB B0 we can disregard the O(B0) terms in Equation (10),

where D Dt ≡ ∂ ∂t + v · ∇ , which no longer has the anisotropic term, B0∂ v . This means that each of the components of δB are affected by the advection and compression equally, causing isotropic fluctuations, without any preferential direction (which we visualised in the top panel of Figure 2) . However, in the opposite case, where B0 δB, the B0∂ v becomes large and we measure δB up to a factor 2 larger than δB ⊥ . This is because the B0∂ v from our quasistatic model causes strong compressions across the field, which increases the magnetic field fluctuations along the mean-field direction, as shown in Figure 3 , and discussed in detail in §2 and §4. The impact that these shocks have on the parallel fluctuations is quite phenomenal, and we see a significant growth in the absolute value of the anisotropy between MA0 = 0.5, which has stronger B ⊥ fluctuations, and MA0 = 0.1.

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

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

One of the key features of Figure 7 is the distinct kink at δB ≈ B0 (indicated by the horizontal blue line in Figure  7) , which happens at MA0 ≈ 2. The same transition can be seen in Figure 1 , where fluctuations are O(B0), and are becoming more isotropic in nature, which is shown in Figure  6 . Clearly this marks a critical transition for the magnetic field evolution, where all of the terms in our quasi-static model play a role, so very little can be deduced using our model in this transition region. Instead, we propose a semianalytical model that encapsulates both the δB/B0 < 1 and δB/B0 > 1 regimes,

where α and β are fit parameters. For δB/B0 < 1 we use the F16 model, but since no analytical model is known in the δB/B0 ≥ 1 regime, and our compressible quasi-static model does not give much insight, we instead directly fit to the data. Using a linear least-squares fit we determine α ≈ 1 and 1 + β = 0.69 ± 0.05. This is somewhat consistent with the qualitative model for the "intermediate regime" in F16

, but here we explicitly show that the mean-field weighted fluctuations are independent of M, regardless of the ratio between δB and B0. Hence the relation is,

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

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

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

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

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

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

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

The remaining structure in the distribution is then most likely due to the turbulent fluctuations, which we showed in §5.1.1 are linear in M. Hence for large M values we mostly see the velocity interactions, and the interaction between δB ⊥ and δB no longer dominates the morphology of the PDFs.

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

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

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

where cs is the sound speed, ρ0 is the mean density and C is the proportionality factor, which changes between the B0 δB and B0 δB regimes. We show in Figure 4 that the linear dependence predicted for all MA0 holds true.

• We plot the velocity divergence and velocity streamlines for one of the MA0 = 0.1 simulations in Figure 3 , revealing the compressive modes coupled with parallel velocity gradients, |∇ · v | = |∂ v |. We show that, even in the absence of self-gravity, shocks that form perpendicular to the mean magnetic field turn into converging flows that are accreted upon by coherent parallel velocity channels, consistent with our compressible quasi-static model.

• Our new fluctuation model reduces to the Federrath (2016a) analytical model for strong-field magnetic field fluctuations by setting C = MA0 /2. We further rewrite the Federrath (2016a) model to incorporate the explicit dependence upon M, δB = cs √ πρ0 MA0 M, and we show in Figure 5 that it is in good agreement with the measured fluctuations for the sub-Alfvénic simulations, across all M.

• We find the ratio between the δB ⊥ and δB components of the magnetic field are similar for different M, but do depend upon MA0 (shown in Figure 6 ). As MA0 gets larger the fluctuations become isotropic, but for MA0 ∼ 0.1 δB ∼ 2δB ⊥ . This is because the strong velocity gradients along the field, when B0 δB, leads to large values of δB . This is a key distinguishing feature from incompressible MHD turbulence with a strong mean magnetic field, where parallel fluctuations exponentially decay and the flow becomes quasi two-dimensional (Alexakis 2011; Verma 2017) . We show that the flow remains very much three-dimensional for compressible MHD turbulence.

• We show that the mean-field normalised fluctuations, δB/B0, are independent of M in Figure 7 . We propose a semi-analytic model using the Federrath (2016a) model in the strong-field regime and a least-squares fit in the weakfield regime, δB B0 = MA0 2 /2, δB/B0 < 1, M 0.69±0.05

, δB/B0 ≥ 1, which is valid until the turbulence reaches the dynamo growth regime, (δB/B0) 1, and independent of M, making it a universal feature of anisotropic, supersonic MHD turbulence.

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

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

The B0/B0 term is just the unit vector of B0,B0, which contributes nothing to the magnitude, but sets direction of the compressive modes. Hence the magnitude of the compression is equal to the magnitude of the velocity gradient along the magnetic field,

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

APPENDIX C: WEAK M DEPENDENCE OF MEAN-FIELD WEIGHTED FLUCTUATIONS Figure C1 shows how the perpendicular (left) and parallel (right) magnetic field fluctuations only very weakly depend upon M. We discuss this result in §5.3. Figure C2 shows the δB /B0 (triangles) and δB ⊥ /B0 fluctuations. We show that they follow a similar trend as the magnitude of the full 3D fluctuations, |δB|/B0, as discussed in §5.3. Figure C2 . The same as Figure 7 , but for δB /B 0 (triangles) and δB ⊥ /B 0 (circles). We see both of the magnetic field components follow a similar trend as |δB|/B 0 , as discussed in §5.3.

High Performance Visualization-Enabling Extreme-Scale Scientific Insight

Numerical Modeling of Space Plasma Flows

The Galaxy Disk in Cosmological Context

The Multi-Messenger Astrophysics of the Galactic Centre

Protostars and Planets VI. p

Spectrum and energy transfer in steady burgers turbulence

Monthly Notices of the Royal Astronomical Society

NumPy: A guide to NumPy

Commmunications of the Konkoly Observatory Hungary

Protostars and Planets VI

We thank the reviewer, Pierre Lesaffre, for providing a thorough and constructive review of our study. We thank the Australian National University and Research School of Astronomy and Astrophysics for the support in place for staff and students during the COVID19 isolation, which

is when most of this study was written. J. R. B. acknowledges funding from the Australian National University, specifically the Deakin PhD and Dean's Higher Degree Research (theoretical physics) Scholarships. C. F. acknowledges funding provided by the Australian Research Council (Discovery Project DP170100603, and Future Fellowship FT180100495) , and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). We further acknowledge high-performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Partnership for Advanced Computing in Europe (PRACE grant pr89mu), the Australian National Computational Infrastructure (grant ek9), and the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, in the framework of the National Computational Merit Allocation Scheme and the ANU Allocation Scheme. The simulation software FLASH was in part developed by the DOEsupported Flash Centre for Computational Science at the University of Chicago.

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

but with a magnetic field decomposition that we show in Equation (9) it is possible to simplify the equation, and gain some physical intuition for the fluctuating component of B.First we write B(t) = B0 + δB(t),since ∂tB0 = 0. Expanding the cross product twice we find,Next we use the identity ∇ × (A × B) = A(∇ · B) − B(∇ · A) + (B · ∇)A − (A · ∇)B to expand each of the two terms. The first term is,which we can simplify because the magnetic field is divergence free, and because ∂x i B0 = 0, hence,Since B0 = B0ẑ the second term can be written in terms of just the velocity gradient in the z direction,Now we turn our attention to the fluctuating term in Equation (A3). Using the same identity as we used for the meanfield term and immediately dropping the ∇ · δB term we find,Now combing the two Equation (A6) and (A7) to construct Equation (A3),which can be simplified to reveal,or more succinctly, using D Dt = ∂t + v · ∇ as the Lagrangian derivative, i.e. the derivative in the frame co-moving with the fluid,The first term, B0∂zv , (we use the notation B0∂ v in the main text) tells us that the fluctuations change when the velocity gradient along the mean magnetic field changes. This term encodes the anisotropy of the fluctuations into the induction equation. Thus, when B0 increases (or equivalently MA0 decreases) the anisotropy in fluctuations increases and the fluctuations become more isotropic. This is also demonstrated via results from our simulations in Figure 6 . The second term, B(∇ · v ), tells us that the change in the fluctuations scale with the compression of the velocity field weighted by the magnetic field. The last term, (δB ·∇)v is the advection of the velocity from the magnetic field fluctuations.

For B0 δB Equation (10) simplifies significantly. The time derivative of δB, on the LHS, is approximately zero, and so are the terms that have strict δB dependence on the RHS. Hence