Measuring magnetization with rotation measures and velocity centroids in supersonic MHD turbulence
DD RAFT VERSION F EBRUARY
11, 2021Typeset using L A TEX preprint2 style in AASTeX63
Measuring magnetization with rotation measures and velocity centroids in supersonic MHD turbulence S IYAO X U AND Y UE H U
2, 31
Institute for Advanced Study, 1 Einstein Drive, Princeton, NJ 08540, USA a2 Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA
ABSTRACTThe interstellar turbulence is magnetized and thus anisotropic. The anisotropy of turbulent magnetic fieldsand velocities is imprinted in the related observables, rotation measures (RMs) and velocity centroids (VCs).This anisotropy provides valuable information on both direction and strength of the magnetic field. However,its measurement is difficult especially in highly supersonic turbulence in cold interstellar phases due to thedistortions by isotropic density fluctuations. By using 3D simulations of supersonic and sub-Alfv´enic magne-tohydrodynamic (MHD) turbulence, we find that the problem can be alleviated when we selectively sample thevolume-filling low-density regions in supersonic MHD turbulence. Our results show that in these low-densityregions, the anisotropy of RM and VC fluctuations depends on the Alfv´enic Mach number as M − / . Thisanisotropy- M A relation is theoretically expected for sub-Alfv´enic MHD turbulence and confirmed by our syn-thetic observations of CO emission. It provides a new method for measuring the plane-of-the-sky magneticfields in cold interstellar phases. INTRODUCTIONThe interstellar medium (ISM) is both turbulent and mag-netized, and turbulence and magnetic fields are dynamicallycoupled together. Understanding their properties is cru-cial for studying many multi-scale physical processes in themulti-phase ISM, including star formation, cosmic ray prop-agation, and turbulent dynamo (Elmegreen & Scalo 2004;McKee & Ostriker 2007).Despite their significance, measuring interstellar turbu-lence and magnetic fields is difficult. One of the main rea-sons for the difficulty is the involvement of densities in ob-servables of turbulent velocities and magnetic fields. Sometechniques aimed at disentangling contributions of turbulentvelocities and densities, such as the Velocity Channel Anal-ysis (VCA) (Lazarian & Pogosyan 2000), the Velocity Co-ordinate Spectrum (VCS) (Lazarian & Pogosyan 2006), havebeen developed and applied to extracting statistical propertiesof turbulent velocities in the ISM from spectroscopic obser-vations (Lazarian 2009; Chepurnov et al. 2010).Turbulence in the presence of magnetic fields isanisotropic. The turbulent energy cascade mainly occurs inthe direction perpendicular to the local magnetic field, re-sulting in larger velocity fluctuations in the perpendiculardirection (Goldreich & Sridhar 1995; Lazarian & Vishniac1999). With the decrease of turbulent velocity along theturbulent energy cascade, the anisotropy becomes more pro- [email protected]; [email protected] a Hubble Fellow The importance of the reference frame of the local magnetic field for MHDturbulence was first discussed in Lazarian & Vishniac (1999). nounced at smaller scales. The scale-dependent anisotropybecomes nearly scale independent when the measurement isperformed in the global reference system with respect to themean magnetic field, as demonstrated by both simulations(Cho & Vishniac 2000; Maron & Goldreich 2001; Cho et al.2002; Beresnyak 2015) and observations of the solar wind(e.g., Matthaeus et al. 1990; Luo & Wu 2010; Wicks et al.2011). Only the anisotropy in the global reference frame isaccessible to observations of the ISM, which is the anisotropyaveraged along the line-of-sight (LOS). While both turbu-lent velocities and magnetic fields are anisotropic (Lazar-ian & Pogosyan 2012), the density fluctuations in highly su-personic MHD turbulence generated by shock compressionare isotropic (Cho & Lazarian 2004; Beresnyak et al. 2005).Consequently, retrieving turbulence anisotropy from observ-ables involving densities in highly supersonic turbulence incold interstellar phases, e.g., molecular clouds, is very chal-lenging.Based on the theoretical understanding on anisotropy ofturbulent velocities, new techniques, such as the Veloc-ity Gradients Technique (VGT, e.g., Gonz´alez-Casanova &Lazarian 2017; Yuen & Lazarian 2017; Hu et al. 2018), andthe Principal Component Analysis (e.g., Heyer et al. 2008)have been introduced for measuring interstellar turbulenceand magnetic fields. Statistical studies of velocity centroids(VCs, e.g., Lazarian & Esquivel 2003; Esquivel & Lazarian2005; Burkhart et al. 2014; Kandel et al. 2017) show that VCscan be used to probe the turbulence anisotropy and the plane-of-the-sky (POS) magnetic fields. However, VC fluctuationsin highly supersonic MHD turbulence are dominated by den-sity fluctuations and are unable to fully reveal the anisotropy a r X i v : . [ a s t r o - ph . GA ] F e b X U & H U of turbulent velocities (Hu et al. 2020a). The specific depen-dence of VC anisotorpy on magnetic field strength is unclear.Similarly, the commonly used observables for tracing mag-netic fields, e.g., Faraday rotation measures (RMs), are alsosubject to the distortions by densities. Statistical analysis ofinterstellar RM fluctuations (Minter & Spangler 1996) sug-gests that they are dominated by density fluctuations (Xu &Zhang 2016), which exhibit a shallow density spectrum dueto the small-scale density enhancements arising from the su-personic turbulence in cold interstellar phases. This findingon a shallow density spectrum is also confirmed by the tem-poral broadening measurements (Xu & Zhang 2017) and dis-persion measures (Xu & Zhang 2020) of pulsars. This find-ing is also consistent with the shallow density spectra mea-sured with gas and dust tracers in the cold neutral mediumand molecular clouds (Hennebelle & Falgarone 2012), show-ing the coupling between ionized and neutral components ofgas in the partially ionized ISM (Xu et al. 2015, 2016).Hu et al. (2020c) investigated the measurement on orienta-tion and strength of magnetic fields from anisotropic turbu-lent velocities, which can be directly obtained by using pointtracers of turbulence, including dense cores (Qian et al. 2018)and young stars (Ha et al. 2021). The dependence of turbu-lence anisotropy on the Alfv´enic Mach number M A , whichis related to the magnetic field strength, was analytically de-rived and numerically tested by Hu et al. (2020c). It leads toa new method for measuring magnetic fields in sub-Alfv´enicturbulence ( M A < ) in the ISM. In this work, we considermore commonly used observables for tracing turbulent mag-netic fields and velocities, i.e., RMs and VCs, which are sub-ject to both projection effect and distortions by density fluc-tuations. By exploring a new analysis method to mitigate theeffect of densities and maximize the information on turbu-lence anisotropy, we will reexamine their applicability to su-personic MHD turbulence in retrieving turbulence anisotropyand measuring magnetization. We will also test our results onVCs by using synthetic observations including optical deptheffects.In Section 2, we present the theoretical formulations forstructure functions (SFs) of turbulent fluctuations, which area basic statistical tool for analyzing turbulence and quantify-ing turbulence anisotropy. In Section 3, we numerically studythe anisotropy of SFs of RMs and VCs by using a set of MHDsimulations of supersonic and sub-Alfv´enic turbulence andproduce synthetic observations for testing our results. Thediscussion and conclusions can be found in Section 4 andSection 5. SF ANALYSIS IN MHD TURBULENCE2.1.
Anisotropy of MHD turbulence
Magnetized turbulence is anisotropic, and the anisotropybecomes more pronounced toward smaller scales with theturbulent energy cascade. Here we focus on the strong turbu-lence regime where the critical balance between the turbulenteddy-turnover time and the Alfv´en wave period is established(Goldreich & Sridhar 1995). Following the turbulent energy cascade, the Kolmogorovscaling for hydrodynamic turbulence still applies to MHDturbulence in the direction perpendicular to the local mag-netic field, as numerically demonstrated with high-resolutionMHD turbulence simulations in, e.g., Beresnyak (2014).Therefore, the turbulent velocity at the length scale l ⊥ per-pendicular to the local magnetic field is given by v ( l ⊥ ) = V st (cid:16) l ⊥ L st (cid:17) , (1)where V st is the turbulent velocity at L st . For sub-Alfv´enicturbulence with the turbulent energy smaller than the mag-netic energy, the Alfv´enic Mach number M A = V L /V A isless than unity, where V L is the injected turbulent velocity atthe injection scale L i , V A = B/ √ πρ is the Alfv´en speed, B is the magnetic field strength, and ρ is the density. Thereis (Lazarian 2006) V st = V L M A , L st = L i M A , (2)and strong MHD turbulence exists on scales smaller than L st .The scale-dependent turbulence anisotropy is described by l (cid:107) = V A V st L st l ⊥ , (3)which can be derived by combining Eq. (1) with the crit-ical balance relation l ⊥ /v l = l (cid:107) /V A , where l (cid:107) is the par-allel length scale measured along the local magnetic field.It shows that the turbulence becomes more anisotropic atsmaller scales.Given the relation in Eq. (3), one can express v in Eq. (1)in terms of l (cid:107) , v ( l (cid:107) ) = V st (cid:16) V st V A (cid:17) (cid:16) l (cid:107) L st (cid:17) . (4)Comparing the above expression with Eq. (1), when v ( l ⊥ ) and v ( l (cid:107) ) are measured at the same length, i.e., l ⊥ = l (cid:107) , wedefine the anisotropy degree (AD) in the system of referencealigned with the local magnetic field asAD loc = v ( l ⊥ ) v ( l (cid:107) ) = V A V st (cid:16) l ⊥ L st (cid:17) − . (5)Inserting Eq. (2) into it yieldsAD loc = M − A (cid:16) l ⊥ L (cid:17) − , ( M A < , (6)which increases with decreasing l ⊥ .The above scale-dependent anisotropy can only be mea-sured in the local reference frame aligned with the locally av-eraged magnetic field direction (Lazarian & Vishniac 1999;Cho & Vishniac 2000; Cho et al. 2003). Observations in theISM are subjected to the projection effect, and only the globalreference frame of the mean magnetic field averaged alongthe LOS is accessible. The measured anisotropy is dominated EASURING MAGNETIZATION WITH ROTATION MEASURES AND VELOCITY CENTROIDS L st is given by(Eq. (3)) L st , (cid:107) ≈ V A V st L st , ⊥ , (7)where we assume L st ≈ L st , ⊥ for anisotropic turbulence. Thescale-independent anisotropy in the global frame is l (cid:107) l ⊥ = L st , (cid:107) L st , ⊥ ≈ V A V st . (8)Then Eq. (1) can be expressed in terms of l (cid:107) , v ( l (cid:107) ) ≈ V st (cid:16) l (cid:107) L st , (cid:107) (cid:17) . (9)The AD in the global frame is the ratio (Eqs. (1) and (9))AD glo = v ( l ⊥ ) v ( l (cid:107) ) ≈ (cid:16) L st , (cid:107) L st (cid:17) ≈ (cid:16) V A V st (cid:17) (10)measured at l ⊥ = l (cid:107) , where l ⊥ and l (cid:107) are the perpendicularand parallel length scales measured with respect to the direc-tion of the mean magnetic field. By using Eq. (2), we furtherget AD glo ≈ M − A , ( M A < . (11)Both AD loc and AD glo have the same dependence on M A andhave larger values for more strongly magnetized turbulencewith a smaller M A . But in contrast to AD loc , AD glo is scale-independent.The above anisotropic scaling for the turbulent velocities inincompressible MHD turbulence is also applicable to com-pressible MHD turbulence, when the turbulent motions aregoverned by Alfv´en modes (Cho & Lazarian 2002, 2003).2.2. Anisotropy of SFs of turbulent fluctuations
For the statistical analysis of turbulent fluctuations carriedout below, we consider the statistics in the global frame ofreference as this is the only available statistics for observa-tions of the ISM. 2.2.1.
3D SF
For the statistical measurement of the fluctuations inducedby turbulence, when the three-dimensional (3D) positions areavailable, we have the 3D correlation function (CF), ξ v ( R, ∆ z ) = (cid:104) v ( X , z ) v ( X , z ) (cid:105) = (cid:104) v (cid:105) L m st L m st + ( R + ∆ z ) m , (12)and 3D SF d v ( R, ∆ z ) = (cid:104) [ v ( X , z ) − v ( X , z )] (cid:105) = 2 (cid:104) v (cid:105) ( R + ∆ z ) m L m st + ( R + ∆ z ) m . (13) We adopt the power-law models for the CF and SF (Lazarian& Pogosyan 2016), with the power-law index m = 2 / cor-responding to the Kolmogorov scaling of turbulence. Theyare applied to turbulent velocities in the above expressionsas an example. Here X and z are the projected positionvector on the plane of sky and the distance along the LOS, R = | X − X | , ∆ z = z − z , and (cid:104) ... (cid:105) denotes the aver-age over all pairs of points with the same 3D separation.The AD of the 3D SF of turbulent velocities is given byEq. (11), d v ( r ⊥ ) d v ( r (cid:107) ) = AD glo ≈ M − A , (14)where we use r ⊥ and r (cid:107) to represent the separations mea-sured in the directions perpendicular and parallel to the meanmagnetic field. The above relation has been numericallytested by Hu et al. (2020c). Their finding suggests that bymeasuring the 3D SF of turbulent velocities along differ-ent directions, one can determine both the orientation of themean magnetic field and the value of M A . The latter leadsto the evaluation of B when ρ is known. This method can beapplied to young stars with full 6D phase-space coordinatesprovided by Gaia, which have been recently found as a newprobe of turbulence in their parent molecular clouds (Ha et al.2021).For some other point tracers of turbulence, such as densecores in molecular clouds (Qian et al. 2012), only the LOScomponent of turbulent velocity v z and the 2D position onthe plane of sky are attainable. In this case, the 3D SF isaveraged over all pairs with the same R but different ∆ z , d v ( R ) = (cid:82) L d v ( R, ∆ z ) d ∆ z (cid:82) L d ∆ z = 2 (cid:104) v z (cid:105) L (cid:90) L ( R + ∆ z ) m L m st + ( R + ∆ z ) m d ∆ z, (15)where we assume that the point sources are distributed uni-formly along the LOS from z = 0 to z = L , and L is thethickness of the sampled volume. As illustrated in Fig. 1(a),when the turbulent volume is thin with L (cid:28) L st , there is d v ( R ) ≈ (cid:104) v z (cid:105) m + 1 (cid:16) LL st (cid:17) m , R < L,d v ( R ) ≈ (cid:104) v z (cid:105) (cid:16) RL st (cid:17) m , L < R < L st ,d v ( R ) ≈ (cid:104) v z (cid:105) , R > L st . (16)Both the scaling and anisotropy properties of turbulence canbe retrieved from d v ( R ) within the range L < R < L st (seeHu et al. 2020c). When the LOS is perpendicular to the meanmagnetic field, we have d v ( R ⊥ ) d v ( R (cid:107) ) = AD glob ≈ M − A , (17)where R ⊥ and R (cid:107) are measured with respect to the directionof the mean magnetic field. X U & H U However, when the turbulent volume has the LOS thick-ness much larger than the correlation length L st , there is al-ways d v ( R ) ≈ (cid:104) v z (cid:105) . (18)The properties of turbulence cannot be retrieved (seeFig. 1(b) and also Hu et al. 2020c). This result can be poten-tially used for constraining the thickness of molecular clouds(Qian et al. 2015). 2.2.2.
2D SF
When using spatially continuous gas tracers of turbulence,we usually need to deal with projected quantities, i.e., inte-grated quantities along the LOS, and their SFs. In the case ofturbulent velocities, we have the 2D SF of projected veloci-ties as (Lazarian & Pogosyan 2016) D v ( R ) = (cid:104) [ v z ( X ) − v z ( X )] (cid:105) = (cid:68)(cid:104) (cid:90) L dzv z ( X , z ) − (cid:90) L dzv z ( X , z ) (cid:105) (cid:69) = 4 (cid:90) L d ∆ z ( L − ∆ z )[ ξ v (0 , ∆ z ) − ξ v ( R, ∆ z )]= 4 (cid:104) v z (cid:105) (cid:90) L d ∆ z ( L − ∆ z ) (cid:34) L m st L m st + ∆ z m − L m st L m st + ( R + ∆ z ) m (cid:35) , (19)where the LOS is in the z direction, and L is the thicknessof the turbulent gas cloud. For a thick cloud with L largerthan L st and for Kolmogorov scaling of turbulent velocities, D v ( R ) can be approximated by D v ( R ) ≈ (cid:104) v z (cid:105) (cid:90) R d ∆ zL R m L m st = 4 (cid:104) v z (cid:105) L − m st LR m +1 (20)within the range R < L st . Compared with 3D SF in Eq. (13), D v ( R ) has a steeper power-law dependence on R because ofthe projection effect, as shown in observations (Lazarian &Pogosyan 2000; Elmegreen et al. 2001; Padoan et al. 2001).The detailed analysis on the anisotropy of D v ( R ) of dif-ferent MHD modes and its relation to the magnetization andthe angle between LOS and mean magnetic field was carriedout by Kandel et al. (2017). Here as a rough estimate, werewrite D v in terms of R ⊥ and R (cid:107) , which are the 2D pro-jected separations measured perpendicular and parallel to themean magnetic field, and have D v ( R ⊥ ) ≈ (cid:104) v z (cid:105) (cid:90) R ⊥ d ∆ zL R m ⊥ L m st = 4 (cid:104) v z (cid:105) L − m st LR m +1 ⊥ , (21)and D v ( R (cid:107) ) ≈ (cid:104) v z (cid:105) (cid:90) R (cid:107) d ∆ zL R m (cid:107) L m st , (cid:107) = 4 (cid:104) v z (cid:105) L − m st , (cid:107) LR m +1 (cid:107) , (22) where we assume that the mean magnetic field is perpendic-ular to the LOS for simplicity. Then the AD measured with2D SFs at R ⊥ = R (cid:107) can be estimated as D v ( R ⊥ ) D v ( R (cid:107) ) ≈ (cid:16) L st , (cid:107) L st (cid:17) m = M − mA , (23)which has the same dependence on M A as AD glo measuredwith 3D SFs (Eq. (14)) for Kolmogorov scaling of turbulentvelocities.The above analysis on anisotropy of SFs is exemplifiedby turbulent velocities as they directly reflect the dynam-ics of MHD turbulence and are the best indicator of tur-bulence anisotropy. The anisotropy of turbulent magneticfields and its application to analyzing the statistics of syn-chrotron fluctuations have been comprehensively studied byLazarian & Pogosyan (2012). Magnetic fields are mainlyperturbed by Alfv´enic turbulent motions even in highly su-personic turbulence and thus follow similar anisotropic prop-erties as turbulent velocities, as shown in numerical simula-tions with solenoidal driving (Beresnyak et al. 2005). Den-sity fluctuations in subsonic and mildly supersonic MHD tur-bulence are passively mixed by Alfv´enic motions and havethe same anisotropy as turbulent velocities (Lithwick & Gol-dreich 2001; Cho & Lazarian 2003; Xu et al. 2019). While inhighly supersonic turbulence, high density contrast is gener-ated due to shock compression, and density fluctuations be-come isotropic (Cho & Lazarian 2004; Beresnyak et al. 2005;Hu et al. 2020a,b).Magnetic fields and velocities are usually coupled withdensities in observational measurements, and the extractionof them is very challenging especially in highly supersonicturbulence. The commonly used observable for studyingmagnetic fields is RM, which is defined asRM = e πm e c (cid:90) L dzn e B z , (24)where e and m e are the charge and mass of an electron, c isthe light speed, n e is the number density of electrons, and B z is the LOS component of magnetic field. Turbulent mag-netic fields and densities both contribute to the fluctuations ofRMs, but the latter can play a dominant role due to the largedensity variation arising from supersonic turbulence. The re-sulting SF of RMs in the ISM measured by, e.g., Minter &Spangler (1996), can be explained by a shallow density spec-trum expected for supersonic turbulence in cold interstellarphases (Xu & Zhang 2016).The VC used for measuring turbulent velocities can be de-fined as VC = (cid:15) (cid:90) L dzρv z (25)in position-position-position space for optically thin emis-sion lines (Kandel et al. 2017), where (cid:15) is the emissivity co-efficient, and ρ is the real space density. Similar to RMs,VCs do not exhibit the scaling properties of the underlyingturbulent velocities in highly supersonic turbulence due to EASURING MAGNETIZATION WITH ROTATION MEASURES AND VELOCITY CENTROIDS R/L -2 -1 d v ( R )2 h v z i -1 R m L st /L (a) L (cid:28) L st R/L st -3 -2 -1 d v ( R )2 h v z i -1 L/L st (b) L (cid:29) L st Figure 1.
Projected 3D SF d v ( R ) (normalized by (cid:104) v z (cid:105) ) of turbulent velocities sampled by point sources in the (a) thin and (b) thick turbulentvolumes. the density distortions (Esquivel & Lazarian 2005; Esquivelet al. 2007).The SFs of RMs and VCs are D RM ( R ) = (cid:104) [ RM ( X ) − RM ( X )] (cid:105) , (26)and D VC ( R ) = (cid:104) [ VC ( X ) − VC ( X )] (cid:105) , (27)respectively. We will next carry out numerical simulationsand investigate whether they can be used to recover the tur-bulence anisotropy in supersonic and sub-Alfv´enic MHD tur-bulence. ANISOTROPY OF SFS OF RMS AND VCS3.1.
Simulations of supersonic and sub-Alfv´enic MHDturbulence
Our scale-free 3D MHD simulations are generated usingthe ZEUS-MP/3D code (Hayes et al. 2006), which solvesthe ideal MHD equations with zero-divergence condition ∇ · (cid:126)B = 0 and an isothermal equation of state in a peri-odic box. Here (cid:126)B is the total magnetic field and the boxis regularly staggered to 792 grid cells. We consider sin-gle fluid and operator-split MHD conditions in the Eulerianframe. The turbulence is solenoidally injected in Fourierspace at wavenumber k i ≈
2, and the dissipation occurs near k d ≈ (see Hu et al. 2020a). The simulations are initial-ized with a uniform density and a uniform magnetic field (cid:126)B .The magnetic field consists of (cid:126)B and a fluctuating compo-nent (cid:126)b , so (cid:126)B = (cid:126)B + (cid:126)b . Initially we set the magnetic fieldalong the x and z axes to be zeros, and the uniform magneticfield is along the y -axis. We choose z -axis to be the LOS sothat the mean magnetic field is perpendicular to the LOS.We consider supersonic and sub-Alfv´enic MHD turbulencefor studying the anisotropy of supersonic turbulence in coldmolecular clouds (MCs) with the sonic Mach number M S = V L /c s ≈ − (Zuckerman & Palmer 1974; Larson 1981),where c s is the sound speed. The values of M S and M A of our numerical models are listed in Table 1. For comparison, wealso include one model (M0) with subsonic and sub-Alfv´enicturbulence.Simulations by Padoan et al. (2016) suggest that the su-pernova diving is able to sustain the turbulence observed inMCs, without invoking other sources of turbulence. The su-pernova energy injection brings both solenoidal and com-pressive modes, but the compressive-to-solenoidal ratio de-creases due to the conversion from compressive to solenoidalmotions along the energy cascade, leading to the dominanceof solenoidal motions at MC scales. This finding is consis-tent with the observations showing the Kolmogorov scalingof turbulent velocities (Qian et al. 2018; Xu 2020) and thedominance of solenoidal motions (Orkisz et al. 2017) at MCscales. Motivated by these numerical and observational re-sults, we adopt the solenoidal driving for studying turbulencein MCs.For our measurements on anisotropic SFs, we select a sam-ple of points for the 3D SFs and lines of sight for the 2D SFsfrom our data cube, and the separations are binned evenlywith the bin width approximately equal to two grid cells onscales smaller than L i ( ≈ /k i = 396 grid units). In the2D case, we use a sample size of . For both the SFs mea-sured perpendicular and parallel to the mean magnetic field,the number of pairs in each separation bin is on the order of − , which is sufficiently large for the SF measurementto be statistically stable (see Hu et al. 2020c). The uncertain-ties in the measured anisotropy of SFs related to the samplesize will be presented in Sections 3.2.2 and 3.3. For the mea-surement on 3D anisotropic SFs, as it requires a larger sam-ple than the 2D case, we use a sample size of to have thenumber of pairs in each bin larger than . To investigatethe effect of density on the anisotropy of RM and VC SFs,we first randomly sample the turbulence volume in Section3.2.1, and then selectively sample the low-density regions inSection 3.2.2. 3.2. Numerical results X U & H U Table 1. M S and M A values of our MHD simulations. B is the initial magnetic field strength corresponding tomean mass density (cid:104) ρ (cid:105) = 300 g cm − and isothermal sound speed c s = 187 m s − .Model M0 M1 M2 M3 M4 M5 M6 M7 M8 M9 M S .
62 7 .
31 6 .
10 6 .
47 6 .
14 6 .
03 10 .
81 10 .
53 10 .
61 10 . A .
56 0 .
22 0 .
42 0 .
61 0 .
82 1 .
01 0 .
26 0 .
51 0 .
68 0 . B [ µG ] 2.5 72 32 22 16 13 90 45 35 25 Comparisons between supersonic and subsonic MHDturbulence
Fig. 2 presents the 3D SFs of turbulent magnetic fields d b ,velocities d v , and densities d ρ in supersonic and subsonicMHD turbulence with similar M A values. d ⊥ and d (cid:107) rep-resent d ( r ⊥ ) and d ( r (cid:107) ) for simplicity. All units used in ourmeasurements here and below are numerical units. From the3D SFs, we see that the numerical dissipation effect sets in at /k d ≈ grid units, below which the SF steepens due todissipation. The difference between the SFs measured in thedirections perpendicular and parallel to the mean magneticfield can be clearly seen in all cases for the subsonic MHDturbulence. The fake dependence of the anisotropy on lengthscales appears due to the isotropic driving and insufficientinertial range (see Cho & Vishniac 2000; Hu et al. 2020c).By contrast, we see in Fig. 2(a) that the SF of density fluctu-ations in highly supersonic turbulence is nearly isotropic asdiscussed in Section 2.2.2. The different density structures inthe supersonic and subsonic MHD turbulence are illustratedin Fig. 3. Small-scale density enhancements with large ρ/ (cid:104) ρ (cid:105) appear in the supersonic MHD turbulence, where (cid:104) ρ (cid:105) is themean density, in comparison with the more uniform densitydistribution in the subsonic case.The anisotropy of 3D SFs was earlier studied by Hu et al.(2020c). In this work we are interested in the anisotropy of2D SFs. As an example, the 2D SFs of the projected mag-netic fields (cid:82) L dzB z ( X , z ) , velocities (cid:82) L dzv z ( X , z ) (Eq.(19)), and densities (cid:82) L dzρ ( X , z ) for the supersonic MHDturbulence (M7) are presented in Fig. 4(a). We use D ⊥ and D (cid:107) to represent D ( R ⊥ ) and D ( R (cid:107) ) for simplicity. Due tothe projection effect, there is D ( R ) ∝ R / (Eq. (20)) forthe Kolmogorov scaling. Due to the solenoidal driving, theKolmogorov scaling for turbulent velocities is expected evenin supersonic MHD turbulence (e.g., Kowal et al. 2007). D ρ distinctively exhibits a shallow slope and isotropy comparedwith D b and D v . Both the shallowness of density spectrumand isotropy of density fluctuations in supersonic MHD tur-bulence were found in earlier studies, e.g., Beresnyak et al.(2005). They account for the shallow slope and suppressedanisotropy seen for D RM and D VC in Fig. 4(b). As expected,the fluctuations of RMs and VCs are dominated by densityfluctuations in highly supersonic turbulence, and thus theanisotropic scalings of turbulent magnetic fields and veloc-ities cannot be fully recovered. For our simulated data, RM and VC are measured asRM = (cid:90) L dzρB z , (28)and VC = (cid:90) L dzρv z . (29)The above approximation for RM is valid when the ioniza-tion fraction can be treated as a constant for the sampled tur-bulence volume.To mitigate the effect of density, we also examine the SFsof normalized RM D RM ( R ) = (cid:68)(cid:104) RM ( X ) DM ( X ) − RM ( X ) DM ( X ) (cid:105) (cid:69) , (30)and normalized VC D VC ( R ) = (cid:68)(cid:104) VC ( X ) DM ( X ) − VC ( X ) DM ( X ) (cid:105) (cid:69) , (31)as shown in Fig. 4(c), whereDM = (cid:90) L dzρ (32)is the column density. DM can be treated as dispersionmeasure (cid:82) L dzn e when the ionization fraction is a constant.We see that the division over DM slightly improves theperformance of RMs and VCs in retrieving the turbulenceanisotropy. To examine the effect of sample size, we alsopresent the measured D RM and D VC with a larger samplesize of in Fig. 4(d). Although increasing the sample sizeyields a smoother SF, the anisotropy is still insignificant.As a comparison, Fig. 5 displays the 2D SFs in subsonicMHD turbulence. For RMs and VCs, the involvement ofdensities does not affect the anisotropic scalings of mag-netic fields and velocities. The applicability of VCs to sub-sonic turbulence has been numerically studied by, e.g., Es-quivel & Lazarian (2005); Esquivel et al. (2007); Burkhartet al. (2014). The above results suggest that the SF of RMscan also be used to obtain turbulence anisotropy in subsonicMHD turbulence.We notice that the effect from isotropic driving and insuf-ficient inertial range is alleviated for 2D SFs. The inertialrange appears to be more extended due to the projection ef-fect, and the anisotropy on smaller scales becomes scale in-dependent. EASURING MAGNETIZATION WITH ROTATION MEASURES AND VELOCITY CENTROIDS r d (r) -2 -1 r (a) Supersonic (M7) r d (r) -6 -5 -4 -3 -2 d b, ⊥ d b,|| d v, ⊥ d v,|| d ρ , ⊥ d ρ ,|| r (b) Subsonic (M0) Figure 2. (a) 3D SFs of turbulent magnetic fields, velocities, and densities measured for run M7. (b) Same as (a) but for M0. The samelinestyles are used for both (a) and (b). The Kolmogorov scaling is indicated by the short solid line. (a) Supersonic (M7) (b) Subsonic (M0)
Figure 3.
Illustrations for the logarithmic density distribution (normalized by the mean density) in the supersonic MHD turbulence (run M7)and subsonic MHD turbulence (run M0). The direction of the mean magnetic field (cid:104) B (cid:105) is indicated by the arrow. Note that different colorscales are used in (a) and (b). Low-density regions in supersonic MHD turbulence
It is clear that the suppressed anisotropy seen from the SFsof RMs and VCs in supersonic MHD turbulence is causedby the isotropic density fluctuations. It was earlier foundin Beresnyak et al. (2005) that by filtering out high-densitypeaks with the logarithm of density used, the density statis-tics even in supersonic MHD turbulence exhibit the turbu-lence anisotropy similar to that in subsonic MHD turbulence.The small-scale density enhancements in supersonic tur-bulence are spatially clustered and concentrated with a smallvolume filling factor. The volume-filling densities are mainlyin the range /M S (cid:46) ρ/ (cid:104) ρ (cid:105) (cid:46) (Robertson & Goldre-ich 2018). To avoid the distortions by high-density peaks,we only select the lines of sight through the relatively low- density regions with DM < (cid:104) ρ (cid:105) L. (33)The resulting SFs of RMs and VCs in low-density regionsare presented in Figs. 6 and 7 for supersonic MHD turbu-lence with M S ∼ − and M S ∼ , respectively (seeTable 1). Compared with the case using random samplingin Fig. 4, the selective sampling adopted here leads to moreanisotropic SFs of both RMs and VCs in supersonic MHDturbulence over an extended range of length scales. The SFson large scales get deformed because of our selective sam-pling method. The incomplete sampling of the data cubecauses loss of information on large-scale turbulent fluctua-tions, but this does not affect the ansitropy measurement onsmall scales (see below). We note that the density effect is not X U & H U R D ( R ) -1 D b, ⊥ D b,|| D v, ⊥ D v,|| D ρ , ⊥ D ρ ,|| R (a) R D ( R ) D RM, ⊥ D RM,|| D VC, ⊥ D VC,|| R (b) R D ( R ) -3 -2 -1 D RM, ⊥ D RM, || D V C, ⊥ D V C, || R (c) R D ( R ) -3 -2 -1 D RM, ⊥ D RM, || D V C, ⊥ D V C, || R (d) Figure 4.
2D SFs for (a) projected magnetic fields, velocities, and densities, (b) RMs and VCs, (c) normalized RMs and VCs in supersonicMHD turbulence for run M7. The Kolmogorov scaling is indicated by the short solid line. (d) Same as (c) but with a larger sample size. fully suppressed, as D RM ( R ) and D VC ( R ) still show devia-tions from the Kolmogorov scaling, especially in the highlysupersonic case in Fig. 7.From Figs. 6 and 7, we easily see that the degree ofanisotropy varies in different numerical runs with differ-ent M A . To quantify the dependence of the anisotropiesof D RM and D VC on M A , we measure the ADs, i.e., D RM ( R ⊥ ) /D RM ( R (cid:107) ) and D VC ( R ⊥ ) /D VC ( R (cid:107) ) , averagedover a range of length scales from one to ten grid cells, wherethe anisotropy remains constant . The SFs on large scalesare subject to the effects of isotropic driving (Cho & Vish-niac 2000; Yuen et al. 2018) and incomplete sampling in thiscase, so they are not used for the AD measurement. The de-pendence of AD on M A is displayed in Fig. 8(a). In most Among our numerical runs, M1 has the smallest M A = 0 . (see Table1). The corresponding L st = L i M A (Eq. (2)) is around grid cells. Therange of length scales for measuring AD is still in the strong turbulenceregime. cases, the AD seen from D RM is slightly smaller than thatfrom D VC . This is expected, as magnetic fluctuations exhibita smaller anisotropy compared with velocity fluctuations insupersonic MHD turbulence, as shown from their 3D and2D SFs (Figs. 2(a) and 4(a)). The numerically measured M A dependence is consistent with M − / , which we ana-lytically derived for the Kolmogorov scaling of turbulencein Section 2.2.2. It shows that despite the deviations fromthe Kolmogorov slope that are still seen in D RM and D VC measured in low-density regions in supersonic MHD turbu-lence, the theoretically expected M A dependence of turbu-lence anisotropy can be recovered. Note that the non-unityAD at M A ≈ is caused by the presence of uniform mag-netic field of our initial setup. The dependence on M S isinsignificant.To further examine the effect of density on the measuredanisotropies of D RM and D VC , we also measure the ra-tio D v ( R ⊥ ) /D v ( R (cid:107) ) in low-density regions in supersonicMHD turbulence (see Fig. 9 as an example) at different M A . EASURING MAGNETIZATION WITH ROTATION MEASURES AND VELOCITY CENTROIDS R D ( R ) D b, ⊥ D b,|| D v, ⊥ D v,|| D ρ , ⊥ D ρ ,|| R (a) R D ( R ) -1 D RM, ⊥ D RM,|| D VC, ⊥ D VC,|| R (b) R D ( R ) -6 -5 -4 -3 D RM, ⊥ D RM, || D V C, ⊥ D V C, || R (c) Figure 5.
Same as Fig. 4 but for subsonic MHD turbulence of run M0.
Obviously, the AD obtained from projected turbulent veloci-ties alone is significant larger than those obtained from D RM and D VC at the same M A (see Fig. 8(a)). Compared withour simple scaling argument in Section 2.2.2, it has the M A -dependence better described by the more rigorous theoreticalresult in Kandel et al. (2017), which can be numerically ap-proximated by − / . It shows that in low-density regionsthe involvement of densities in D RM and D VC affects the ob-served AD, but not its dependence on M A . Since usually D v is not directly accessible to observations, our numericalresult can be used to estimate M A with RMs and VCs mea-sured in low-density regions of supersonic turbulence in coldinterstellar phases.As a comparison, we also present D RM ( R ⊥ ) /D RM ( R (cid:107) ) and D VC ( R ⊥ ) /D VC ( R (cid:107) ) measured with random samplingin supersonic MHD turbulence at different M A in Fig. 8(b).Compared with the measurements in low-density regions,the ADs are considerably smaller and have a much weakerdependence on M A . We use the standard errors of thebin averages to calculate the uncertainties in SF measure- ments related to the sample size. The error of the ratio D ( R ⊥ ) /D ( R (cid:107) ) is further determined with error propagation,represented by the error bars in Fig. 8(b). We note that theerror bars are not shown in Fig. 8(a) as they are smaller thanthe symbol size. For the selectively sampled low-density re-gions in smaller sub-volumes of the data cube, the numberof pairs at small separations is larger than that with randomsampling, leading to smaller uncertainties.3.3. Synthetic observations of VCs
To test the applicability of VCs to measuring the mag-netization in cold interstellar phases, we produce syntheticobservations of CO emission by performing radiative trans-fer calculations based on our supersonic MHD simulationsM1-M9. The synthetic emission lines of CO isotopologs,i.e., CO (1-0) and C O (1-0), are generated through theSPARX radiative transfer code (Hsieh et al. 2019). TheSPARX solves the radiative transfer equation through theAccelerated Lambda Iteration (ALI). The radiative interac-tion accounts for the molecular spontaneous emission, stim-0 X U & H U R D ( R ) -5 -4 -3 -2 -1 D V C, || D V C, ⊥ D RM, || R D RM, ⊥ (a) M1 R D ( R ) -4 -3 -2 -1 D RM, ⊥ D V C, ⊥ D RM, || D V C, || R (b) M2 R D ( R ) -5 -4 -3 -2 -1 D RM, ⊥ D RM, || D V C, ⊥ D V C, || R (c) M3 R D ( R ) -4 -3 -2 -1 D RM, ⊥ D RM, || D V C, ⊥ D V C, || R (d) M4 R D ( R ) -4 -3 -2 -1 D RM, || D V C, ⊥ D RM, ⊥ D V C, || R (e) M5 Figure 6.
SFs of selectively sampled RMs and VCs in low-density regions in supersonic MHD turbulence for runs M1-M5 with M S ∼ − . EASURING MAGNETIZATION WITH ROTATION MEASURES AND VELOCITY CENTROIDS R D ( R ) -4 -3 -2 -1 D RM, || D V C, ⊥ D V C, || R D RM, ⊥ (a) M6 R D ( R ) -3 -2 -1 D RM, ⊥ D RM, || D V C, ⊥ D V C, || R (b) M7 R D ( R ) -2 -1 D RM, ⊥ D RM, || D V C, || D V C, ⊥ R (c) M8 R D ( R ) -2 -1 R D RM, ⊥ D RM, || D V C, ⊥ D V C, || (d) M9 Figure 7.
Same as Fig. 6 but for runs M6-M9 with M S ∼ . U & H U M A D ( R ⊥ ) / D ( R k ) vKandel et al. 2017 RMV C (a) Selective sampling of low-density regions M A D ( R ⊥ ) / D ( R k ) RMV C (b) Random sampling
Figure 8. (a) D v ( R ⊥ ) /D v ( R (cid:107) ) , D RM ( R ⊥ ) /D RM ( R (cid:107) ) , and D VC ( R ⊥ ) /D VC ( R (cid:107) ) as a function of M A obtained in low-density regions ofsupersonic MHD turbulence for runs M1-M9. The open symbols correspond to runs M1-M5, and the filled symbols correspond to runs M6-M9. The dashed lines indicate M − / and − / . The solid line shows the analytical prediction by Kandel et al. (2017). (b) Same as (a) butfor D RM and D VC with random sampling. The dashed line corresponds to M − / . R D v ( R ) D v, ⊥ D v,|| R Figure 9.
SFs of projected velocities in low-density regions of su-personic MHD turbulence for run M7. ulated emission, and the collision with gas particles in the Lo-cal Thermodynamic Equilibrium (LTE) condition. FollowingHsieh et al. (2019), the fractional abundance of the CO iso-topolog CO (1-0) is set to × − , which comes fromthe cosmic value of C/H = 3 × − and the assumption that15% of C is in the molecular form. By adopting CO/C O= 588, we have the fractional abundance of C O (1-0) as . × − .The information about molecular gas density and velocityis taken from our MHD simulations (see Section 3.1). Foreach scale-free MHD simulation, we choose a cube size of pc, a gas temperature of K (which corresponds to soundspeed c s = 187 m s − ), and a mean mass density of g cm − . The initial magnetic field strength is then determinedbased on the M A value (see Table 1). The mean optical depth, which mainly depends on the molecular abundance and meandensity, is approximately 200 for optically thick CO (1-0)and 0.33 for optically thin C O (1-0), which is close to theobservational value (Lin et al. 2016).The synthetic VC map is calculated using:VC = (cid:82) T R · v z dv z (cid:82) T R dv z , (34)where T R is the radiation temperature, and v z is the LOSvelocity. We note that for optically thin emission lines, theabove definition of VC is equivalent to that used in Eq. (31)(Kandel et al. 2017). As an example, the VC maps of COand C O for run M7 are displayed in Fig. 10. The small-scale dense regions with small velocity fluctuations are pref-erentially traced by C O, while the relatively diffuse gason large scales with large velocity fluctuations is primarilytraced by CO, resulting in larger VCs of CO. In addi-tion, we can easily see that the structures in the CO mapare more anisotropic than those in the C O map, which arealigned with the mean magnetic field.For the quantitative anisotropy analysis, we measure theVC SFs of both CO and C O in the directions perpen-dicular and parallel to the mean magnetic field, as shown inFig. 11(a) as an example. Different from the 2D SFs with-out radiative transfer, the VC SFs of CO and C O steepentoward smaller scales, due to the stronger self-absorptionin small-scale dense regions. This is different from thecase with homogeneous density distribution, where the self-absorption effect becomes more important toward larger LOSseparations, and the turbulence scaling can only be seen onsufficiently small scales (Kandel et al. 2017). Obviously,the difference between D VC ( R ⊥ ) and D VC ( R (cid:107) ) for COis significantly larger than that for C O. The AD is mea-sured as the ratio D VC ( R ⊥ ) /D VC ( R (cid:107) ) at around ten grid EASURING MAGNETIZATION WITH ROTATION MEASURES AND VELOCITY CENTROIDS (a) CO (b) C O Figure 10.
VC maps of (a) CO and (b) C O emission for run M7. Note that the color scales in (a) and (b) are different. cells, i.e., . pc, where both the effects from isotropic driv-ing and self-absorption are insignificant. In realistic situa-tions with a more extended inertial range of turbulence, oneshould measure the AD over a range of length scales wherethe anisotropy remains constant. The dependence of AD on M A is presented in Fig. 11(b). The error bars correspond tothe standard errors with error propagation (see Section 3.2.2).When using the low-density tracer CO, we obtain the resultsimilar to that with selective sampling of low-density regionsin Section 3.2.2 (see Fig. 8(a)). D VC recovers the theoreti-cally expected dependence of anisotropy on M A as M − / .However, when the higher-density tracer C O is used, theanisotropy is suppressed and has a much weaker dependenceon M A . This is similar to the result with random sampling inSection 3.2.1 (see Fig. 8(b)). The above tests with syntheticobservations demonstrate that VC s of both CO and C Ocan be used to measure the orientation of the plane-of-skymagnetic field, and the former can also be used to determineits strength.We further display the M A values obtained from the ADsmeasured with CO, i.e., AD − / , in Fig. 12. They agreewell with the real values for small M A . The discrepancy seenat M A ∼ is due to the presence of initial uniform mag-netic field in our simulations as mentioned before. We alsonote that for runs M1-M5 with smaller M S values than thoseof M6-M9, the measured ADs with CO tend to be largerthan the theoretical expectation, leading to the underestimateof M A . The additional uncertainties due to the M S depen-dence of AD should be taken into account when measuring M A with CO. DISCUSSIONNote that even in low-density regions in supersonic MHDturbulence, the density effect cannot be eliminated from theSFs of RMs and VCs. As a result, the SFs have shallowerslopes than the Kolmogorov one and smaller anisotropiesthan that of turbulent velocities. In addition, as argued in Leeet al. (2016), it is not easy to reveal the true spectral slope ofturbulence by using SFs even with a high numerical resolu-tion. However, our results suggest that the M A dependenceof anisotropy is not very sensitive to the exact slope of SF.When dealing with real observations, the orientation of thePOS mean magnetic field is unknown. One should mea-sure the SFs of RMs or VCs at different position angles onthe POS. The direction corresponding to the minimum SFis the direction of the POS mean magnetic field. Then theanisotropy can be measured by using the ratio between themaximum and minimum SFs (measured in directions orthog-onal to each other) averaged over a range of length scaleswith a constant anisotropy. The anisotropy of our interestoriginates from the magnetization of turbulence, but the ob-served anisotropy can be affected by various effects. Be-sides the driving effect and density effect, in dense regionsof molecular clouds, both the scaling and anisotropy of tur-bulence can be affected by gravity. Once the gas traced by CO undergoes gravitational collapse, the dependence ofAD on M A would get changed.In this work we only considered the RMs and VCs of aturbulence layer with a given thickness, which is applicablewhen studying, e.g., the interstellar RMs for extragalactic ra-dio sources, the VCs of a molecular cloud. In the case ofRMs of Galactic pulsars with different distances from the ob-4 X U & H U R [pc] -1 D V C ( R ) [ m s − ] D V C, ⊥ , CO D V C, || , CO D V C, ⊥ ,C O D V C, || ,C O R (a) M A D V C ( R ⊥ ) / D V C ( R k ) COC O (b) Figure 11. (a) SFs of VCs of CO and C O for run M7. (b) D VC ( R ⊥ ) /D VC ( R (cid:107) ) as a function of M A measured with CO and C O. Theopen symbols correspond to runs M1-M5, and the filled symbols correspond to runs M6-M9. The dashed line shows M − / . M A M ea s u r ed M A Figure 12. M A measured with CO vs. the real value. The opencircles correspond to runs M1-M5, and the filled circles correspondto runs M6-M9. The dashed line represents the equality betweenthe measured and real values. server, there are extra terms in their SF accounting for the RMfluctuations induced by distance differences, which are inde-pendent of the separation R between a pair of lines of sight(Lazarian & Pogosyan 2016; Xu & Zhang 2020). The scalingof turbulence can still be revealed as long as the RM fluctua-tions induced by interstellar turbulence, which increase with R , can dominate over those induced by distance differencesover a range of length scales.There are other methods for measuring the anisotropyof turbulent magnetic fields and velocities. For instance,the anisotropy of turbulent magnetic fields can be extractedfrom synchrotron intensity and polarization data (Lazar-ian & Pogosyan 2012, 2016). The sturcture function and quadrupole ratio modulus can be used to quantify theanisotropy of the synchrotron polarization intensity (Leeet al. 2019; Wang et al. 2020). In addition, in spatially-coincident synchrotron emission and Faraday rotation re-gions, the statistical properties of the plane-of-sky and LOScomponents of magnetic fields can be separately obtained atdifferent wavelengths (Lazarian & Pogosyan 2016; Zhanget al. 2018). The application of our method to studyingthe anisotropy of RM fluctuations in interstellar synchrotron-emitting media deserves further investigation. We note thatalthough RMs are related to the LOS magnetic fields, theanisotropy obtained from RM fluctuations reveals the prop-erties of the plane-of-sky mean magnetic field. Kandel et al.(2016) extended the VCA to study the anisotropy of turbulentvelocities from velocity channel maps, which can be appliedto supersonic MHD turbulence. The synergy of the abovedifferent methods, as well as other techniques for measur-ing magnetization, e.g., the VGT (Lazarian et al. 2018), thesynchrotron intensity and polarization gradients (Lazarian &Yuen 2018; Carmo et al. 2020), is necessary for achieving acomprehensive picture of the interstellar magnetic fields andturbulence. SUMMARYWe have studied the anisotropy of SFs of RMs and VCsin the global reference frame of the mean magnetic field byusing a set of 3D simulations of supersonic and sub-Alfv´enicMHD turbulence. Due to the overwhelmingly large contribu-tions from isotropic density fluctuations in supersonic MHDturbulence, the turbulence anisotropy measured with RMsand VCs is significantly suppressed.This result is consistentwith earlier findings in observations and simulations (e.g., Xu& Zhang 2016; Hu et al. 2020a).By selectively sampling the turbulence volume using thelines of sight with relatively low column densities, we findthat the SFs of RMs and VCs normalized by column densities
EASURING MAGNETIZATION WITH ROTATION MEASURES AND VELOCITY CENTROIDS M A as M − / . This is consis-tent with our theoretical expectation for anisotropic turbulentmagnetic fields and velocities in sub-Alfv´enic MHD turbu-lence. The numerical result also shows that the anisotropy of2D SFs has the same dependence on M A as the 3D SF of tur-bulent velocities studied in Hu et al. (2020c). We see that byselectively sampling the relatively diffuse regions in super-sonic turbulence in cold interstellar phases, the anisotropicfluctuations of RMs and VCs can be used to determine theorientation of the POS magnetic field and its strength if thedensity is known.By applying the radiative transfer calculations to our super-sonic MHD simulations, we produce synthetic observationsof CO and C O emission. In the presence of small-scalehigh-density structures in supersonic turbulence, the SF ofnormalized VCs steepens toward smaller scales due to the op-tical depth effect. As expected, the anisotropy obtained with the low-density tracer CO has the M A dependence consis-tent with M − / . It can be used to measure both the directionand strength of the POS magnetic field. The anisotropy mea-sured with the higher-density tracer C O is smaller and hasa weaker dependence on M A , which can still indicate the di-rection of the POS magnetic field.ACKNOWLEDGMENTSS.X. acknowledges the support for this work provided byNASA through the NASA Hubble Fellowship grant Software:
MATLAB (MATLAB 2018), ZEUS-MP/3Dcode (Hayes et al. 2006), Paraview (Ahrens et al. 2005)REFERENCES