NONIDEAL MAGNETOHYDRODYNAMIC TURBULENT DECAY IN MOLECULAR CLOUDS

T. P. Downes; S. O'Sullivan

doi:10.1088/0004-637X/701/2/1258

1. INTRODUCTION

The role of turbulence in molecular cloud evolution has been a subject of much study in the literature (see, for example, the excellent reviews of Mac Low & Klessen 2004; Elmegreen & Scalo 2004). Observations of the properties of gas and dust motions in molecular clouds (Larson 1981) suggest that, indeed, turbulence is present. It is clear that turbulent motion could influence the star formation rate and efficiency as well as the initial mass function (Elmegreen 1993; Klein et al. 2003). Given all this there are several interesting questions which arise:

1.
What is the source of molecular cloud turbulence?
2.
How fast does it decay?
3.
How does it affect star formation?
4.
How does it affect the evolution of molecular clouds?

The first two questions are clearly inter-related: if the turbulence decays very quick, then we need a lot of energy from its source in order to maintain it. Indeed, to address the final two questions we must first gain insight into the first two.

In order to study turbulence in molecular clouds we must resort to numerical simulations. Quite simply, there are no satisfactory analytic techniques for addressing compressible magnetohydrodynamic (MHD) turbulence (e.g., Elmegreen & Scalo 2004).

Many authors have performed sophisticated numerical simulations in order to investigate both the qualitative nature of MHD turbulence and its decay (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). Most of this work has been carried out for the case where ideal MHD is valid i.e., on relatively large length scales. When smaller length scales are considered (e.g., lengths of significantly less than a parsec), ambipolar diffusion becomes non-negligible in molecular clouds (Oishi & Mac Low 2006). Some authors (Li et al. 2008; Kudoh & Basu 2008; Oishi & Mac Low 2006) 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.

It has been suggested (Wardle 2004) that although the Hall resistivity is generally at least an order of magnitude lower than the ambipolar resistivity in molecular clouds, its effect should not be ignored due to the qualitative change it induces in the behavior of the magnetic field. 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 (Servidio et al. 2007; Mininni et al. 2006; Matthaeus et al. 2003). 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 low resolution simulations by Downes & O'Sullivan 2008). In particular, to our knowledge, the work presented here represents the first systematic study of molecular cloud turbulence incorporating both the Hall effect and ambipolar diffusion simultaneously.

The main aim of this work is to examine in detail the differences between the decay of ideal MHD turbulence and that of more realistic nonideal MHD turbulence with a full tensor resistivity, incorporating the effects of ambipolar diffusion, the Hall effect, and Ohmic resistivity. This work is new in two respects: no previous work has focused on decaying (i.e., un-driven) turbulence in the presence of nonideal terms and, in addition, no previous work has addressed the issue of turbulence in the presence of both ambipolar diffusion and the Hall effect simultaneously. This is the first of a short series of papers describing a comprehensive study of multifluid MHD turbulence in a parameter regime appropriate to molecular clouds.

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

In this work, we use the code HYDRA (O'Sullivan & Downes 2006, 2007) to integrate the equations of nonideal 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. For this work, we assume spatially and temporally constant Ohmic, ambipolar, and Hall resistivities (see Section 2.1).

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 both the full abilities of the HYDRA code and the assumptions underlying the equations used.

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). The equations solved in this work are then

$\begin{eqnarray} \frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho \mbox{\boldmath $q$}\right) & = & 0, \end{eqnarray} \tag{ 1 }$

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

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

$\begin{eqnarray} \nabla \cdot \mbox{\boldmath $B$} & = & 0, \end{eqnarray} \tag{ 4 }$

$\begin{eqnarray} \nabla \times \mbox{\boldmath $B$} & = & \mbox{\boldmath $J$}, \end{eqnarray} \tag{ 5 }$

where ρ, q, a, $\mbox{\boldmath $\mathsf {I}$}$ , B, and J are the neutral mass density, neutral velocity, sound speed, identity matrix, magnetic field, and current density, respectively. The electric field in the frame of the fluid, E', 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{ 6 }$

where

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

$\begin{eqnarray} \mbox{\boldmath $E$}_{\rm H} & = & \mbox{\boldmath $J$}\times \mbox{\boldmath $a$}_{\rm H}, \end{eqnarray} \tag{ 8 }$

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

using the definitions a_O ≡ f_OB, a_H ≡ f_HB, a_A ≡ f_AB, where $f_{\rm O}\equiv \sqrt{r_{\rm O}}/B$ , f_H ≡ r_H/B, $f_{\rm A}\equiv \sqrt{r_{\rm A}}/B$ . r_O, r_H, and r_A are the Ohmic, Hall, and ambipolar resistivities, respectively.

In this work, these resistivities are kept constant in both space and time. We note that, physically, they depend on both the magnetic field and the density of the various charged species in the fluid (e.g., Ciolek & Roberge 2002; Falle 2003; O'Sullivan & Downes 2006, 2007) and hence, in reality, do vary in both space and time. Treating the resistivities in such a simple way allows us to gain a deeper understanding of their influence on turbulence without having to consider the complicating effects of dynamically varying resistivities at the same time. As such, and as a first step away from the approximation of ideal MHD and a single form of magnetic diffusion, we believe this is an interesting study. Having gained some insight into this simplified model a follow-up paper will address multifluid MHD turbulence under the influence of self-consistently calculated resistivities.

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. Equation (4) is enforced using the method of Dedner et al. (2002). The effects of the diffusive terms are then incorporated in an operator split fashion.

2.2. Initial Conditions

We examine the decay of supersonic MHD turbulence in conditions suitable for dense regions of molecular clouds. While the simulations presented here are of course scale-free, we present the initial conditions used in standard astrophysical units for ease in reading.

The simulations are carried out in a cube of side L = 0.2 pc with periodic boundary conditions being enforced on all faces. The sound speed is set to 0.55 km s⁻¹ and the initial density is chosen to be uniform with a value of 10⁶ cm⁻³. 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⁻¹ (see Figure 1, Wardle & Ng 1999). 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 can use our simulations to find whether the Hall effect is ever likely to be important in molecular cloud turbulence.

The initial velocity field is used to instigate the turbulence in these simulations. Each component of the velocity field is defined to be the sum of waves with 16 wave-vectors, each with random amplitude and phase, i.e.,

$\begin{equation} q_i(x,y,z) = \sum _{l,m,n=1}^{4} A^{lmn}_i \cos \big(k^l_i x + k^m_i y + k^n_i z + \phi ^{lmn}_i\big), \end{equation} \tag{ 10 }$

where i = 0, 1, 2 defines the component (x, y or z respectively) of the velocity, A^lmn_i and ϕ^lmn_i are the random amplitudes and phases and

$\begin{equation} k^l_i \equiv \frac{2 \pi l }{L} (1-\delta _{li}), \end{equation} \tag{ 11 }$

where δ_li is the usual Dirac delta function. The inclusion of the term in parenthesis in the definition of k^l_i restricts the initial velocity field to be solenoidal. Note that we perform all the analysis of these simulations at t ⩾ 0.2 t_c (1 flow crossing time) at which time the effects of the precise initial conditions used should be negligible.

The nomenclature for the simulations is xx-ab-c where xx denotes the type of physics (e.g., a standard molecular cloud run is "mc," ideal MHD is "mhd" etc.,), ab is the initial rms Mach number, and c is the resolution used. The initial root-mean-square (rms) of the field is chosen to be either Mach 2.5, 5, or 10 depending on the simulation in question. These correspond to Alfvénic Mach numbers of approximately 0.96, 1.9, and 3.85, respectively. In addition to the nonideal MHD simulations described we also run four further simulations. The first is an ideal MHD simulation (mhd-5-512) which we use for comparison purposes, another is a pure hydrodynamic simulation (hd-256-0.5), and the other two (ambi-5-512 and hall-5-512) only incorporate one of ambipolar diffusion or the Hall effect, respectively. We use these latter simulations to separate out the effects of each of these diffusions to better understand the physics occurring. Table 1 contains definitions of the various simulations used in this work.

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

Simulation	Mach Number^a	Resolution	Comment
mc-5-64	5	64³	⋅⋅⋅
mc-5-128	5	128³	⋅⋅⋅
mc-5-256	5	256³	⋅⋅⋅
mc-5-512	5	512³	⋅⋅⋅
mc-2.5-512	2.5	512³	⋅⋅⋅
mc-10-512	10	512³	⋅⋅⋅
ambi-5-512	5	512³	r_H = 0
hall-5-512	5	512³	r_A = 0
mhd-5-512	5	512³	Ideal MHD
hd-5-256	5	256³	Hydrodynamic

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 main aim of this paper is to investigate the decay rate of supersonic turbulence in molecular clouds. Hence, the main analysis carried out on the simulation results is the calculation of the volume-averaged kinetic, magnetic, and total energy in the simulation as a function of time. These quantities are defined as

$\begin{eqnarray} e_{\rm k} & = & \langle \rho |\mbox{\boldmath $q$}|^2 \rangle, \end{eqnarray} \tag{ 12a }$

$\begin{eqnarray} e_{\rm b} & = & \left\langle \frac{|\mbox{\boldmath $B$}|^2}{2}\right\rangle - \frac{\langle \mbox{\boldmath $B$}\rangle ^2}{2}, \end{eqnarray} \tag{ 12b }$

$\begin{eqnarray} e_{\rm tot} & = & e_b + e_k, \end{eqnarray} \tag{ 12c }$

where the angle brackets denote averaging over the computational domain. Note that e_b is therefore the difference in the magnetic energy at the current time and the magnetic energy at t = 0 (see, for example, Vestuto et al. 2003; Lemaster & Stone 2009), given the assumption that no external electromotive force is applied and that our boundary conditions are periodic. 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{ 13 }$

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

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

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 spectrum in the x-, y-, and z-directions and then integrating the power for all k satisfying k ⩽ |k| < k + dk for each k with dk = 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{ 15 }$

where

$\begin{equation} \sigma _\alpha (l) = \left\lbrace \big\langle q_\alpha ^2\big\rangle _l-\langle q_\alpha \rangle _l^2\right\rbrace ^{\frac{1}{2}}, \end{equation} \tag{ 16 }$

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 nonoverlapping cubes within the domain.

4. RESULTS

We now present the results of the simulations carried out. Each simulation was run for one sound crossing time, t_c = 3.56 × 10⁴ yr, of the simulation domain.

Figure 1 shows the density distributions in a slice at x = 0.1 pc (i.e., the midplane) and times t = 0.2, 0.5 and 1 t_c for the mc-5-512 and the mhd-5-512 simulations. Both the nonideal and ideal MHD simulations show some anisotropy with respect to the projected initial magnetic field direction, i.e., the (1, 1) direction with filaments both perpendicular and parallel to this direction. This is what would be expected as the material will flow preferentially along the (1, 1, 1) direction and hence we expect the shock fronts to be normal to it. This anisotropy becomes more evident as time progresses because the kinetic energy decays to below the energy of the mean magnetic field (which is conserved with time), hence the flow evolves from being dominated by kinetic energy to dominated by magnetic energy over the lifetime of the simulations. Figure 2 contains a plot of the same density distribution for hydro-5-256 at t = 0.2 t_c for comparison. It can be seen that this is qualitatively different to both the ideal and nonideal simulations shown in Figure 1 and does not display any signs of anisotropy.

**Figure 2.** log-scale plot of a slice in the mass density at time t = 0.2 t_c for the hydrodynamic simulation (hd-5-256).
Download figure:
Standard image High-resolution image

There are obvious qualitative differences between the mc-5-512 and mhd-5-512 simulations at all times with the density distribution in the ideal case containing much more small-scale structure than the nonideal case. This results from the ability of the magnetic field to move with respect to the fluid. This means that it is not compressed as much as in the ideal case, and hence the same level of small structure is not formed. Since the coupling between the magnetic field and the flow field is strong, even though it is imperfect, this lack of small-scale structure in the magnetic field becomes mirrored in the density field. The implications of this are discussed in more detail in Section 4.2. Clearly, in the nonideal simulation, energy is not as efficiently transported from large scales to smaller scales as in the ideal case.

We now continue our discussion of the results by considering a resolution study to demonstrate convergence of our numerical results. We then go on to discuss the properties of the density, velocity, and magnetic field distribution resulting from the turbulence, and finally the properties of the energy decay.

4.1. Resolution Study

In order to be reasonably confident that the results we present in this work are converged we have performed a resolution study for the conditions used in the simulations. We have run simulations which are identical in all respects except for the resolution, which varies from 64³ up to 512³ (see Table 1). We then analyze these in an identical fashion and consider the differences between our results for the different resolutions.

Figure 3 shows plots of the kinetic energy normalized to its initial value as a function of time for each of the simulations in our resolution study. The 64³ and 128³ plots do appear to be significantly different from the higher resolution simulations. However, the 256³ and 512³ simulations are much more similar with a maximum difference in the total kinetic energy in the simulation at any one time being less than 10%. Table 2 contains least squares fits of the decay over the time interval [0.2t_c, t_c] assuming it to be of the form t^−β for each of the kinetic energy (β_K), the magnetic energy (β_B), and their total (β_Tot). It can be seen that β_K varies by approximately 4% over the entire range of the resolution study (simulations mc-5-64 through mc-5-512). We feel, therefore, that we can be fairly confident of the value of this exponent when comparing it with the other simulations presented here.

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

Simulation	β_K	β_B	β_Tot
mc-5-64	1.34	1.40	1.35
mc-5-128	1.33	1.35	1.34
mc-5-256	1.37	1.36	1.36
mc-5-512	1.40	1.37	1.39
ambi-5-512	1.40	1.35	1.38
hall-5-512	1.25	1.18	1.22
mhd-5-512	1.26	1.19	1.23
hydro-5-256	1.10	⋅⋅⋅	1.10
mc-2.5-512	1.21	1.29	1.23
mc-10-512	1.42	1.39	1.41

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

Download table as: ASCII Typeset image

The decay of the energy in the magnetic field is more rapid than that in the kinetic energy with a difference between β_K and β_B of about 4%. At low resolution (mc-5-64) β_B is at its highest which is a sign that the low resolution is introducing sufficient numerical viscosity to induce large amounts of numerical reconnection. As the resolution is increased, β_B reduces to approximately 1.36 and stays at about this value even up to the maximum resolution of 512³. We can be reasonably confident then that resolution is not affecting our estimate of the decay rate of the magnetic energy at our maximum resolution.

The decay of the total energy, being derived from the decay of kinetic and magnetic energy, is also reasonably well converged with a total change over the entire range of the resolution study of around 4%.

Note that, while we can be confident from our results that the decay rate is converged, we have to be more careful when considering our power-spectra results presented in Sections 4.3 and 4.5.2. For example, Kritsuk et al. (2007), Schmidt et al. (2009) and Lemaster & Stone (2009) present power spectra in ideal MHD which suggest that at resolutions of 1024³ the turbulent inertial range is established over at most a decade in k, while for resolutions of 512³ this falls to around half a decade.

4.2. Energy Decay

We now discuss the behavior of the kinetic and magnetic energy in our nonideal simulations and compare them with those for our ideal simulation.

4.2.1. Kinetic Energy Decay

Figure 4 contains plots of the kinetic energy as a function of time for the mc-5-512, ambi-5-512, hall-5-512, and mhd-5-512 simulations. It is apparent that the behavior of mhd-5-512 is rather similar to hall-5-512 and the behavior of ambi-5-512 is similar to that of mc-5-512. This indicates that, at least for the kinetic energy decay, the Hall effect has little impact. However, there is a marked difference between those simulations incorporating ambipolar diffusion (mc-5-512 and ambi-5-512) and those which do not (hall-5-512 and mhd-5-512). Ambipolar diffusion clearly increases the decay rate of the turbulence. This can also be seen from the data presented in Table 2 where the exponents of the simulations containing ambipolar diffusion are greater by about 10% than those without.

In general, we expect that turbulence will decay more rapidly in systems which have higher viscosity than those without. Clearly the ambipolar diffusion, although it corresponds to a "viscosity" in the magnetic field, acts in a similar way to the usual viscous forces in a fluid when considering this system. This can be seen in Figure 1, where the density features are more smeared out in the nonideal case, as would be expected if we simply introduced a viscous term into the momentum equations.

It is important to note here that the nonideal diffusive terms in the induction equation will lead to enhanced reconnection in the magnetic field. Since we assume an isothermal equation of state there is no path by which the energy released by reconnection can find its way to the kinetic energy of the system. There may be significant differences between the effects of the nonideal terms on the decay rate of turbulence in the isothermal and nonisothermal regimes in molecular clouds.

4.2.2. Magnetic Energy Growth and Decay

Figure 5 contains plots of the energy of the magnetic perturbations induced by the turbulence (normalized to the initial kinetic energy) as a function of time for the mc-5-512, ambi-5-512, hall-5-512, and mhd-5-512 simulations. A similar trend is evident in the results for magnetic energy decay as that already noted in Section 4.2.1. Again, the main result is that ambipolar diffusion increases the energy decay rate while the Hall effect does little to influence it. This result is borne out by the data in Table 2 in which the decay exponents are greater for the simulations incorporating ambipolar diffusion by about 14% than those which do not.

To gain a little more insight into the interplay between the magnetic and kinetic energy in the system we now focus on simulations mc-5-512 and mhd-5-512 at early times. Figure 6 contains plots of the magnetic and kinetic energies of these two simulations as a function of time. The magnetic energy grows initially as the flow converges and compresses the magnetic field in regions throughout the computational domain. The kinetic energy gradually decays during this time. Once shocks form, equipartition between the total magnetic energy and kinetic energy is reached and the magnetic field begins to decay. At this time the decay of the kinetic energy accelerates due to the dissipative effect of the shocks which have just formed. This dependence on shocks forming is confirmed when the growth and decay of the magnetic energy is examined for simulations mc-2.5-512 and mc-10-512 and compared with mc-5-512 (see Figure 7). Here, we can see that the peak in magnetic energy is reached at later times as the initial rms Mach number is decreased. The time at which shocks begin to form is determined by the initial rms Mach number (giving us a characteristic velocity) and half the shortest wavelength in the initial velocity field. The approximate shock formation times are 0.0125 t_c, 0.025 t_c, and 0.05 t_c for the mc-10-512, mc-5-512, and mc-2.5-512 simulations, respectively. These times match up reasonably well, particularly given the temporal granularity of the simulation data, with the peak of the fluctuating part of the magnetic energy in these simulations.

**Figure 6.** log-scale plot of the magnetic and kinetic energies (normalized to the initial kinetic energy) for simulations mc-5-512 and mhd-5-512 for a larger range of times. See text.
Download figure:
Standard image High-resolution image

**Figure 7.** Plot of the growth and decay of magnetic energy (normalized to the initial kinetic energy) for simulations mc-10-512, mc-5-512, and mc-2.5-512.
Download figure:
Standard image High-resolution image

Just after equipartition is reached the decay of the kinetic energies of mc-5-512 and mhd-5-512 begin to behave differently. This is unsurprising as it is only when the magnetic field perturbations are reasonably strong that nonideal effects can play a dynamically significant role. Just after the magnetic energy peaks the kinetic energy of simulation mhd-5-512 goes through a short period during which it does not decay particularly rapidly—presumably a result of the transfer of energy from the compressed magnetic field back to the kinetic energy as the fluid expands after the initial compressions in the simulation. This does not happen in simulation mc-5-512 since the magnetic field does not become so compressed in the first place due to nonideal effects. This can be seen from the fact that the magnetic energy peaks at a lower level in this simulation. In addition, energy from this field will not be transferred efficiently back to the kinetic energy because of the imperfect coupling between the magnetic and velocity fields resulting from the nonideal terms in the induction equation.

The picture which emerges is the following: initially kinetic energy is transferred into magnetic energy through compressions. If there are diffusive terms present in the induction equation then this extra magnetic energy will be dissipated and there will be less energy available to transfer back to the kinetic energy as the compressed regions expand. This process will be repeated throughout the simulation as shocks form and dissipate. Hence, we expect a more rapid decay of both magnetic and kinetic energy when diffusive terms are present. Recalling that the Hall effect does not actually diffuse magnetic energy, it is no surprise that simulation hall-5-512 is so similar in terms of energy decay to mhd-5-512.

Finally, we note that the peak of the fluctuating part of the magnetic energy occurs at a slightly later time in the mhd-5-512 simulation than the mc-5-512 one. We explain this as follows. As noted in Section 4.3.1, the ambipolar diffusion appears to set a dissipation length scale corresponding to k = 10. This is significantly larger than the dissipation scale in the ideal mhd simulation. Hence, energy in the magnetic field need not cascade to such short length scales in order to be dissipated in the mc-5-512 and ambi-5-512 simulations as in the mhd-5-512 and hall-5-512 simulations. Since the energy cascade to shorter length scales takes time (particularly at very early times when shocks are only beginning to form), the energy in the magnetic field can be removed from the former simulations at earlier times than from the latter simulations as it need not cascade so far.

4.3. Power Spectra

We now turn to power spectra of the simulations. These tell us about the scale of the structures formed and are useful in gaining some insight into the effects of the nonideal physics incorporated. It is important to stress that since these simulations are of decaying turbulence, rather than driven turbulence, we do not necessarily expect to get the often-quoted power-law dependence of power on wavenumber. All the power spectra presented have been calculated at t = t_c when we can be fairly confident that the turbulence is fully developed and the initial conditions have been effectively forgotten. In addition, the concerns noted in Section 4.1 should be borne in mind, i.e., the turbulent inertial range is likely to be of order half a decade in k for our 512³ simulations.

We discuss the power spectra of density, magnetic field strength, and velocity in turn.

4.3.1. Density Power Spectra

Figure 8 contains plots of the spherically integrated power spectrum of the density for the mc-5-512, hall-5-512, ambi-5-512, and mhd-5-512 simulations at t = t_c. It is clear that there is little difference between the mc-5-512 and ambi-5-512 simulations—again strongly indicating that the Hall effect has little influence on the behavior of the simulations. The mhd-5-512 and hall-5-512 simulations are rather similar, although there is less power in small-scale structures in the hall-5-512 simulation than the mhd-5-512 one. We leave discussion of this until Section 4.3.3.

Even at large length scales there is a significant difference in the power spectra, indicating that dissipation on short length scales due to the nonideal terms in the induction equation indeed affects the behavior of the density distribution at large scales. The simulations with ambipolar diffusion have considerably steeper slopes than those without (see Table 3) being in qualitative agreement with the results of Li et al. (2008). That the mc-5-512 simulation has less power at shorter length scales than mhd-5-512 can be seen by a cursory examination of Figure 1—the mc-5-512 images appear more "blurry" than those from mhd-5-512—so this result is not a surprise.

Table 3. The 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-5-512	2.09^a, 4.06^b	1.47	2.17^a, 4.96^b
ambi-5-512	2.04^a, 4.03^b	1.49	2.21^a, 4.75^b
hall-5-512	1.41	1.20	1.65
mhd-5-512	1.45	1.17	1.59
mc-2.5-512	2.37^a, 3.32^b	2.14	2.86^a, 3.80^b
mc-10-512	2.20^a, 4.20^b	1.44	2.59^a, 4.70^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

Further structure in the power spectrum for mc-5-512 is apparent which is absent in the mhd-5-512 simulation: at k ≈ 10 there is a break to a steeper slope in the simulations containing ambipolar diffusion. For the mhd-5-512 simulation, at about k ≈ 35 numerical viscosity begins to affect the power spectrum, as evidenced by the roll over of the spectrum from a power law to a steeper slope. The same effect does not appear in the mc-5-512 simulation because the physical viscosity due to the nonideal terms in the induction equation dominates the numerical viscosity up to much higher values of k. While the power spectrum for the latter simulation does change slope at k ≈ 10 it then maintains a strong power law up to k ≈ 100.

Hence, the break at k ≈ 10 for mc-5-512 and ambi-5-512 appears to be a physical result, rather than a numerical one. This suggests that, apparently in contradiction with a result of Oishi & Mac Low (2006), ambipolar diffusion can set a length scale in turbulence—in these simulations that length scale is about 0.02 pc. We discuss this further in Section 4.3.3 but note that we should be cautious about drawing general inferences given that the resistivities used in these simulations are constant in space and time.

Ultimately, at very short length scales (k ⩾ 100 or l ⩽ 0.002 pc), all four simulations steepen significantly. This is expected since at these values of k, which correspond to lengths of less than about 5 grid zones in these simulations, numerical viscosity will certainly dominate structure generation/dissipation.

4.3.2. Velocity Power Spectra

We now turn to the velocity power spectra. Figure 9, contains plots of the spherically integrated power spectrum of the velocity for the mc-5-512, hall-5-512, ambi-5-512, and mhd-5-512 simulations at t = t_c. In broad terms of the differences between the simulations we obtain the same result as for the density power spectrum (compare Figure 8). It is, however, obvious that the density and velocity power spectra differ quite substantially in qualitative terms.

At very low k, the velocity spectra have a rather shallow slope with little difference between any of the simulations up to k ≈ 4—the maximum value of k at which the initial velocity field was nonzero. In the range 5 ⩽ k ⩽ 20, both mhd-5-512 and hall-5-512 follow a similar power law, while mc-5-512 and ambi-5-512 attain a steeper slope in this range. In general, for k ⩾ 5 the simulations with ambipolar diffusion have significantly steeper power spectra, as observed in the density power spectra and in qualitative agreement with the driven turbulence results of Li et al. (2008).

There is a further apparent break in the power spectrum at around k ≈ 100 in the simulations with ambipolar diffusion occurring where numerical viscosity begins to dominate the ambipolar diffusion (see Section 4.3.1). It is interesting to note that the series of breaks in the velocity power spectra in mc-5-512 and ambi-5-512 are not mirrored in the density power spectra indicating a level of decoupling between these two fields.

Kritsuk et al. (2007) noted that under certain circumstances a power spectrum of ρ^1/3q could follow a power law with the classical Kolmogorov slope of −5/3 even for compressible (hydrodynamic) turbulence. Figure 10, contains power spectra of this variable for the mc-5-512 and mhd-5-512 simulations. It is apparent that the power spectra for this variable for these simulations is shallower than the −5/3 law, at least up to the value of k where dissipative effects may be important. In the range 4 ⩽ k ⩽ 10, the exponents are −1.289 and −1.14 for the mc-5-512 and mhd-5-512 simulations, respectively. The theory under which these variables could be expected to have the Kolmogorov slope has the assumption that the system has reached a statistical steady state—since we investigate decaying turbulence here this is unlikely to be the case. It might be expected that if we do not continually supply energy at low k then the slope after a turbulent crossing time when the turbulence is well established and has also decayed significantly would be too shallow as the power spectra would be "too low" at low k.

**Figure 10.** log-scale plot of the normalized spherically integrated power spectrum of ρ^1/3q for the mc-5-512 and mhd-5-512 simulations at t = t_c. Also shown is the k^−5/3 Kolmogorov power-law.
Download figure:
Standard image High-resolution image

4.3.3. Magnetic Field Power Spectra

Figure 11 contains plots of the power spectrum of the magnetic field at time t = t_c. Once again these plots indicate that ambipolar diffusion is the dominant nonideal effect in molecular cloud turbulence. We can see that there are some differences between each of the simulations. The mc-5-512 and ambi-5-512 simulations are rather similar to the density power spectra (Figure 8). However, at high k the ambi-5-512 simulation does have noticeably more power than mc-5-512. This indicates that the Hall effect, present in the mc-5-512 simulation, is having some impact in the structure of the magnetic field at short length scales. This effect is very small in the power spectra of both velocity and density and is not apparent in the energy decay rates either. It would appear that when the ambipolar diffusion is strong in comparison to the Hall effect, as is the case for molecular clouds, the Hall effect may change the structure of the magnetic field but this change will not propagate into the rest of the fluid variables.

The power spectra for mhd-5-512 and hall-5-512 are rather similar to each other, again with the exception of high k where the hall-5-512 simulation has significantly less power. We explain this as follows. While the Hall effect does not actually diffuse the magnetic field it does reorient it. This process of reorientation can give rise to a topology favoring some magnetic reconnection—particularly on the small scales at which the Hall effect operates. We therefore attribute the lower power at high k in hall-5-512 to the interplay between the Hall effect and numerical viscosity at these length scales. Recall that since hall-5-512 contains no ambipolar diffusion and negligible parallel resistivity any reconnection which occurs must be almost entirely due to numerical viscosity. The same argument holds for the difference between mc-5-512 and ambi-5-512 at high k, but in this case the interplay is between the Hall effect and ambipolar diffusion. The marked turn up in power for k ⩾ 200 for the mhd-5-512 simulation is due to numerical effects.

Overall the plots are rather similar to the density power spectra (Figure 8): the density and magnetic field distributions appear to be fairly well coupled and different in nature to the velocity distribution.

4.4. Velocity Dispersion

It has been generally accepted that observations of line-of-sight velocity dispersion in molecular clouds exhibits a power law with the size of the field of view (e.g., Larson 1981). It is worth noting, however, that recent observational results call the so-called Larson's law into question (Heyer et al. 2009). However, for completeness we feel it is worthwhile to examine how our nonideal simulations behave in this regard. In this section, we examine the velocity dispersion as a function of characteristic length for each of the simulations.

Figure 12 contains plots of the velocity dispersions at t = t_c for each of the 512³ simulations. Once again, the results follow the general trend of indicating that the Hall effect has almost no impact with the data for mhd-5-512 and hall-5-512 being almost indistinguishable. At all length scales the velocity dispersion for the simulations containing ambipolar diffusion (mc-5-512 and ambi-5-512) is lower, with the difference being larger at short length scales.

From the power-spectra study in Section 4.3 we already know that the power in short-scale variations of the velocity is decreased by the presence of ambipolar diffusion. It is therefore no surprise that comparison of the results for the velocity dispersion shows that the dispersion is decreased by the presence of this diffusion.

There is no obvious power law in these results, in agreement with the simulations of driven, ideal MHD turbulence of Lemaster & Stone (2009). As noted by these authors, the varying strength and direction of the magnetic field give rise to a large number of signal speeds within the domain, rather than the simple sound speed frequently assumed in hydrodynamic turbulence and hence there is no reason to expect a power law in the velocity dispersion.

4.5. The Effect of Initial rms Mach Number

Finally, we turn to the issue of the initial Mach number chosen in the simulations. To study this we consider simulations mc-2.5-512, mc-5-512, and mc-10-512. These correspond to initial sonic Mach numbers of 2.5, 5, and 10, respectively and Alfvénic Mach numbers of 0.96, 1.9, and 3.85.

4.5.1. Energy Decay

Figure 13 contains plots of the decay of kinetic energy for each simulation. The energy in the plots has been normalized to the starting energy for each simulation. It is clear that the higher Mach number flows lose their energy more rapidly than the lower Mach number ones. This behavior is repeated in the magnetic energy and the total energy. Table 2 contains the power-law indices of decay for the kinetic, magnetic, and total energy in these simulations. The exponents for the kinetic energy decay for the mc-2.5-512, mc-5-512, and mc-10-512 simulations are 1.21, 1.40, and 1.42 respectively, while those for magnetic energy are 1.29, 1.37, and 1.39, respectively. It is clear that, indeed, all the indices increase with initial rms Mach number. We explain this by noting that high Mach number flows tend to have strong shocks which will dissipate energy more effectively than lower Mach number flows.

4.5.2. Power Spectra

Figures 14–16 contain the power spectra for the mc-10-512, mc-5-512, and mc-2.5-512 simulations for each of the density, velocity, and magnetic field, respectively taken at t = t_c.

**Figure 14.** log-scale plot of the power spectrum of the density for the mc-10-512, mc-5-512, and mc-2.5-512 simulations.
Download figure:
Standard image High-resolution image

**Figure 15.** log-scale plot of the power spectrum of the velocity for the mc-10-512, mc-5-512, and mc-2.5-512 simulations.
Download figure:
Standard image High-resolution image

**Figure 16.** log-scale plot of the power spectrum of the magnetic field for the mc-10-512, mc-5-512, and mc-2.5-512 simulations.
Download figure:
Standard image High-resolution image

The density power spectrum does not appear to change much between the mc-10-512 and mc-5-512 simulations, although it is markedly shallower for the mc-2.5-512 simulation. It is worth noting that the mc-2.5-512 simulation is initially sub-Alfvénic and so is qualitatively different in nature to the other two. In addition, the rms sonic Mach number drops below 1 at t ≈ 0.11 t_c and hence we have rather well evolved subsonic, decaying turbulence in this case.

The velocity power spectrum for the mc-10-512 and mc-5-512 simulations are also remarkably similar, while the mc-2.5-512 simulation is steeper at low k and shallower for 20 ⩽ k ⩽ 100. Once again, we attribute this difference in behavior to the fact that the initial rms velocity of the mc-2.5-512 simulation is sub-Alfvénic.

The magnetic power spectrum follows a similar pattern in terms of differences between the mc-10-512, mc-5-512, and mc-2.5-512 simulations as for the other two sets of power spectra.

4.5.3. Velocity Dispersion

Figure 17, shows plots of the velocity dispersion at t = t_c for the mc-10-512, mc-5-512, and mc-2.5-512 simulations. Again, as noted in Section 4.4, no overall power law is observed. At short length scales (less than about 0.01 pc), the slope of the relations are all approximately the same. Above this scale the slope is somewhat lower for mc-10-512 and mc-5-512 than for mc-2.5-512. This is to be expected since large velocity variations will be preferentially suppressed by strong shocks during the early evolution of the system, leading to each of these simulations having more similar velocity dispersion at large length scales. Hence, we expect higher Mach number simulations to retain somewhat higher velocity dispersions at all length scales than their lower Mach number counterparts, but that the fractional differences in these dispersions will be lower at large length scales than at shorter ones.

5. CONCLUSIONS

We have presented the first study of turbulent decay in the presence of nonideal terms in the induction equation and also the first simulations to incorporate both the Hall effect and ambipolar diffusion simultaneously in simulations of molecular cloud turbulence. This is the first stage in a comprehensive study of nonideal MHD turbulence. The nonideal effects are therefore included in a simplistic way in this work with a view to developing a complete, intuitive understanding of their impact on turbulence as we continue the study and gradually add in more realistic components of these effects. The results of this further stage of the study are the subject of a forthcoming paper currently in preparation.

We have used a resolution study to determine that our simulations are well resolved for the purposes of the results which we present. In particular, the turbulent energy decay is well resolved from a resolution of about 256³ and the results we present here are derived from 512³ simulations.

We have found that the Hall effect has little influence on kinetic energy decay when present at the levels found in molecular clouds. However, ambipolar diffusion increases the rate of energy decay on length scales of 0.2 pc or less. Similar results are found for the behavior of magnetic energy. Ambipolar diffusion increases the rate of energy decay through diffusing away energy while it is stored in the magnetic field. Viewed in this way, it is unsurprising that the Hall effect does not influence kinetic energy decay strongly since it does not actually diffuse magnetic energy. The kinetic and magnetic energy decay is faster in higher Mach number flows as would be expected both from this argument and from the fact that we are simulating isothermal flows which will lose energy more efficiently when stronger shocks are present.

The power spectra of the density for these simulations again suggests that the Hall effect has little impact on the flows, while ambipolar diffusion cannot be ignored. As might be expected from a diffusion term, the power spectrum is softer (i.e., steeper) when it is present. At a resolution of 512³ and an assumed length scale of 0.2 pc we appear to have resolved the length at which ambipolar diffusion begins to influence the flow. Ambipolar diffusion does appear to impose a characteristic length-scale on the turbulence. However, we must be cautious about this interpretation since our assumption of spatially constant resistivities will have some impact on this result.

When comparing the velocity power spectra with the density power spectra we find that there appears to be a decoupling between the two fields with breaks in the power laws not mirrored between the two sets of spectra when nonideal effects are included.

The Hall effect does have some impact on the magnetic power spectra at high k. It decreases the energy at high k, probably due to reorientation of the magnetic field at small scales in such a way as to favor reconnection (either numerical or physical) and hence destruction of structure on these scales.

Calculations of the velocity dispersion as a function of length scale show that, again, ambipolar diffusion is the dominant diffusive term and that it has a significant impact. It preferentially reduces the velocity dispersion at small scales. We do not find a power law dependence between length scale and velocity dispersion, in agreement with Lemaster & Stone (2009).

Comparisons of decaying turbulence with varying initial rms Mach number show that higher Mach number flows decay more quickly than their low Mach number counterparts. There are also differences in the power spectra with the Mach 2.5 flow being significantly different to the Mach 5 and Mach 10 flows. The Mach 2.5 flow is slightly sub-Alfvénic initially and this may explain the qualitative difference seen. The velocity dispersions are lower at short scales for lower Mach number flows, but are similar at large scales. This is due to the tendency of strong shocks to decay very quickly and hence we do not expect even high Mach number flows to maintain high velocity dispersion at large length scales for long.

The next step in this work is to study turbulent decay in the presence of resistivities determined consistently from the magnetic field and the density of charged species throughout the computational domain. This will allow us to understand precisely the impact of spatially varying resistivities on turbulent decay. We can expect, for example, that the behavior of the power spectra and velocity dispersions will be strongly affected. In addition, such simulations will then incorporate all the nonideal effects likely to be of importance in the study of molecular cloud turbulence.

The authors thank the anonymous referee for useful suggestions during the refereeing process. This material is based upon works supported by the Science Foundation Ireland under Grant No. 07/RFP/PHYF586. The authors wish to 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.

NONIDEAL MAGNETOHYDRODYNAMIC TURBULENT DECAY IN MOLECULAR CLOUDS

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

Dates

ABSTRACT

1. INTRODUCTION