A multi-shock model for the density variance of anisotropic, highly-magnetised, supersonic turbulence
James R. Beattie, Philip Mocz, Christoph Federrath, Ralf S. Klessen
MMNRAS , 1–16 (2020) Preprint 2 February 2021 Compiled using MNRAS L A TEX style file v3.0
A multi-shock model for the density variance of anisotropic,highly-magnetised, supersonic turbulence
James R. Beattie (cid:63) , Philip Mocz , † , Christoph Federrath , andRalf S. Klessen , Research School of Astronomy and Astrophysics, Australian National University, Canberra, ACT 2611, Australia Australian Research Council Centre of Excellence in All Sky Astrophysics (ASTRO3D), Canberra, ACT 2611, Australia Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA Universit¨at Heidelberg, Zentrum f¨ur Astronomie, Institut f¨ur Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany Universit¨at Heidelberg, Interdisziplin¨ares Zentrum f¨ur Wissenschaftliches Rechnen, Im Neuenheimer Feld 205, 69120 Heidelberg, Germany † Einstein Fellow
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Shocks form the basis of our understanding for the density and velocity statistics ofsupersonic turbulent flows, such as those found in the cool interstellar medium (ISM).The variance of the density field, σ ρ/ρ , is of particular interest for molecular clouds(MCs), the birthplaces of stars in the Universe. The density variance may be used toinfer underlying physical processes in an MC, and parameterises the star formation(SF) rate of a cloud. However, models for σ ρ/ρ all share a common feature – thevariance is assumed to be isotropic. This assumption does not hold when a trans/sub-Alfv´enic mean magnetic field, B , is present in the cloud, which observations suggestis relevant for some MCs. We develop an anisotropic model for σ ρ/ρ , using contribu-tions from hydrodynamical and fast magnetosonic shocks that propagate orthogonalto each other. Our model predicts an upper bound for σ ρ/ρ in the high Mach number( M ) limit as small-scale density fluctuations become suppressed by the strong B .The model reduces to the isotropic σ ρ/ρ − M relation in the hydrodynamical limit.To validate our model, we calculate σ ρ/ρ from 12 high-resolution, three-dimensional,supersonic, sub-Alfv´enic magnetohydrodynamical (MHD) turbulence simulations andfind good agreement with our theory. We discuss how the two MHD shocks may bethe bimodally oriented over-densities observed in some MCs and the implications forSF theory in the presence of a sub-Alfv´enic B . By creating an anisotropic, super-sonic density fluctuation model, this study paves the way for SF theory in the highlyanisotropic regime of interstellar turbulence. Key words:
MHD – turbulence – ISM: kinematics and dynamics – ISM: magneticfields – ISM: structure
Shocks are the fundamental building blocks of supersonicturbulence and are key to understanding both the densityand velocity statistics of these flows which are ubiquitousin many astrophysical phenomena (Burgers 1948; Vazquez-Semadeni 1994; Padoan et al. 1997; Klessen et al. 2000;Padoan & Nordlund 2011; Federrath 2013; Lehmann et al.2016; Robertson & Goldreich 2018; Mocz & Burkhart 2018;Park & Ryu 2019; Abe et al. 2020; Federrath et al. 2021). (cid:63)
E-mail: [email protected]
Density statistics are of particular interest for understand-ing the structure and physical processes that shape and gov-ern the dynamics of the compressible, supersonic interstel-lar medium (ISM; Mac Low & Klessen 2004; Krumholz &McKee 2005; Federrath et al. 2008; Federrath et al. 2009;Burkhart et al. 2009; Padoan & Nordlund 2011; Hennebelleet al. 2011; Federrath & Klessen 2012; Burkhart & Lazarian2012; Schneider et al. 2012; Konstandin et al. 2012a; Schnei-der et al. 2013; Burkhart et al. 2014; Klessen & Glover 2016;Burkhart 2018; Mocz & Burkhart 2018; Mocz & Burkhart2019; Beattie & Federrath 2020; Menon et al. 2021). Onesuch important statistic is the density probability density © a r X i v : . [ a s t r o - ph . GA ] F e b Beattie, Mocz, Federrath & Klessen function, the ρ -PDF. For example, simple models of molec-ular clouds (MCs) in the ISM, which are supersonic andmagnetised, and have not yet started to collapse under theirown self-gravity have an approximately log-normal volume-weighted ρ -PDF, p s ( s ; σ s ) d s = 1 √ πσ s exp (cid:18) − ( s − s ) σ s (cid:19) d s, (1) s ≡ ln( ρ/ρ ) , (2) s = − σ s , (3) σ s = f ( M , M A , b, γ, Γ , (cid:96) D ) , (4)where ρ is the cloud density, with ρ being its volume-weighted mean value. The log-density variance, σ s , is thekey parameter for this model. It is a function of (i) the tur-bulent Mach number, M = σ V /c s , (5)where σ V is the root-mean-square (rms) velocity and c s is the sound speed (Vazquez-Semadeni 1994; Padoan et al.1997; Passot & V´azquez-Semadeni 1998; Price et al. 2011;Konstandin et al. 2012b), (ii) the Alfv´en Mach number, M A = σ V /V A = c s M /V A , (6)where V A = B/ √ πρ , is the Alfv´en wave velocity and B isthe magnetic field (Padoan & Nordlund 2011; Molina et al.2012), (iii) the turbulent driving parameter b (Federrathet al. 2008; Federrath et al. 2010), (iv) the adiabatic index γ (Nolan et al. 2015), (v) the polytropic index Γ (Feder-rath & Banerjee 2015), and (vi) the turbulence driving scale (cid:96) D , i.e., the scale on which energy is injected into the sys-tem (Bialy & Burkhart 2020). According to this model, thedensity variance plays two important roles: First, it encapsu-lates all of the physical processes that influence the densityfluctuations of a cloud, fully parameterising the log-normaldensity fluctuation theory, and second, the density varianceis the key ingredient for turbulence-regulated star formationtheories, including for both determining the star formationrate of an MC, directly through integrating the ρ -PDF andfor setting the width of the Gaussian component of the ini-tial core and stellar mass function, (e.g. Padoan & Nordlund2002; Hennebelle & Chabrier 2008, 2009; Krumholz & Mc-Kee 2005; Padoan & Nordlund 2011; Hennebelle et al. 2011;Federrath & Klessen 2012; Hopkins 2013; Kainulainen et al.2014; Burkhart 2018; Krumholz & Federrath 2019).However, all of these studies explicitly (or implicitly)assume that the variance is distributed isotropically in space,i.e., that the magnetic, density or velocity fluctuations donot have a preferential direction. This assumption becomesespecially violated in the presence of strong magnetic meanfields creates significant anisotropy in the flow which changeshow the momentum and thus the energy is transported alongand across the mean magnetic field, B (for incompressible,sub-Alfv´enic turbulence see Goldreich & Sridhar 1995; Choet al. 2002; Boldyrev 2006; Skalidis & Tassis 2020).Nature seems to be more than content with exploringthis anisotropic regime in the ISM. In a detailed analysisof the velocity gradients Hu et al. (2019) finds a numberof star-forming MCs that are trans- to sub-Alfv´enic, i.e. M A0 (cid:46) ⇒ ρ σ V (cid:46) B , where M A0 is the Alfv´en Mach number of the mean-field . Hence in the sub-Alfv´enic mean-field regime the magnetic energy is greater than or withinequipartition of the turbulent energies. Hu et al. (2019)reported upon Taurus ( M A0 = 1 . ± . . ± . . ± . . ± . . ± . . Another extremely magnetisedcloud is the central molecular zone cloud, G0.253+0.016,studied in Federrath et al. (2016). Pillai et al. (2015) mea-sures a mean field of B = (2 . ± .
95) mG for this cloud,and with the density and velocity quantities measured inFederrath et al. (2016) we find a corresponding Alfv´enMach number, M A0 = 0 . ± .
2, placing it well within thisanisotropic regime. For more details of this calculation, seeBeattie et al. (2020).Furthermore, recent analyses of both CO and columndensity maps of the Taurus cloud hint that sub-Alfv´enicflows are present, simply based upon the observed densitystructures in the clouds. For example, in the CO map ofthe Taurus cloud Palmeirim et al. (2013) and Heyer et al.(2020) find bimodal distributions of density orientationswith respect to the plane-of-sky magnetic field, which havebeen associated with sub-Alfv´enic turbulence (Li et al. 2013;Soler et al. 2013; Burkhart et al. 2014; Soler et al. 2017;Soler & Hennebelle 2017; Tritsis & Tassis 2018; Beattie &Federrath 2020; Seifried et al. 2020; K¨ortgen & Soler 2020).We discuss these oriented density structures further, and inmuch more detail in §
7. Hence, to understand star formationand the physical processes that create density fluctuationsin this extremely magnetised, supersonic turbulence regime,one must first construct a model for σ ρ/ρ , which is the pri-mary goal of this study.This study is organised as follows: in §
2, we considerthe underlying motivation for the ρ -PDF model in the con-text of turbulent fluctuations and the central limit theoremfrom the pioneering work of Vazquez-Semadeni (1994). In § § § § § § We use the same definition as Equation 6, but using the meanmagnetic field instead of the total field, hence, M A0 = σ V /V A0 ,where V A0 = B / √ πρ . Not all MCs are sub-Alfv´enic in nature. See for example Lunt-tila et al. (2008) for quite convincing results that clouds observedin Troland & Crutcher (2008) are super-Alfv´enic based on themagnetic and column density field correlations. The clouds we listas sub-Alfv´enic are not part of the survey reported in Troland &Crutcher (2008). MNRAS000
2, we considerthe underlying motivation for the ρ -PDF model in the con-text of turbulent fluctuations and the central limit theoremfrom the pioneering work of Vazquez-Semadeni (1994). In § § § § § § We use the same definition as Equation 6, but using the meanmagnetic field instead of the total field, hence, M A0 = σ V /V A0 ,where V A0 = B / √ πρ . Not all MCs are sub-Alfv´enic in nature. See for example Lunt-tila et al. (2008) for quite convincing results that clouds observedin Troland & Crutcher (2008) are super-Alfv´enic based on themagnetic and column density field correlations. The clouds we listas sub-Alfv´enic are not part of the survey reported in Troland &Crutcher (2008). MNRAS000 , 1–16 (2020) he anisotropic density variance Log-normal models for the density can be traced backto Vazquez-Semadeni (1994), where they explore differentfunctional forms for the ρ -PDF in order to model theself-similar, hierarchical structure of density in the cool,pressureless ( M (cid:29) δ = ρ/ρ , and motivates the log-normalPDF by assuming that for some time, t n , the density can beexpressed as a multiplicative interaction through indepen-dent density fluctuations, hence, ρ ( t n ) = δ n δ n − . . . δ δ ρ ( t ) = (cid:32) n (cid:89) i =0 δ i (cid:33) ρ ( t ) , (7)where ρ ( t ) is the initial density. Under the log-transform,the product becomes a sum,ln ρ ( t n ) = (cid:32) n (cid:88) i =0 ln δ i (cid:33) + ln ρ ( t ) . (8)If each δ follows the same underlying distribution and areindependent from each other (i.e., not correlated) then thecentral limit theorem states that the distribution of the log-density affected by this additive random process is approxi-mately normal, specifically: √ n (cid:2) ln ρ ( t n ) − (cid:104) ln ρ (cid:105) t (cid:3) = √ ns n → N (0 , σ s ) , (9)where N (0 , σ s ) is the normal distribution with mean zero,as n , the number of density fluctuations, approaches infin-ity. This lets us understand the nature of a single densityfluctuation changing in time. However, Vazquez-Semadeni(1994) further argues that since the hydrodynamical equa-tions are self-similar in space (i.e., invariant under arbitrarylength scalings), the fluctuations should be log-normal glob-ally. Passot & V´azquez-Semadeni (1998) extend this idea,recasting δ as a perturbation in the density from a transientshock, rather than just a generic turbulent fluctuation. Thismeans that one could relate the well-known shock-jump re-lations to the σ ρ variable , formulating the relation, σ ρ/ρ ∝ M . (10)As described in §
1, the density σ ρ/ρ − M relation has beenstudied intensely over the last two decades. The main ideais that σ ρ/ρ , or the spread of the ρ -PDF is a function of theunderlying physical processes that govern the fluid, whetherit be motivated by the supersonic plasmas of the ISM (e.gMac Low & Klessen 2004; Federrath et al. 2008; Federrathet al. 2010; Price et al. 2011; Padoan & Nordlund 2011; Fed-errath & Klessen 2012; Ginsburg et al. 2013; Federrath &Banerjee 2015; Klessen & Glover 2016) or the subsonic den-sity fluctuations in the hot, stratified, intracluster medium(ICM; e.g. Mohapatra et al. 2020a,b). More specifically, foran isothermal, turbulent, hydrodynamic medium, σ s = ln(1 + b M ) , (11) For a log-normally distributed variable, ρ/ρ , the variance ofthe log-variable obeys σ ρ/ρ = ln(1 + σ ρ/ρ ) and is the valuethat is commonly reported in the literature. Throughout thisstudy we however consider σ ρ/ρ , without making any assump-tions of the underlying distribution. where b is controlled by the amount of solenoidal ( ∇ × F )and compressive ( ∇ · F ) modes injected by the turbulencedriving field F (Federrath et al. 2008; Federrath et al. 2010;Price et al. 2011). Introducing a non-isothermal equation ofstate gives the relation, σ s = ln(1 + b M (5 γ +1) / ) , b M ≤ , ln (cid:18) γ + 1) b M ( γ − b M + 2 (cid:19) , b M > , (12)where γ is the adiabatic index of the medium (Nolan et al.2015). For isothermal and magnetised turbulence, the vari-ance changes with the correlation between B and ρ , σ s = ln (cid:18) b M β β + 1 (cid:19) , B ∝ ρ / , ln (cid:18) (cid:20)(cid:113) (1 + β ) + 4 b β (cid:21) − − β (cid:19) , B ∝ ρ, (13)where β = 2 M / M is the plasma beta with respect tothe mean magnetic field, which was explored in Molina et al.(2012) and will be discussed in more detail in § The relation between density contrasts and the volume-weighted variance of the underlying field is given by σ ρ/ρ = 1 V (cid:90) V d V (cid:18) ρρ − (cid:28) ρρ (cid:29)(cid:19) = 1 V (cid:90) V d V (cid:18) ρρ − (cid:19) , (14)where V is the volume for the fluid region of interest (Padoan& Nordlund 2011). Molina et al. (2012) turn Equation 14into an integral over density contrasts by constructing thefunction d V = f ( ρ /ρ ) d( ρ /ρ ) based on the shock geom-etry, where ρ /ρ is the density contrast for a single shock.Molina et al. (2012) use this to derive an analytical modelfor σ ρ/ρ for different B − ρ correlations, assuming that theshock geometry is the same for all shocks in the magnetisedplasma, and that there is no preferential direction for thedensity contrast produced by the shocks, i.e., that they areisotropic. This is a good approximation for M A0 ≥
2, whenthe turbulent component of the magnetic field is larger thanor equal to the mean-field component (Beattie et al. 2020).However, it breaks down for M A0 (cid:46)
1, i.e., when the meanfield is very strong, relevant to this study.Anisotropic, sub-Alfv´enic, compressive density struc-tures were studied in detail in Beattie & Federrath (2020),showing that the anisotropy of the density fluctuationsis a function of both M and M A0 . They attributed theanisotropy to hydrodynamical shocks that form perpendic-ular to B , causing parallel fluctuations, and fast magne-tosonic waves that form parallel to B , causing perpendic-ular fluctuations, illustrated schematically in Figure 1. Thisis consistent with findings from Tritsis & Tassis (2016), whoshowed that observations of striations (density perturbationsthat form parallel to B in sub-Alfv´enic flows) are repro-ducible when one considers fast magnetosonic wave pertur-bations. In this study, we take this phenomenology furtherand model the variance of density structures arising from asupersonic velocity field with a strong magnetic guide field,considering these two types of shocks. MNRAS , 1–16 (2020)
Beattie, Mocz, Federrath & Klessen
The general form of this variance can be constructed asfollows , σ ρ/ρ = σ ρ/ρ (cid:107) + σ ρ/ρ ⊥ + 2 σ ρ/ρ (cid:107) σ ρ/ρ ⊥ , (15)which decomposes the total variance into components of thedensity variance parallel to B , termed σ ρ/ρ (cid:107) , and perpen-dicular to B , termed σ ρ/ρ ⊥ . By assuming that these fluc-tuations are independent from one another we can furthersimplify the equation to σ ρ/ρ ≈ σ ρ/ρ (cid:107) + σ ρ/ρ ⊥ , (16)where the correlation term, 2 σ ρ/ρ (cid:107) σ ρ/ρ ⊥ , between the fluc-tuations becomes zero in the statistical average . For a shortdiscussion on how the velocity and magnetic field is arrangedin this turbulence regime we refer the reader to Appendix A.The main purpose of this study is therefore to explore andmodel each of these components, extending the work ofMolina et al. (2012) and Beattie & Federrath (2020) andindeed the numerous works on σ ρ/ρ , into the sub-Alfv´enic,anisotropic, supersonic turbulent regime. Furthermore, withthe development of a density fluctuation model, this is thefirst step towards an anisotropic, magneto-turbulent star for-mation theory. In this section we construct an anisotropic variance model.We first explore some geometrical features of shocks, sum-marising the key works of Padoan & Nordlund (2011) andMolina et al. (2012), and derive Equation 14. Next we cre-ate our two-shock model for the anisotropic density field,based on hydrodynamical shocks that travel parallel to themean magnetic field, which we call type I shocks, and fastmagnetosonic shocks that travel perpendicular to the meanmagnetic field, which we call type II shocks. Finally, we dis-cuss the volume-filling fractions of the two shock types andthe limiting behaviour of the model.
Consider a planar shock with surface area L and shockwidth λ , propagating in a system with volume L . The vol-ume of the shock is then V shock = λL . (17) By expressing the variance as we have below we are assumingthat the density fluctuations are symmetric around the magneticfield. This was shown explicitly through the elliptic symmetrywith respect to the magnetic field of the 2D density power spectrain Figure 2 of Beattie & Federrath (2020) and structure functionin Hu et al. (2020). This is a justified assumption because the statistically-averagedelliptic power spectra in Figure 2 of Beattie & Federrath (2020)do not show any rotations in the k ⊥ − k (cid:107) plane, i.e. the principalaxes of the ellipses fitted to the power spectra are aligned withthe k ⊥ − k (cid:107) coordinate axis, where k ⊥ corresponds to wave num-bers perpendicular to B and k (cid:107) for parallel wave numbers. Thismeans that the density fluctuations along and across the meanfield are statistically independent from each other, and hence2 σ ρ/ρ (cid:107) σ ρ/ρ ⊥ ≈ The shock width is proportional to the density jump betweenthe pre-shock ( ρ ) and post-shock ( ρ ) densities, multipliedby the integral scale of the turbulence, θL , i.e., the scalewhere velocity structure is no longer correlated, λ ≈ θL ρ ρ , (18)which is derived in Padoan & Nordlund (2011) by balancingthe ram and thermal pressures for a shock in hydrodynami-cal turbulence. We substitute Equation 18 into 17 to revealthe volume of the shock in terms of the density contrast, V shock ≈ θL ρ ρ . (19)Now assuming that d V ≈ d V shock we can construct the vol-ume differential, substitute it into Equation 14 and integrateit to construct the variance as a function of density contrasts.Hence, d V ≈ θL (cid:18) ρ ρ (cid:19) d (cid:18) ρ ρ (cid:19) , (20)and therefore we can rewrite Equation 14 in terms of densitycontrasts between 1 ≤ ρ /ρ ≤ ρ/ρ , which is the densitydomain for which shocks are well-defined in, σ ρ/ρ ≈ V (cid:90) ρ/ρ d (cid:18) ρ ρ (cid:19) θ V (cid:18) ρ ρ (cid:19) (cid:18) ρ ρ − (cid:19) , (21) ≈ θ (cid:90) ρ/ρ d (cid:18) ρ ρ (cid:19) (cid:18) − ρ ρ (cid:19) , (22) σ ρ/ρ θ ≈ ρρ (cid:124)(cid:123)(cid:122)(cid:125) over-densities − under-densities (cid:122)(cid:125)(cid:124)(cid:123) ρ ρ − (cid:18) ρρ (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) logarithmicover-densities , (23)where V = L cancels between the differential element ofthe shock and volume-weighted integral. What this meansis that for a single-shock model we average over the shock-jump conditions through the whole volume of interest. Thiswill be important for generalising this model to multipleshocks in § ρ/ρ , the over-densities, dominates the total variance, since the amplitudeof the density fluctuations is large in supersonic turbulence.Both ρ /ρ , the under-densities, and ln( ρ/ρ ), the logarith-mic over-densities are small in comparison, allowing us toapproximate the total variance as solely dependent upon theover-densities (Padoan & Nordlund 2011). Assuming the in-tegral scale factor is θ ≈
1, i.e., of order the system scale(Federrath & Klessen 2012), and making the above approx-imation we find the relation, σ ρ/ρ ≈ ρρ , (24)which is the key result that Padoan & Nordlund (2011) andMolina et al. (2012) use to relate the variance to the shock-jump conditions. Before generalising this result, we first con-sider the details of the two shock types that we use to con-struct our model of the anisotropic density variance. We consider the transformation d V = | d V d( ρ /ρ ) | d( ρ /ρ ) tosimplify the treatment of the limits in Equation 21. Note that we integrate ρ /ρ from 1, i.e., when the pre- andpost-shock densities are the same, up to an arbitrarily large den-sity contrast, ρ/ρ . This is the range of density contrasts wherea shock is well defined. MNRAS000
1, i.e., of order the system scale(Federrath & Klessen 2012), and making the above approx-imation we find the relation, σ ρ/ρ ≈ ρρ , (24)which is the key result that Padoan & Nordlund (2011) andMolina et al. (2012) use to relate the variance to the shock-jump conditions. Before generalising this result, we first con-sider the details of the two shock types that we use to con-struct our model of the anisotropic density variance. We consider the transformation d V = | d V d( ρ /ρ ) | d( ρ /ρ ) tosimplify the treatment of the limits in Equation 21. Note that we integrate ρ /ρ from 1, i.e., when the pre- andpost-shock densities are the same, up to an arbitrarily large den-sity contrast, ρ/ρ . This is the range of density contrasts wherea shock is well defined. MNRAS000 , 1–16 (2020) he anisotropic density variance Figure 1.
A schematic of the orientation, propagation directionand shock jump, ρ/ρ , of the two shocks discussed in § Left: type I, hydrodynamical shocks form from streaming material upand down B . Right: type II, fast magnetosonic shocks that formfrom longitudinal compressions of the magnetic field lines. Thesecause longitudinal compressions in the density field because thedensity is flux-frozen to the magnetic field.
We have now seen that the variance of a stochastic den-sity field full of shocks can be related to the shock-jumpconditions. The shock-jump conditions can be derived inthe regular Rankine-Hugoniot fashion, where we equate theupstream and downstream B , ρ v and ρ (for an isothermalshock) across the shock boundary, in the local frame of theshock, to conserve energy, momentum and mass (Landau &Lifshitz 1959), and we follow Molina et al. (2012) to includethe magnetic pressure contribution in the derivation. In thefollowing two subsections we will describe two different typesof shocks that are derived using this method. In isotropic, supersonic turbulence shocks are randomly ori-ented in space. However, in the presence of a strong meanmagnetic field, where δB (cid:28) B , shocks form from velocitygradients along the magnetic field (Beattie et al. 2020). Wecall these shocks type I shocks. The shocks create large den-sity contrasts perpendicular to the mean field that propagateparallel to it (Beattie & Federrath 2020). Since the propaga-tion is parallel to B the shocks do not feel the Lorentz force(nor does the magnetic field feel the shock, with upstreamand downstream B remaining unchanged), and hence theshock-jump conditions are hydrodynamical in nature (Mocz& Burkhart 2018). This is the supersonic (non-linear) equiv-alent of slow or acoustic modes from linearised MHD turbu-lence theory. For an isothermal fluid, with turbulence driving parameter b , the shock jump is (cid:18) ρρ (cid:19) (cid:107) = b M , (25)where b M is the compressive component of the Mach num-ber (Federrath et al. 2008; Konstandin et al. 2012b), whichis the Mach number associated with the compressive modesin the flow. It follows from Equation 24 that σ ρ/ρ (cid:107) ≈ b M . (26) Fast magnetosonic waves are the only MHD waves (in lin-earised MHD theory) that can propagate perpendicular tothe magnetic field and compress the gas (Landau & Lifshitz1959; Lehmann et al. 2016; Tritsis & Tassis 2016). The per-pendicular waves propagate with velocity v = (cid:113) v + c s , (27)where v A0 = B / √ πρ is the mean-field Alfv´en velocity,and c s is the sound speed. The non-linear counterpart isthe fast magnetosonic shock, which too propagates orthog-onal to B . We call this a type II shock. These shocks formthrough longitudinal compressions of the field lines, visu-alised in Figure 2 of Beattie et al. (2020). The field linescompress together, causing longitudinal, striated compres-sions in the density field. Through the flux-freezing condi-tion, the magnetic field is locked perpendicular to the shockwith B ∝ ρ , hence these shocks are equivalent to the B ∝ ρ shocks described in Molina et al. (2012). In an isothermalfluid the gas pressure scales directly with ρ and hence themagnetic field becomes larger at higher ρ/ρ in the fluid(Mocz & Burkhart 2018). The shock jump is given by (cid:18) ρρ (cid:19) ⊥ = 12 (cid:115)(cid:18) M b M (cid:19) + 8 M − (cid:18) M b M (cid:19) , (28)and hence, σ ρ/ρ ⊥ ≈ (cid:115)(cid:18) M b M (cid:19) + 8 M − (cid:18) M b M (cid:19) . (29)Unlike the hydrodynamical shock-jump condition in Equa-tion 25, the density jump for fast magnetosonic waves isasymptotic to ( ρ/ρ ) ⊥ = ( (cid:112) M − / M → ∞ .The limit is controlled entirely by the strength of the mag-netic field. Physically, this corresponds to the strong mag-netic field suppressing small-scale fluctuations that are in-troduced to the density as M increases (i.e., steepening ofthe ρ power spectra; Kim & Ryu 2005; Beattie & Federrath2020). If M A0 → ρ/ρ ) ⊥ →
0. This is because thedivergence of the velocity field only has a parallel to B component as B (cid:29) δB , where δB is the fluctuating com-ponent of the magnetic field. Hence one cannot create anyperpendicular density contrasts for B → ∞ (Beattie et al.2020). MNRAS , 1–16 (2020)
Beattie, Mocz, Federrath & Klessen
Now we combine Equation 26 and Equation 29 to modelthe total variance of the anisotropic density field, i.e., weassume the stochastic medium of interest is composed outof the type I and type II shocks discussed above. Thus, wemodel the total fluctuations as the sum of integrals over theshock-jump conditions in the parallel and perpendicular di-rections, as discussed in § N -shockmodel would sum over N types of uncorrelated shocks withshock-jump conditions ( ρ/ρ ) i and shock volumes V i , givenby σ ρ/ρ = 1 V N (cid:88) i (cid:90) V i d V i (cid:20)(cid:18) ρρ (cid:19) i − (cid:21) . (30)For our two-shock model ( N = 2, for type I and type II),this leads to σ ρ/ρ = σ ρ/ρ (cid:107) + σ ρ/ρ ⊥ = V V (cid:90) ρ/ρ d (cid:18) ρ ρ (cid:19) (cid:107) (cid:18) ρ ρ (cid:19) − (cid:107) (cid:18) ρ ρ − (cid:19) (cid:107) (cid:124) (cid:123)(cid:122) (cid:125) σ (cid:107) + V V (cid:90) ρ/ρ d (cid:18) ρ ρ (cid:19) ⊥ (cid:18) ρ ρ (cid:19) − ⊥ (cid:18) ρ ρ − (cid:19) ⊥ (cid:124) (cid:123)(cid:122) (cid:125) σ ⊥ , (31) σ ρ/ρ = f (cid:107) b M (cid:124) (cid:123)(cid:122) (cid:125) σ (cid:107) + f ⊥ (cid:115)(cid:18) M b M (cid:19) + 8 M − (cid:18) M b M (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) σ ⊥ , (32)where the parameters f (cid:107) = V / V and f ⊥ = V / V controlthe weighting of the contributions from each of the shocktypes, and come directly from integrating Equation 14 overtwo different sub-volumes that the fluctuations occupy. Wecall these the volume-fractions of the parallel and perpen-dicular fluctuations, respectively, and discuss and determine f (cid:107) and f ⊥ in § Here we describe the numerical data that we use to testour density variance model. These are high-resolution, three-dimensional, ideal magnetohydrodynamical (MHD) simula-tions of supersonic turbulence. The ideal, isothermal MHDequations are ∂ρ∂t + ∇ · ( ρ v ) = 0 , (33) (cid:18) ∂∂t + v · ∇ (cid:19) ρ v = ( B · ∇ ) B π − ∇ (cid:18) c s ρ + | B | π (cid:19) + ρ F , (34) ∂ B ∂t = ∇ × ( v × B ) , (35) ∇ · B = 0 , (36)where v is the fluid velocity, ρ the density, B the magneticfield, c s the sound speed and F , a stochastic function thatdrives the turbulence. We use a modified version of flash based on version 4.0.1 (Fryxell et al. 2000; Dubey et al.2008) to solve the MHD equations in a periodic box withdimensions L , on a uniform grid with resolution 512 , us-ing the multi-wave, approximate Riemann solver frameworkdescribed in Bouchut et al. (2010) and implemented in flash Table 1.
Main simulation parameters.Sim. ID M ( ± σ ) M A0 ( ± σ ) σ ρ/ρ ( ± σ ) N (1) (2) (3) (4) (5) M2Ma0.1 ± ± ± M4Ma0.1 ± ± ± M10Ma0.1 ± ± ± M20Ma0.1 ± ± ± M2Ma0.5 ± ± ± M4Ma0.5 ± ± ± M10Ma0.5 ± ± ± M20Ma0.5 ± ± ± M2Ma1.0 ± ± ± M4Ma1.0 ± ± ± M10Ma1.0 ± ± ± M20Ma1.0 ± ± ± Notes:
For each simulation we extract 51 realisations at times t/T ∈ { . , . , . . . , . } , where T is the turbulent turnovertime, and all 1 σ fluctuations listed are from the time-averagingover the 5 T time span. Column (1): the simulation ID. Col-umn (2): the rms turbulent Mach number, M = σ V /c s . Column(3): the Alfv´en Mach number for the mean- B component, B , M A0 = (2 c s M√ πρ ) / | B | , where ρ is the mean density, c s isthe sound speed. Column (4) the volume-weighted density vari-ance, computed as shown in Equation 14. Column (5): the numberof computational cells for the discretisation of the spatial domainof size L . by Waagan et al. (2011). For a detailed discussion of the sim-ulations we refer the reader to Beattie & Federrath (2020)and Beattie et al. (2020). Table 1 provides a summary of theimportant time-averaged parameter values. Here we brieflysummarise the key methods and properties of the simula-tions. MNRAS000
For each simulation we extract 51 realisations at times t/T ∈ { . , . , . . . , . } , where T is the turbulent turnovertime, and all 1 σ fluctuations listed are from the time-averagingover the 5 T time span. Column (1): the simulation ID. Col-umn (2): the rms turbulent Mach number, M = σ V /c s . Column(3): the Alfv´en Mach number for the mean- B component, B , M A0 = (2 c s M√ πρ ) / | B | , where ρ is the mean density, c s isthe sound speed. Column (4) the volume-weighted density vari-ance, computed as shown in Equation 14. Column (5): the numberof computational cells for the discretisation of the spatial domainof size L . by Waagan et al. (2011). For a detailed discussion of the sim-ulations we refer the reader to Beattie & Federrath (2020)and Beattie et al. (2020). Table 1 provides a summary of theimportant time-averaged parameter values. Here we brieflysummarise the key methods and properties of the simula-tions. MNRAS000 , 1–16 (2020) he anisotropic density variance M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . B B M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . M ≈ M ≈ M A0 ≈ M A0 ≈ log ( ρ/ρ ): [ − . , . ( ρ/ρ ): [ − . , . Figure 2.
Slices through the y = 0 plane of log ( ρ/ρ ) at t = 6 T for the 12 simulations listed in Table 1. The mean magnetic field, B ,is oriented up the page as shown in the top-left panel. The plots are ordered such that the simulations with the highest Mach numbersare on the right column ( M ≈
20) and weakest (
M ≈
2) on the left. Likewise, the strongest B simulations ( M A0 ≈ .
1) are on the toprow and weakest ( M A0 ≈ .
0) on the bottom. The top row reveals slices of the density structures for the most anisotropic simulations.For these simulations, fast magnetosonic waves (type II shocks, § B and leaving striations from the compressions and rarefactions in the field. The largest over-densities however form along the B , whichare space-filling for low- M , and narrow and highly-compressed for high- M (type I shocks, § M A0 = 1 . The initial velocity field is set to v ( x, y, z, t = 0) = , withunits c s = 1, and the density field ρ ( x, y, z, t = 0) = ρ ,with units ρ = 1. The turbulent acceleration field, F , inEquation 34 follows an Ornstein-Uhlenbeck process in timeand is constructed such that we can control the mixture ofsolenoidal and compressive modes in F (see Federrath et al.2008; Federrath et al. 2009, 2010 for a detailed discussion ofthe turbulence driving). We choose to drive with a naturalmixture of the two modes (Federrath et al. 2010), which cor-responds to b ≈ .
4. We isotropically drive in wavenumberspace, centred on k = 2, corresponding to real-space scalesof (cid:96) D = L/
2, and falling off to zero with a parabolic spec-trum between k = 1 and k = 3. Thus, energy is injectedonly on large scales and the turbulence on smaller scales develops self-consistently through the MHD equations. Theauto-correlation timescale of F is equal to T = (cid:96) D σ V = L c s M . (37)We vary the sonic Mach number between M = 2 and 20,encapsulating the range of observed M values for turbu-lent MCs (e.g., Schneider et al. 2013; Federrath et al. 2016;Orkisz et al. 2017; Beattie et al. 2019b). We run the simu-lations from 0 ≤ t/T ≤
10 and extract 51 time realisationsof ρ ( x, y, z ) over 5 ≤ t/T ≤
10 to gather data only when theturbulence is in a statistically stationary state (Federrathet al. 2009; Price & Federrath 2010). For simulations with astrong guide field, a statistically steady state is reached afterabout 5 T , while for purely hydrodynamical turbulence andsuper-Afv´enic turbulence, stationarity is already reached af- MNRAS , 1–16 (2020)
Beattie, Mocz, Federrath & Klessen
Figure 3.
The structure of the 3D density field, ρ/ρ , for the highly-magnetised M2MA0.1 (left) and
M20MA0.1 (right) simulations (seeTable 1). The direction of the mean magnetic field, B , is shown at the bottom-right corner of the simulations. Large-scale vortices formin the density plane perpendicular to B . Density structures are tightly wound and stretched into ribbons and tubes as the supersonicvortices advect the flux-frozen density field around the field lines. The vortices are found on similar scales for both simulations, roughly L/
2, the driving scale of the turbulence, however the
M20MA0.1 simulation has more small-scale density fluctuations than the
M2MA0.1 simulation (note the different scaling of the colour bars). ter about 2 T . All our results are based on time-averagesacross these 51 realisations, unless explicitly indicated oth-erwise.In Figure 2 we show density slices through a plane per-pendicular to the direction of B , indicated in the top-leftplot. The plots are organised such that the top-left plothas the weakest (but still supersonic) turbulence ( M = 2)and strongest B ( M A0 = 0 . M = 20) and weakest B ( M A0 = 1 . M increases, extending across ∝ L/ M ∝ L/ M ≈
2, to a tiny fraction of the boxat
M ≈
20, qualitatively consistent with the shock-widthmodel discussed in § ospray ray-tracer in visit (Childs et al. 2012)we show examples of the volumetric density field for the M2MA0.1 and
M20MA0.1 simulations in Figure 3. Large-scalevortices are revealed in the
M2MA0.1 simulation, with sheetsof density wrapped around in the direction of the B . Thevortices form at roughly the driving scale ( L/
2) of the sim-ulation. The volumetric rendering for
M20MA0.1 shows muchmore small-scale structure compared to
M2MA0.1 , both alongand perpendicular to the direction of B . The large-scalevortices found in M2MA0.1 are roughly at the same scale asthe vortices in
M20MA0.1 , suggesting that it is the drivingscale that sets the size of the vortices. Dense, shocked mate-rial is compressed along the magnetic field lines, creating the filamentary structures that were found in the slices, orientedperpendicular to the field.
The initial magnetic field, B ( x, y, z ) = B ˆ z at t = 0, inEquations (34–36) is a uniform field with field lines threadedthrough the ˆ z direction of the simulations. B is set to en-sure the desired M A0 , using the definition of the Alfv´envelocity (see footnote 1) and the turbulent Mach number(see Equation 5), M A0 = σ V /V A0 = 2 c s √ πρ M /B , (38)with B constant in space and time. We set M A0 = 0 . , . . B ( t ) = B ˆz + δ B ( t ) , (39)where the fluctuating component of the field, δ B , evolvesself-consistently from the MHD equations, with (cid:104) δ B (cid:105) t = 0.From previous experiments in the anisotropic regime dis-cussed throughout this study, | δ B | / | B | ≈ − − − (Fed-errath 2016a; Beattie et al. 2020). For this reason M A ≈M A0 in this regime, however, we explicitly formulate ourmodels around M A0 , which is an invariant across different M simulations. In this section we describe and model the fluctuationvolume-filling fractions that we introduced in Equation 32.
MNRAS000
MNRAS000 , 1–16 (2020) he anisotropic density variance The fluctuation volume fractions, f (cid:107) and f ⊥ in Equation 32,determine the contributions of the parallel and perpendic-ular fluctuations in our model. These parameters naturallycome about from considering multiple terms of a volume-weighted statistic, where the total volume is partitioned intodistinct sub-volumes for the parallel and perpendicular fluc-tuations. The filling fractions can be derived by returning toEquation 17, but instead of using a planar shock, consideringa non-planar shock of the form V shock = λL D , (40)where D is the fractal dimension, or space-filling factor of theshock, e.g., for D = 1 the shock is contained within a line,and D = 2 is the planar shock from Equation 17. Consideringfractal shocks is an important consideration to make in ourtwo-shock model, because the amount of volume that thefluctuations fill encodes how much they contribute to thevariance through the density variance integral, Equation 14.Propagating Equation 40 through Equations 20-24, makingthe same assumptions but for non-planar shocks, shows that σ ρ/ρ ≈ L D (cid:107) +1 L (cid:34)(cid:18) ρρ (cid:19) (cid:107) − (cid:18) ρ ρ (cid:19) (cid:107) − (cid:18) ρρ (cid:19) (cid:107) (cid:35) (41)+ L D ⊥ +1 L (cid:20)(cid:18) ρρ (cid:19) ⊥ − (cid:18) ρ ρ (cid:19) ⊥ − (cid:18) ρρ (cid:19) ⊥ (cid:21) , (42)where D (cid:107) is the fractal dimension of the type I shocks and D ⊥ is the fractal dimension of the type II shocks. In Equa-tion 31 we show that ( L D (cid:107) +1 ) /L = V / V = f (cid:107) is thenthe volume-filling fraction for the parallel fluctuations, andlikewise for f ⊥ . Note that the geometry of the fractal shockpropagates through to the total volume fraction of the par-allel and perpendicular fluctuations, including the rarefac-tions. This is because the rarefaction region in the flow willshare some of the same geometrical properties as the shock,since the gas density is compressed from the rarefaction intothe over-density. However, as we discussed in § σ ρ/ρ ≈ f (cid:107) (cid:18) ρρ (cid:19) (cid:107) + f ⊥ (cid:18) ρρ (cid:19) ⊥ , (43)and Equation 32 immediately follows by substituting in theappropriate shock-jump conditions. Since the volume-fillingfractions are for the total parallel and perpendicular fluctu-ations it makes sense that they add to unity,1 = N (cid:88) i V i V ⇐⇒ f (cid:107) + f ⊥ . (44)This means our model is an N − N is the number of shock types, and more explicitly, a one-parameter model for the two-shock model we describe here.The volume fractions need not be constants, and in-deed, may depend upon both M and M A0 . For example,Konstandin et al. (2016), Beattie et al. (2019a) and Beattieet al. (2019b) demonstrate that the global D for supersonicturbulence does depend upon M , varying between 3 and 1.This is because the flow is more compressible for higher M ,reducing D , which corresponds to the emergence of highly-compressed structures like one-dimensional filaments. In-creasing M also directly affects the geometry of the shocks in the flow, which can be shown by expressing Equation 18,the shock thickness, in terms of the shock-jump conditionsfor an isothermal, hydrodynamical shock, λ ≈ L M . (45)This establishes that the shocks become more compressed(thinner), filling less space, as M increases, which meansthat a reasonable model for the volume fraction of the par-allel fluctuations, which is dominated by type I shocks is f (cid:107) ∝ M − . However, this only holds for M (cid:29)
1, but not for
M ∼
1. Thus, we require a phenomenological function forthe volume fraction that describes its dependence on M forboth the trans to supersonic flow regime. We consider theform f (cid:107) = f (cid:34) (cid:18) MM c (cid:19) (cid:35) − / , (46)which describes a function that tends towards ∼ / M when M > M c (cid:29) f , for M < M c , e.g., in the sub- trans-sonic regime, where thetype I shock thickness becomes ill-defined, and linear MHDwaves (slow, parallel and fast, perpendicular modes) weaklycompress the flow. Here, M c defines a critical Mach num-ber for the transition between the constant and the ∼ M − regime. Thus, the parameter f defines the volume-fillingfraction of the parallel fluctuations, in the presence of over-densities formed by slow modes, as discussed in § f ⊥ = 1 − f (cid:107) , as a direct consequence of Equation 44, whichmeans we only need to consider a model for f (cid:107) . The first step to fitting our model is to develop a dataset,( M i , f (cid:107) ,i ), for the different M A0 simulation listed in Table 1.Since, as discussed in § f (cid:107) , we can rearrange Equation 32 and solve analyt-ically , f (cid:107) = σ ρ/ρ − σ ρ/ρ ⊥ σ ρ/ρ (cid:107) − σ ρ/ρ ⊥ , (47)with the single constraint that f (cid:107) ∈ [0 , f (cid:107) ( M ) for each simulation. We show thisdata in Figure 4. We then fit the functional form for f (cid:107) ( M ),Equation 46, with parameters f and M c using a non-linearleast squares fitting routine weighted by 1 / ∆ f (cid:107) , where ∆ f (cid:107) is the uncertainty in f (cid:107) , propagated through Equation 47.We show the fits using the solid lines, coloured by M A0 .For each of the datasets, the critical sonic Mach numberis M c ≈
10 and the f parameter varies monotonically be-tween ≈ . .
7, listed in Table 2. There are two importantconclusions to make: (1) magnetised turbulence is signifi-cantly saturated with shocks above
M ≈
10, and (2) the This follows from σ ρ/ρ = f (cid:107) σ ρ/ρ (cid:107) + f ⊥ σ ρ/ρ ⊥ = f (cid:107) σ ρ/ρ (cid:107) +(1 − f (cid:107) ) σ ρ/ρ ⊥ and then solving for f (cid:107) . MNRAS , 1–16 (2020) Beattie, Mocz, Federrath & Klessen M − f k = − f ⊥ ∝ M − M A0 = 0 . M A0 = 0 . M A0 = 1 . Figure 4.
The volume-filling fraction of fluctuations that run par-allel to B , f (cid:107) , as a function of M , as discussed in § σ uncertainties propagated from Equation 47. We fit ourphenomenological models in Equation 46 to each of the M A0 datasets, shown the with solid lines. This models suggests that f (cid:107) ∼ λ ∼ M − in the high- M limit, where λ is the shock-width,and f (cid:107) ∼ const. in the low- M limit, which is supported well bythe data. This means f ⊥ (cid:29) f (cid:107) for high- M flows, i.e., theperpendicular fluctuations, which are dominated by fastmagnetosonic shocks, contribute most to volume fraction athigh- M , which is visualised in Figure 5. Table 2.
Fit parameters for f (cid:107) ( M ), as per Equation 46.Parameter M A0 = 0 . M A0 = 0 . M A0 = 1 . f ± ± ± M c ± ± ± parallel fluctuations become less confined by the magneticfield, and occupy more of the total volume when the mag-netic field is weakened. Accordingly, one should expect that f ∼ . M A0 → ∞ , i.e., when there is no confinementfrom the magnetic field. This reclaims the hydrodynamical σ ρ/ρ −M relation. Regardless of the field strength, the high- M behaviour is the same between the simulations, which isexpected since the shock width, λ ∼ M − , encodes howthe ( b M ) term contributes to the total variance once thefluid is sufficiently shocked. This transition is consistent withwhat Beattie & Federrath (2020) found, where M ≈ − B . Since f ⊥ = 1 − f (cid:107) (Equation 44)as M grows and the parallel fluctuations occupy less and lessof the volume until the perpendicular fluctuations contributethe most to the total volume budget for the fluid.We show some of the shocked regions for both the par-allel and perpendicular fluctuations in the top panels of Fig-ure 5, for M ≈
20 turbulence by taking the divergence alongand across B . We also show the corresponding density fluc- tuations in the bottom two panels, of the Figure, where wehighlight the −∇ · v > −∇· v > Now that we have constrained f (cid:107) we put this back into Equa-tion 48 and generate our final variance model. We depictour models with solid lines in Figure 6. The hydrodynami-cal limit is drawn in red, which provides an upper bound ofthe variance ( f (cid:107) = 1, f ⊥ = 0), and the three sets of simula-tion data at different M A0 in blue. The models fit the datawell, never deviating from the data by more than a smallfraction of the total 1 σ fluctuations in σ ρ/ρ .Our results establish the following picture of the densityfield in this regime: at low M , hydrodynamical, parallel to B , fluctuations are large-scale, occupying a significant frac-tion of the fluid volume with a relatively large volume-fillingfraction. The type I shocks dominate the parallel fluctuationcontribution to the volume-weighted variance of the densityfield. However, as M increases, the small-scale fluctuationsgrow in the field, with and the type I shocks shrink in sizewith their shock width ∼ M − , and thus the parallel fluctu-ations begin to occupy only very small fractions of the fluidvolume, decreasing the contribution from type I shocks tothe variance. When M (cid:29)
1, the total volume-weighted vari-ance becomes by type II shocks, which are highly magnetisedand do not support small-scale fluctuations. This flattens σ ρ/ρ ( M ) out, because the small-scale fluctuations that areadded to the field with increasing M are suppressed. Even-tually this leads to an upper bound on σ ρ/ρ , as the type IIshock density contrasts dominate the total variance. Let us now consider the limiting behaviour of the densityvariance in our model. The density variance is summarisedas σ ρ/ρ = f (cid:107) σ ρ/ρ (cid:107) + (1 − f (cid:107) ) σ ρ/ρ ⊥ . (48)In the high- M limit we havelim M→∞ σ ρ/ρ = f b M + 12 (cid:18)(cid:113) M − (cid:19) , (49)which, like the density contrast caused by type II shocksshown in § B . This tells us that even inthe high- M limit the turbulent driving parameter influencesthe spread of the PDF, frozen into the variance at high- M . MNRAS000
1, the total volume-weighted vari-ance becomes by type II shocks, which are highly magnetisedand do not support small-scale fluctuations. This flattens σ ρ/ρ ( M ) out, because the small-scale fluctuations that areadded to the field with increasing M are suppressed. Even-tually this leads to an upper bound on σ ρ/ρ , as the type IIshock density contrasts dominate the total variance. Let us now consider the limiting behaviour of the densityvariance in our model. The density variance is summarisedas σ ρ/ρ = f (cid:107) σ ρ/ρ (cid:107) + (1 − f (cid:107) ) σ ρ/ρ ⊥ . (48)In the high- M limit we havelim M→∞ σ ρ/ρ = f b M + 12 (cid:18)(cid:113) M − (cid:19) , (49)which, like the density contrast caused by type II shocksshown in § B . This tells us that even inthe high- M limit the turbulent driving parameter influencesthe spread of the PDF, frozen into the variance at high- M . MNRAS000 , 1–16 (2020) he anisotropic density variance L L − L B B − L − L L L L L − L − L − L − L L L − − − ∇ k · ( v / c s ) − − − ∇ ⊥ · ( v / c s ) ρ / ρ ρ / ρ Figure 5. Top left:
A slice of the parallel convergence of the velocity field, with respect to the direction parallel to B . The convergenceis in units of c s /L . Here we show a single time realisation from the M = 20, M A0 = 0 . B , reminiscent in morphology of the hydrodynamical bow shocks studied in Robertson& Goldreich (2018). Top right:
The same time realisation as the left plot, but for the perpendicular convergence of the velocity field,showing fast magnetosonic compressions of the velocity field, propagating orthogonal to B . The amplitude of the convergence is an orderof magnitude larger for the hydrodynamical shocks compared to the magnetosonic fast shocks due to the strong magnetic cushioning effectperpendicular to B (Molina et al. 2012). Bottom left: the parallel shocks in the density field, corresponding to the positive convergentstructures in the upper-left plot, with density contrasts up to ( ρ/ρ ) (cid:107) ∼ M = 400, and very low volume-filling fractions. These arethe type I shocks discussed in § σ ρ/ρ (cid:107) inEquation 32. Bottom right: the same as the bottom left plot, but for the perpendicular shocks in the density field. The volume-fillingfraction of the perpendicular fluctuations is much larger than the parallel component, however the density contrast is orders of magnitudelower. These are the type II shocks outlined in § σ ρ/ρ ⊥ . Interestingly, the limiting behaviour is independent of M ,as the effect of increasingly sharp density contrasts ( M ) ofsupersonic hydrodynamical shocks is cancelled out by thereduced volume-filling factor of the parallel fluctuations.The next limit of interest is the low- M A0 limit,lim M A0 → σ ρ/ρ = f (cid:107) b M . (50)Thus, the parallel component of the variance is retained inthis limit. We can interpret this as hydrodynamical fluctu-ations are able to survive along the magnetic field but are confined in a small volume, f (cid:107) V , and hence are reduced bythe volume-filling fraction, f (cid:107) . This value will be significantlyless than 1 for M (cid:29)
10, as measured in § σ ρ/ρ = 0 in this limit.Finally, in the high- M A0 limit, as the turbulence tran-sitions from magnetised to hydrodynamic,lim M A0 →∞ σ ρ/ρ = f (cid:107) b M + (1 − f (cid:107) ) b M = b M , (51) MNRAS , 1–16 (2020) Beattie, Mocz, Federrath & Klessen
102 3 4 5 6 7 8 9 20 30 M σ ρ / ρ b M σ ρ/ρ ( M ; M A0 = 0 . σ ρ/ρ ( M ; M A0 = 0 . σ ρ/ρ ( M ; M A0 = 1 . M A0 = 0 . M A0 = 0 . M A0 = 1 . Figure 6.
Two-shock model (Equation 32) for the density variance, σ ρ/ρ , as a function of M , for each set of M A0 simulations. Blue pointsare numerical data from the 12 simulations (Table 1), averaged over a time interval of 5 T . The red line indicates the hydrodynamicallimit for the density variance, b M . We find that the variance grows from low M to M ≈ M values) are suppressed by the strong mean magnetic field. Weplot our model, Equation 32, on each M A0 dataset, shown with the solid lines. We find good agreement between the data and our models,which suggests there is an upper bound for σ ρ/ρ in this highly-magnetised, anisotropic turbulent regime. We discuss the implicationsfor this result in § which is the well-known, isotropic, hydrodynamic M− σ ρ/ρ relation. We motivated in § N H ∝ cm − , where N H is the number density of the molecular gas column density) tend to be oriented perpen-dicularly to, and low-density ( N H ∝ cm − ) structuresparallel to, the mean magnetic field (Planck Collaborationet al. 2016a; Soler et al. 2017; Tritsis et al. 2018; Heyer et al.2020; Pillai et al. 2020). Through the analysis of numer-ical MHD turbulence simulations we have identified someof the physics that causes the bimodal alignment of gasdensity structures during this density transition. To sum-marise, the alignment in simulations like ours comes aboutthrough an interplay between the divergence (compressibil-ity) and the strain (most likely associated with vortex cre-ation; such as those visualised in Figure 3) in the flow. Thishas been found in multiple studies of highly-magnetised,compressible plasmas, including those simulations presentedin Figure 5, (Soler & Hennebelle 2017; Beattie & Federrath2020; Seifried et al. 2020; K¨ortgen & Soler 2020). A differ-ent model is developed in Xu et al. (2019). They use theGoldreich & Sridhar (1995) anisotropy theory to model theextent of low-density filaments in the warm, diffuse and sub-sonic ISM. This model is most likely not applicable for thehighly-compressible flows in the supersonic, cool and densemolecular clouds, which need not conform to the incompress-ible Goldreich & Sridhar (1995) theory. The transition oc-curs with or without the presence of self-gravity, and hence MNRAS000
Two-shock model (Equation 32) for the density variance, σ ρ/ρ , as a function of M , for each set of M A0 simulations. Blue pointsare numerical data from the 12 simulations (Table 1), averaged over a time interval of 5 T . The red line indicates the hydrodynamicallimit for the density variance, b M . We find that the variance grows from low M to M ≈ M values) are suppressed by the strong mean magnetic field. Weplot our model, Equation 32, on each M A0 dataset, shown with the solid lines. We find good agreement between the data and our models,which suggests there is an upper bound for σ ρ/ρ in this highly-magnetised, anisotropic turbulent regime. We discuss the implicationsfor this result in § which is the well-known, isotropic, hydrodynamic M− σ ρ/ρ relation. We motivated in § N H ∝ cm − , where N H is the number density of the molecular gas column density) tend to be oriented perpen-dicularly to, and low-density ( N H ∝ cm − ) structuresparallel to, the mean magnetic field (Planck Collaborationet al. 2016a; Soler et al. 2017; Tritsis et al. 2018; Heyer et al.2020; Pillai et al. 2020). Through the analysis of numer-ical MHD turbulence simulations we have identified someof the physics that causes the bimodal alignment of gasdensity structures during this density transition. To sum-marise, the alignment in simulations like ours comes aboutthrough an interplay between the divergence (compressibil-ity) and the strain (most likely associated with vortex cre-ation; such as those visualised in Figure 3) in the flow. Thishas been found in multiple studies of highly-magnetised,compressible plasmas, including those simulations presentedin Figure 5, (Soler & Hennebelle 2017; Beattie & Federrath2020; Seifried et al. 2020; K¨ortgen & Soler 2020). A differ-ent model is developed in Xu et al. (2019). They use theGoldreich & Sridhar (1995) anisotropy theory to model theextent of low-density filaments in the warm, diffuse and sub-sonic ISM. This model is most likely not applicable for thehighly-compressible flows in the supersonic, cool and densemolecular clouds, which need not conform to the incompress-ible Goldreich & Sridhar (1995) theory. The transition oc-curs with or without the presence of self-gravity, and hence MNRAS000 , 1–16 (2020) he anisotropic density variance seems to be a MHD phenomena (however, self-gravitatingMHD simulations are shown to enhance the density struc-tures that are perpendicular to B , e.g. Soler & Hennebelle2017; G´omez et al. 2018; Barreto-Mota et al. 2021).Our two-shock model suggests a simple explanation forthis density transition described above, irrespective of thephysical processes that create the aligned structures, whichis beyond the scope of this study. The shock-jump condi-tions predict that perpendicularly oriented hydrodynami-cal shocks formed from material flowing along the magneticfield are at least an order of magnitude larger in densitycontrast compared to parallel oriented fast magnetosonicshocks, formed through shuffling magnetic field lines. Weshow this qualitatively in Figure 5, with the density con-trasts produced by the shocks visualised in the two bottompanels. Hence, we suggest that the observed transition inthe density arises from observing these two distinct types ofMHD modes.A further transition is also now reported at even highercolumn densities, within the filaments ( < φ , between ∇ ρ and B . They showed that cos φ = ± φ = 0 constitute equilibrium points that an idealMHD system tends towards. In our work we demonstratedthat the physical realisation of this insight corresponds tohydrodynamical shocks that form along B (cos φ = 0) andfast magnetosonic shocks that form perpendicular to B (cos φ = ± ∇ · v < Bottom-up star formation theories treat MCs as the funda-mental building blocks that determine galactic star forma-tion rates, and are parameterised by σ ρ/ρ , as discussed in § σ ρ/ρ , and ( ρ/ρ ) crit , the critical density atwhich the cloud becomes Jeans unstable and collapses (Fed-errath 2018). The magnetic field plays a role in these modelsby reducing the total density variance (as shown in Equa-tion 13), and by introducing additional support throughmagnetic pressure in the critical density, preventing collapse(Krumholz & McKee 2005; Federrath & Klessen 2012; Fed-errath & Klessen 2013). However, all of these models treatthe magnetic field only as an isotropic contribution via the magnetic pressure. The effects of magnetic tension or theanisotropic effects introduced by a strong guide field arenot included in the current theories of star formation. Herewe show that the magnetic field, specifically, a sub-Alfv´enicmean field, which encodes tension into the theory, acts pref-erentially to suppress small-scale density fluctuations (andhence turbulent fragmentation), preventing σ ρ/ρ from grow-ing beyond a specific limit set by the strength of the field,shown in Equation 49. This is an extra form of suppressionthat B has on the star formation rate.A complete theory for anisotropic star formation wouldtake into account that fluctuations perpendicular to themagnetic field are suppressed significantly, whereas along thefield they are not. As such, the theory would predict thatstar formation may become bimodal with respect to the di-rection of the large-scale, coherent B -field. There are someobservational signatures that this may be the case with thestar formation rate of some MCs seeming to depend uponthe large-scale orientation of the magnetic field (Law et al.2019, 2020).In this study we describe a model for the variance appli-cable to these types of MCs, but to properly predict the starformation rate in a strongly magnetised environment onealso must create an anisotropic model for ( ρ/ρ ) crit , whichcontains both information about the scale in which the tur-bulence transitions from supersonic to subsonic in rms ve-locities, i.e., the sonic scale (Federrath & Klessen 2012; Fed-errath et al. 2021), and the Jeans scale. The morphologyof the sonic scale in the presence of a strong magnetic fieldis unknown, but one can speculate that it most likely willbecome stretched along the field lines, changing the natureof the critical density and how cloud collapse happens inthis regime. What we emphasise here is that there is muchwork to do in this supersonic, anisotropic, highly-magnetisedregime, much of which we intend to pursue in future studies. In this study we explore the density variance, σ ρ/ρ , ofhighly-magnetised, anisotropic, supersonic turbulence acrossand along the mean magnetic field, B . In § σ ρ/ρ − M relation and highlight the fun-damental connection between shock-jump conditions for thedensity field and σ ρ/ρ . In §
3, we describe in detail how wedefine the geometry of shocks, and finally, we derive our two-shock model for σ ρ/ρ . In § M >
1; where M is thesonic Mach number), trans-/sub-Alfv´enic ( M A0 ≤
1; where M A0 is the mean-field Alfv´enic Mach number) MHD simula-tions. We show examples of the 2D density slices and the full3D density fields in Figure 2 and Figure 3, respectively. In § § σ ρ/ρ = f (cid:107) σ ρ/ρ (cid:107) + f ⊥ σ ρ/ρ ⊥ , where σ ρ/ρ (cid:107) comes from type I shocks and σ ρ/ρ ⊥ comes from type IIshocks.(ii) The entire volume of the turbulence must contribute MNRAS , 1–16 (2020) Beattie, Mocz, Federrath & Klessen to the total variance, hence f (cid:107) = 1 − f ⊥ . This defines a singleparameter model σ ρ/ρ = f (cid:107) σ ρ/ρ (cid:107) + (1 − f (cid:107) ) σ ρ/ρ ⊥ (iii) We propose a phenomenological model for f (cid:107) , f (cid:107) = f (cid:20) (cid:16) MM c (cid:17) (cid:21) − / , based on the shock thickness of typeI shocks.(iv) We use numerical simulations parameterised by( M , M A0 ) to directly measure σ ρ/ρ , and calculate σ ρ/ρ (cid:107) and σ ρ/ρ ⊥ .(v) Using the numerical data we fit for the two parametersin the f (cid:107) model, f , which is associated with the volume-filling fraction of the parallel fluctuations in the subsoniclimit, and M c , is the M when the supersonic flow is signif-icantly saturated with shocks.Finally in § • We derive a model for the density variance, σ ρ/ρ , ina sub-Alfv´enic, supersonic, anisotropic flow regime, wherea strong mean magnetic field, B , creates dynamicallydifferent fluctuations parallel and perpendicular to field,characterised by Equation 32. To do this, we generalisethe shock-variance relations from Padoan & Nordlund(2011) and Molina et al. (2012), discussed in §
3, partition-ing the total volume into two sub-volumes that containhydrodynamical shocks moving along (parallel to) themagnetic field, and fast magnetosonic shocks moving across(perpendicular to) the field, with details of the orientationand compression mechanism shown in Figure 1. • Our density variance model relies upon the volume-filling fraction, f (cid:107) – a measure of how much relative volumethe parallel fluctuations occupy along B . We propose aphenomenological model for f (cid:107) , Equation 46, using theshock width described in Padoan & Nordlund (2011),discussed in detail in § M flows, whilst theperpendicular fluctuations dominate in high- M flows,consistent with the compressible structures visualised inFigure 2 and Figure 5. • Our new model predicts a finite value of σ ρ/ρ inthe high- M limit, shown in Equation 49, which is setby the strength of B . This is because a strong B fieldacts preferentially to suppress the small-scale fluctuationsintroduced in high- M flows. The new model also predictsa finite value of σ ρ/ρ as B → ∞ , as shown in Equa-tion 50, corresponding to density fluctuations that canpersist along B . In the hydrodynamical limit, as shown inEquation 51, our model reduces to the well-known relation, σ ρ/ρ = b M , where b is the turbulent driving parameter.We demonstrate that our variance model provides a goodfit to the simulation data in Figure 6. • In § B may addi-tionally suppress star formation by limiting the small-scaledensity fluctuations, and hence turbulent fragmentation inMCs with M (cid:38) ACKNOWLEDGEMENTS
We thank the anonymous reviewer for their detailed review,which increased the clarity and presentation of the study.J. R. B. thanks Shyam Menon for the many productivediscussions and acknowledges financial support fromthe Australian National University, via the Deakin PhDand Dean’s Higher Degree Research (theoretical physics)Scholarships, the Research School of Astronomy andAstrophysics, via the Joan Duffield Research Scholarshipand the Australian Government via the Australian Govern-ment Research Training Program Fee-Offset Scholarship.P. M. acknowledges support for this work provided by NASAthrough Einstein Postdoctoral Fellowship grant numberPF7-180164 awarded by the
Chandra
X-ray Centre, whichis operated by the Smithsonian Astrophysical Observatoryfor NASA under contract NAS8-03060. C. F. acknowledgesfunding provided by the Australian Research Council(Discovery Project DP170100603 and Future FellowshipFT180100495), and the Australia-Germany Joint ResearchCooperation Scheme (UA-DAAD). R. S. K. acknowledgesfinancial support from the German Research Foundation(DFG) via the Collaborative Research Center (SFB 881,Project-ID 138713538) ’The Milky Way System’ (subpro-jects A1, B1, B2, and B8). He also thanks for funding fromthe Heidelberg Cluster of Excellence STRUCTURES in theframework of Germany’s Excellence Strategy (grant EXC-2181/1 - 390900948) and for funding from the EuropeanResearch Council via the ERC Synergy Grant ECOGAL(grant 855130). We further acknowledge high-performancecomputing resources provided by the Australian NationalComputational Infrastructure (grant ek9) in the frameworkof the National Computational Merit Allocation Schemeand the ANU Merit Allocation Scheme, and by the LeibnizRechenzentrum and the Gauss Centre for Supercomputing(grants pr32lo, pr48pi, pr74nu). The simulation software flash was in part developed by the DOE-supported FlashCentre for Computational Science at the University ofChicago.Data analysis and visualisation software used in thisstudy:
C++ (Stroustrup 2013), numpy (Oliphant 2006; Har-ris et al. 2020), matplotlib (Hunter 2007), cython (Behnelet al. 2011), visit (Childs et al. 2012), scipy (Virtanen et al.2020), scikit-image (van der Walt et al. 2014).
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author.
MNRAS000
MNRAS000 , 1–16 (2020) he anisotropic density variance REFERENCES
Abe D., Inoue T., Inutsuka S.-i., Matsumoto T., 2020,arXiv e-prints, p. arXiv:2012.02205Barreto-Mota L., de Gouveia Dal Pino E. M., BurkhartB., Melioli C., Santos-Lima R., Kadowaki L. H. S., 2021,arXiv e-prints, p. arXiv:2101.03246Beattie J. R., Federrath C., 2020, MNRAS, 492, 668Beattie J. R., Federrath C., Klessen R. S., 2019a, MNRAS,487, 2070Beattie J. R., Federrath C., Klessen R. S., Schneider N.,2019b, MNRAS, 488, 2493Beattie J. R., Federrath C., Seta A., 2020, MNRAS, 498,1593Behnel S., Bradshaw R., Citro C., Dalcin L., SeljebotnD. S., Smith K., 2011, Computing in Science & Engineer-ing, 13, 31Bialy S., Burkhart B., 2020, ApJ, 894, L2Boldyrev S., 2006, Phys. Rev. Lett., 96, 115002Bouchut F., Klingenberg C., Waagan K., 2010, NumerischeMathematik, 115, 647Burgers J., 1948, Advances in Applied Mechanics, 1, 171Burkhart B., 2018, ApJ, 863, 118Burkhart B., Lazarian A., 2012, ApJ, 755, L19Burkhart B., Falceta-Gon¸calves D., Kowal G., Lazarian A.,2009, ApJ, 693, 250Burkhart B., Lazarian A., Le˜ao I. C., de Medeiros J. R.,Esquivel A., 2014, ApJ, 790, 130Busquet G., 2020, Nature Astronomy, 4, 1126Childs H., et al., 2012, in , High Performance Visualization–Enabling Extreme-Scale Scientific Insight. Taylor & Fran-cis, pp 357–372Cho J., Lazarian A., Vishniac E. T., 2002, ApJ, 564, 291Cox N. L. J., et al., 2016, A&A, 590, A110Dubey A., et al., 2008, in Pogorelov N. V., Audit E., ZankG. P., eds, Astronomical Society of the Pacific Confer-ence Series Vol. 385, Numerical Modeling of Space PlasmaFlows. p. 145Federrath C., 2013, MNRAS, 436, 1245Federrath C., 2016a, Journal of Plasma Physics, 82,535820601Federrath C., 2016b, MNRAS, 457, 375Federrath C., 2018, Physics Today, 71, 38Federrath C., Banerjee S., 2015, MNRAS, 448, 3297Federrath C., Klessen R. S., 2012, ApJ, 761Federrath C., Klessen R. S., 2013, ApJ, 763, 51Federrath C., Klessen R. S., Schmidt W., 2008, ApJ, 688,L79Federrath C., Klessen R. S., Schmidt W., 2009, ApJ, 692,364Federrath C., Roman-Duval J., Klessen R., Schmidt W.,Mac Low M. M., 2010, A&A, 512Federrath C., et al., 2016, ApJ, 832, 143Federrath C., Klessen R. S., Iapichino L., Beattie J. R.,2021, Nature AstronomyFryxell B., et al., 2000, ApJS, 131, 273Ginsburg A., Federrath C., Darling J., 2013, ApJ, 779Goldreich P., Sridhar S., 1995, ApJ, 438, 763G´omez G. C., V´azquez-Semadeni E., Zamora-Avil´es M.,2018, MNRAS, 480, 2939Harris C. R., et al., 2020, Nature, 585, 357Hennebelle P., Chabrier G., 2008, ApJ, 684, 395 Hennebelle P., Chabrier G., 2009, ApJ, 702, 1428Hennebelle P., Inutsuka S.-i., 2019, Frontiers in Astronomyand Space Sciences, 6, 5Hennebelle P., Commer¸con B., Joos M., Klessen R. S.,Krumholz M., Tan J. C., Teyssier R., 2011, A&A, 528,A72Heyer M., Soler J. D., Burkhart B., 2020, MNRAS, 496,4546Hopkins P. F., 2013, MNRAS, 430, 1653Hu Y., et al., 2019, Nature Astronomy, 3, 776Hu Y., Lazarian A., Bialy S., 2020, ApJ, 905, 129Hunter J. D., 2007, Computing in Science & Engineering,9, 90Kainulainen J., Federrath C., Henning T., 2014, Science,344, 183Kim J., Ryu D., 2005, ApJ, 630, L45Klessen R. S., Glover S. C. O., 2016, Saas-Fee AdvancedCourse, 43, 85Klessen R. S., Heitsch F., Mac Low M.-M., 2000, ApJ, 535,887Konstandin L., Federrath C., Klessen R. S., Schmidt W.,2012a, Journal of Fluid Mechanics, 692, 183Konstandin L., Girichidis P., Federrath C., Klessen R. S.,2012b, ApJ, 761, 149Konstandin L., Schmidt W., Girichidis P., Peters T., ShettyR., Klessen R. S., 2016, MNRAS, 460, 4483K¨ortgen B., Soler J. D., 2020, MNRAS, 499, 4785Krumholz M. R., Federrath C., 2019, Frontiers in Astron-omy and Space Sciences, 6, 7Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250Landau L., Lifshitz E., 1959, Fluid Mechanics: Landauand Lifshitz: Course of Theoretical Physics. Butterworth-HeinemannLaw C. Y., Li H. B., Leung P. K., 2019, MNRAS, 484, 3604Law C. Y., Li H. B., Cao Z., Ng C. Y., 2020, MNRAS, 498,850Lee Y.-N., Hennebelle P., 2019, A&A, 622, A125Lee Y.-N., Offner S. S. R., Hennebelle P., Andr´e P., Zin-necker H., Ballesteros-Paredes J., Inutsuka S.-i., KruijssenJ. M. D., 2020, Space Science Reviews, 216, 70Lehmann A., Federrath C., Wardle M., 2016, MNRAS, 463,1026Li H.-B., Henning T., 2011, Nature, 479, 499Li H.-b., Fang M., Henning T., Kainulainen J., 2013, MN-RAS, 436, 3707Lunttila T., Padoan P., Juvela M., Nordlund ˚A., 2008, ApJ,686, L91Mac Low M. M., Klessen R. S., 2004, Reviews of ModernPhysics, 76, 125Malinen J., et al., 2016, MNRAS, 460, 1934Menon S. H., Federrath C., Klaassen P., Kuiper R., ReiterM., 2021, MNRAS, 500, 1721Mocz P., Burkhart B., 2018, MNRAS, 480, 3916Mocz P., Burkhart B., 2019, ApJ, 884, L35Mocz P., Burkhart B., Hernquist L., McKee C. F., SpringelV., 2017, ApJ, 838, 40Mohapatra R., Federrath C., Sharma P., 2020a, arXiv e-prints, p. arXiv:2010.12602Mohapatra R., Federrath C., Sharma P., 2020b, MNRAS,493, 5838Molina F. Z., Glover S. C. O., Federrath C., Klessen R. S.,2012, MNRAS, 423, 2680
MNRAS , 1–16 (2020) Beattie, Mocz, Federrath & Klessen
Nolan C. A., Federrath C., Sutherland R. S., 2015, MN-RAS, 451, 1380Oliphant T., 2006, NumPy: A guide to NumPy, USA: Trel-gol Publishing,
Orkisz J. H., et al., 2017, A&A, 599, A99Padoan P., Nordlund ˚A., 2002, ApJ, 576, 870Padoan P., Nordlund ˚A., 2011, ApJ, 730, 40Padoan P., Nordlund P., Jones B. J. T., 1997, Commmu-nications of the Konkoly Observatory Hungary, 100, 341Palmeirim P., et al., 2013, A&A, 550, A38Park J., Ryu D., 2019, ApJ, 875, 2Passot T., V´azquez-Semadeni E., 1998, Phys. Rev. E, 58,4501Pillai T., Kauffmann J., Tan J. C., Goldsmith P. F., CareyS. J., Menten K. M., 2015, The Astrophysical Journal,799, 74Pillai T. G. S., et al., 2020, Nature AstronomyPlanck Collaboration et al., 2016a, A&A, 586, A137Planck Collaboration et al., 2016b, A&A, 586, A138Price D. J., Federrath C., 2010, MNRAS, 406, 1659Price D. J., Federrath C., Brunt C. M., 2011, ApJ, 727,1380Robertson B., Goldreich P., 2018, ApJ, 854, 88Schneider N., et al., 2012, A&A, 540, L11Schneider N., et al., 2013, ApJ, 766, L17Seifried D., Walch S., Weis M., Reissl S., Soler J. D.,Klessen R. S., Joshi P. R., 2020, arXiv e-prints, p.arXiv:2003.00017Skalidis R., Tassis K., 2020, arXiv e-prints, p.arXiv:2010.15141Soler J. D., Hennebelle P., 2017, A&A, 607, A2Soler J. D., Hennebelle P., Martin P. G., Miville-DeschˆenesM. A., Netterfield C. B., Fissel L. M., 2013, ApJ, 774, 128Soler J. D., et al., 2017, A&A, 603, A64Stroustrup B., 2013, The C++ Programming Language,4th edn. Addison-Wesley ProfessionalTritsis A., Tassis K., 2016, MNRAS, 462, 3602Tritsis A., Tassis K., 2018, Science, 360, 635Tritsis A., Federrath C., Schneider N., Tassis K., 2018, MN-RAS, 481, 5275Troland T. H., Crutcher R. M., 2008, ApJ, 680, 457Vazquez-Semadeni E., 1994, ApJ, 423, 681Virtanen P., et al., 2020, Nature Methods, 17, 261Waagan K., Federrath C., Klingenberg C., 2011, Journal ofComputational Physics, 230, 3331Xu S., Ji S., Lazarian A., 2019, ApJ, 878, 157van der Walt S., et al., 2014, PeerJ, 2, e453
APPENDIX A: FLOW ORIENTATION IN THESUB-ALFV´ENIC, SUPERSONIC REGIME
In the sub-Alfv´enic regime large-scale vortices aligned with B form (Beattie et al. 2020), which arranges the magnetic π π π πθ . . . . . p ( θ | M ) M = 2 M = 4 M = 10 M = 20 Figure A1.
The distribution of the angle, θ , between the local mag-netic and velocity field for the ensemble of M A0 = 0 . ≤ t T ≤ and velocity fields such that on average v ⊥ B , i.e., (cid:104) θ (cid:105) = (cid:104) arccos [( v · B ) / ( (cid:107) v (cid:107)(cid:107) B (cid:107) )] (cid:105) ≈ π/
2, as opposed to isotropicturbulence where θ follows a uniform distribution. We showthe full distribution of θ averaged over 5 ≤ t T ≤
7, forour M A0 = 0 . θ = π/ B intertwined with intermittent events perturbing θ awayfrom the mean. The density variance model that we pro-pose in this study applies to this limiting case, where type Ishocks form the intermittent events that perturb θ , and thetype II shocks form by the magnetic field lines being draggedthrough the supersonic, vortical flow. This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000