Curvature of Magnetic Field Lines in Compressible Magnetized Turbulence: Statistics, Magnetization Predictions, Gradient Curvature, Modes, and Self-gravitating Media

The magnetic field line can be considered to be "the path" of an imaginary particle with the particle speed B at each point in a small neighborhood. Such characterization of "magnetic field lines" globally does not exist since the concept of magnetic field lines becomes ambiguous when we face regions with magnetic reconnections or magnetic field crossing (Newcomb 1958). The consideration of the geometry of the local field lines provides an immediate advantage: a natural curvilinear frame called the Frenet–Serret frame could be defined locally with two geometric properties about magnetic field called curvature κ and torsion τ, which characterize how the magnetic field lines deviate from a line and a plane, respectively. While their velocity counterparts (velocity curvature, torsion) are well studied (see, e.g., Braun et al. 2006; Kadoch et al. 2011), the study of magnetic field curvature has not been popular until recently (Yang et al. 2019).

In the studies of magnetic field topology, the concept of curvature is important since magnetic field lines are expected to have a smaller curvature as the strength of the magnetic field increases. A recent numerical work by Yang et al. (2019) shows the probability density function (pdf) of the magnetic field curvature has a power-law tail of κ⁻² in 2D and κ^−2.5 in 3D for incompressible simulation with initially fluctuation energy equipartitioned between the kinetic and magnetic ones. They also show that magnetic field lines follow a proportionality relation of κ ∝ f_B/B², where f_b represents the normal force component. Observationally Li et al. (2015) use the curvature of magnetic field lines as an estimate of magnetic field strength in NGC 6334. Together, curvature of magnetic field lines becomes an important physical quantity in characterizing the strength of magnetic field.

Mathematically, we would parameterize the magnetic field lines by the line variable s assuming that the magnetic field forms a vector field for an imaginary particle. Intuitively, we would imagine a particle placed at at an arbitrary initial position r₀ and allow it to evolve following the vector integral $\tfrac{d{{\boldsymbol{L}}}_{B}}{{dt}}={\boldsymbol{B}}({\boldsymbol{r}}-{{\boldsymbol{r}}}_{0})$, where the magnetic field B is a velocity field of the particle. The path length of the imaginary particle is then

$$\begin{eqnarray}&&s(l)={\int }_{0}^{l}{dl}^{\prime} \sqrt{{B}_{x}^{2}(l^{\prime} )+{B}_{y}^{2}(l^{\prime} )+{B}_{z}^{2}(l^{\prime} )},\end{eqnarray} \tag{ 1 }$$

and the Frenet–Serret frame of the magnetic fields lines would be (Callen et al. 2003)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d\hat{t}}{{ds}} & = & +\kappa \hat{n}\\ \displaystyle \frac{d\hat{n}}{{ds}} & = & -\kappa \hat{t}+\tau \hat{b}\\ \displaystyle \frac{d\hat{b}}{{ds}} & = & -\tau \hat{n},\end{array}\end{eqnarray} \tag{ 2 }$$

where $\hat{t}=\hat{B}$ is the unit vector of the magnetic field, $\hat{n},\hat{b}$ are the normal and binormal vectors, respectively, $\tfrac{{dx}}{{ds}}=\hat{t}\cdot {\rm{\nabla }}x$ is the line derivative along a vector x, and τ is the torsion of the magnetic field line (Section 8.1). Under this formulation, the signed curvature could be given by

$$\begin{eqnarray}&&\kappa =\hat{n}\cdot \displaystyle \frac{d\hat{t}}{{ds}}.\end{eqnarray} \tag{ 3 }$$

The Frenet–Serret frame is a natural frame in studying the local geometry of magnetic field. We expect that in the presence of magnetized turbulence, there should be a power law κ ∝ B^−γ, where γ is a constant yet to be defined. From Equation (3), apparently the magnetic field line curvature is inversely proportional to the squared amplitude of magnetic field. Indeed, as argued in Yang et al. (2019), if the force term is constant, then κ ∝ B⁻², but in the case of MHD turbulence the interaction of velocity and magnetic fields would tend to reduce the dependencies of κ to B. Therefore, we would expect the properties of turbulence are characterized not only by magnetic field strength but also by the velocities and densities of the fluid elements, a more appropriate estimate would be

$$\begin{eqnarray}&&\kappa \propto {M}_{{\rm{A}}}^{\gamma },\end{eqnarray} \tag{ 4 }$$

where M_A = v_inj/v_A is the ratio of the injection and Alfvén velocities, ${v}_{{\rm{A}}}={B}_{0}/\sqrt{4\pi {\rho }_{0}}$.

We can make a similar estimation on the dependence of κ ∝ ${M}_{A}^{\gamma }$ from the theory of MHD turbulence (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999). As Goldreich & Sridhar (1995) formulated for trans-Alfvénic, i.e., M_A = 1 turbulence, in what follows we are mostly using the Lazarian & Vishniac (1999) expressions obtained for M_A < 1. For M_A > 1 the Goldreich & Sridhar (1995) approach can be easily generalized (see Lazarian 2006) as we discuss below.

In the incompressible limit, the parallel and perpendicular length scales of the turbulent eddies would be related by the following relations based on the theory of MHD turbulence. In the following we shall only consider cases where we have strong magnetic field (M_A = v_inj/v_A < 1, where v_inj is the injection velocity and v_A is the Alfvén velocity) or we are in the regime of dynamically important magnetic field (M_A > 1 but the length scale $l\lt {l}_{{\rm{A}}}={L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{-3}$, where L_inj is the injection scale). Here we consider that the magnetic field eddies that are coherent to that of the velocity eddies (which is the case when we are considering Goldreich & Sridhar turbulence; see Goldreich & Sridhar 1995), so that one could simply use the same scaling law for magnetic field structures without further approximations.

In the cases of strong magnetization, there exists a scale $l={L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{2}$ such that scales smaller than that would be the strong turbulence. We shall omit discussion of length scales larger than that since those usually have a limited spatial range. We stress that the parallel length scale λ_∥ and perpendicular length scale λ_⊥ are defined in terms of the local direction of magnetic field. While in earlier works (Shebalin et al. 1983; Higdon 1984; Matthaeus et al. 1996) the anisotropy is observed and theoretically tested, these anisotropies are generally measured along the mean magnetic field. In fact, the concept of the local magnetic field is absent in Goldreich & Sridhar (1995), where all the closure relations used for the derivation are formulated in the reference frame of the mean field (Lazarian & Vishniac 1999). The concept of local reference frame was shown to be valid numerically (see Cho & Vishniac 2000; Maron & Goldreich 2001; Cho et al. 2001, or the Appendix of Lazarian et al. 2018 for a comprehensive discussion). It is important to stress that the relations between λ_∥ and λ_⊥ are not valid if $\parallel$ and $\perp$ distances are measured with respect to the mean field (Cho & Vishniac 2000).

Notice that the Goldreich-Sridhar turbulence has to be understood under the reference frame defined by the local magnetic field directions (Lazarian & Vishniac 1999; Lazarian et al. 2018). That means that the anisotropy should be computed with respect to magnetic field at the scale of the eddies. In the local system of magnetic field of the eddies, the parallel and perpendicular scales of the eddies are related for M_A < 1 as (Lazarian & Vishniac 1999)

$$\begin{eqnarray}&&{\lambda }_{\parallel }\sim {L}_{\mathrm{inj}}{\left(\displaystyle \frac{{\lambda }_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{2/3}{M}_{{\rm{A}}}^{-4/3},\end{eqnarray} \tag{ 5 }$$

where λ_∥,⊥ are the parallel and perpendicular length scales of the turbulent eddies under Lazarian & Vishniac (1999) formalism. Equation (5) has two differences from the original expression by Goldreich & Sridhar (1995). First of all, it relates the physical scales of the eddies rather than wavenumbers. The latter are given in the frame of the mean field and do not exhibit the scale-dependent anisotropy of Equation (5). Second, Goldreich & Sridhar theory is formulated for M_A = 1.

Readers should be careful that the parallel and perpendicular scales that we are discussing here are referring to the local scales argument in Lazarian & Vishniac (1999). Numerically the 2/3 scaling law is tested in Cho & Vishniac (2000), and the energy spectrum is tested in Jungyeon & Lazarian (2003). Those studies confirmed that the aspect ratio is larger for smaller eddies and, in fact, follows the predicted "critical balance" relation. The justification of such a procedure follows from the eddy description of MHD turbulence based on the Lazarian & Vishniac (1999) study. The subsequent studies (see Cho et al. 2002; Beresnyak et al. 2005) provided the numerical support for these scalings. Recently, Yuen et al. (2018, see their Appendix) also revisited this issue and stated again clearly that only the local computation would yield the desired 2/3 scaling law, in agreement with the expectations based on Lazarian & Vishniac (1999) reconnection theory.

For super-Alfvénic turbulence M_A > 1 and at large scales magnetic fields do not change the Kolmogorov picture. In the case of weak magnetization and $l\lt {l}_{{\rm{A}}}={L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{-3}$, the parallel length scale λ_∥ and perpendicular length scale λ_⊥ are related by

$$\begin{eqnarray}&&{\lambda }_{\parallel }\sim {L}_{\mathrm{inj}}\,{\left(\displaystyle \frac{{\lambda }_{\perp }}{{L}_{\mathrm{inj}}}\right)}^{2/3}.\end{eqnarray} \tag{ 6 }$$

In either case, the minimal curvature of the eddies measured for the eddy of size λ_⊥ could be estimated by ${\lambda }_{\perp }/{\lambda }_{\parallel }^{2}$. Therefore, we would obtain the expected curvature for eddies as a variable of the minor length λ_⊥:

$$\begin{eqnarray}&&\begin{array}{l}\kappa ({\lambda }_{\perp })\sim {L}_{\mathrm{inj}}^{-2/3}{\lambda }_{\perp }^{-1/3}\\ \ \ \cdot \ \left\{\begin{array}{ll}{M}_{{\rm{A}}}^{8/3} & ({M}_{{\rm{A}}}\,\lt \,1,l\lt {L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{2})\\ 1 & ({M}_{{\rm{A}}}\,\gt \,1,l\lt {L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{-3})\end{array}\right.,\end{array}\end{eqnarray} \tag{ 7 }$$

for which the curvature at length scale λ_⊥ is proportional to the Alfvénic Mach number. Here the curvature that we measured, κ(λ_⊥), corresponds to the curvature we measure in real space when we only discuss eddies of size λ_⊥. For eddies larger than the scales in the two individual cases, i.e., for $l\gt {L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{2}$ if M_A < 1 and $l\lt {L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{-3}$ for M_A > 1, the curvature of these isotropic eddies with scale λ is simply λ⁻¹ with no relation between λ and M_A.

The measured curvature value in real space would be mainly contributed by two factors: (1) the largest scale that has turbulent anisotropy and (2) contributions from eddies with scales larger than the largest scale that has turbulent anisotropy. Since the curvature value of real space is basically a statistical sum of all curvature values from the eddies with different sizes, it is natural to consider for what size the eddy dominates over the measurements. In the case of sub-Alfvénic turbulence, the largest scale measurable with anisotropy is $l={L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{2}$. Therefore, the curvature of eddies at scale $l={L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{2}$ would be

$$\begin{eqnarray}&&\kappa ({\lambda }_{\perp }={L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{2})={L}_{\mathrm{inj}}^{-1}{M}_{{\rm{A}}}^{2}\quad ({M}_{{\rm{A}}}\,\lt \,1),\end{eqnarray} \tag{ 8 }$$

while for super-Alfvénic turbulence, the respective transition scale is $l={L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{-3}$,

$$\begin{eqnarray}&&\kappa ({\lambda }_{\perp }={L}_{\mathrm{inj}}{M}_{{\rm{A}}}^{-3})={L}_{\mathrm{inj}}^{-1}{M}_{{\rm{A}}}^{1}\quad ({M}_{{\rm{A}}}\,\gt \,1).\end{eqnarray} \tag{ 9 }$$

The contribution of curvature from eddies with scales larger than the scale having turbulence anisotropy has no dependence on M_A. The observed curvature would then be a statistical sum of curvature values acting on different size of eddies. Moreover, notice that the discussion above is based on the turbulence scaling laws from incompressible turbulence theory. Since compressible modes in compressible magnetized turbulence (see Jungyeon & Lazarian 2003) rarely have nonzero curvature and the observed curvature is the statistical average of the curvature values obtained from the three MHD modes, the larger weight of compressible modes would decrease the observed curvature (see discussion in Section 6). Observationally we cannot perform an analysis of curvature as a function of λ_⊥ unless precise magnetic field orientations are given, which requires the knowledge of the 3D magnetic field distribution. Statistically one would only obtain a value of curvature contributed by all three modes of MHD turbulence and all scales. As a result, we expect that the measured curvature value has a weaker dependence on M_A for both cases, i.e., we expect

$$\begin{eqnarray}\kappa \sim {L}_{\mathrm{inj}}^{-1}\cdot \left\{\begin{array}{ll}{M}_{{\rm{A}}}^{2-\alpha } & ({M}_{{\rm{A}}}\lt 1)\\ {M}_{{\rm{A}}}^{1-\alpha } & ({M}_{{\rm{A}}}\gt 1)\end{array}\right.\end{eqnarray} \tag{ 10 }$$

for some constant α dependent on the properties of injection and the composition of modes. As we will see in Section 4, the inclusion of α is necessary in explaining the behavior of curvature in compressible turbulence.

We also expect the index $\gamma =2-\alpha ({M}_{{\rm{A}}}\lt 1)$  or $1\,-\alpha ({M}_{{\rm{A}}}\gt 1)$ to be shallower in 2D, i.e γ_2D < γ_3D since the projection effect of magnetic field would tend to cancel out the magnetic field deviations, which are not a straight line. As a result, the curvature observed in the projected space should be systematically smaller than that in 3D.

The numerical data cubes are obtained by 3D MHD simulations that are from a single-fluid, operator-split, staggered grid MHD Eulerian code ZEUS-MP/3D to set up a 3D, uniform turbulent medium. Our simulations are isothermal with T = 10 K. To simulate the part of the interstellar cloud, periodic boundary conditions are applied. We inject turbulence solenoidally.⁶

For our controlling simulation parameters, various Alfvénic Mach numbers M_A = V_inj/V_A and sonic Mach numbers M_s = V_inj/V_s are employed,⁷ where V_inj is the injection velocity, while V_A and V_s are the Alfvén and sonic velocities, respectively, which are listed in Table 1. For the case of M_A < M_s, it corresponds to the simulations of turbulent plasma with thermal pressure smaller than the magnetic pressure, i.e., plasma with low confinement coefficient $\beta /2={V}_{s}^{2}/{V}_{{\rm{A}}}^{2}\lt 1$. In contrast, the case M_A > M_s corresponds to the magnetic-pressure-dominated plasma with high confinement coefficient β/2 > 1.

Table 1.  Description of MHD Simulation Cubes, Some of Which Have Been Used in the Series of Papers about VGT (Yuen & Lazarian 2017a, 2017b; Lazarian & Yuen 2018a, 2018b)

Model	MS	MA	$\beta =2{M}_{{\rm{A}}}^{2}/{M}_{S}^{2}$	Resolution
huge-0	6.17	0.22	0.0025	7923
huge-1	5.65	0.42	0.011	7923
huge-2	5.81	0.61	0.022	7923
huge-3	5.66	0.82	0.042	7923
huge-4	5.62	1.01	0.065	7923
huge-5	5.63	1.19	0.089	7923
huge-6	5.70	1.38	0.12	7923
huge-7	5.56	1.55	0.16	7923
huge-8	5.50	1.67	0.18	7923
huge-9	5.39	1.71	0.20	7923

Note. M_s and M_A are the rms values at each snapshot.

Download table as:  ASCIITypeset image

From now on we refer to the simulations in Table 1 by their model name. For example, the figures with model name indicate which data cube was used to plot the corresponding figure. Each simulation name follows the rule that is the name is with respect to the varied M_s and M_A in ascending order of confinement coefficient β. The selected ranges of M_s, M_A, and β are determined by possible scenarios of astrophysical turbulence from very subsonic to supersonic cases.

The raw data from simulation cubes are converted to synthetic maps for studies of curvature. Since we would investigate both polarization and magnetic field probed by the velocity gradient technique (VGT; Yuen & Lazarian 2017a; Lazarian & Yuen 2018a, see Section 3.3), we shall deliver the method of synthesizing the Stokes parameters, intensities, centroids, and velocity channels here.

The Stokes parameters for synthetic observation are given by

$$\begin{eqnarray}\begin{array}{rcl}Q & \propto & \displaystyle \int {dzn}\cos (2\theta ){\sin }^{2}{\gamma }_{\mathrm{inc}}\\ U & \propto & \displaystyle \int {dzn}\sin (2\theta ){\sin }^{2}{\gamma }_{\mathrm{inc}}\\ {\theta }_{\mathrm{pol}} & = & \displaystyle \frac{1}{2}\tan {2}^{-1}(U/Q),\end{array}\end{eqnarray} \tag{ 11 }$$

where n is the number density and θ and γ_inc are the 3D planer angle and inclination angle of magnetic field vectors with respect to the line of sight (LOS), respectively. The dispersion of the polarization angle θ_pol is directly proportional to the perpendicular Alfvénic Mach number M_A,⊥.

The normalized velocity centroid C(R) in the simplest case⁸ is defined as

$$\begin{eqnarray}\begin{array}{rcl}C({\boldsymbol{R}}) & = & {I}^{-1}\displaystyle \int {\rho }_{v}({\boldsymbol{R}},v){vdv},\\ I({\boldsymbol{R}}) & = & \displaystyle \int {\rho }_{v}({\boldsymbol{R}},v){dv},\end{array}\end{eqnarray} \tag{ 12 }$$

where ρ_v is density of the emitters in the position–position–velocity (PPV) space, v is the velocity component along the LOS, and R is the 2D vector in the pictorial plane. The integration is assumed to be over the entire range of v. Naturally, I(R) is the emission intensity. C(R) is also an integral of the product of velocity and LOS density, which follows from a simple transformation of variables (see Lazarian & Esquivel 2003). For constant density, C(R) is just the LOS velocity averaged over the LOS.

We also consider the velocity channel at v = 0 here as a case study. Mathematically, the density in PPV space of emitters with local sonic speed ${c}_{s}({\boldsymbol{x}})=\sqrt{\gamma {k}_{B}T/{\mu }_{\mathrm{MMW}}}$, where μ_MMW is the mean molecular weight of the emitter, moving along the LOS with stochastic turbulent velocity u(x) and regular coherent velocity, e.g., the galactic shear velocity, v_g(x), is (Lazarian & Pogosyan 2004)

$$\begin{eqnarray}&&{\rho }_{s}({\boldsymbol{X}},v)={\int }_{0}^{S}{dz}\displaystyle \frac{\rho ({\boldsymbol{x}})}{\sqrt{2\pi {\beta }_{T}}}\exp \left[-\displaystyle \frac{{\left(v-{v}_{{\rm{g}}}({\boldsymbol{x}})-u({\boldsymbol{x}})\right)}^{2}}{2{c}_{s}^{2}({\boldsymbol{X}},z)}\right],\end{eqnarray} \tag{ 13 }$$

where sky position is described by 2D vector X = (x, y), z is the LOS coordinate, γ is the adiabatic index, and S is the cloud depth. Notice that c_s would be a function of distance if the emitter is not isothermal. Equation (13) is exact, including the case when the temperature of emitters varies in space. The observed velocity channel at velocity position v₀ and channel width Δv is then, assuming a constant velocity window W(v) = 1 and v_g(x) = 0,

$$\begin{eqnarray}&&\begin{array}{l}{Ch}({\boldsymbol{X}};{v}_{0},{\rm{\Delta }}v)={\displaystyle \int }_{{v}_{0}-{\rm{\Delta }}v/2}^{{v}_{0}+{\rm{\Delta }}v/2}{dv}{\rho }_{s}({\boldsymbol{X}},v)\\ \ =\ {\displaystyle \int }_{0}^{S}{dz}\displaystyle \frac{\rho ({\boldsymbol{x}})}{\sqrt{2\pi {c}_{s}^{2}}}{\displaystyle \int }_{{v}_{0}-{\rm{\Delta }}v/2}^{{v}_{0}+{\rm{\Delta }}v/2}{{dve}}^{-\tfrac{{\left(v-u({\boldsymbol{x}})\right)}^{2}}{2{c}_{s}^{2}}}.\end{array}\end{eqnarray} \tag{ 14 }$$

We shall deliver the method of computing the gradients of velocity channel and use it as a probe of tracing magnetic fields.

The VGT (see González-Casanova & Lazarian 2017) is an innovative method that uses the properties of turbulence anisotropy in MHD turbulence to probe the direction of magnetic field. The basic idea is to obtain the magnetic field predictions by rotating the output of the sub-block averaging of the gradients of observables by 90° (Yuen & Lazarian 2017a), which is supported by a number of theoretical and numerical works (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999; Maron & Goldreich 2001; Cho & Lazarian 2002; Jungyeon & Lazarian 2003) The VGT is applicable to velocity centroids (González-Casanova & Lazarian 2017; Yuen & Lazarian 2017a), intensities (Yuen & Lazarian 2017a; Hu et al. 2019), and also velocity channel maps (Lazarian & Yuen 2018a). The same method is also migrated to synchrotron studies and applicable to both synchrotron intensities (Lazarian et al. 2017) and multifrequency synchrotron polarization (Lazarian & Yuen 2018b).

The gradients of the three observables (intensity, centroid, velocity channel) we introduced in Section 3.2 would be computed as follows. We would first compute the Sobel kernel of the observables, which we would call the raw gradients. The distribution histogram peak of raw gradients would provide us with the predicted direction of magnetic field directions probed by the gradients of observables, provided that the Gaussian fitting requirement stated in Yuen & Lazarian (2017a) is satisfied, which is called sub-block averaging in Yuen & Lazarian (2017a). The 90° rotated gradients are our predicted magnetic field directions by the gradients of observables. We would not apply further improvements of the technique (e.g., Hu et al. 2018; Lazarian & Yuen 2018a) since these techniques tend to straighten the estimated magnetic field lines.

Computing curvature by Equation (2) is difficult since the computation of B · ∇B does not often yield a vector parallel to $\hat{n}$ numerically. The reason is that the spatial derivative of magnetic field is not guaranteed to have (B · ∇B) · B = 0. In view of that, Yang et al. (2019) use the expression $\kappa \,=| \hat{{\boldsymbol{B}}}\times (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}\hat{{\boldsymbol{B}}})|$ to extract the curvature-only part of magnetic field. Below we discuss an algorithm that could possibly bypass the problem with the cost of interpolation accuracy: we treat the magnetic field as the velocity field of an imaginary particle (see Section 2) and find the path integral in a sufficiently small area, so that one could have a representation of the "magnetic field line function L_B(t)" that has

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{L}}}_{B}(t=0) & = & {{\boldsymbol{r}}}_{0}=\mathrm{Spatial}\ \mathrm{position}\ \mathrm{of}\ \mathrm{the}\ \mathrm{pixel}\\ \displaystyle \frac{d{{\boldsymbol{L}}}_{B}}{{dt}} & = & {\boldsymbol{B}}({\boldsymbol{r}}-{{\boldsymbol{r}}}_{0}).\end{array}\end{eqnarray} \tag{ 15 }$$

Readers should be reminded that if we put $\tfrac{d{{\boldsymbol{L}}}_{B}}{{dt}}={\boldsymbol{v}}$, then we are effectively converting the Eulerian hydrodynamic variables to the Lagrangian one, which is simply the standard Lagrangian particle-tracking algorithm (see Ouellette et al. 2006; Xu et al. 2007).

We perform an integrator method that integrates Equation (15) by the Runge-Kutta method (order 2/4/4.5 depending on the accuracy requirement). See the Appendix for the convergence test. We first impose an interpolation field so that the magnetic vector field B(r) is well defined in a local patch of the 3D real space ${\boldsymbol{r}}\in {\rm{\Delta }}{L}^{3}\ \subset {{\mathbb{R}}}^{3}$. The interpolation is usually performed using the family of splines to estimate values between the grid points. Cubic spline is a popular option. The interpolation could be done very easily through Julia's Interpolation package.⁹ Then, one could obtain the tangent vector:

$$\begin{eqnarray}&&{\boldsymbol{T}}(t)=\displaystyle \frac{{\dot{{\boldsymbol{L}}}}_{B}(t)}{| {\dot{{\boldsymbol{L}}}}_{B}(t)| }.\end{eqnarray} \tag{ 16 }$$

Following Equation (3), one could obtain both κ and τ very easily.

The number of steps required for the integrator method is explicitly dependent on how one expresses the differential operator $\tfrac{d}{{dt}}$. For instance, if one uses the 1D five-point stencil in estimating the differentiation on the position vector r(t) at time t₀ with some custom step size dt,

$$\begin{eqnarray}&&\begin{array}{l}{\boldsymbol{r}}^{\prime} ({t}_{0})\\ \ =\ \displaystyle \frac{-{\boldsymbol{r}}({t}_{0}+2{dt})+8{\boldsymbol{r}}({t}_{0}+{dt})-8{\boldsymbol{r}}({t}_{0}-{dt})+{\boldsymbol{r}}({t}_{0}-2{dt})}{12{dt}},\end{array}\end{eqnarray} \tag{ 17 }$$

then the number of points required in obtaining κ(t = t₀) and τ(t = t₀) (see Section 8.1) at position t₀ is 2 × (5 − 1) + 1 = 9 and 3 × (5 − 1) + 1 = 13 points, respectively. The advantage of the current method is that we are sticking to the definition of the Frenet–Serret frame based on the tangent of the (magnetic field) line. The small dt will guarantee that the integrated line would only be locally defined. The method is valid for both 2D and 3D. However, readers should be reminded that torsion is nonzero only when we are tackling a 3D line.

This method is very versatile since one could also compute curvature of scalar structures, say, intensity map I, by replacing B in Equation (15) with the gradients of I rotated by 90°. This application is especially useful for the studies based on the VGT since curvature is an important quantity aside from orientation and amplitudes of gradients that characterize the underlying magnetic field properties.

Readers should be reminded that, to compute curvature, one must make sure the respective circle that has radius 1/κ could be represented by the grid resolution. One caveat here is that it is impossible to obtain curvature values larger than 0.5 (in units of 1 pixel⁻¹) since the respective circle with radius smaller than 2 pixels would not be resolved in the numerical grid.

In 3D we do not need to employ the algorithm as listed in Section 3.4 since it is straightforward to show $\kappa \,=| \hat{{\boldsymbol{B}}}\,\times (\hat{{\boldsymbol{B}}}\cdot {\rm{\nabla }}\hat{{\boldsymbol{B}}})|$. Figure 1 shows how the 3D statistics would behave as a function of the Alfvénic Mach number. Here we use the mean and standard deviation as ways to extract simple statistics (top panel and middle panel of Figure 1). From our theoretical discussion (Section 2), we expect that statistically the magnetic field curvature should be proportional to ${M}_{{\rm{A}}}^{\gamma }$ for some constant γ. From the top panel of Figure 1, we see a two-section power law for $\langle \kappa \rangle$ (in 1 pixel⁻¹):

$$\begin{eqnarray}\langle \kappa \rangle \propto \left\{\begin{array}{ll}{M}_{{\rm{A}}}^{1.67} & {M}_{{\rm{A}}}\lt 1\\ {M}_{{\rm{A}}}^{0.71} & {M}_{{\rm{A}}}\geqslant 1\end{array}\right.,\end{eqnarray} \tag{ 18 }$$

with a turning point $\langle \kappa \rangle \sim 0.15\,{\mathrm{pixel}}^{-1}$, while for standard deviation of curvature, σ_κ (middle panel of Figure 1, in 1 pixel⁻¹),

$$\begin{eqnarray}{\sigma }_{\kappa }\propto \left\{\begin{array}{ll}{M}_{{\rm{A}}}^{1.67} & {M}_{{\rm{A}}}\lt 1\\ {M}_{{\rm{A}}}^{0.78} & {M}_{{\rm{A}}}\geqslant 1\end{array}\right.,\end{eqnarray} \tag{ 19 }$$

with a turning point $\langle \kappa \rangle \sim 0.25\,{\mathrm{pixel}}^{-1}$. The fitting lines in Equations (18) and (19) look surprisingly similar, with a relatively high coefficient of determination of 0.9 or above. These measurements are consistent with the theoretical prediction in Section 2.2, with the coefficient α defined in Section 2.2 to be

$$\begin{eqnarray}\alpha =\left\{\begin{array}{ll}0.33 & ({M}_{{\rm{A}}}\,\lt \,1,\langle \kappa \rangle )\\ 0.29 & ({M}_{{\rm{A}}}\,\gt \,1,\langle \kappa \rangle )\\ 0.33 & ({M}_{{\rm{A}}}\,\lt \,1,{\sigma }_{\kappa })\\ 0.22 & ({M}_{{\rm{A}}}\,\gt \,1,{\sigma }_{\kappa })\end{array}\right.,\end{eqnarray} \tag{ 20 }$$

which are fairly similar for these four cases. The split of the power laws between M_A ≥ 1 and M_A < 1 is also expected and consistent with our findings in one of our previous works using velocity gradients in probing magnetization (Lazarian et al. 2018) since there are different scaling laws below and above M_A.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also employ the method of histogram tails since it was suggested that f(κ) ∼ κ^−2.5 in Yang et al. (2019). The bottom panel of Figure 1 shows the power-law tail slope as a function of M_A. One could see that the slope power-law tail is actually a function of M_A. Our result (see Figure 2) should not be directly compared to that in Yang et al. (2019), who have a general slope of −2.5 since we are performing simulations for different settings: (1) We are performing full 3D, compressible, driven simulations, while Yang et al. perform 3D, incompressible, decaying simulations. The compressible simulation allows one to store energy in the two other compressible modes, effectively reducing the curvature of magnetic field lines driven by turbulence driving. (2) We are testing our theory using the turbulence statistics, while the work of Yang et al. is based on the formulation of the magnetic force term in the Cauchy momentum equation (3). We tested the dependencies of the histogram tail as a function of global Alfvénic Mach number with multiple simulations, while Yang et al. (2019) characterize the curvature as a distribution of local magnetic field strength. Moreover, it is expected that the pdf of the curvature should be a function of Alfvénic Mach number since in a strongly magnetized medium we do not expect a large variation of curvature values owing to the constraints of a strong field, resulting in a relatively steep power-law tail in the low M_A limit. In the case of weak magnetic field, the allowed values of curvature are less bounded by the magnetic field itself compared to the strong-field cases, which would lead to an extended pdf tail and a shallower slope.

We would like to see whether the dependencies would be seen in observation when we view the cloud parallel or perpendicular to the LOS. Figure 3 shows how $\langle \kappa \rangle$, σ_κ, and the slope of the pdf tail would respond as a function of M_A when we project the Stokes parameter perpendicular and parallel to the mean magnetic field using the procedure laid out in Section 3.

Zoom In Zoom Out Reset image size

We immediately recognize that there are different behaviors for the relation of $\langle \kappa \rangle$ to M_A in 2D and 3D. In 3D we have a two-section power law that has a cutoff of M_A = 1. However, such a power law is not seen in both cases of B⊥ LOS and B∥ LOS. In fact, in the case of B⊥ LOS we see that a uniform power law describes the data better, which has $\langle \kappa \rangle \propto {M}_{{\rm{A}}}^{1.17}$. Nevertheless, it is shallower than the section of the power law that we have in M_A < 1 in Section 4.1 yet steeper than the other section of the power law that we have in M_A > 1. A very similar effect also happens to σ_κ, that only a single power law would be sufficient to describe the data here (${\sigma }_{\kappa }\propto {M}_{{\rm{A}}}^{1.39}$).

We also see that the previously seen power law disappears when we project the magnetic field data along its mean field direction, i.e., B∥ LOS. In this scenario both $\langle \kappa \rangle$ and σ_κ are more or less a constant of M_A. This is expected since we should only see the hydrodynamic nature of the magnetized turbulence if we are observing it along its mean field direction. We see that the value of $\langle \kappa \rangle$ arrives at the maximum value allowed in the curvature algorithm (see Section 3.4). This suggests that the measured Alfvénic Mach number is actually the perpendicular Mach number M_A,perp = M_A cos θ, where θ is the angle between the mean magnetic field and the LOS.

In the case of the power-law tail slope of the pdf, we see that there still exists a power law between the pdf slope and M_A. One interesting thing to note here is that the values of the power-law tail slope actually become more disperse when the synthetic maps are obtained by being projected perpendicular to the mean magnetic field directions than those in 3D. When we view the numerical cube with its mean field parallel to the LOS, we also see that the pdf tail slope has different responses as a function of M_A compared to that when B⊥ LOS. By using these tools together, it is possible to extract the Alfvénic Mach number in the range of M_A ∈ (0.2, 1.7).

Aside from curvature, torsion is also a very important geometric quantity for the magnetic field. While the use of torsion is less popular in literature, we should not underestimate the importance of torsion since it records how sharply the magnetic field lines twist out of the plane of curvature. The concept of torsion is especially useful when we are dealing with the natural oscillation modes in magnetized turbulence, namely, the fast, slow, and Alfvénic modes (see Biskamp 2003; Jungyeon & Lazarian 2003). In fact, when a magnetic field mode is propagating along some directions, the speed of rotation of the eddy that is along the azimuthal direction of the propagating directions records the torsion of the mode. The signature of magnetic field rotates when propagating is especially important when we are studying modes in MHD turbulence.

Curvature and torsion could be easily visualized when the magnetic field lines have some special geometry. For instance, if the magnetic field line could be parameterized by a circular helix ${L}_{B}(t)=\left(a\cos (t),a\sin (t),{bt}\right)$, then the curvature and torsion of this magnetic field line are simply $\kappa =a/\sqrt{{a}^{2}+{b}^{2}}$ and $\tau =b/\sqrt{{a}^{2}+{b}^{2}}$, respectively. From Equation (2) we could have a formula for the signed torsion:

$$\begin{eqnarray}&&\tau =\hat{b}\cdot \displaystyle \frac{d\hat{n}}{{ds}}.\end{eqnarray} \tag{ 23 }$$

The estimation of curvature and torsion would also be very important in advancing the VGT since essentially Equation (2) describes the behavior of first three derivatives of the orientation of magnetic fields. In the VGT we approximated $\hat{{\boldsymbol{t}}}$ by block averaging (Yuen & Lazarian 2017a). The curvature could be possibly approximated by the magnetization technique by Lazarian et al. (2018) with Figure 1 in 3D or Figure 3 in 2D. Numerically if we know both the tangent and curvature, we could then solve the magnetic field line equation by Equation (2) and Equation (15). By studying the properties of the torsion field, we would know what the second derivative of the orientation of magnetic field would look like, and that would provide us with a geometrically more plausible way of reconstructing the magnetic field directions.

One of the very important applications of the Frenet–Serret frame derived from the magnetic field line is to recognize how the three MHD modes are related to the curvature and torsion of magnetic lines. Below we shall discuss how the three MHD modes would be related to the local Frenet–Serret frame of magnetic field by the theory of MHD turbulence (Goldreich & Sridhar 1995; Lazarian & Vishniac 1999). In the case of Jungyeon & Lazarian (2003) decomposition the mean field has to be locally a straight line.¹⁰ In this scenario there is an ambiguity in defining the directions of both $\hat{n}$ and $\hat{b}$. One could simply assign the Alfvén wave unit vector to be along one of these directions. However, in the case of a mean magnetic field with nonzero curvature, the effect of frame changes in the local regions has to be taken into account, and we have to drop some assumptions that were valid in usual compressible 3D magnetized turbulence. For instance, the assumption that Alfvén mode is divergence free is not correct when the local mean field has nonzero curvature. Southwood & Saunders (1985) propose that the Alfvén mode has a nonzero divergence related to the curvature of the mean magnetic field:

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\hat{\zeta }}_{{\rm{A}}}\sim \kappa {\hat{\zeta }}_{{\rm{A}}}\cdot \hat{n},\end{eqnarray} \tag{ 24 }$$

where ${\hat{\zeta }}_{{\rm{A}}}$ is the Alfvén mode unit displacement vector. Following Southwood & Saunders (1985), there is a source term for the slow wave equation to be proportional to $\kappa {\hat{\zeta }}_{{\rm{A}}}\cdot \hat{n}$, while the Alfvén mode wave equation has a similar term that is $\propto \kappa \hat{n}$. Physically the bent magnetic field would induce a change of fluid element volume toward the center of curvature, whose size is proportional to κ. As a result, the plasma pressure due to the compression of volume is changed under this magnetic field geometry, which would drive the slow mode to oscillate along with the azimuthal direction of Alfvén mode. However, the Alfvén mode in this situation is no longer solenoidal (see Equation (23) or Figure 3 of Southwood & Saunders 1985), which is different from the case when we perform mode decomposition with a mean field with κ = 0 (see Jungyeon & Lazarian 2003 for a discussion of MHD modes). Unlike Jungyeon & Lazarian (2003), where the mean magnetic field was approximated by a straight line, the study in Southwood & Saunders (1985) points out that in the presence of uniform nonzero curvature magnetic field the transfer of energy between Alfvén and slow modes increases. This effect is, however, reduced by the Alfvénic cascade happening in one eddy turnover time, which makes the effect important only when the curvature is comparable to the wavenumber of the perturbations under consideration.

Equation (24) suggests that the Alfvén mode wavevector, which lies on the plane spanned by the normal and the binormal vector, would rotate as a function of curvature with its rotation axis aligned with the tangent vector t = B₀. Equation (24) suggests that if λ_⊥ is the perpendicular wavelength and ${\phi }_{k}={\cos }^{-1}(\hat{k}\cdot \hat{b})$, then the angle of rotation ϕ_A of the Alfvén mode with respect to the tangent vector is given by

$$\begin{eqnarray}&&\cos {\phi }_{{\rm{A}}}\sim \displaystyle \frac{\cos {\phi }_{k}}{{\lambda }_{\perp }\kappa }.\end{eqnarray} \tag{ 25 }$$

Therefore, for a magnetized turbulent system with a nonzero mean magnetic field curvature, one could (1) compute the components of wavevector with respect to the Frenet–Serret frame, (2) compute the orientation of MHD modes according to Jungyeon & Lazarian (2003) and represent these eigenvectors by the Frenet–Serret frame, and (3) rotate the three MHD mode vectors with the rotation matrix $R(\hat{t},{\phi }_{{\rm{A}}})$, where $\hat{t}$ is the rotation axis. This method allows one to use the formulation of Jungyeon & Lazarian (2003) to compute the MHD modes.

The important insight of increasing coherent coupling between Alfvén and slow modes in the presence of nonzero curvature mean field suggests that using the concept of modes in systems that have nonzero mean magnetic field curvature should be done with caution. In Section 2, we see that the curvature naturally introduces a length scale related to the radius of curvature r_c = κ⁻¹. Above the aforementioned length scale Alfvén and slow modes are coupled, while below that they act independently. In the case that there are no external forces holding the magnetic field, the magnetic field curvature would restore to infinity, so that the magnetic field will become a straight line when the magnetic energy from the curvature is all transferred to the dissipation of Alfvén and slow modes.

However, in real astrophysical scenarios, there are external forces that can keep magnetic fields bent. One of the best examples would be self-gravity. It is shown observationally by Li et al. (2015) that the curvature of magnetic field could be used in estimating the magnetic field strength on a self-similar molecular cloud, suggesting that the curvature of magnetic field does not dissipate to Alfvén and slow mode driving in a self-gravitating system. In fact, it is easy to imagine that the gravitational field provides support for the curvature of magnetic field. Following Li et al. (2015), for a spherical self-gravitating object of constant density ρ, radius R, and mean magnetic field strength B, along the normal direction of magnetic field,

$$\begin{eqnarray}&&\displaystyle \frac{{B}^{2}}{R}\sim \displaystyle \frac{4\pi }{3}G{\rho }^{2}R,\end{eqnarray} \tag{ 26 }$$

for which one would have κ (v_At_ff) ∼ 1, where t_ff ∼ (Gρ)^−1/2 is the freefall time. When r_c = κ⁻¹ ≫ v_At_ff, magnetic field dominates over gravity, and vice versa. This suggests that curvature could be supported by gravity if Equation (26) is satisfied.

Under the scenario that κ(v_At_ff) ∼ 1, the gravitational energy would become the primary energy source in driving both (non-incompressible) Alfvén and slow waves. When the density of the self-gravitating object increases, it is expected to have a runoff effect by having a larger curvature and a much stronger coupling between the Alfvén and slow modes. Moreover, the dependence of curvature on magnetic field changes from κ ∝ B⁻² in nongravitating systems to κ ∝ ρB⁻¹ in self-gravitating systems.