MULTIFLUID MAGNETOHYDRODYNAMIC TURBULENT DECAY

T. P. Downes; S. O'Sullivan

doi:10.1088/0004-637X/730/1/12

1. INTRODUCTION

Turbulence is recognized as a possible source of support against gravitational collapse for molecular clouds. The precise role and source of the observed motions interpreted as evidence of turbulence in these clouds has been studied extensively by many researchers (see the reviews of Mac Low & Klessen 2004; Elmegreen & Scalo 2004). Clearly, if turbulence can support molecular clouds then it can influence star formation in terms of rate, efficiency and initial mass function (Elmegreen 1993; Klein et al. 2003).

Many studies of turbulence in molecular clouds have focused on ideal magnetohydrodynamics (MHD) as an approximation of the physics governing this system (Mac Low et al. 1998; Mac Low 1999; Ostriker et al. 2001; Vestuto et al. 2003; Gustaffson et al. 2006; Glover & Mac Low 2007; Lemaster & Stone 2008, 2009). The assumption of ideal MHD, while desirable for technical reasons, is perhaps risky in the context of turbulence. The reason for this is that while ideal MHD is valid in molecular clouds on fairly large length scales, on shorter length scales non-ideal effects are thought to become significant (Wardle 2004; Oishi & Mac Low 2006). Given that turbulence in three dimensions involves the transfer of energy from large scales to ever smaller scales, the assumption of ideal MHD will be invalid below some critical spatial scale and the correct nature of the energy cascade may not be observed at this range.

The most important of the non-ideal effects for molecular cloud dynamics is ambipolar diffusion. Some authors (Oishi & Mac Low 2006; Li et al. 2008; Kudoh & Basu 2008) have studied driven MHD turbulence in the presence of ambipolar diffusion. All these authors find that ambipolar diffusion produces significant differences in the properties of the turbulence.

While most likely of lesser significance, it has been suggested that although the Hall resistivity is generally at least an order of magnitude lower than the ambipolar resistivity in molecular clouds (Wardle 2004), its effect should not be ignored. Although relatively weak, it is capable of inducing topological changes in the magnetic field which are quite distinct to any influence caused by ambipolar diffusion. In support of this assertion, we note that researchers working on reconnection and the solar wind have studied the Hall effect in the context of turbulence and found that, although the overall decay rate appears not to be affected, the usual coincidence of the magnetic and velocity fields seen in MHD does not occur at small scales (Matthaeus et al. 2003; Mininni et al. 2006; Servidio et al. 2007). Almost no work has been done on comparing the influences of this effect coupled with that of ambipolar diffusion on turbulence with the exception of Downes & O'Sullivan (2009, hereafter Paper I).

In Paper I, a series of simulations of decaying supersonic non-ideal MHD turbulence incorporating both ambipolar diffusion and the Hall effect were performed. These simulations, however, were constrained in that the resistivities associated with each of ambipolar diffusion, the Hall effect, and the Pederson resistivity were kept fixed in both space and time. The authors found that, at length scales of 0.2 pc, ambipolar diffusion has a significant impact on the decay of the turbulence. The Hall effect was less significant in this respect but does have an influence on the magnetic field at short length scales. Here, we present simulations in which the resistivities are self-consistently calculated from the evolution of both the magnetic field and the densities of all of the component species of the fluid. Using these dynamically evolving resistivities we study the decay of fully multifluid MHD turbulence. This is the first such study presented in the literature, with the exception of the low resolution simulations presented by Downes & O'Sullivan (2008).

The aim of this work is to examine in detail the differences between the decay of ideal MHD turbulence and that of multifluid MHD turbulence with a full tensor resistivity incorporating the effects of ambipolar diffusion, the Hall effect, and Pederson resistivity. We will use the results of Paper I in our discussion of these differences as it represents an intermediate stage between the calculations presented here and those of ideal MHD. This work is new in two respects: notwithstanding Paper I, no previous work has focused on decaying (e.g., undriven) multifluid MHD turbulence and, in addition, no previous work has addressed the issue of turbulence in the presence of both ambipolar diffusion and the Hall effect simultaneously.

In Section 2, we outline the numerical techniques used in this work, as well as the initial conditions and general setup for the simulations while in Section 3 we describe the methods used to analyze the simulation data. In Section 4, we present and discuss the results of our simulations of turbulent decay. Finally, Section 5 contains a summary of our results.

2. NUMERICAL METHOD

As in Paper I, we use the code HYDRA (O'Sullivan & Downes 2006, 2007) to integrate the equations of weakly ionized multifluid MHD (see Section 2.1). We assume that the molecular cloud material we are simulating can be treated as isothermal and that initially the density and magnetic field are uniform. We use the capabilities of HYDRA to extend the physics incorporated in the simulations here beyond those presented in Paper I so that the turbulence here is fully multifluid MHD.

2.1. Equations and Algorithm

We briefly outline the equations and assumptions in our model here but refer the reader to O'Sullivan & Downes (2006, 2007) for a comprehensive description of the underlying assumptions for the weakly ionized model of multifluid MHD.

We assume that the cloud material can be treated as weakly ionized. This is clearly valid for molecular clouds and allows us to ignore the inertia of the charged species (Ciolek & Roberge 2002; Falle 2003). For a system composed of N fluids, one of which is neutral, the equations to be solved are then

$\begin{equation} \frac{\partial \rho _i}{\partial t} + \mbox{\boldmath $\nabla $} \cdot (\rho _i \mbox{\boldmath $q$}_i) = 0,\quad \mbox{ ($1 \le i \le N$), } \end{equation} \tag{ 1 }$

$\begin{equation} \frac{\partial \rho _1 \mbox{\boldmath $q$}_1}{\partial t} + \nabla \cdot (\rho \mbox{\boldmath $q$}_1 \mbox{\boldmath $q$}_1 + a^2\rho \mbox{\boldmath $\mathsf {I}$}) = \mbox{\boldmath $J$}\times \mbox{\boldmath $B$}, \end{equation} \tag{ 2 }$

$\begin{equation} \frac{\partial \mbox{\boldmath $B$}}{\partial t} + \nabla \cdot (\mbox{\boldmath $q$}_1\mbox{\boldmath $B$}-\mbox{\boldmath $B$}\mbox{\boldmath $q$}_1) = -\nabla \times \mbox{\boldmath $E$}^{\prime }, \end{equation} \tag{ 3 }$

$\begin{equation} \alpha _i \rho _i (\mbox{\boldmath $E$} + \mbox{\boldmath $q$}_i \times \mbox{\boldmath $B$}) + \rho _i \rho _1 K_{i\,1}(\mbox{\boldmath $q$}_1-\mbox{\boldmath $q$}_i) = 0, \quad \mbox{ $2 < i \le N$, } \end{equation} \tag{ 4 }$

$\begin{equation} \nabla \cdot \mbox{\boldmath $B$} = 0, \end{equation} \tag{ 5 }$

$\begin{equation} \nabla \times \mbox{\boldmath $B$} = \mbox{\boldmath $J$}, \end{equation} \tag{ 6 }$

$\begin{equation} \sum _{i=2}^N \alpha _i \rho _i = 0, \end{equation} \tag{ 7 }$

where ρ₁, $\mbox{\boldmath $q$}_1$ , a, $\mbox{\boldmath $B$}$ , and $\mbox{\boldmath $J$}$ are the neutral mass density, neutral velocity, sound speed, magnetic field, and current density, respectively. K_i 1, α_i, and ρ_i (i>1) are the collision coefficients between species i and the neutrals, the charged fluid charge-to-mass ratios, and mass densities, respectively. Equations (1)–(7) express conservation of mass for each fluid, conservation of neutral momentum, the induction equation, force balance for the charged species, the inadmissibility of magnetic monopoles, Faraday's law, and charge neutrality, respectively.

The electric field in the frame of the fluid, $\mbox{\boldmath $E$}^{\prime }$ , is calculated from the generalized Ohm's law for weakly ionized fluids (e.g., Falle 2003; O'Sullivan & Downes 2006) and is given by

$\begin{equation} \mbox{\boldmath $E$}^{\prime } = \mbox{\boldmath $E$}_{\rm O} +\mbox{\boldmath $E$}_{\rm H} +\mbox{\boldmath $E$}_{\rm A}, \end{equation} \tag{ 8 }$

where

$\begin{equation} \mbox{\boldmath $E$}_{\rm O} = (\mbox{\boldmath $J$}\cdot \mbox{\boldmath $a$}_{\rm O})\mbox{\boldmath $a$}_{\rm O}, \end{equation} \tag{ 9 }$

$\begin{equation} \mbox{\boldmath $E$}_{\rm H} = \mbox{\boldmath $J$}\times \mbox{\boldmath $a$}_{\rm H}, \end{equation} \tag{ 10 }$

$\begin{equation} \mbox{\boldmath $E$}_{\rm A} = -(\mbox{\boldmath $J$}\times \mbox{\boldmath $a$}_{\rm A})\times \mbox{\boldmath $a$}_{\rm A}, \end{equation} \tag{ 11 }$

using the definitions $\mbox{\boldmath $a$}_{\rm O}\equiv f_{\rm O} \mbox{\boldmath $B$}$ , $\mbox{\boldmath $a$}_{\rm H}\equiv f_{\rm H} \mbox{\boldmath $B$}$ , $\mbox{\boldmath $a$}_{\rm A}\equiv f_{\rm A} \mbox{\boldmath $B$}$ , where $f_{\rm O}\equiv \sqrt{r_{\rm O}}/B$ , f_H ≡ r_H/B, and $f_{\rm A}\equiv \sqrt{r_{\rm A}}/B$ . r_O, r_H, and r_A are the Ohmic, Hall, and ambipolar resistivities, respectively, and are given by

$\begin{equation} r_{\rm O} = \frac{1}{\sigma _{\rm O}}, \end{equation} \tag{ 12 }$

$\begin{equation} r_{\rm H} = \frac{\sigma _{\rm H}}{\sigma _{\rm H}^2 + \sigma _{\rm A}^2}, \end{equation} \tag{ 13 }$

$\begin{equation} r_{\rm A} = \frac{\sigma _{\rm A}}{\sigma _{\rm H}^2 + \sigma _{\rm A}^2}, \end{equation} \tag{ 14 }$

with the conductivities given by

$\begin{equation} \sigma _{\rm O} = \frac{1}{B}\sum _{i=2}^N \alpha _i \rho _i \beta _i, \end{equation} \tag{ 15 }$

$\begin{equation} \sigma _{\rm H} = \frac{1}{B}\sum _{i=2}^N \frac{\alpha _i \rho _i}{1+\beta _i^2}, \end{equation} \tag{ 16 }$

$\begin{equation} \sigma _{\rm A} = \frac{1}{B}\sum _{i=2}^N \frac{\alpha _i \rho _i \beta _i}{1+\beta _i^2}, \end{equation} \tag{ 17 }$

where β_i is the Hall parameter for species i and is given by

$\begin{equation} \beta _i = \frac{\alpha _i B}{K_{1\,i}\rho _1}. \end{equation} \tag{ 18 }$

As noted by Falle (2003) and O'Sullivan & Downes (2006), the main difficulty with standard numerical techniques for integrating Equation (3) lies with the Hall term. As this term becomes dominant the stable time step goes to zero. However, O'Sullivan & Downes (2006, 2007) presented a novel, explicit numerical method for integrating this term such that the limit on the stable time step is not overly restrictive. We use this "Hall Diffusion Scheme" in this work. Of course, all explicitly differenced diffusion terms give rise to a stable time step which is proportional to Δx², where Δx is the resolution of the simulation. To ameliorate this we use standard subcycling of the Hall terms and super time stepping to accelerate the ambipolar diffusion terms (see Alexiades et al. 1996; O'Sullivan & Downes 2006, 2007).

Equations (1)–(3) are solved using a standard shock-capturing, second-order, finite volume, conservative scheme. The numerical techniques employed in this work are slightly different to those used in Paper I in one respect: we have altered the calculation of the advective fluxes in Equation (3) to use the method suggested by Falle (2003). In Paper I, these fluxes were derived from interface values of the neutral gas velocity and the magnetic field. We find that at high resolutions with variable resistivities the latter method is prone to introducing grid scale features in the solution while the former is not. This undesirable effect was not an issue for the investigations carried out in Paper I since resistivities were fixed. The downside of the described variation between the numerical approaches is that it must be considered as a possible source of discrepancy in comparisons between the results of Paper I and this work. However, in order to provide evidence of the small influence, we have also run a fixed resistivity simulation for this work (see Section 2.2).

Equation (5) is enforced using the method of Dedner et al. (2002). The effects of the diffusive terms in Equation (3) are then incorporated in an operator split fashion.

2.2. Initial Conditions

We examine the decay of MHD turbulence in conditions suitable for dense regions of molecular clouds. The conditions we use are similar to those used in Paper I. We briefly review them here for completeness.

The computational domain is set up as a cube of side L = 0.2 pc. Periodic boundary conditions are enforced on all faces of the simulation domain. The sound speed is set to 0.55 km s⁻¹, the initial density is chosen to be uniform with a value of 10⁶ cm⁻³, and the magnetic field is also initially uniform in the (1, 1, 1) direction with a magnitude of 1 mG. For these conditions, suitable conductivities are σ_O = 1 × 10¹⁰ s⁻¹, σ_H = 10⁻² s⁻¹, and σ_A = 10⁻¹ s⁻¹ (Wardle & Ng 1999). We choose a three-fluid setup for our multifluid simulation: one neutral species and two charged species. The densities, charge-to-mass ratios, and collisional coefficients of the charged species are chosen in order to achieve these conductivities. We choose these particular physical conditions with a view to maximizing the influence of the Hall effect in our simulations (Wardle & Ng 1999). In this way, we hope to determine whether the Hall effect is ever likely to be important in molecular cloud turbulence.

The initial velocity field is defined to be the sum of waves with 64 wavevectors each with random amplitude and phase, i.e.,

$\begin{equation} q_{\alpha } = \sum _{j=0}^{64} A_{\alpha, j} \cos (\mbox{\boldmath $k$}_j\cdot \mbox{\boldmath $x$}+ \phi _{\alpha, j}), \end{equation} \tag{ 19 }$

where α defines the component (x, y, or z) of the appropriate quantity, A_α,j and ϕ_α,j are the random amplitudes and phases, and $\mbox{\boldmath $x$}$ is the position vector. We restrict the velocity field to be solenoidal (e.g., non-compressional). By construction the mean velocity over the domain is zero.

Table 1 presents a complete list of the various simulations carried out in this work. The nomenclature we employ in referencing the simulations is xx-c where xx denotes the type of physics (e.g., a standard molecular cloud run is "mc," ideal MHD is "mhd," etc.) and c is the resolution used. The initial root mean square (rms) of the field is chosen to be 5 with a corresponding Alfvénic Mach number of approximately 1.9. In addition to the four multifluid MHD simulations run at different resolutions, we also run four further simulations. The first is an ideal MHD simulation (mhd-512) which we use for comparison purposes and the other two (ambi-512 and hall-512) only incorporate one of either ambipolar diffusion or the Hall effect, respectively. The final case (fr-512) is a fixed resistivity simulation used to make contact with the simulations of Paper I. We use these four latter simulations to investigate separately the influence of each non-ideal effect.

Table 1. Definition of the Initial Conditions Used in the Simulations in this Work

Simulation	Mach Number^a	Resolution	Comment
mc-64	5	64³	⋅⋅⋅
mc-128	5	128³	⋅⋅⋅
mc-256	5	256³	⋅⋅⋅
mc-512	5	512³	⋅⋅⋅
ambi-512	5	512³	r_H = 0
hall-512	5	512³	r_A = 0
mhd-512	5	512³	Ideal MHD
fr-512	5	512³	As mc-512 but for fixed resistivities

Note. ^aInitial rms Mach number of the flow.

Download table as: ASCII Typeset image

3. ANALYSIS

In this section, we discuss the method of analysis of the output of the simulations described in Section 2.2. The aim of this paper is to investigate the decay rate of supersonic turbulence in molecular clouds. Hence, the main analysis carried out is of the kinetic, magnetic, and total energy as functions of time. These quantities are defined respectively as

$\begin{equation} e_{\rm k} = \int _{\rm domain} \rho |\mbox{\boldmath $q$}|^2\,dV \end{equation} \tag{ 20a }$

$\begin{equation} e_{\rm b} = \int _{\rm domain} \frac{|\mbox{\boldmath $B$}|^2}{2}\,dV - \frac{\langle \mbox{\boldmath $B$}\rangle ^2}{2}V \end{equation} \tag{ 20b }$

$\begin{equation} e_{\rm tot} = e_b + e_k, \end{equation} \tag{ 20c }$

where V is the volume of the computational domain.

We also calculate the mass-weighted average Mach number, defined by

$\begin{equation} M = \frac{1}{a}\left\lbrace \sigma _x^2 + \sigma _y^2 + \sigma _z^2\right\rbrace ^{1/2} \end{equation} \tag{ 21 }$

where a is the sound speed and the velocity dispersions, σ_α, are defined by

$\begin{equation} \sigma _\alpha = \left\lbrace \frac{\left<\rho q_\alpha ^2\right>}{\left<\rho \right>}\right\rbrace ^{1/2}, \end{equation} \tag{ 22 }$

where α is either x, y, or z and the angle brackets denote averaging over the computational domain (see Lemaster & Stone 2009).

In Section 4.3, we present the power spectra for the velocity, density, and magnetic field for each of the 512³ simulations. These spectra are calculated by taking the power spectra in the x-, y-, and z-directions separately and summing the power over the interval $k \le |\mbox{\boldmath $k$}| < k + \Delta k$ (where we take Δk = 1). This gives us some insight into the scale of structures being formed by the turbulence for the various initial conditions and range of physics examined.

Finally, in Section 4.4 we calculate the velocity dispersion as a function of length scale, l. For these purposes we define the velocity dispersion to be

$\begin{equation} \sigma (l) = \left\lbrace \big\langle \sigma _x^2(l)\big\rangle _{\rm domain} + \big\langle \sigma _y^2(l)\big\rangle _{\rm domain} + \big\langle \sigma _z^2(l)\big\rangle _{\rm domain}\right\rbrace ^{\frac{1}{2}} \end{equation} \tag{ 23 }$

where

$\begin{equation} \sigma _\alpha (l) = \left\lbrace \left<q_\alpha ^2\right>_l-\left< q_\alpha \right>_l^2\right\rbrace ^{\frac{1}{2}}, \end{equation} \tag{ 24 }$

where 〈 · 〉_l indicates an average taken over a cube of side l in the simulation domain and 〈 · 〉_domain indicates averaging of the quantity over all such non-overlapping cubes within the domain.

4. RESULTS

Each of the simulations detailed in Table 1 was run for one sound crossing time, t_c = 3.56 × 10⁴ yr, of the simulation domain. All analyses were carried out for t ⩾ 0.2 t_c (i.e., one flow crossing time) at which point we expect significant turbulent mixing to have taken place and the system's memory of the initial state to be largely forgotten.

As an illustration of the differences between an ideal MHD turbulence simulation and a multifluid MHD simulation, Figure 1 contains plots of the density distribution at t = t_c in a slice through the computational domain for simulations mhd-512 and mc-512. It is clear that there is much less fine structure in mc-512. Also shown is the same slice for simulation fr-512 (i.e., fixed resistivities). While the similarities in terms of the levels of structure are relatively small between mc-512 and fr-512, there are clear differences between the distributions indicating that calculating the resistivities self-consistently has some impact on the dynamics of the system.

In Figure 2, we present plots of the ambipolar and Hall resistivities for simulation mc-512 for the same times and slices as Figure 1. It is clear that the resistivities vary considerably throughout the computational domain with the features strongly correlated with the features in the density distribution. We also show $\eta \equiv \frac{r_{\rm A}}{|r_{\rm H}|}$ to give an indication of the relative importance of each of the resistivities. This, as we shall see, is an important parameter. Finally, Figure 3 is the same as Figure 2 except that the data are taken at time t = 0.2t_c—i.e., after one flow crossing time. Here we can see that the variation of the resistivities in space is dramatic with, for example, the ambipolar resistivity varying by almost four orders of magnitude with η varying by around two orders of magnitude.

**Figure 2.** Plot of the normalized ambipolar resistivity (top panel), normalized Hall resistivity (middle panel), and $\eta \equiv \frac{r_{\rm A}}{|r_{\rm H}|}$ . Note that the scale is logarithmic in the top and middle panels while it is linear for the bottom panel.
Download figure:
Standard image High-resolution image

**Figure 3.** Same as Figure 2 but the snapshots are taken at t = 0.2 t_c. Note that the scale here is logarithmic in each case including the bottom panel.
Download figure:
Standard image High-resolution image

4.1. Resolution Study

Four simulations identical in every way except for the resolution were run. Specifically, the resolutions used were 64³, 128³, 256³, and 512³. We now focus our attention on how the energy decay behaves with resolution. Figure 4 contains plots of the kinetic energy as a function of time for each of the simulations in the resolution study. It is clear that the lower the resolution, the faster the decay—this is what one would expect since lower resolution results in a higher numerical viscosity and hence one expects faster dissipation of energy.

**Figure 4.** Log-scale plot of the kinetic energy (normalized to its initial value) for each of the simulations in the resolution study.
Download figure:
Standard image High-resolution image

Simulations mc-256 and mc-512 are, however, quite similar in terms of the energy decay with a maximum relative difference of around 10% between the kinetic energies in the simulations at any one time—an almost identical result to that obtained from the resolution study in Paper I. This is notable since in Paper I the resistivities were kept constant in space and time, whereas here the resistivities locally increase significantly during the course of the simulations. A reasonable inference is that the influence of local variations in resistivities averages out in some sense on the global scale.

As we shall see later, however, the effect of spatially varying resistivities is noticeable in properties such as the power spectrum of the density and magnetic field.

The various energy decay rates can be modeled approximately as power laws in time, t^−β. Fitting the kinetic energy, β_K, magnetic energy, β_B, and total energy, β_Tot, as functions of time in this way we obtain the values given in Table 2. These data confirm quantitatively what can be observed in Figure 5 and extend it to the decay of the energy in magnetic perturbations: increasing resolution decreases the rate of energy decay, but the difference between the 256³ and 512³ simulations is relatively minor. We note in passing that the decay in the energy in magnetic perturbations is considerably more sensitive to resolution than kinetic energy: β_B varies between 1.43 and 1.28 while β_K only varies between 1.38 and 1.34. We know that the resistivities are a critical factor in determining the decay of the magnetic energy since they facilitate loss of magnetic energy through reconnection. We attribute the extra sensitivity to the resolution of β_B to the necessity of properly resolving the diffusive effects in the induction equation, including the variation of the resistivities themselves. The variation of the resistivities throughout the computational domain is significant (see Figure 2) and it is interesting to note that β_B in simulations with fixed resistivities is not so sensitive to resolution (see Paper I).

**Figure 5.** Log-scale plot of the kinetic energy (normalized to its initial value) for each of the 512³ simulations.
Download figure:
Standard image High-resolution image

4.2. Energy Decay

We now discuss the behavior of the kinetic and magnetic energy in our multifluid simulations and compare with those in Paper I.

4.2.1. Kinetic Energy Decay

Figure 5 contains plots of the decay of kinetic energy with time for the 512³ resolution simulations outlined in Table 1. Note that the kinetic energy decay in all of the simulations is very similar until around t ≈ 0.02 t_c. This is because at such early times compressions are only just starting to form and so the non-ideal terms in the induction equation have had almost no effect on the dynamics. The subsequent energy decay of simulations mhd-512 and hall-512 is almost identical to each other. The energy decay of simulations mc-512 and ambi-512 is virtually identical over the full plotted range, while the data plotted for fr-512 coincide with the former simulations for times in the range 0.3 t_c ⩽ t ⩽ t_c. We have fitted the kinetic energy decay as a power law in the range 0.2t_c ⩽ t ⩽ t_c (i.e., after approximately one initial flow crossing time) and the exponents are given in the first column of Table 2.

Table 2. Values of the Exponent for the Kinetic, Magnetic, and Total Energy Decay for the Simulations Presented in this Work

Simulation	β_K	β_B	β_Tot
mc-64	1.38	1.43	1.39
mc-128	1.37	1.35	1.36
mc-256	1.35	1.30	1.33
mc-512	1.34	1.28	1.32
hall-512	1.12	1.05	1.09
ambi-512	1.35	1.30	1.33
mhd-512	1.12	1.06	1.10
fr-512	1.30	1.28	1.30

Note. These exponents are calculated by fitting the data over the time interval [0.2t_c, t_c].

Download table as: ASCII Typeset image

It is clear that the presence of ambipolar diffusion has a significant impact on the behavior of the kinetic energy in the turbulent system. This is a result of the exchange of energy between kinetic and magnetic energies as will be discussed in Section 4.2.2.

From Figure 5 and Table 2, it is evident that the Hall effect has almost no impact on the kinetic energy decay in turbulence in molecular clouds. In order to emphasize any possible impact of the Hall effect we have plotted the time evolution of the ratio of the kinetic energy in each of our simulations to that in mc-5-512 in Figure 6 on a linear scale. It is clear even in this figure that the Hall effect has little impact on the evolution of the turbulence. This result is also reproduced if we examine the evolution of the magnetic energy (see Section 4.2.2). This supports our conclusion from Paper I in which the simulations were run using fixed resistivities (see also fr-512 in this work).

**Figure 6.** Plot of the ratio of the kinetic energy in our 512³ simulations to that in mc-5-512.
Download figure:
Standard image High-resolution image

Also shown in Figure 5 is the energy decay for simulation fr-512. Given the wide variation of the resistivities in both space and time (see Figures 2 and 3) it is somewhat surprising that the energy decay is so similar to that of mc-512. In fact, the volume average of the ambipolar resistivity at t = 0.5 t_c in mc-512 is approximately 60% higher than that in the fr-512 simulation. It would appear that, while ambipolar diffusion enhances energy loss, the expected spatial and temporal variation of it does not have that much influence.

One possible reason for the small divergence between mc-512 and fr-512 would be if the locations in which the resistivity is high are regions in which the magnetic field, B, is varying weakly. To explore this we define a scalar, δB, by

$\begin{equation} \delta B \equiv |{\bf grad}(B_x)| + |{\bf grad}(B_y)| + |{\bf grad}(B_z)|, \end{equation} \tag{ 25 }$

where grad(B_x), for example, is a normalized gradient defined by

$\begin{equation} {\bf grad}(B_x) \equiv \frac{1}{B_0} \left(\begin{array}{@{}c@{}}\delta _x B_x\\ \delta _y B_x\\ \delta _z B_x\end{array} \right), \end{equation} \tag{ 26 }$

where δ means centered differencing in the indicated direction without normalizing by the zone spacing and B₀ is the magnitude of the magnetic field throughout the domain at t = 0. δB is then a dimensionless measure of the variation of B at any point in space or time: if δB is large it means there is a high gradient in one or more of the components of B and therefore resistivity will have an important influence here. Figure 7 contains snapshots of δB and |B| at t = 0.5 t_c. There is rich structure in δB which is not apparent in |B|.

**Figure 7.** Slices showing the z = 0.5 plane at t = 0.5 t_c for both |B| (left panel) and δB (right panel).
Download figure:
Standard image High-resolution image

Figure 8 contains a histogram plot of the two-dimensional probability density function for δB and r_A (the ambipolar resistivity). What is striking about this plot is that there is a notable lack of high δB with corresponding high resistivity. Of course, if we have any system in which there are regions of high and low resistivity we expect that, over time, the regions with high diffusion will have lower variation in B so this, in itself, does not tell us much. However, it does prompt us to look a little more closely at the behavior of r_A.

Comparison of the middle panel of Figure 1 with the top panel of Figure 2 shows that the ambipolar resistivity is higher in regions of low density, as would be expected. Now, in supersonic/super-Alfvénic turbulence kinetic energy is dissipated most strongly at strong shocks. However, shocks propagating into regions of very low density will not dissipate kinetic energy effectively. Since it is precisely these regions in which our resistivities are high we must conclude that, in fact, the regions of enhanced resistivity do not contribute significantly to energy dissipation and hence we would not expect the introduction of spatially varying resistivities to increase the rate of energy decay in our simulations. Further, since the volume average of, for example, the ambipolar resistivity is actually higher than that used in fr-512 we would not expect the spatial variation of this resistivity to reduce the rate of energy decay either.

4.2.2. Magnetic Energy Decay

We now move on to discuss the decay in magnetic energy. Initially, as outlined in Paper I, the magnetic energy increases as the flow compresses and stretches the magnetic field throughout the computational domain. Once this initial increase in the energy has occurred it is gradually lost through two main avenues: magnetic reconnection and transfer of magnetic energy to kinetic energy which can then be dissipated in shocks and other viscous processes.

Figure 9 contains plots of the decay of magnetic energy with time. It is clear that, as is common in the case of kinetic energy, the hall-512 and mhd-512 simulations are almost identical, while the ambi-512 and mc-512 simulations are also well matched. This supports our inference from Section 4.2.1 that the Hall effect has almost no impact on energy decay in molecular clouds on the global scale.

Ambipolar diffusion does, however, have a significant impact on the behavior of both the kinetic and magnetic energies. This was also noted in Paper I and we explain it in the same way, recapitulated briefly here for completeness. Consider a region of the flow undergoing compression. During this compression kinetic energy will be converted into both magnetic and internal energy through the increase in the magnetic pressure and the thermal pressure. Given that this is a turbulent flow we expect that, after some time, this region will begin to expand again. However, during the compression ambipolar diffusion will have diffused away some of the magnetic energy thereby leaving less to be converted back to kinetic energy. In this way, the presence of ambipolar diffusion creates a new path through which energy can be lost from the system and reduces the level of all forms of energy in the system.

4.3. Power Spectra

We now move on to a study of the power spectra obtained from the multifluid MHD turbulence simulations. These spectra are important from the point of view of understanding the types of structures formed by the turbulence and are a more discerning tool for exploring any structural differences caused by multifluid effects. Table 3 contains the exponents of the power spectra assuming a power-law relationship between power and wave number. All the analyses presented in this section are performed on data taken at t = t_c.

Table 3. Values of the Exponent for the Power Spectra of Density, Velocity, and Magnetic Field Measured at t = t_c

Simulation	Density	Velocity	Magnetic Field
mc-512	1.82,^a 4.33^b	1.34	2.14,^a 5.38^b
ambi-512	1.79,^a 4.31^b	1.34	2.15,^a 5.43^b
hall-512	1.25	1.01	1.51
mhd-512	1.27	1.00	1.55
fr-512	2.11,^a 4.21^b	1.28	2.20,^a 5.53^b

Notes. All fits are over the range 5 ⩽ k ⩽ 20 unless otherwise noted. ^aFitted over 4 ⩽ k ⩽ 10. ^bFitted over 10 ⩽ k ⩽ 100.

Download table as: ASCII Typeset image

4.3.1. Density Power Spectra

We turn first to the scales of the structures formed in the density distributions for our various simulations. Figure 10 contains plots of the power spectra of the neutral density (or density, in the case of simulations mhd-512 and fr-512) for all the 512³ resolution simulations.

The power spectra for the simulations including the effects of ambipolar diffusion are approximately broken power laws made up of three distinct power laws: 1 ≲ k ≲ 10, 10 ≲ k ≲ 100, and 100 ≲ k. Below the low k break, the spectrum is dependent on the scale of the computational domain. At high k approaching the grid scale, numerical viscosity will begin to dominate. Similar to the results presented so far we see that there is little difference between simulations mc-512 and ambi-512—in fact it is difficult to distinguish between the two spectra without careful examination of Figure 10. The fr-512 power spectrum is also similar to mc-512 and ambi-512 although it has slightly less power at intermediate values of k with the difference here being at most 10%.

There is almost no detectable difference between simulations hall-512 and mhd-512. Evidently, the Hall effect has a much weaker influence on the density structure in molecular clouds than ambipolar diffusion. Ambipolar diffusion, on the other hand, has a very significant impact with pronounced damping of density structures at scales less than one-tenth of the domain size (corresponding to a physical scale of approximately 0.02 pc). This damping is evident in Figure 1 where the density structures in simulations with ambipolar diffusion are more smeared than in mhd-512.

From a quantitative perspective, the data in Table 3 show that the inclusion of ambipolar diffusion significantly steepens the power spectra, increasing the exponent by more than 0.5 over the cases which do not include the effect. The data again suggest that the Hall effect has minimal impact.

The results for fr-512 show a softer spectrum at large length scales and a harder spectrum at scales shorter than mc-512 and ambi-512, indicating that the self-consistent calculation of the resistivities reduces the level of fine structure. This result must, however, be confirmed by higher resolution simulations before it can be regarded as reliable.

Figure 11 contains plots of the power spectra of the neutral density and the density of the negatively charged species at t = t_c for comparison. It can be seen that the neutral mass density has more power for k ⩾ 10 although the qualitative shape of the power spectra is the same in each case. This is in qualitative agreement with the results presented in Li et al. (2008) for driven turbulence simulations. In comparison with the results in Table 3, the exponent for the charged species mass density in the range 4 ⩽ k ⩽ 10 is −2.21 while in the range 10 ⩽ k ⩽ 100 it is −4.84. The power spectrum for the neutral mass density is harder than that for the magnetic field and, in the range 10 ⩽ k ⩽ 100, the spectrum for the charged species mass density lies somewhere between the two. This is not surprising as it is a function of the neutral density and velocity (through drag) and the magnetic field (through the Lorentz force).

4.3.2. Velocity Power Spectra

Figure 12 contains plots of the velocity power spectrum for the neutral velocity. Similar to Section 4.3.1, those simulations incorporating ambipolar diffusion have significantly less power at almost all scales than those without. This is to be expected given the increased rate of loss of turbulent energy in the presence of ambipolar diffusion (see Section 4.2).

As found in Paper I and again here, there are clear differences of a qualitative nature between the density power spectra and the velocity power spectra for the multifluid simulations. Those simulations incorporating ambipolar diffusion (mc-512 and ambi-512) exhibit a strong power law in the range 10 ⩽ k ⩽ 100 in the density power spectrum. This is not true of the velocity power spectra. For the latter spectra there is a break at roughly k ∼ 20 and again at k ∼ 100. The latter break we interpret as being at the scale where numerical diffusive effects begin to dominate the non-ideal effects in the induction equation. The lower break can reasonably be interpreted as the scale at which the non-ideal effects become important. The marked qualitative differences between the density and velocity power spectra indicate a considerable decoupling between the two variables. It is worth recalling here that the power spectra are calculated at t = t_c so it is reasonable to expect that the turbulence is well developed at this stage.

As with the density power spectra, the presence of ambipolar diffusion produces steeper velocity power spectra (see Table 3) in qualitative agreement with the results of Li et al. (2008).

The results at high k are, of course, dominated by numerical diffusive effects. However, it is interesting to note that simulation fr-512 actually has very slightly less power at short scales than mc-512 despite the volume average of r_A being roughly twice as high as the resistivity used in fr-512. This adds weight to the inference from Section 4.2.1 that regions of high resistivity tend not to be coincident with regions of high gradients in the magnetic field and therefore do not have the level of influence one would naively expect.

Figure 13 contains plots of the power spectra for the velocity of the neutral and negatively charged species. It can be seen that, except at very high k, they are virtually identical. Given that the charged velocity is defined by balance between drag with the neutrals and the Lorentz force this is, perhaps, unsurprising.

4.3.3. Magnetic Field Power Spectra

Figure 14 contains plots of the spherically integrated power spectra for the magnetic field. Once again, ambipolar diffusion has a much bigger impact on these spectra than the Hall effect. The magnetic field power spectra are considerably steeper at all k in its presence. The absolute power at any scale is also considerably lowered by ambipolar diffusion, as would be expected from the discussion in Section 4.2.2. There is a qualitative similarity between the density power spectra (Figure 10) and the magnetic field power spectra which is absent when comparing the latter with the velocity power spectra (see Figure 12).

Interestingly, we see the phenomenon that fr-512 has very slightly less power at k ⩾ 20 than mc-512 indicating that fixing the resistivities at t = 0 actually results in slightly more dissipation than allowing it to vary in time and space.

We find that the Hall effect has a slightly more noticeable effect on the magnetic field power spectra than on the spectra of velocity or density: it gives rise to a little more structure on short scales which is absent in its absence. This is expected as the Hall effect is a dispersive effect acting directly on the magnetic field. The results for the density and velocity power spectra, however, demonstrate that this influence over the magnetic power spectra does not translate into an influence over the other variables.

4.4. Velocity Dispersion

It has been widely reported that the observed velocity dispersion in molecular clouds behaves as a power law in the size of the field of view. Figure 15 contains plots of the velocity dispersion for each of the 512³ resolution simulations presented in this work. As noted in Paper I and Lemaster & Stone (2009), no power law is observed. This may be due to the fact that in (multifluid) MHD turbulence there are many relevant signal speeds due to the variation in the orientation and intensity of the magnetic field, in contrast to the situation with hydrodynamic turbulence. Hence, there is no reason to expect to see a power law. Finally, it should be noted that, at this point in the simulation, the turbulence has decayed such that it is largely subsonic or transonic. This may also explain the lack of a power law. Indeed, it should also be noted that our results here are not necessarily in contradiction with observations since the velocity dispersion–size correlation can only be accurately measured in the supersonic regime.

**Figure 15.** Log-scale plot of the velocity dispersion as a function of length scale for each of the 512³ simulations. There is no power law apparent in these data.
Download figure:
Standard image High-resolution image

It is clear that ambipolar diffusion reduces the velocity dispersion at all length scales—this is what would be expected given the results from the velocity power spectra presented in Section 4.3.2. Once again the results here indicate that the Hall effect has little impact. Again, it is interesting to note that the spatial variation of the resistivity in mc-512 appears to have almost no impact as the velocity dispersion seen in fr-512 is almost identical to that in mc-512.

5. CONCLUSIONS

We have presented results from a suite of 512³ resolution simulations of fully multifluid MHD decaying turbulence. The effects incorporated include the Hall effect and ambipolar diffusion. We have performed a resolution study to ensure that the energy decay rate, being the main result presented here, is reliable. We have confirmed the results of the simplified calculations in Paper I that the Hall effect has little impact on the nature and behavior of turbulence in molecular clouds under the well-motivated physical parameters assumed in this work. Further, the presence of ambipolar diffusion increases the rate of energy decay at length scales of 0.2 pc and less. The same conclusion is drawn for the behavior of the energy in the magnetic field.

The power spectra for these simulations again suggest that the Hall effect has little impact on the flows with the exception of the spectrum of magnetic field variations. We must keep in mind that the maximum resolution used here (512³) is only enough to resolve about half a decade in k-space and it is therefore difficult to be confident of the details of the power spectra. Notwithstanding this consideration, it does appear clear that ambipolar diffusion steepens the power spectra of the neutral velocity, density, and the magnetic field. As noted in Paper I, it appears that at a resolution of 512³ and an assumed length scale of 0.2 pc we have resolved the length at which ambipolar diffusion begins to influence the turbulent cascade. In Paper I, only constant resistivities were implemented and hence it was unclear whether this latter result would survive the inclusion of more realistic fully multifluid MHD in which the resistivities vary strongly in space and time. The results presented here imply that it does.

The power spectra of the neutral velocity and the magnetic field differ qualitatively from that of the density with breaks occurring in the former which are not seen in the latter. This suggests a decoupling between these fields.

The velocity dispersion as a function of length does not behave as a power law. This is not unexpected as the nature of MHD turbulence implies a wide range of applicable signal speeds which can, when combined, remove the power-law behavior which might be expected if only one signal speed were relevant.

The authors acknowledge the SFI/HEA Irish Centre for High-End Computing (ICHEC) for the provision of computational facilities and support. The work described in this paper was carried out using resources provided to ICHEC through the Irish National Capability Computing Initiative, a partnership between all the major third level research institutions and IBM coordinated by the Dublin Institute for Advanced Studies and supported by the HEA under PRTLI cycles 3 and 4 with funding from the ERDF and the NDP. This work was partly funded by SFI through the Research Frontiers Programme, grant No. 07 PHYF586. Finally, the authors thank the anonymous referee for helpful suggestions and comments on this work.

MULTIFLUID MAGNETOHYDRODYNAMIC TURBULENT DECAY

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION