ALFVÉNIC TURBULENCE BEYOND THE AMBIPOLAR DIFFUSION SCALE

Turbulence and magnetic fields are critical components of the interstellar medium (ISM) of galaxies from kiloparsec to sub-astronomical unit scales. Magnetohydrodynamic (MHD) turbulence is a key element in the study of star formation and molecular cloud structure, magnetic reconnection, heat transport, and cosmic ray propagation (see Elmegreen & Scalo 2004; Tilley & Balsara 2006; Tilley et al. 2006; Ballesteros-Paredes et al. 2007; McKee & Ostriker 2007; Balsara et al. 2008; Brandenburg & Lazarian 2013). Despite the importance of turbulence for ISM studies, many mysteries remain, including the nature of turbulence driving and damping scales.

The smallest scales of the turbulence cascade, including the damping scale, may play a pivotal role in the dynamics of giant molecular clouds (GMCs) and star formation. GMCs are partially ionized, with neutrals coupled to the magnetic field through collisions with ions. The drift of the neutrals toward the central gravitational potential through the ionized particles tied to the magnetic field, known as ambipolar diffusion, is often invoked as a source of the dissipation of the MHD cascade (Zweibel & Josafatsson 1983; Ciolek & Basu 2000; Tassis & Mouschovias 2004). The ambipolar diffusion scale (L_AD), or the scale at which neutrals and ions decouple, has been thought to set the dissipation scale of turbulence in molecular clouds and to set a fundamental characteristic scale for gravitational collapse and star formation (Balsara 1996; Houde et al. 2000a, 2000b; Klessen et al. 2000; Li et al. 2006, 2008a; Li & Houde 2008; Hezareh et al. 2010; Tilley & Balsara 2010; Meyer et al. 2014).

The ambipolar diffusion scale can be estimated as the scale at which the Reynolds number, with diffusivity given by the ambipolar diffusivity, is equal to unity (Brandenburg & Zweibel 1994, 1995; Balsara 1996; Oishi & Mac Low 2006). The ambipolar diffusivity is given by

$$\begin{eqnarray}&&{{\nu }_{{\rm AD}}}=\frac{{{B}^{2}}}{4\pi {{\rho }_{i}}{{\rho }_{n}}\alpha },\end{eqnarray} \tag{ 1 }$$

where ${{\rho }_{i}}$ and ${{\rho }_{n}}$ are the density of the ions and neutrals, B is the magnetic field strength, and α is the frictional coupling coefficient⁵ between the ions and neutrals.

The Reynolds number for ion-neutral drift is defined as:

$$\begin{eqnarray}&&{{R}_{{\rm AD}}}=\frac{LV}{{{\nu }_{{\rm AD}}}},\end{eqnarray} \tag{ 2 }$$

where V is a characteristic velocity (e.g., for trans-Alfvénic turbulence it is the Alfvén speed, ${{V}_{{\rm A}}}=\frac{B}{\sqrt{4\pi {{\rho }_{n}}}}$) and L = L_AD when R_AD = 1. This gives the form of the ambipolar diffusion scale as often found in the literature:

$$\begin{eqnarray}&&{{L}_{{\rm AD}}}={{V}_{{\rm A}}}/\alpha {{\rho }_{i}}.\end{eqnarray} \tag{ 3 }$$

More details about the derivation and significance of L_AD can be found in a number of works including Langer (1978), Zweibel & Josafatsson (1983), Balsara (1996), Klessen et al. (2000), Oishi & Mac Low (2006), Li & Houde (2008), Hezareh et al. (2010), and Meyer et al. (2014).

The application of ambipolar diffusion extends beyond direct studies of star formation to include general studies of magnetic fields. For example, Houde et al. (2000a, 2000b), Li & Houde (2008), and Hezareh et al. (2010) have proposed that the magnetic field in the plane of the sky may be obtained from observations via calculation of the ambipolar diffusion length scale. In light of these diverse interpretations of the meaning of L_AD it is important to understand to what extent L_AD is relevant to turbulence damping in partially ionized gasses.

More generally, ambipolar diffusion has also been proposed to damp particular families of linear MHD waves (see Balsara 1996 and reference therein). On scales larger than L_AD it was predicted that two-fluid turbulence acts like single-fluid MHD turbulence. In particular, on the basis of 1D dispersion analysis, Balsara (1996) showed that two separate mode damping situations can occur at scales at or smaller than L_AD, which are based on the value of plasma beta (i.e., the ratio of the gas to magnetic pressure). When the Alfvén speed is greater than the sound speed the fast and Alfvén wave families are damped at or below L_AD. When the Alfvén speed is smaller than the sound speed, the slow and Alfvén wave families are damped. For either high or low plasma beta, Balsara (1996) predicts that the Alfvénic waves should damp at L_AD and thus the MHD cascade should damp past the ambipolar diffusion scale. To put Balsara (1996) in context, large scale three-dimensional (3D) MHD simulations were not possible at the time. Tilley & Balsara (2011) have extended the study of Balsara (1996) to include flows with radiative effects. They find that the general analysis of the MHD wave modes that are damped remained unchanged.

Does MHD turbulence, specifically the Alfvén wave modes, damp at the decoupling scale L_AD? MHD turbulence is known to be different from a collection of linear Alfvénic waves. Cascading rates and the anisotropy of turbulence should be accounted for carefully before we can make a definitive conclusion about turbulent damping in the partially ionized media. The purpose of this paper is to numerically address the question above, as well as to numerically test the validity of the Alfvén mode scaling relations of the Goldreich & Sridhar (1995, henceforth GS95) theory for ion-neutral turbulence. The GS95 theory was first extended to a partially ionized compressible medium in subsequent works by Lithwick & Goldreich (2001, henceforth LG01) and Lazarian et al. (2004, henceforth LVC04). In this paper we will apply the Cho & Lazarian (2002, 2003; henceforth CL02, CL03, respectively) MHD mode decomposition technique⁶ to two-fluid MHD simulations, first presented in Meyer et al. (2014), in order to investigate the behavior of the Alfvén modes in the ions and neutrals. Our paper is organized as follows: in Section 2 we review the basic scaling relations that are predicted by the GS95 model, in Section 3 we describe the numerical simulations and relevant scales, in Section 4 we present the structure function scalings for the full data cubes and the Alfvénic modes of the ions and neutrals for our simulations, and finally, in Section 5 we discuss our results, followed by our conclusions in Section 6.

MHD turbulence is a subject with an extended history (see the book by Biskamp 2003). However, it has recently been given a boost with the advent of 3D MHD simulations that allow for testing theoretical predictions (see Beresnyak et al. 2005; Brandenburg & Lazarian 2013). The modern theory of MHD turbulence is based on the Goldreich-Sridhar (1995) idea, which was extended and tested in subsequent publications (Lazarian & Vishniac 1999; Maron & Goldreich 2001; CL02; CL03). The applicability of GS95 theory to the partially ionized gas was discussed in LG01 and LVC04.

It has been numerically shown that in the presence of dynamically important magnetic fields, eddies become elongated along the magnetic field lines. GS95's approach to Alfvénic modes can be easily understood: for the eddies that are perpendicular to the magnetic field, the original Kolmogorov energy scaling applies, i.e., ${{V}_{l}}\sim _{\bot }^{1/3}$, where ${{l}_{\bot }}$ denotes scales measured perpendicular to the local magnetic field. Mixing motions induce Alfvénic perturbations that determine the parallel size of the magnetized eddy. This is the concept of critical balance, i.e., the equality of the eddy turnover time (${{l}_{\bot }}/{{V}_{l}}$) and the period of the corresponding Alfvén wave $\sim {{l}_{\parallel }}/{{V}_{{\rm A}}}$, where ${{l}_{\parallel }}$ is the parallel eddy scale and V_A is the Alfvén velocity. Making use of ${{V}_{l}}\sim _{\bot }^{1/3}$, one finds the scaling relation for the parallel and perpendicular eddies as ${{l}_{\parallel }}\sim l_{\bot }^{2/3}$. This represents the scale-dependent tendency of eddies to become more elongated along the magnetic field lines as the energy cascades proceed to smaller scales, and corresponds to the scaling of the slow and Alfvén wave anisotropy. The power spectrum for the slow and Alfvén waves scales as $E\sim k_{\bot }^{-5/3}$ (see review by Brandenburg & Lazarian 2013). Numerical studies of the scaling of compressible MHD turbulence based on the decomposition into Alfvén, Slow , and Fast modes were first performed in CL02 and CL03 and these results were later reconfirmed in Kowal & Lazarian (2010) with a wavelet approach.

GS95 theory assumes the isotropic injection of the energy at scale L and the injection velocity equal to the Alfvén velocity in the fluid V_A, i.e., the Alfvén Mach number ${{M}_{{\rm A}}}\equiv ({{V}_{L}}/{{V}_{{\rm A}}})=1$. The GS95 model was later generalized for both sub-Alfvénic, i.e., ${{M}_{{\rm A}}}\lt 1$, and super-Alfvénic, i.e., ${{M}_{{\rm A}}}\gt 1$, cases (see Lazarian & Vishniac 1999 and Lazarian 2006) and thus the simulations used in this paper (which are sub- and super-Alfvenic) should be understood in this context. In the next two paragraphs we provide a brief synopsis of the differences of single-fluid turbulence with ${{M}_{{\rm A}}}\gt 1$ and ${{M}_{{\rm A}}}\lt 1$.

For ${{M}_{{\rm A}}}\gt 1$ magnetic fields are not dynamically important at the largest scales and hence turbulence from the driving scale (L) to a transition scale L_A follows an isotropic cascade. At scales smaller than L_A, critical balance occurs and scale-dependent anisotropy proceeds down to the dissipation range. Scale L_A is given by

$$\begin{eqnarray}&&{{L}_{{\rm A}}}=L{{({{V}_{{\rm A}}}/{{V}_{L}})}^{3}}=LM_{{\rm A}}^{-3}\end{eqnarray} \tag{ 4 }$$

and the relationship between parallel and perpendicular scales that occurs at scales smaller than L_A is

$$\begin{eqnarray}&&{{l}_{\parallel }}\sim L{{({{l}_{\bot }}/L)}^{2/3}}M_{{\rm A}}^{-1},\;\;\;{{M}_{{\rm A}}}\gt 1,\end{eqnarray} \tag{ 5 }$$

where $\parallel$ and $\bot$ are relative to the direction of the local magnetic field.

Similarly, for ${{M}_{{\rm A}}}\lt 1$, turbulence does not obey the GS95 scaling starting at scale L, but from a smaller scale L_trans, as it transitions from weak to strong turbulence

$$\begin{eqnarray}&&{{L}_{{\rm trans}}}\sim L{{({{V}_{L}}/{{V}_{{\rm A}}})}^{2}}\equiv LM_{{\rm A}}^{2}.\end{eqnarray} \tag{ 6 }$$

In the range $[L,{{L}_{{\rm trans}}}]$ the turbulence is "weak," meaning that it is dominated by interacting MHD waves, rather than eddies.

For scales less than L_trans the turbulence is eddy-like (i.e., strong) and it follows a GS95-type scaling

$$\begin{eqnarray}&&{{l}_{\parallel }}\sim L{{({{l}_{\bot }}/L)}^{2/3}}M_{{\rm A}}^{-4/3},\;\;\;{{M}_{{\rm A}}}\lt 1.\end{eqnarray} \tag{ 7 }$$

CL02, CL03, and Kowal & Lazarian (2010) focused their analysis and comparisons with theoretical predictions of single-fluid MHD turbulence simulations. Ions and neutrals are generally expected to behave as a single MHD fluid at length scales from L_trans or L_A to scale ${{L}_{{\rm AD}}}={{V}_{{\rm A}}}/\alpha {{\rho }_{i}}$. At scales smaller than L_AD, the ion and neutral energetics separate and MHD turbulence may dissipate (Zweibel & Josafatsson 1983; Balara 1996; Klessen et al. 2000).

All of these abovementioned studies assumed that the diffusivities arising from viscosity and resistivity are the same, i.e., that the Prandtl number of turbulence is in unity.⁷ However, viscosity induced by neutrals can result in a very different regime of turbulence. The high Prandtl number MHD turbulence was described theoretically in LVC04 with numerical simulations published in Cho et al. (2002). In this regime, magnetic and kinetic energy spectra are different from theoretical predictions for the kinetic energy ${{E}_{K}}\sim {{k}^{-4}}$ and magnetic energy ${{E}_{M}}\sim {{k}^{-1}}$ (LVC04). These predictions agree with single-fluid MHD numerical simulations and we therefore use these studies as a touchstone for our analysis. At the same time, it would be an oversimplification to assume that turbulence in a partially ionized gas is equivalent to high Prandtl number turbulence. It can have some of the same features, but ion-neutral damping of turbulent motions and the decoupling of ion and neutrals, change the picture in a significant way (LG01; LVC04; Xu et al. 2014). Therefore in what follows we numerically investigate how MHD turbulence evolves in the partially ionized gas.

We use the simulations first presented in Meyer et al. (2014) and reference that work for the details of the numerical setup, here providing only the essential points that are relevant for this work.

The MHD capabilities of the RIEMANN code (Balsara 1998a, 1998b, 2004, 2010, 2012; Balsara & Spicer 1999a, 1999b) have recently been upgraded to treat two-fluid MHD. In two-fluid MHD, the ambipolar effects are modeled with a neutral fluid that obeys the isothermal Euler equations and an ionized fluid that obeys the isothermal MHD equations (Inoue & Inutsuka 2008; Tilley & Balsara 2008; Tilley et al. 2012). The large Alfvén velocities are typically dealt with by invoking the "heavy ion approximation" (HIA, Oishi & Mac Low 2006; Li et al. 2008a), which increases the ion mass in order to reduce the Alfvén speed. However, it was shown and discussed in Tilley & Balsara (2010) that the HIA can alter the dispersion relationship and hence may be missing key physics.

The two-fluid version of RIEMANN, which does not use the HIA, has been applied to astrophysical turbulence in prior papers (Tilley & Balsara 2010; Meyer et al. 2014). A variety of higher order reconstruction methods are available in our code and we used r = 3 WENO reconstruction (Jiang & Shu 1996; Balsara & Shu 2000) because it represents a good compromise between accuracy and speed. The density of the neutrals was scaled to unity and the mean molecular weights of the neutrals and ions were given by ${{\mu }_{N}}=2.3$ amu and ${{\mu }_{i}}=29$ amu (corresponding to HCO+) respectively. Consequently, for a particular ionization fraction χ , the ion density is given by ${{\rho }_{i}}=\chi {{\rho }_{N}}{{\mu }_{i}}/{{\mu }_{N}}$ . The magnetic field was initialized along the x-direction. The simulations were started with uniform density in the ions and neutrals, and they were forced via random Gaussian fluctuations (peaked at k = 2 and spanning $1\leqslant k\leqslant 4$) into a turbulent state so that the desired Mach numbers were represented in the velocity field. Once a steady state was reached, the simulations were continued for a few further turnover times so that  the statistics of the two-fluid turbulence could be gathered.

In this paper, we investigate three forced turbulence simulations at resolution 512³. We set the Alfvén speed for the combined fluid to be one, three, and six times the sound speed. We refer to these runs as A1, A3, and A6, respectively. Three obvious scales of interest exist in these simulations: the turbulence driving scale (L), the ambipolar diffusion scale (L_AD), and the numerical dissipation scale, which begins around 10–20 grid units. We note that in the regime of strong MHD turbulence, one must think of the ambipolar diffusion scale in both the direction parallel and perpendicular to the mean magnetic field, i.e., L$_{{\rm AD}\parallel }$ and L$_{{\rm AD}\bot }$, respectively. L$_{{\rm AD}\parallel }$ as set in the code is given by L$_{{\rm AD}\parallel }={{v}_{{\rm A}}}/\alpha {{\rho }_{i}}$. L$_{{\rm AD}\bot }$ is given by either Equation (5) or (7) depending on if the turbulence is super-Alfvénic or sub-Alfvénic, respectively. We list values of L$_{{\rm AD}\bot }$ in Table 1 for the range of driving scales, i.e., from L = 512 to 128. Table 1 also provides the values for L, L$_{{\rm AD}\parallel }$, and the transition scale between weak and strong turbulence (l_trans) for models A3 and A6 as given by Equation (2) and, for the super-Alfvénic model A1, the transition scale between hydrodynamic and MHD turbulence (L_A) given by Equation (4). In columns two, three, and four of Table 1, we list the sonic Mach number and the volume averaged neutral-ion Alfvén Mach number.

Table 1.  Description of the Simulation Parameters

Run	Ms	$\langle {{M}_{{\rm A}}}\rangle$	χ	L$_{{\rm AD}\bot }$	L$_{{\rm AD}\parallel }$	LA	ltrans	LDrive
A1	3.0	1.8	10−4	58–9	40	77	N/A	512–128
A3	2.5	0.7	10−4	33–16.4	80	N/A	256	512–128
A6	2.5	0.4	2x10−4	11–5.3	80	N/A	86	512–128

Download table as:  ASCIITypeset image

We perform mode decomposition of the velocity field of the ions and neutrals as described in CL03. We briefly summarize this procedure and direct the reader to the original work by CL03 for a more detailed presentation. The slow and fast velocity components can be obtained by projecting the velocity Fourier component v_k onto ${{\xi }_{s}}$ and ${{\xi }_{f}}$, where ${{\xi }_{f}}$ and ${{\xi }_{s}}$ are the basis vectors for the slow and fast modes, respectively, which lie in the plane defined by the mean magnetic field direction ${{{\boldsymbol{B}} }_{0}}$ and ${\boldsymbol{k}}$. They are defined as

$$\begin{eqnarray}&&\widehat{{{\xi }_{s}}}\propto (-1+{{\alpha }_{a}}-\sqrt{D}){{k}_{\parallel }}\widehat{{{{\boldsymbol{k}} }_{\parallel }}}+(1+\alpha -\sqrt{{\boldsymbol{D}} }){{{\boldsymbol{k}} }_{\bot }}{{{\boldsymbol{\hat{k}}} }_{\bot }},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{{\hat{\xi }}_{f}}\propto (-1+{{\alpha }_{a}}+\sqrt{D}){{k}_{\parallel }}\widehat{{{{\boldsymbol{k}} }_{\parallel }}}+(1+\alpha +\sqrt{{\boldsymbol{D}} }){{{\boldsymbol{k}} }_{\bot }}{{{\boldsymbol{\hat{k}}} }_{\bot }},\end{eqnarray} \tag{ 9 }$$

and the basis vector for the Alfvén mode is

$$\begin{eqnarray}&&\widehat{{{\xi }_{{\rm A}}}}=\widehat{{{{\boldsymbol{k}} }_{\bot }}}\times \widehat{{{{\boldsymbol{k}} }_{\parallel }}}\end{eqnarray} \tag{ 10 }$$

where $D={{(1+{{\alpha }_{a}})}^{2}}$ − $4{{\alpha }_{a}}{\rm cos} (\theta )$, ${{\alpha }_{a}}=c_{s}^{2}/V_{{\rm A}}^{2}=\beta (\gamma /2)$. θ is the angle between ${\boldsymbol{k}}$ and ${{{\boldsymbol{B}} }_{0}}$, c_s is the sound speed, and γ is the adiabatic index. These expressions for the mode decomposition basis are rigorously derived in Appendix A of CL03 and illustrated in their Figure 2.

4.1. Velocity Scaling without Mode Decomposition

First we investigate the structure functions of the velocity field in the local frame relative to the mean magnetic field. In this subsection we present our results without  mode decomposition.

We calculate the second-order structure functions, SF₂:

$$\begin{eqnarray}&&{\rm S}{{{\rm F}}_{2}}({{l}_{\bot }},{{l}_{\parallel }})=\langle {{[{\boldsymbol{v}} ({\boldsymbol{x}} +{\boldsymbol{l}} )-{\boldsymbol{v}} ({\boldsymbol{x}} )]}^{2}}\rangle ,\end{eqnarray} \tag{ 11 }$$

where the angle brackets denote the spatial average over ${\boldsymbol{x}}$ and ${\boldsymbol{l}} ={{l}_{\bot }}{{{\boldsymbol{e}} }_{\bot }}+{{l}_{\parallel }}{{{\boldsymbol{e}} }_{\parallel }}$. The vectors ${{{\boldsymbol{e}} }_{\bot }}$ and ${{{\boldsymbol{e}} }_{\parallel }}$ are unit vectors that are perpendicular and parallel to ${{{\boldsymbol{B}} }_{L}}$ respectively. The vector ${{{\boldsymbol{B}} }_{L}}$ denotes the local mean magnetic field.

We stress that we separately obtain the structure functions for cases in which the axis is aligned parallel and perpendicular with the local mean field in a similar manner to past studies such as Cho et al. (2002), CL02, CL03, and Kowal & Lazarian (2010). The GS95 theory of Alfvénic turbulence summarized in the previous section predicts an inertial range scaling of ${{l}_{\parallel }}\propto l_{\bot }^{2/3}$, which implies that scalings in the second-order structure function as ${\rm S}{{{\rm F}}_{2}}(0,{{l}_{\parallel }})\propto {{l}_{\parallel }}$ and ${\rm S}{{{\rm F}}_{2}}({{l}_{\bot }},0)\propto l_{\bot }^{2/3}$. Deviations from these scaling relationships at small scales indicate the damping of the Alfvénic cascade.

We present the structure function analysis for the ion velocity field in Figure 1 and for the neutral velocity field in Figure 2. In the ions for all simulations, the GS95 scaling is observed until L$_{{\rm AD},\parallel }$. Model A1 (top panel) is super-Alfvénic and thus the transition to strong turbulence and the GS95 scaling does not begin until scale L$_{{\rm A}}\approx 77$ grid units; however, past this scale a scaling of ${{l}^{2/3}}$ is clearly seen down to L$_{{\rm AD},\parallel }$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Models A3 and A6 (middle and bottom panels, respectively), are sub-Alfvénic and thus are in the regime of weak or wave-like turbulence until scale L_trans given in Table 1. Both A3 and A6 have L$_{{\rm AD},\parallel }$ and L$_{{\rm AD},\bot }$ separated by a larger dynamic range of scales than model A1 due to sub-Alfvénic turbulence. At scales smaller than the driving scale and L_trans, Alfvénic turbulence in the ions develops and is clearly seen in the scaling relation ${{l}_{\parallel }}\sim {{l}^{2/3}}$. For models A3 and A6, the range of scales for GS95 is limited, as either the driving scale or L_trans is close to L$_{{\rm AD},\parallel }$. However, for model A1 the range of L_AD is separated from both L and L_A , with a sufficient dynamic range to generate an Alfvénic cascade with the GS95 scaling in the ions.

The neutral structure functions shown in Figure 2 also have GS95 scalings about L$_{{\rm AD},\parallel }$ in all three simulations studied. This suggests that the neutrals and ions behave as a single MHD fluid until they decouple. At $L={{L}_{{\rm AD},\parallel }}$ neutrals decouple from the ions and no longer exhibit the scaling relations of an MHD cascade. The ${{l}^{2/3}}$ scaling is present in the neutrals until L = 40 for model A1. Models A3 and A6 transition to an isotropic l¹ scaling at larger scales than model A1, with model A6 transitioning to l¹ at slightly larger scales than model A3. Since models A3 and A6 have identical values for L$_{{\rm AD},\parallel }$ this suggests that the driving scale and/or Alfvén Mach number also play a role in the dissipation scale of partially ionized fluid turbulence. It is also interesting to note that in the A1 model the ions exhibit a greater dynamic range than the neutrals over which the ${{l}^{2/3}}$ is observed. This suggests that the ions may continue the MHD cascade even after they have decoupled from the neutrals. As most of the energy in the MHD cascade lies in the Alfvén modes (see Kowal & Lazarian 2010) it is important to determine if the damping of the cascade at these scales is due to damping in the Alfvén modes.

4.2. Velocity Scaling with Mode Decomposition

In order to address the nature of the Alfénic cascade in ion-neutral turbulence we must separate the Alfvén modes from the fast and slow modes. We perform the mode decomposition on the the ion and neutral velocity fields to extract the Alfvén modes and then apply the second-order structure function analysis as shown in Figures 1 and 2.

We show the ion Alfvén mode structure functions in Figure 3. The Alfvén modes of model A1 (top panel) obey the GS95 scalings to scales smaller than the ambipolar diffusion scale, which indicates that the Alfvénic cascade can continue to scales smaller than L_AD. This suggests that turbulence does not necessarily damp at the ambipolar diffusion scale but can continue to the viscous dissipation scale and would not necessarily provide a characteristic scale for star formation, in agreement with Oishi & Mac Low (2006).

Zoom In Zoom Out Reset image size

Unlike model A1, Alfvénic turbulence in the ions begins to damp at larger scales for models A3 and A6. Alfvénic turbulence in models A3 and A6 damps at or before the perpendicular ambipolar diffusion scale The larger dynamic range of scales between L$_{{\rm AD},\parallel }$ and L$_{{\rm AD},\bot }$ coupled with the limited dynamic range between the driving scale and L$_{{\rm AD},\parallel }$, increases the effect of damping due to ion-neutral collision. In addition, the strength of the damping due to ion-neutral collisions corresponds to V_A² (Klessen et al. 2000), which explains why the damping is the strongest for model A6.

Figure 4 plots the structure functions of the Alfvén mode neutral velocity field of each model in Table 1. The organization of the figure is identical to Figure 3.

Zoom In Zoom Out Reset image size

In all of the simulations, the Alfvén mode scaling in the neutrals is seen at large scales but breaks off from the GS95 scaling at scales smaller than ${{L}_{{\rm AD},\parallel }}$ and/or ${{L}_{{\rm AD},\bot }}$. In all three models the perpendicular motions do not follow ${{l}_{\parallel }}\sim {{l}^{2/3}}$ at scales smaller than ${{L}_{{\rm AD},\parallel }}$ but we see another effect, which is that both ${{l}_{\parallel }}$ and ${{l}_{\bot }}$ begin to steepen to slopes greater than l¹. This is close to the numerical dissipation scale, but also most likely an unphysical effect due to the fact that once the neutrals decouple from the ions we can no long apply mode decomposition, as Alfvén modes will not exist in the neutral cascade.

4.3. Power Spectrum

Finally, we investigate the kinetic energy spectra for simulations with and without mode decomposition. We calculate the energy power spectra as:

$$\begin{eqnarray}&&P({\boldsymbol{k}} )=\mathop{\sum }\limits_{{\boldsymbol{k}} ={\rm const}.}\tilde{F}({\boldsymbol{k}} )\cdot {{\tilde{F}}^{*}}({\boldsymbol{k}} ),\end{eqnarray} \tag{ 12 }$$

where $\tilde{F}$ is the Fourier transform of the velocity field.

In hydrodynamic turbulence the viscous damping scale sets a minimal scale for turbulent motions and the kinetic energy power spectrum is dissipated exponentially. This marks the end of the hydrodynamic cascade, but in MHD turbulence the viscous damping scale is not the end of magnetic structure evolution. On these scales magnetic field structures will be created by shear and magnetic tension. As a result, LVC04 predicted a power-law tail in the energy distribution, rather than an exponential cutoff.

Figure 5 shows the neutral kinetic energy power spectrum (black lines) and the ion kinetic energy power spectrum (red lines) for models A1, A3, and A6. For model A1, the neutral Alfvén modes damp out well before the ion Alfvén modes. The ions retain the ${{k}^{-5/3}}$ scaling even after the neutrals have decoupled. For models A3 and A6, there is a small portion of the spectrum with the ${{k}^{-5/3}}$ scaling for the Alfvén mode in the ions that is not noticeable in the neutrals due to the limited dynamic range over which the Alfvén modes propagate. In the high k damped regime, the slopes approach ${{k}^{-4}}$ rather than behaving as  exponential decays of energy. This result should be confirmed with higher resolution simulations in the future.

Zoom In Zoom Out Reset image size

What are the possible damping mechanisms that can explain the dampening behavior seen in these simulations? The role of neutral-ion damping in MHD turbulence has previously been discussed in the context of the ISM (in particular, see Spangler 1991; Minter & Spangler 1997; LG01; LVC04). In a partially ionized medium a combination of neutral particle viscosity and ion-neutral collisional coupling drives damping. Neutral-ion friction will compete with and eventually dominate the Alfvén wave restoring force, which will damp oscillations in the magnetic field. As soon as the neutral-ion collisional rate is approximately equal to the eddy turnover time, the neutrals will begin to form a hydrodynamic cascade. At this scale, and all smaller scales, the ionic fluid motions will damp at the rate of ion-neutral collisions. The damping mechanism in A3/A6 is most likely the result of ion-neutral collisions; however, an additional damping mechanism will be provided by neutral viscosity, which may be non-negligible in dissipating the MHD cascade in partially ionized gases (LVC04). A detailed investigation of the dominance of one damping mechanism over another is beyond the scope of this work; however, comparisons with Equations (37) and (46) of LG01 suggest that the neutral mean free path for our models A3 and A6 is smaller than A1 and hence damping by neutral collisions may play an important role.

Supersonic MHD turbulence is known to play a critical role in both the support of GMCs and their subsequent small-scale collapse. Ambipolar diffusion has been thought to be the most important dissipation mechanism for turbulence in GMCs and could set characteristic mass and length scales for star formation. However this study shows that turbulence, in particular the Alfvénic cascade, in a partially ionized medium does not necessarily damp at the ambipolar diffusion scale L_AD. Other authors, such as Oishi & Mac Low (2006), for example, investigated two-fluid simulations and found that L_AD did not set a characteristic length scale for mass sizes or dissipation and that turbulence may proceed to smaller scales. Previous studies, however, often used the HIA which dramatically changes the mode propagation characteristics of two-fluid turbulence. We show in this work that Alfvénic modes can survive to scales much smaller than L_AD in some cases. This finding also implies that an important diffusion process called reconnection diffusion, which is mediated by turbulence (see Lazarian & Vishniac 1999; LVC04) may persist at small scales.

While this work is a theoretical and numerical study, our results have very important direct implications for observations. We have shown that the Alfvén Mach number is a critical parameter for obtaining the damping and scaling characteristics of the turbulence that may play a role in the interpretation of velocity dispersion methods for obtaining the ambipolar diffusion scale (Hezareh et al. 2014). Our results work well in the context of interpreting the observational signatures of decoupling found in Tilley & Balsara (2010) and Meyer et al. (2014). Meyer et al. (2014) showed that the PDFs of sub-Alfvénic turbulence (i.e., models A3 and A6) were different in the ions and neutrals when looking at a line of sight (LOS) parallel to B, but similar with a sightline perpendicular to B. As we have shown, the perpendicular decoupling scale for these models occurs at much smaller scales along the cascade, hence density fluctuations may be coupled in this direction even after the cascade has begun to damp. Contrary to the sub-Alfvénic behavior, Meyer et al. (2014) showed that the PDFs of ions and neutrals for model A1 are very similar and are generally LOS independent, which is expected from our results as turbulence persists in this model to scales smaller than the decoupling scale.

One of our most important findings suggests that the MHD turbulence cascade does not simply damp at the decoupling scale L_AD but rather depends on other parameters such as the Alfvén Mach number, which determines the level of anisotropy in the cascade. Furthermore, we show that neutrals and ions behave as a single MHD fluid at certain scales, which can explain the correlation between 21 centimeter filaments and the magnetic field (Clark et al. 2014). There are several methods for obtaining the Alfvén Mach number that extend beyond direct observational techniques, such as Zeeman spiting (Crutcher et al. 2009). These include anisotropy in the structure function of velocity centroids (Esquivel & Lazarian 2005, 2011; Burkhart et al. 2014), Principle Component Analysis (Correia et al. 2014), and the bispectrum (Burkhart et al. 2009). The power spectrum of partially ionized fluids can also be an indication of the Alfvénic state of the gas, as suggested by Hezareh et al. (2010) and Meyer et al. (2014) and shown in Figure 5.

Our work shows that the use of scale L_AD as the scale for the suppression of the MHD cascade in partially ionized gasses may not be the final story. It is clear that a more detailed theoretical study of turbulence dissipation in partially ionized gas is required in order to understand the precise damping mechanism for a given parameter space (see more in Xu et al. 2014).

We studied the structure function scaling relations and power spectra of the velocity field Alfvén modes of two-fluid (ion-neutral) MHD turbulence. We investigated at what scales MHD turbulence is damped in the partially ionized gas, treating ions and neutrals separately. We showed that the GS95 scalings and anisotropy are present in the velocity structure functions of the ions and neutrals in all models studied; however, the dynamic range of scales over which GS95-like turbulence is seen varies based on the Alfvén Mach number and the range of the parallel and perpendicular ambipolar diffusion scales bounded by the viscous dissipation scale and the driving scale. We found that Alfvénic turbulence does not necessarily damp at the ambipolar diffusion scale in super-Alfvénic turbulence and that the damping depends on other turbulence parameters such as the Alfvén Mach number and driving scale. We also showed that the power spectrum of the kinetic energy in the damped regime is consistent with a ${{k}^{-4}}$ scaling, as predicted by LVC04.

B.B. acknowledges support from the NASA Einstein Fellowship. A.L. and B.B. thank the Center for Magnetic Self-Organization in Astrophysical and Laboratory Plasmas for financial support and acknowledge financial support from the Physics Graduate Program/UFRN, at Natal, for hospitality. A.L. is supported by the NSF grant AST 1212096. D.S.B. acknowledges support via NSF grants NSF-AST-1009091, NSF-ACI-1307369, and NSF-DMS-1361197. D.S.B. also acknowledges support via NASA grants from the Fermi program as well as NASA-NNX 12A088G. Several simulations were performed on a cluster at UND that is run by the Center for Research Computing. Computer support for NSF's XSEDE computing resources is also acknowledged. J.C.'s work is supported by the National R & D Program through the National Research Foundation of Korea (NRF), funded by the Ministry of Education (No. 2011-0012081)