Magnetic field fluctuations in anisotropic, supersonic turbulence

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

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

Turbulent magnetic fields can be separated into two components,

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

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

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

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

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

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

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

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

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

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

For sub-Alfvénic molecular clouds, such as G0.253+0.016, $B_{0}\gg\delta B$ (in this case, $B_{0}>\delta B$ by an order of magnitude; Federrath et al. 2016b), and the time derivative of $\delta B$ is very small, i.e., much smaller than the time scale of compression. Hence, by assuming $B_{0}\gg\delta B$ and $(D/Dt)\delta B\approx 0$,¹¹1$\frac{D}{Dt}\equiv\frac{\partial}{\partial t}+\mathrm{{\mn@boldsymbol{\mathit{v}}}}\cdot\nabla$ is the material derivative. Equation (2.1) 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²²2Note that these are in addition to the mixture of various modes caused by the turbulent driving (Federrath et al., 2008, 2010).. We show the full derivation of Equation (12) in Appendix B. This is consistent with other studies, where hydrodynamical shocks appear perpendicular to the mean magnetic field, with the compression happening in the parallel direction to the mean magnetic field (Mocz & Burkhart, 2018; Beattie & Federrath, 2020b; Seifried et al., 2020). We will discuss this further in §4.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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