DIFFUSION OF MAGNETIC FIELD AND REMOVAL OF MAGNETIC FLUX FROM CLOUDS VIA TURBULENT RECONNECTION

Astrophysical flows are known to be turbulent and magnetized. The specific role played by MHD turbulence in different branches of astrophysics is still highly debated, but it is generally regarded as important. In particular, for the interstellar medium (ISM) and star formation, the role of turbulence has been discussed in many reviews (see Elmegreen & Scalo 2004; McKee & Ostriker 2007). The opinion on the role of magnetic field in these environments vary from magnetic field being regarded as absolutely dominant in the processes (see Tassis & Mouschovias 2005; Galli et al. 2006) to moderately important, as in super-Alfvénic models of star formation (see Padoan et al. 2004).

The vital question that frequently permeates these debates is the diffusion of the magnetic field in astrophysical fluids. The conductivity of most of the astrophysical fluids is high enough to make the Ohmic diffusion negligible on the scales involved, which means that the "frozen-in" approximation is a good one for many astrophysical environments. However, without considering diffusive mechanisms that can violate the flux freezing, one faces problems attempting to explain many observational facts. For example, simple estimates show that if all the magnetic flux is brought together with the material that collapses to form a star in molecular clouds, then the magnetic field in a protostar should be several orders of magnitude higher than the one observed in T Tauri stars (this is the "magnetic flux problem;" see Galli et al. 2006 and references therein, for example).

To address the problem of the magnetic field diffusion both in the partially ionized ISM and in molecular clouds, researchers usually appeal to the ambipolar diffusion concept (see Mestel & Spitzer 1956; Shu 1983). The idea of the ambipolar diffusion is very simple and may be easily exemplified in the case of gas collapsing to form a protostar. As the magnetic field is acting on charged particles only, it does not directly affect neutrals. Neutrals move under the gravitational pull but are scattered by collisions with ions and charged dust grains which are coupled with the magnetic field. The resulting flow dominated by the neutrals will be unable to drag the magnetic field lines and these will diffuse away through the infalling matter. This process of ambipolar diffusion becomes faster as the ionization ratio decreases and therefore, becomes more important in poorly ionized cloud cores.

Shu et al. (2006) have explored the accretion phase in low-mass star formation and concluded that there should exist an effective diffusivity about 4 orders of magnitude larger than Ohmic diffusivity in order for an efficient magnetic flux transport to occur. They have argued that ambipolar diffusion could work, but only under rather special circumstances like, for instance, considering particular dust grain sizes. In other words, at the moment it is unclear if ambipolar diffusion is really high enough to solve the magnetic flux transport problem in collapsing flows.

Does magnetic field remain absolutely frozen-in within highly ionized astrophysical fluids? The answer to this question affects the description of numerous essential processes in the interstellar and intergalactic gas.

Magnetic reconnection was appealed to in Lazarian (2005) as a way of removing magnetic flux from gravitating clouds, e.g., from star-forming clouds. That work referred to the reconnection model of Lazarian & Vishniac (1999) and Lazarian et al. (2004) for the justification of the concept of fast magnetic reconnection in the presence of turbulence. The advantage of the scheme proposed by Lazarian (2005) was that robust removal of magnetic flux can be accomplished both in partially and fully ionized plasma, with only marginal dependence on the ionization state of the gas.⁴ The concept of "reconnection diffusion" introduced in Lazarian (2005) is relevant to our understanding of many basic astrophysical processes. In particular, it suggests that the classical textbook description of molecular clouds supported both by hourglass magnetic field and turbulence is not self-consistent. Indeed, turbulence is expected to induce "reconnection diffusion" which should enable fast magnetic field removal from the cloud. However, in the absence of numerical confirmation of the fast reconnection, the scheme of magnetic flux removal through "reconnection diffusion" as opposed to ambipolar diffusion stayed somewhat speculative.

Fortunately, it has been recently shown numerically (see Kowal et al. 2009) that three-dimensional magnetic reconnection in turbulent fluid follows the predictions of the reconnection model of Lazarian & Vishniac (1999) and therefore, is fast. This naturally increased the interest to the "reconnection diffusion" (see Lazarian & Vishniac 2009). Motivated by this fact, here we perform simulations aiming to gain understanding of the diffusion of magnetic field induced by turbulence. As Kowal et al. (2009) tested magnetic reconnection in fully ionized gas in the present paper we shall focus our efforts on one fluid MHD simulation.

We shall compare "reconnection diffusion" with the ambipolar diffusion and discuss the effects of ambipolar diffusion qualitatively focusing on the comparison of our results on "reconnection diffusion" with those on "ambipolar turbulent diffusion" described in Heitsch et al. (2004). The latter study reported the enhancement of ambipolar diffusion in the presence of turbulence, which raised the question of how important is the simultaneous action of turbulence and ambipolar diffusion and whether turbulence alone, i.e., without any effect from ambipolar diffusion, can equally well induce de-correlation of magnetic field and density.

Our work on "reconnection diffusion" should also be distinguished from the research on the de-correlation of magnetic field and density within compressible turbulent fluctuations. Cho & Lazarian (2002, 2003) performed three-dimensional MHD simulations and reported the existence of separate turbulent cascades of Alfvén and fast modes in strongly driven turbulence as well as a cascade of slow modes driven by Alfvénic cascade. Slow modes in magnetically dominated plasma are associated with density perturbations with marginal perturbation of magnetic fields, while the same is true for fast modes in weakly magnetized or high beta plasmas. Naturally, these two modes de-correlate magnetic fields and density on the crossing time of the wave. This was the effect studied in more detail in one-dimension both analytically and numerically by Passot & Vázquez-Semadeni (2003), who stressed that the enhancements of magnetic field strength and density may correlate and anti-correlate in turbulent interstellar gas within the fluctuations and this can introduce the dispersion of the mass-to-flux ratios within the turbulent volume. Each of the fluctuations provides a transient change of the pointwise magnetization. In the absence of other effects, e.g., related to the thermal instability, the de-correlation is reversible. In comparison, the "reconnection diffusion," similar to the ambipolar diffusion, deals with the permanent de-correlation of magnetic field and density. Acting alone, the "reconnection diffusion" increases entropy, making magnetic field-density de-correlation irreversible.

What are the laws that govern magnetic field diffusion in turbulent magnetized fluids? Could those affect our understanding of basic interstellar and star formation processes? These are the questions that we address in this paper.

In this study, we try to understand the diffusion of magnetic field in a couple of idealized models in the presence of turbulence. We explore setups both with and without gravity and compare the diffusion of magnetic field with that of a passive scalar. In the context of star formation an important issue that we will address is an alternative way of decreasing the magnetic flux-to-mass ratio without appealing to ambipolar diffusion. We claim that since turbulence in astrophysics is really ubiquitous, our results should be widely applicable. We also perform simulations including turbulent heat diffusion, which allow us to compare results obtained with our code with those in the literature.

This paper is organized as follows. In Section 2, we draw the theoretical grounds about fast magnetic reconnection. In Section 3, we describe the numerical code employed. In Section 4, we present the results concerning the diffusion of magnetic field in a setup without external gravitational forces. In Section 5, we present the results of our experiments of diffusion of magnetic field in the presence of a gravitational field. In Section 6, we discuss our results and compare them with previous works. In Section 7, we discuss the accomplishments and limitations of our present study. In Section 8, we discuss our findings in the context of strong turbulence theory, and finally in Section 9, we summarize our conclusions. While our work is focused on the diffusion of magnetic fields, we address in Appendices A and B the heat transport in magnetized turbulent plasma. We confirm with higher resolution the results in Cho et al. (2003a) that the heat advection within turbulent flows in the presence of magnetic fields is very similar to that induced by hydrodynamic turbulence.

The dynamical response of magnetic fields in turbulent fluids, as we discussed above, depends on the ability of magnetic fields to change their topology via reconnection. We know from observations that magnetic field reconnection may be both fast and slow. Indeed, a slow phase of reconnection is necessary in order to explain the accumulation of free energy associated with the magnetic flux that precedes eruptive flares in magnetized coronae. Thus, it is important to identify the conditions for the reconnection to be fast. Different mechanisms prescribe different necessary requirements for this to happen.

The problem of magnetic reconnection is most frequently discussed in terms of solar flares. However, this is a general basic process underlying the dynamics of magnetized fluids in general. If the magnetic field lines in a turbulent fluid do not easily reconnect, the properties of the fluid should be dominated by intersecting magnetic flux tubes which are unable to pass through each other. Such fluids cannot be simulated with the existing codes as magnetic flux tubes readily reconnect in the numerical simulations which are currently very diffusive compared to the actual astrophysical flows.

The famous Sweet–Parker model of reconnection (Sweet 1958; Parker 1958; see Figure 1, upper panel) produces reconnection rates which are smaller than the Alfvén velocity by a square root of the Lundquist number, i.e., by S^−1/2 ≡ (L_xV_A/η)^−1/2, where L_x in this case is the length of the current sheet. Astrophysical values of S can be as large as 10¹⁵  or  10²⁰; thus, this scheme produces reconnection at a rate which is negligible for most of the astrophysical circumstances. If the Sweet–Parker were proven to be the only possible model of reconnection, it would have been possible to show that MHD numerical simulations do not have anything to do with real astrophysical fluids. Fortunately, faster schemes of reconnection are available.

The first model of fast reconnection proposed by Petschek (1964) assumed that magnetic fluxes get into contact not along the astrophysically large scales of L_x, but instead over a scale comparable to the resistive thickness δ, forming a distinct X-point, where magnetic field lines of the interacting fluxes converge at a sharp point to the reconnection spot. The stability of such a reconnection geometry in astrophysical situations is an open issue. At least for uniform resistivity, this configuration was proven to be unstable and to revert to a Sweet–Parker configuration (Biskamp 1986; Uzdensky & Kulsrud 2000).

Recent years have been marked by the progress in understanding some of the key processes of reconnection in astrophysical plasmas. In particular, substantial progress has been obtained by considering reconnection in the presence of the Hall effect (Shay et al. 1998, 2004). The condition for which the Hall–MHD term becomes important for reconnection is that the ion skin depth δ_ion becomes comparable with the Sweet–Parker diffusion scale δ_SP. The ion skin depth is a microscopic characteristic and it can be viewed as the gyroradius of an ion moving at the Alfvén speed, i.e., δ_ion = V_A/ω_ci, where ω_ci is the cyclotron frequency of an icon. For the parameters of the ISM (see Table 1 in Draine & Lazarian 1998), the reconnection is collisional (see further discussion in Yamada et al. 2006).

To deal with both collisional and collisionless plasma, Lazarian & Vishniac (1999, henceforth LV99) proposed a model of fast reconnection in the presence of weak turbulence when magnetic field back-reaction is extremely important. We stress that the LV99 model does not assume that the magnetic field can be easily bent by fluid motions, which is a usual assumption to "justify" the erroneous concept of turbulent diffusivity in the mean field dynamo.

The middle and bottom panels of Figure 1 illustrate the key components of the LV99 model.⁵ The reconnection events happen on small scales λ_|| where magnetic field lines get into contact. As the number of independent reconnection events that take place simultaneously is L_x/λ_|| ≫ 1 the resulting reconnection speed is not limited by the speed of individual events on the scale λ_||. Instead, the constraint on the reconnection speed comes from the thickness of the outflow reconnection region Δ, which is determined by the magnetic field wandering in a turbulent fluid. The model is intrinsically three dimensional⁶ as both field wandering and simultaneous entry of many independent field patches, as shown in Figure 1, are three-dimensional effects. In the LV99 model, the magnetic reconnection speed becomes comparable with V_A when the scale of magnetic field wandering Δ becomes comparable with L_x.⁷

Zoom In Zoom Out Reset image size

For a quantitative description of the reconnection, one should adopt a model of MHD turbulence (see Iroshnikov 1963; Kraichnan 1965; Dobrowolny et al. 1980; Shebalin et al. 1983; Montgomery & Turner 1981; Higdon 1984). The most important for magnetic field wandering is the Alfvénic component. Adopting the Goldreich & Sridhar (1995, henceforth GS95) scaling of the Alfvénic component of MHD turbulence extended to include the case of weak turbulence, LV99 predicted that the reconnection speed in a weakly turbulent magnetic field is

$$\begin{equation} V_{R}=V_A (l/L_x)^{1/2} (v_{l}/V_A)^2, \end{equation} \tag{ 1 }$$

where the level of turbulence is parameterized by the injection velocity, the combination V_A(v_l/V_A)² is the velocity of the largest strong turbulence eddies V_strong, i.e., the velocity at the scale at which the Alfvénic turbulence transfers from the weak to the strong regimes. Thus Equation (1) can also be rewritten as V_R = V_strong(l/L_x)^1/2; v_l < V_A and l is the turbulence injection scale.

The scaling predictions given by Equation (1) have been tested successfully by three-dimensional MHD numerical simulations in Kowal et al. (2009). This stimulates us to adopt the LV99 model as a starting point for our discussion of magnetic reconnection.

How can λ_|| be determined? In the LV99 model, as many as L²_x/λ_⊥λ_|| localized reconnection events take place, each of which reconnects the flux at the rate V_rec,local/λ_⊥, where V_rec,local is the velocity of local reconnection events at the scale λ_||. The individual reconnection events contribute to the global reconnection rate, which in three dimensions becomes a factor of L_x/λ_|| larger, i.e.,

$$\begin{equation} V_{\rm rec,global}\approx L_x/\lambda _{\Vert } V_{\rm rec, local}. \end{equation} \tag{ 2 }$$

The local reconnection speed, conservatively assuming that the local events are happening at the Sweet–Parker rate, can be easily obtained by identifying the local resistive region δ_SP with λ_⊥ and using the relations between λ_|| and λ_⊥ that follow from the MHD turbulence model. The corresponding calculations in LV99 provided the local reconnection rate v_lS^−1/4. Substituting this local reconnection rate in Equation (2) one estimates the global reconnection speed, if this speed were limited by Ohmic resistivity, which is larger than V_A by a factor S^1/4. As a result, one has to conclude that the reconnection does not depend on resistivity.

The LV99 model is applicable to situations when plasma effects are included, e.g., the Hall–MHD effect, which increases effective resistivities for local reconnection events. While the latter point is difficult to test directly with existing plasma codes, e.g., with PIC codes, due to the necessity of simulating both plasma microphysics effects as well as macrophysical effects of magnetic turbulence, Kowal et al. (2009) simulated the action of plasma effects by parameterizing them via anomalous resistivities. The values of such resistivities are a steep function of the separation between oppositely directed magnetic field lines, which also determines the current separating magnetic fluxes. With anomalous resistivities, the structure of the fractal current sheet of the turbulent reconnection changed substantially, but no significant changes of the reconnection rate were reported, which agrees well with the theoretical expectations of LV99. Within the model the explanation of this stems from the fact that the reconnection is already fast (i.e., independent of resistivity) even when small-scale reconnection events are mediated by the Ohmic resistivity, while the bottleneck for the reconnection process is provided by magnetic field wandering. Thus, the increase of the local reconnection rate does not increase the global reconnection speed.

In this work, we address the problems that are relevant to the reconnection in a partially ionized, weakly turbulent gas. The corresponding model of reconnection was proposed in Lazarian et al. (2004, henceforth LVC04). The extensive calculations summarized in Table 1 in LVC04 show that the reconnection for realistic circumstances varies from 0.1V_A to 0.03V_A, i.e., is also fast, which should enable fast diffusion arising from turbulent motions.

The fact that magnetic fields reconnect fast in turbulent fluids ensures that the large-scale dynamics that we can reproduce well with numerical codes is not compromised by the difference in reconnection processes in the computer and in astrophysical flows. This motivates our present study in which we investigate diffusion processes in turbulent magnetized fluids via three-dimensional simulations. We feel that the issue of the dimension of the simulations is crucial for the turbulence dynamics and the reconnection and, as a result, for our final results.

Visualization of turbulent diffusion of heat or any passive scalar field is easy within the GS95 model, which can be interpreted as a model of Kolmogorov cascade perpendicular to the local direction of the magnetic field (LV99; more discussion in Section 8). The corresponding eddies are expected to advect heat similarly to the case of the hydrodynamic heat advection (see Appendix A). The corresponding visualization of the magnetic field diffusion is more involved. Every time that magnetic field lines intersect each other, they change their configuration draining free energy from the system. In the presence of self-gravity this may mean the escape of magnetic field, which is a "light fluid" from the self-gravitating gaseous "heavy fluid." Naturally, if the turbulence gets very strong, the system gets unbound and then the mixing of magnetic field and gas, rather than their segregation, is expected. In what follows we test these ideas numerically.

The systems studied numerically in this work are described by the resistive MHD equations, assuming an isothermal equation of state:

$$\begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot \left(\rho \mathbf {u} \right) = 0 \end{equation} \tag{ 3 }$$

$$\begin{equation} \rho \left(\frac{\partial }{\partial t} + \mathbf {u} \cdot \nabla \right) \mathbf {u} = - c_{s}^{2} \nabla \rho + (\nabla \times \mathbf {B}) \times \mathbf {B} - \rho \nabla \Psi + \mathbf {f} \end{equation} \tag{ 4 }$$

$$\begin{equation} \frac{\partial \mathbf {B}}{\partial t} = \nabla \times \left(\mathbf {u} \times \mathbf {B} \right) + \eta _{\rm Ohm} \nabla ^{2} \mathbf {B} \end{equation} \tag{ 5 }$$

plus the divergenceless condition for the magnetic field ∇ · B = 0. The spatial coordinates are given in units of a typical length L_*. The density ρ is normalized by a reference density ρ_*, and the velocity field u by a reference velocity U_*. The constant sound speed c_s is also given in units of U_*, and the magnetic field B is measured in units of $U_{*} \sqrt{4 \pi \rho _{*}}$. The uniform Ohmic resistivity η_Ohm is given in units of U_*L_*. In our numerical calculations we will use both non-zero and zero values of η_Ohm. In the latter case, the calculations will include only numerical resistivity. Time t is measured in units of L_*/U_*. The external gravitational potential Ψ is given in units of U²_*. The source term f is a random force term responsible for injection of turbulence.

The above equations are solved inside a three-dimensional box with periodic boundary conditions. We use a shock-capturing Godunov-type scheme with an HLL solver (see, for example, Kowal et al. 2007; Falceta-Gonçalves et al. 2008). Time integration is performed with the Runge–Kutta method of second order. Unless we say explicitly the opposite, we assume η_Ohm = 0.

We employ an isotropic, non-helical, solenoidal, delta correlated in time forcing f. This forcing acts in a thin shell around the wave number k_f = 2.5(2π/L), that is, the scale of turbulence injection l_inj that is about 2.5 times smaller than the box size L (in all the simulations, we choose L = 1 in code units). In most of the experiments, the rms velocity V_rms induced by turbulence in the box is close to unity (in code units). Therefore, in these cases, the turnover time of the energy-carrying eddies (or the turbulent timescale t_turb) is t_turb ∼ l_inj/v_turb ∼ (L/2.5)/V_rms ≈ 0.4 units of time in code units.

We note that for our present purposes we use a one-fluid approximation, which does not include ambipolar diffusion. This choice is appropriate for approximating a fully ionized gas. One may argue that the code also describes the dynamics of partially ionized gas, but on the scales where ions and neutrals are strongly coupled, i.e., on scales larger than the scale of ambipolar diffusion (see discussion in Lazarian et al. 2004).

Observations of different regions of the diffuse ISM compiled by Troland & Heiles (1986) indicate that magnetic fields and density are not straightforwardly correlated. These observations motivated Heitsch et al. (2004) to perform 2.5 dimensional numerical calculations in the presence of both ambipolar diffusion and turbulence. The results in Heitsch et al. (2004) also indicated de-correlation of magnetic field and density⁸ and one may wonder whether ambipolar diffusion is always required to de-correlate magnetic field and density or if, otherwise, to what extent the concept of "turbulent ambipolar diffusion" introduced in Heitsch et al. (2004) is useful (see also Zweibel 2002). To address these issues we performed three-dimensional simulations of magnetic diffusion in the absence of ambipolar diffusion effects.

4.1. Numerical Setup

The magnetic field is assumed to have initially only the component in the z-direction. The initial configuration of the magnetic and density fields are

$$\begin{equation} B_{z} (x,y) = B_{0} + B_{1} \cos \left(\frac{2\pi }{L} x \right) \cos \left(\frac{2\pi }{L} y \right) \end{equation} \tag{ 6 }$$

$$\begin{eqnarray} \rho (x,y) = & \rho _{0} - \frac{1}{c_{s}^{2}} \left\lbrace B_{0} B_{1} \cos \left(\frac{2\pi }{L} x \right) \cos \left(\frac{2\pi }{L} y \right) \right. \nonumber \\ & + 0.5 \left.\left[ B_{1} \cos \left(\frac{2\pi }{L} x \right)\cos \left(\frac{2\pi }{L} y \right) \right] ^{2} \right\rbrace, \end{eqnarray} \tag{ 7 }$$

where (x, y) = (0, 0) is the center of the x, y-plane. Boundary conditions are periodic.

This initial magnetic field configuration has an uniform component B₀ plus a harmonic perturbation of amplitude B₁. The density field is distributed in such a way that the gas pressure, given by the isothermal equation of state p = c²_sρ, equilibrates exactly the magnetic pressure, giving a magneto-hydrostatic solution. We choose the parameters B₁ = 0.3, ρ₀ = 1, and c_s = 1 in all our simulations. Figure 2 illustrates these initial fields when B₀ = 1.0.

Zoom In Zoom Out Reset image size

We should remark that the simulations presented in Heitsch et al. (2004) do not start at the equilibrium, like ours. There is no pressure term in their equation for the evolution of the momentum of the ions to counterbalance the magnetic pressure. In addition, in their work the ion-density field is kept constant in time and space.

Another difference between our setup and that in Heitsch et al. (2004) is that our parameter B₁ is assumed to be the same for all the models studied, and not a fraction of B₀. Also, the amplitude of the perturbation of the homogeneous component of the density field is a free parameter in Heitsch et al. (2004), while here it is constrained by the imposed equilibrium between gas and magnetic pressures.

In addition, we introduce a passive scalar field Φ initially identical to B_z. The parameters of our relevant simulations are presented in Table 1.

Table 1. Parameters of the Simulations in the Study of Turbulent Diffusion of Magnetic Flux without Gravity

Model	B0	linj	Vrms	tturb	ηturb − Bz	ηturb − Φ	Resolution
B1	0.5	0.4	0.8	0.5	0.18(0.03)	0.19(0.02)	2563
B2	1.0	0.4	0.8	0.5	0.09(0.02)	0.14(0.02)	2563
B3	1.5	0.4	0.8	0.5	0.10(0.02)	0.08(0.02)	2563
B4	2.0	0.4	0.8	0.5	0.13(0.02)	0.11(0.02)	2563
B2l	1.0	0.4	0.8	0.5	0.15(0.02)	0.11(0.02)	1283
B2h	1.0	0.4	0.8	0.5	0.10(0.01)	0.13(0.02)	5123

Download table as:  ASCIITypeset image

In these simulations we keep the random velocity approximately constant, i.e., V_rms ≈ 0.8 for all the models after one time step. Therefore, all these models are slightly subsonic.

4.2. Notation

4.3. Results

Figure 3 shows the evolution of the amplitude of the mode that is identical to the initial harmonic perturbation of the magnetic field (i.e., the rms of the amplitude of the Fourier modes (k_x, k_y) = (±1, ± 1)), for 〈B_z〉_z and 〈B_z〉_z/〈ρ〉_z. Most right plot in Figure 14 (Appendix B) shows the evolution of the amplitude of the same mode for 〈Φ〉_z and 〈Φ〉_z/〈ρ〉_z. All the curves presented were smoothed in order to make the visualization clearer. We see that the decay of the magnetic field occurs at a similar rate to that of the passive field. The mode decays nearly exponentially at roughly the same rate for most of the models. Only the model B1 (B₀ = 0.5) exhibits a higher decay rate. This may be due to the large-scale field reversals that are common in super-Alfvénic turbulence. Table 1 shows the fitted values (and the uncertainty) of η_turb in the curves corresponding to the evolution of the amplitude of the modes for 〈B_z〉_z. The fitted curve is exp{−k²η_turbt}, where k² = k²_x + k²_y = 2 is the square of the module of the corresponding wave-vector. We observe that the decay of the amplitude of the modes of 〈B_z〉_z is not continuous but saturates at a value that is naturally maintained by the turbulence (in Figure 3 it occurs after about t = 6).

Zoom In Zoom Out Reset image size

The diffusion of B/ρ on large scales was also observed in Heitsch et al. (2004) for two-fluid simulations and there it was associated with the difference between the velocity field of the ions and neutrals, at small scales. However, here we observe a similar effect, but in one-fluid simulations, which is suggestive that turbulence rather than the details of the microphysics are responsible for the diffusion.

The left and center panels of Figure 4 shows the distribution of 〈ρ〉_z versus 〈B_z〉_z for the model B2 (B₀ = 1.0) at the initial configuration (t = 0) and after 10 time steps. We see in this projected view that the initial magnetic flux-to-mass relation is quickly spread and, in contrast with the Φ–ρ distribution (see Figure 15 in Appendix B), we do not see any tendency for the magnetic field and density to become correlated.

Zoom In Zoom Out Reset image size

To give a quantitative measure of the evolution of the flux-to-mass relation in the models, let us consider 〈δB, δρ〉, the correlation between fluctuations of the magnetic field δB and density δρ, defined by

$$\begin{equation} \left\langle \delta B, \delta \rho \right\rangle = \frac{\int (B - \bar{B}) (\rho - \bar{\rho }) d^{3}x}{ \sqrt{\int (B - \bar{B})^{2} d^{3}x} \sqrt{\int (\rho - \bar{\rho })^{2} d^{3}x} }. \end{equation} \tag{ 8 }$$

The right panel of Figure 4 shows the evolution of 〈δB, δρ〉 (see right side of Figure 15 for the evolution of 〈δΦ, δρ〉 which is similarly defined). Different from the passive scalar field that quickly becomes correlated to the density field, the magnetic field keeps a residual anti-correlation with it.

A more careful analysis of our results indicates that the correlation between magnetic field intensity and density depends on the Mach number M_s. For example, when we calculate the correlation 〈δB, δρ〉 using the turbulent models of Table 5 (for study of diffusion of passive scalar fields, see Appendix A), we find weak positive correlations for the supersonic models and negative correlations for the subsonic ones. These correlations increase with M_s. Thus, the anti-correlation detected in Figure 4 can be due to the slightly subsonic regime of the turbulence. These correlations and anti-correlations at this level cannot be excluded by the observational data as discussed, e.g., in Troland & Heiles (1986). We shall address this issue in more detail elsewhere.

To summarize, the results of Figures 3 and 4 suggest that the turbulence can substantially change the flux-to-mass ratio B/ρ without any effect of ambipolar diffusion. The diffusion of the magnetic flux occurs in a rate similar to the rate of the turbulent diffusion of heat (passive scalar), even for sub-Alfvénic turbulence.

As remarked in Section 2, we should emphasize that the efficient turbulent diffusion of magnetic field that we are observing in the simulations above is due to fast magnetic reconnection because otherwise, if the tangled magnetic lines by turbulence were not reconnecting, then they would be behaving like a Jello-type substance and this would make the diffusive transport of magnetic flux very inefficient (contrary to what is observed in the simulations). The issue of magnetic reconnection was avoided in Heitsch et al. (2004) due to the settings in which magnetic field was assumed perpendicular to the plane of the fluid motions. Magnetic reconnection, however, is an effect that is present within any realistic three-dimensional setup.

4.4. Effects of Resolution on the Results

To convince the reader that the above results are not being affected by numerical effects, we ran one of the models (model B2) with increased and decreased resolutions (models B2h and B2l, respectively; see Table 1). Figure 5 compares the same quantities presented in Figure 3 for these models. We do not observe significant difference between them. Thus, we can expect that the results presented for the models with resolution of 256³ are robust.

Zoom In Zoom Out Reset image size

The "reconnection diffusion" of magnetic field in the absence of gravity can represent the magnetic field dynamics in diffuse interstellar gas where cloud self-gravity is not important. In molecular clouds, clumps, and accretion disks, the diffusion of magnetic field should happen in the presence of gravity. We study this process below.

5.1. Numerical Approach

In order to get an insight into the magnetic field diffusion in a turbulent fluid immersed in a gravitational potential, we have performed experiments in the presence of a gravitational potential with cylindric symmetry Ψ, given in cylindrical coordinates (R, ϕ, z) by

$$\begin{equation} \Psi (R \le R_{\rm max}) = - \frac{A}{R + R_{*}}, \end{equation} \tag{ 9 }$$

$$\begin{equation} \Psi (R > R_{\rm max}) = - \frac{A}{R_{\rm max} + R_{*}}, \end{equation} \tag{ 10 }$$

where R = 0 is the center of the (x, y)-plane, and we fixed R_* = 0.1L and R_max = 0.45L (L = 1 is the size of the computational box, as remarked in Section 3). We assume a relatively high value of R_* in order to limit the values of the gravitational force and prevent the system to be initially Parker–Rayleigh–Taylor unstable. We assume an outer cut-off R_max on the gravitational force to ensure the cylindrical symmetry while using periodic boundary conditions.

In one class of experiments, we start the simulation with a magneto-hydrostatic equilibrium with β = P_gas/P_mag = c²_sρ/(B²/8π) constant. The initial density and magnetic fields are, respectively,

$$\begin{equation} \rho (R) = \rho _0 \exp \left\lbrace (\Psi (R_{\rm max}) - \Psi (R)) / c_{s}^{2} (1 + \beta ^{-1}) \right\rbrace \end{equation} \tag{ 11 }$$

$$\begin{equation} B_{z} (R) = c_{s} \sqrt{2 \beta ^{-1} \rho (R)}. \end{equation} \tag{ 12 }$$

Figure 6 illustrates this initial configuration for one of the studied models (model C2, see Table 2).

Zoom In Zoom Out Reset image size

Table 2. Parameters for the Models with Gravity Starting at Magneto-hydrostatic Equilibrium with Initial Constant β

Model	βa	VA,i	Ab	$\bar{\rho }$	$\bar{B_{z}}$	Vrms	tturb	ηturb	ηturb/Vrmslinj	Resolution
C1	1.0	1.4	0.6	1.26	1.59	0.8	0.5	⪅0.005	⪅0.015	2563
C2	1.0	1.4	0.9	1.52	1.74	0.8	0.5	≈0.01	≈0.03	2563
C3	1.0	1.4	1.2	1.95	1.98	0.8	0.5	≈0.03	≈0.09	2563
C4	1.0	1.4	0.9	1.52	1.74	1.4	0.3	≈0.10 − 0.20	≈0.18 − 0.36	2563
C5	1.0	1.4	0.9	1.52	1.74	2.0	0.2	⪆0.30	⪆0.37	2563
C6	3.3	0.8	0.9	2.40	1.20	0.8	0.5	≈0.02	≈0.06	2563
C7	0.3	2.4	0.9	1.18	2.66	0.8	0.5	≈0.01	≈0.03	2563
C2l	1.0	1.4	0.9	1.52	1.74	0.8	0.5	⋅⋅⋅	⋅⋅⋅	1283
C2h	1.0	1.4	0.9	1.52	1.74	0.8	0.5	⋅⋅⋅	⋅⋅⋅	5123

Notes. ^aInitial β parameter for the plasma: β = P_gas/P_mag. ^bThe parameter A for the strength of gravity (see Equations (9) and (10)) is given in units of c²_sL.

Download table as:  ASCIITypeset image

We restricted our experiments to the trans-sonic case c_s = 1 (in most of the experiments, we keep V_rms ≈ 0.8, see Table 2). We also fixed ρ₀ = 1. The turbulence is injected at t = 0 and we follow the evolution of 〈B_z〉_R and 〈ρ〉_R for eight time steps. Table 2 lists the parameters used for these experiments. $\bar{\rho }$ and $\bar{B_{z}}$ represent the average of the density and magnetic field over the entire box. V_A,i refers to the initial Alfvén speed of the system. The rms velocity of the system V_rms is measured after the turbulence is well-developed.

We have also performed experiments starting out of equilibrium, with homogeneous fields: the system starts in free fall. We leave the system to evolve for eight time steps, applying turbulence from the very beginning. For comparison, we have also performed these experiments without turbulence. The initial uniform magnetic field is parallel to the z-direction for these models. Table 3 lists the parameters for these runs. The listed values of β refer to the initial conditions.

Table 3. Parameters for the Models with Gravity Starting Out-of-equilibrium, with Initially Homogeneous Fields

Model	β	VA,i	A	$\bar{\rho }$	$\bar{B_{z}}$	Vturb	Resolution
D1	1.0	1.4	0.9	1.0	1.41	0.8	2563
D1a	1.0	1.4	0.9	1.0	1.41	0.0	2563
D2	3.3	0.8	0.9	1.0	0.77	0.8	2563
D2a	3.3	0.8	0.9	1.0	0.77	0.0	2563
D3	0.3	2.4	0.9	1.0	2.45	0.8	2563
D3a	0.3	2.4	0.9	1.0	2.45	0.0	2563
D1l	1.0	1.4	0.9	1.0	1.41	0.8	1283
D1h	1.0	1.4	0.9	1.0	1.41	0.8	5123

Download table as:  ASCIITypeset image

Concerning the diffusion of the magnetic field, in order to provide a quantitative comparison between the models, we have also performed simulations with similar initial conditions to the models presented in Table 2, but without turbulence and with the explicit presence of Ohmic diffusivity η_Ohm in the induction equation. As these models have perfect symmetry in the z-direction, we simulated only a plane cutting the z-axis, that is, they are 2.5 dimensional simulations. We use a resolution comparable to the turbulent three-dimensional models. Table 4 lists the parameters for these runs. We simulated a three-dimensional model equivalent to the model E7r1, and we found exact agreement in the time evolution of the magnetic flux distribution (not shown). Therefore, we can believe that in this case these 2.5 dimensional simulations give results which are equivalent to three-dimensional simulations.

Table 4. Parameters for the 2.5 dimensional Resistive Models with Gravity Starting with Magneto-hydrostatic Equilibrium and Constant β

Model	β	A	ηOhm	Resolution
E1r0	1.0	0.6	0.005	2562
E1r1	1.0	0.6	0.01	2562
E2r1	1.0	0.9	0.01	2562
E2r2	1.0	0.9	0.02	2562
E2r3	1.0	0.9	0.03	2562
E2r4	1.0	0.9	0.05	2562
E3r1	1.0	1.2	0.01	2562
E3r2	1.0	1.2	0.02	2562
E3r3	1.0	1.2	0.03	2562
E4r2	1.0	0.9	0.10	2562
E4r3	1.0	0.9	0.20	2562
E5r3	1.0	0.9	0.30	2562
E6r1	3.3	0.9	0.01	2562
E6r2	3.3	0.9	0.02	2562
E6r3	3.3	0.9	0.03	2562
E7r1	0.3	0.9	0.01	2562

Download table as:  ASCIITypeset image

5.2. Results

5.2.1. Evolution of the Equilibrium Distribution

Top row of Figure 7 shows the evolution of 〈B_z〉_0.25 (left), 〈ρ〉_0.25 (middle), and 〈B_z〉_0.25/〈ρ〉_0.25 (right), normalized by the respective characteristic values inside the box ($\bar{B_{z}}$, $\bar{\rho }$ and $\bar{B_{z}}/\bar{\rho }$), for the models C1, C2, and C3 (β = 1). We compare the evolution of these quantities for different strengths of gravity A, maintaining the other parameters identical. The central magnetic flux reduces faster the higher the value of A. The flux-to-mass ratio has similar behavior. The other plots of Figure 7 show the profile of the quantities 〈B_z〉_z (upper panels), 〈ρ〉_z (middle panels), and 〈B_z〉_z/〈ρ〉_z (bottom panels) along the radius R, each column corresponding to a different value of A both for t = 0 (in magneto-hydrostatic equilibrium and constant β) and for t = 8. We see the deepest decay of the magnetic flux toward the central region for the highest value of A at t = 8.

Zoom In Zoom Out Reset image size

In Figure 8, we compare the rate of the "reconnection diffusion" when we change the turbulent velocity and maintain the other parameters identical as in models C2, C4, and C5. An inspection of the left panel shows that the central magnetic flux 〈B_z〉_0.25 decreases faster for the two highest turbulent velocities. Fluctuations are higher in the cases with higher velocity. The central density, however, gets smaller for the highest turbulent velocity. This is explained by the fact that the dynamic pressure is higher for the largest velocities. The central flux-to-mass ratio 〈B_z〉_0.25/〈ρ〉_0.25 decreases for the two smallest velocities. However, for the largest velocity, it is not clear if this ratio decays or not. Looking at the middle graph of Figure 8 (bottom row), we see that the central density decreases for the highest forcing. This indicates that the turbulence driving overcomes the gravitational potential making the system less bound.

Zoom In Zoom Out Reset image size

Top row of Figure 9 compares the evolution of 〈B_z〉_0.25 (left), 〈ρ〉_0.25 (middle), and 〈B_z〉_0.25/〈ρ〉_0.25 (right), normalized by the characteristic average values inside the box ($\bar{B_{z}}$, $\bar{\rho }$, and $\bar{B_{z}}/\bar{\rho }$), for models C2, C6, and C7 with different β. Both the central magnetic flux and the flux-to-mass ratio decrease faster for the less magnetized model (β = 3.3). The other plots of Figure 9 show the radial profile of the quantities 〈B_z〉_z (upper panels), 〈ρ〉_z (middle panels), and 〈B_z〉_z/〈ρ〉_z (bottom panels) for each model. We can again observe a lower value of the flux in the central region (relative to the external regions) for the highest values of β at the time step t = 8. The contrast between the central and the more external values for the flux-to-mass ratio is quite different for the three models, being higher for the more magnetized models. This is expected, as turbulence brings the system in the state of minimal energy. The effect of varying magnetization in some sense is analogous to the effect of varying gravity. The equilibrium flux-to-mass ratio is larger in both the case of higher gravity and higher magnetization. The physics is simple, the lighter fluid (magnetic field) gets segregated from the heavier fluid (gas).

Zoom In Zoom Out Reset image size

All in all, we clearly see that turbulence substantially influences the quasi-static evolution of magnetized gas in the gravitational potential. The system in the presence of turbulence relaxes faster to its minimum potential energy state. This explains the change of the flux-to-mass ratio, which for years was a problem to deal with invoking ambipolar diffusion.

5.2.2. Equilibrium Models: Comparison of Magnetic Diffusivity and Resistivity Effects

In terms of the removal of the magnetic field from quasi-static clouds, does the effect of "reconnection diffusion" act similar to the effect of diffusion induced by resistivity? To address this question, we have performed a series of simulations with enhanced Ohmic resistivity (see models of Table 4).

In Figure 10, we compare the evolution of 〈B_z〉_R (at different radius) for model C2 of Table 2 with similar resistive models without turbulence of Table 4, with different values of Ohmic diffusivity η_Ohm. The decay seems initially faster and comparable with the highest value of η_Ohm (η_Ohm = 0.05). But after this initial phase, the turbulent model (C2) seems to have a behavior similar to the resistive models with η_Ohm between 0.01 and 0.02.

Zoom In Zoom Out Reset image size

Figure 16 (Appendix C) compares the turbulent models C1, C3, C4, C5, C6, and C7 of Table 2 with similar resistive models of Table 4. After roughly one time step, the model C1 (weaker gravitational field) seems to be consistent with a value of η_Ohm between 0.005 or lesser, while the model C3 (stronger gravitational field) seems to be consistent with η_Ohm ≈ 0.03. These results show that the effective turbulent magnetic diffusivity is sensitive to the strength of the gravitational field.

The resistive simulations with increasing η_Ohm values are more comparable with models with increasing turbulent velocity. The model C4 (with smaller turbulent velocity) seems to be consistent with the resistive model with η_Ohm = 0.10. The model C5 (with larger turbulent velocity) seems to be more comparable with the model with η_Ohm ∼ 0.30 (or higher); however, in this case we cannot associate a representative value of η_Ohm due to the very quick diffusion which occurs even before t = 1, when the turbulence becomes well developed.

The turbulent curve for the less magnetized model (C6) seems to follow the resistive curve with η_Ohm = 0.02, while the more magnetized model (C7) is comparable to the resistive model with η_Ohm ≈ 0.01. This result indicates that the effective turbulent diffusivity is also sensitive to the strength of the magnetic field.

In summary, the results above indicate a correspondence between the two different effects. In other words, the turbulent magnetic diffusion may mimic the effects of Ohmic diffusion of magnetic fields in gravitating clouds. However, we should keep in mind that the physics of turbulent diffusion and Ohmic resistivity is different. Thus this analogy should not be overstated.

5.2.3. Evolution of Non-equilibrium Models

Figure 11 shows the same set of comparisons as in Figure 9 for the models D1, D2, and D3 of Table 3—these models have started out of the equilibrium with a homogeneous density and magnetic field in a free fall system. Besides the runs with turbulence, we also present, for comparison, the evolution for the systems without turbulence (models D1a, D2a, and D3a). The strong oscillations seen in the evolution of the central magnetic flux and density for these models (which are more pronounced in the models without turbulence) are acoustic oscillations, since the time for the virialization of these systems is larger than the simulated period. We note that the initial flux-to-mass ratio does not change in the cases without turbulence. We also observe similar trends as in Figure 9: the higher the value of β, the faster the decrease of the central magnetic flux relative to the mean flux into the box. We also note that the radial profile of the flux-to-mass ratio for the turbulent models crosses the mean value for the models without turbulence at nearly the same radius. This is due to the fact that the effective gravity potential in all these simulations acts up to this radius approximately.

Zoom In Zoom Out Reset image size

This set of simulations shows that the change of mass-to-flux ratio can happen at the timescale of the gravitational collapse of the system and therefore, turbulent diffusion of magnetic field is applicable also to dynamic situations, e.g., to the formation of supercritical cores.

5.3. Effects of Resolution on the Results

As in the case of the study presented in Section 4, we would like to know how the results shown in this section are sensitive to changes in resolution. Again, we ran some models employing higher resolution and we inspected the changes in the results concerned.

Figure 12 compares the evolution of some of the quantities studied through this section for models C2l and C2h, which are identical to C2, except by the lower (C2l) and higher resolution (C2h; see Table 2). It shows no significant difference between these models. Figure 12 also depicts the evolution of the same quantities for the models D1, D1l, and D1h. Both models have the same parameters as in model D1, but model D1h (D1l) has higher (lower) resolution (see Table 3). Again, we see no disagreement between the models.

Zoom In Zoom Out Reset image size

Therefore, the results presented in this section are not expected to change with an increase in resolution.

5.4. Magnetic Field Expulsion Revealed

Through this work we have performed the comparison of our results with the study by Heitsch et al. (2004). Below, we provide yet another outlook of the connection of that study with the present paper. We also discuss the work by Shu et al. (2006), which was the initial motivation of our study of the diffusion of magnetic field in the presence of gravity.

6.1. Comparison with Heitsch et al. (2004): Ambipolar Diffusion versus Turbulence and 2.5 dimensional versus Three-dimensional

6.2. Transient De-correlation of Density and Magnetic Field

6.3. Relation to Shu et al. (2006): Fast Removal of Magnetic Flux During Star Formation

The concept of "reconnection diffusion" is related to the LV99 model of fast reconnection which makes use of the model of strong turbulence proposed by Goldreich & Sridhar (1995, henceforth GS95) the turbulence is being injected at the large scales with the injection velocity V_inj equal to the Alfvén velocity V_A (see Cho et al. 2003b for a review). The turbulent eddies mix up magnetic field mostly in the direction perpendicular to the local magnetic field thus forming a Kolmogorov-type picture in terms of perpendicular motions. Naturally, these eddies are as efficient as hydrodynamic eddies are expected in terms of heat advection. One also can visualize how such eddies can induce magnetic field diffusion.

It is important to note that the GS95 model deals with motions with respect to the local rather than the mean magnetic field. Indeed, it is natural that the motions of the parcel of fluid are affected only by the magnetic field of the parcel and of the near vicinity, i.e., by local fields. At the same time, in the reference frame of the mean field, the local magnetic fields of different parcels vary substantially. Thus we do not expect to see a substantial anisotropy of the heat advection when V_inj ∼ V_A.

It was noted in LV99 that one can talk about turbulent eddies perpendicular to the magnetic field only if the magnetic field can reconnect fast. The rates of reconnection predicted in LV99 ensured that the magnetic field changes topology over one eddy turnover period. If the reconnection were slow, the magnetic fields would form progressively complex structures consisting of unresolved knots, which would invalidate the GS95 model. The response of such a fluid to mechanical perturbations would be similar to "Jello," making the turbulence-sponsored diffusion of magnetic field and heat impossible.

What happens when V_inj < V_A? In this case the turbulence at large scales is weak and therefore magnetic field mixing is reduced. Thus one may expect a partial suppression of magnetic diffusivity. However, as turbulence cascades the strength of interactions increases and at a scale L_inj(V_inj/V_A)² the turbulence gets strong. According to Lazarian (2006), the diffusivity in this regime decreases by the ratio of (V_inj/V_A)³, with the eddies of strong turbulence playing a critical role in the process. When we compare the turbulent diffusivity η_turb estimated for the sub-Alfvénic models described in Table 2 (see Section 5) with L_injV_turb(V_turb/V_A)³, we find that the values are roughly consistent with the predictions of Lazarian (2006), although a more detailed study is required in this regard. For instance, we know that Lazarian (2006) theory was not intended for high Mach number turbulence.

All in all, we believe that the high diffusivity that we observe is related to the properties of strong magnetic turbulence. While the latter is still a theory which is subject to intensive study (see Boldyrev 2005, 2006; Beresnyak & Lazarian 2006, 2009a, 2009b; Gogoberidze 2007), we believe that for the purpose of describing magnetic and heat diffusion the existing theory and the present model catch all the essential phenomena.

8.1. Major Findings

8.2. Applicability of the Results

8.3. Magnetic Field Reconnection and Different Stages of Star Formation

8.4. Unsolved Problems and Future Studies

Motivated by a vital problem of the dynamics of magnetic fields in astrophysical fluids, in particular, by the magnetic flux removal in star formation, in this paper we have numerically studied the diffusion of magnetic field both in the absence and in the presence of gravitational potential. Recent work on validating the idea of the LV99 model of reconnection supports our assertion that our results obtained at moderate resolution represent the dynamics of turbulent magnetic field lines in astrophysics. Our findings obtained on the basis of three-dimensional MHD numerical simulations can be briefly summarized as follows.

• 1.  
  In the absence of gravitational potential the "reconnection diffusion" removes strong anti-correlations of magnetic field and density that we impose at the start of our simulations. The system after several turbulent eddy turnover times relaxes to a state with no clear correlation between magnetic field and density, reminiscent of the observations of the diffuse ISM by Troland & Heiles (1986).
• 2.  
  Our simulations that started with a quasi-static equilibrium in the presence of a gravitational potential, revealed that the turbulent diffusivity induces gas to concentrate at the center of the gravitational potential, while the magnetic field is efficiently pushed to the periphery. Thus the effect of the magnetic flux removal from collapsing clouds and cores, which is usually attributed to ambipolar diffusion effect, can be successfully accomplished without ambipolar diffusion, but in the presence of turbulence.
• 3.  
  Our simulations that started in a state of dynamical collapse induced by an external gravitational potential showed that in the absence of turbulence, the flux-to-mass ratio is preserved for the collapsing gas. On the other hand, in the presence of turbulence, fast removal of magnetic field from the center of the gravitational potential occurs. This may explain the low magnetic flux-to-mass ratio observed in stars compared to the corresponding ratio of the interstellar gas.
• 4.  
  As an enhanced Ohmic resistivity to remove magnetic flux from cores and accretion disks has been appealed to in the literature, e.g., by Shu et al. (2006), we have also compared models with a turbulent fluid and models without turbulence but with substantially enhanced Ohmic diffusivity. We have shown that, in terms of the magnetic flux removal, the reconnection diffusion can mimic the effect of an enhanced Ohmic resistivity.
• 5.  
  In addition, our results extend earlier findings of Cho et al. (2003a) for heat advection by magnetized turbulence. We show that heat advection can be parameterized by the product of the turbulence injection scale and the turbulent velocity for a range of Alfvénic and sonic Mach numbers.

R.S.L. and E.M.G.D.P. acknowledge partial support from grants of the Brazilian Agencies FAPESP (2006/50654-3 and 2007/04551-0), CAPES (PDEE 3979-08-3), and CNPq. A.L. acknowledges the NSF grant AST 0808118 and the NSF-sponsored Center for Magnetic Self-Organization. J.C. was supported in part by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST; No. 2009-0077372). Helpful comments by Ellen Zweibel and discussions with Enrique Vazquez-Semadeni and Fabian Heitsch are acknowledged.

We present here some studies of three-dimensional transport properties of MHD turbulence. Our setup is identical to that adopted in Cho et al. (2003a). Our goal is to provide simulations with higher numerical resolution (the resolution employed in Cho et al. 2003a is mostly 192³) and for a larger parameter space. In the present work, we study how compressibility and magnetization affect the results. We note that at the time when the results in Cho et al. (2003a) were obtained the issues of reconnection in turbulent media were more speculative. Thus, after numerical testing of the model of reconnection in LV99 by Kowal et al. (2009), we feel it is appropriate to revisit that domain.

Besides the ideal MHD set of equations, we also evolve four independent fields of passive scalars Φ_α(x) (α = 1, 2, 3, 4) by the continuity equation:

$$\begin{equation} \frac{\partial \Phi _{\alpha }}{\partial t} + \nabla \cdot \left(\Phi _{\alpha } \mathbf {u} \right) = 0. \end{equation} \tag{ A1 }$$

These passive scalar fields can trace, for example, the volume concentration of some physical property of the gas, like metallicity or heat,¹² as long as the time-scale associated with the molecular (microscopic) diffusion is larger than the typical dynamical time-scale of the flow. In our simulations this is ensured by fluids being turbulent.¹³

These passive scalar fields are assumed to have initial spherical symmetry with a Gaussian radial profile:

$$\begin{equation} \Phi _{\alpha } (\mathbf {x}, t=0) = \exp \left\lbrace \frac{3}{2} \frac{(\mathbf {x - x_{0}})^{2}}{\sigma _{0}^{2}} \right\rbrace, \end{equation} \tag{ A2 }$$

where x₀ is the center of the box and, as in Cho et al. (2003a), we choose the initial dispersion $\sigma _{0} = L\sqrt{3} / 16\sqrt{2}$ (≈19 grid cells, for the resolution employed of 256³), to ensure that the characteristic width σ₀ is inside the inertial range of the turbulence. This initial distribution of Φ_α is a natural choice to study its diffusion. If we had, for instance, a microscopic diffusion coefficient κ, the Gaussian shape of the distribution would remain unaltered, with the dispersion increasing linearly in time at a rate κ.

The initial magnetic field B₀ is uniform and parallel to the x-direction, and the initial density field is also uniform with ρ = 1. The passive scalar fields are injected after the turbulence is fully developed. The first scalar-field Φ₁ is injected at t = 3 (approximately 6 turnover times), and the other fields (α = 2, 3, 4) at each time step after.

Table 5 shows the set of parameters studied. The magnetic field is given in terms of the initial Alfvén speed. Turbulence is injected with the same power = 1 in all the models. These runs cover four combinations of regimes of sonic and Alfvénic Mach numbers (M_s and M_A, respectively). The quantities shown in Table 5 (M_s, M_A, V_rms, B_rms) are averages taken over the time after the injection of the first passive field (these time averages are over the available computed cubic domains taken at every 0.25 time step). The standard deviation of these averages are shown within parentheses. For each data cube, the Mach numbers correspond to the average of the local values computed over the entire box. The resolution employed in all the simulations is 256³.

A.1. Results

Compressibility can change the properties of flows substantially. For instance, high Mach number turbulence is known to create regions of enhanced density which will coexist with the expanses of low density. At the same time, one should expect to see similarity between the properties of fluid in low-Mach number flows and incompressible flows.

We first consider a case of heat transport, which can be considered as a test case, as we can compare our results with the earlier heat transport simulations in Cho et al. (2003a). This is the case where the back-reaction can be ignored and we can use a passive scalar diffusion to represent the diffusion of heat.

Let us relate the evolution of the dispersion σ with the turbulent diffusion coefficient η_turb, that is, the coefficient that gives the rate at which the passive scalar field diffuses in scales larger than the turbulent scale. The dispersion σ is calculated through the definition: $\sigma ^{2} = \frac{\int (\mathbf {x} - \bar{\mathbf {x}})^2 \Phi (\mathbf {x}) d^{3}x}{ \int \Phi (\mathbf {x}) d^{3}x}$, where $\bar{\mathbf {x}} = \frac{\int \mathbf {x} \Phi (\mathbf {x}) d^{3}x}{ \int \Phi (\mathbf {x}) d^{3}x}$. If λ is the rms distance between two fluid elements being advected, for λ higher than the injection scale of the turbulence l_inj, its evolution is related to η_turb by δλ² ∼ η_turbδt.

Considering ordinary hydrodynamic turbulence, we can write, for the rms distance l between two fluid elements, within the inertial range,

$$\begin{equation} \delta l^2 \sim v_{l} l \delta t, \end{equation} \tag{ A3 }$$

where v_l is the velocity at scale l. From the assumption that both relations above should be valid at the scale l_inj, it must be true that η_turb ∼ l_injv_turb, and we could define a constant C_dyn so that

$$\begin{equation} \eta _{\rm turb} = C_{\rm dyn} l_{\rm inj} v_{\rm turb}. \end{equation} \tag{ A4 }$$

Now, one wonders if C_dyn is reduced by the magnetic field in magnetized turbulence (due to suppression of turbulent mixture). C_dyn is related to the constant of proportionality in Equation (A3). Again, let l be the distance between two fluid elements initially separated by a distance l₀ within the inertial range in ordinary hydrodynamic turbulence. Repeating Equation (A3) and using Kolmogorov's phenomenology,

$$\begin{equation} \delta l \sim v_{l} \delta t \sim (\epsilon l)^{1/3} \delta t, \end{equation} \tag{ A5 }$$

where is the power of injection (or dissipation) of the kinetic energy. Integrating the last expression,

$$\begin{equation} l^{2/3} - l_{0}^{2/3} = C_{R} \epsilon ^{1/3} (t - t_{0}), \end{equation} \tag{ A6 }$$

where the dimensionless constant C_R is the Richardson constant (Richardson 1926; Lesieur 1990). Therefore, C_dyn is related to C_R.

Using the dispersion σ of our experiments as the rms distance l in the Equation (A6), and supposing that this equation remains valid in MHD turbulence,

$$\begin{equation} \sigma ^{2/3} - \sigma _{0}^{2/3} = C_{*} \epsilon ^{1/3} (t - t_{0}), \end{equation} \tag{ A7 }$$

where C_* has not necessarily the same value as that of C_R.

Figure 13 shows the evolution of (σ^2/3 − σ^2/3₀) for our models. We have made an offset of the values for each model in order to make the visualization easier. Table 5 lists the approximate values of C_* for the models, obtained from a linear fit of the data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All in all, our results presented in Figure 13 are consistent with the scalings of diffusion shown in Cho et al. (2003a). This indicates that the relation (A7) continues to be approximately valid in the MHD case. The slopes do not seem to be very sensitive to M_s, but we can note a small difference for the trans-Alfvénic and super-Alfvénic cases, the constant C_* being a little higher in the super-Alfvénic cases. This finding is consistent with the prediction of changes of the diffusion efficiencies of Lazarian (2006). There it is shown that for higher magnetization, corresponding to an Alfvénic Mach number M_A = V/V_A < 1, where V is the injection velocity and V_A is the Alfvénic velocity, the turbulence at the large scales is weak, which produces weak mixing. In our parameterization this corresponds to the decrease of C_*, which is the case of A2 simulation. Interestingly enough, in the case of M_A < 1 and M_s ≫ 1, which is the case of A4 simulation, we see the increase of C_*. The case of high Mach number turbulence is not covered by Lazarian (2006), but we expect that the concentration of gas in dense regions plays an important role for the effect. In regions of gas concentration the Alfvénic velocity drops and the corresponding M_A increases (see Burkhart et al. 2009). As a result, the mixing should be less constrained by the magnetic field as we see in the A4 model. Comparing cases A1 and A3 we observe that for super-Alfvénic turbulence and different Mach number M_s the diffusivity is rather similar.

Note that Cho et al. (2003a) obtained the values of C_*/V_rms = 0.4 (V_A = 0.1 and M_s = 2.3), 0.3 (V_A = 1 and M_s = 2.3), 0.4 (V_A = 1 and M_s = 0.3), and 0.5 (hydrodynamic case and M_s = 0.3), which are consistent with our findings. In Cho et al. (2003a), V_rms ≈ 0.8 in all the cases (to make a comparison we had to convert the timescale to the present code units, i.e., 1 time step here = 2π time steps there). We conclude that, for a wide range of magnetizations and Mach numbers, turbulence is efficient in inducing turbulent diffusivity of heat.

Our results show that the turbulent diffusivity of passive scalar fields η_turb is well described by η_turb = C_*v_turbl_inj, similar to the hydrodynamic case, and that the coefficient C_* is not very sensitive to the strength of the magnetic field. Therefore, the magnetic field does not impose a strong suppression of the turbulent diffusion perpendicular to it as far as heat transfer is concerned.

As already pointed in Cho et al. (2003a; see also Lazarian 2006), the effectiveness of the turbulent diffusion has important consequences on the thermal diffusion in the ISM and intracluster medium (ICM). The nearly isotropic turbulent diffusivity can be 1 order of magnitude higher in the gas of central regions of galaxy clusters, like the Hydra A cluster, compared to the Spitzer value, based on the kinetic theory. In the ISM, in mixing layers, where turbulence is very strong, the turbulent diffusivity can be 2 orders of magnitude higher than the laminar values (see Esquivel et al. 2006). Our results extend the range of the applicability of the turbulent heat advection model of Cho et al. (2003a).

In our simulations described in Section 4 (see Table 1) we also keep track of passive scalar field representative either of heat or passive impurity.

Figure 14 (left and center panels) shows the time evolution of $\left\langle B_{z} \right\rangle _{0.25} / \left\langle \rho \right\rangle _{0.25} - \bar{B_{z}} / \bar{\rho }$ and $\left\langle \Phi \right\rangle _{0.25} / \left\langle \rho \right\rangle _{0.25} - \bar{\Phi } / \bar{\rho }$. Both quantities have a similar behavior, and all models seem to achieve the characteristic average values ($\bar{B_{z}}/\bar{\rho }$ and $\bar{\Phi }/\bar{\rho }$) roughly at the same time.

The right plot of Figure 14 is an analog of Figure 3. A comparison between these figures reinforces that the quantity 〈B_z〉_z/〈ρ〉_z evolves similar to 〈Φ〉_z/〈ρ〉_z. How can we understand this behavior?

Using the expression $\frac{\partial {\bf B}}{\partial t}=\nabla \times ({\bf v}\times {\bf B})$ one gets $\frac{\partial {\bf B}}{\partial t}= ({\bf B}\cdot \nabla){\bf v}-({\bf v}\cdot \nabla){\bf B}$, when ∇ · B = 0 is accounted for. Using the continuity equation in the form $\nabla \cdot {\bf v}= -\frac{1}{\rho }\frac{\partial \rho }{\partial t}-\frac{{\bf v}}{\rho }\cdot \nabla \rho$ one can get

$$\begin{equation} \left(\frac{\partial }{\partial t} + {\bf v}\cdot \nabla \right) \frac{\bf B}{\rho }=\frac{\bf B}{\rho } \cdot \nabla {\bf v} \end{equation} \tag{ B1 }$$

which coincides with Equation (A1) if Φ_α is substituted by B/ρ.

Figure 15 (left and center panels) shows the distribution of 〈ρ〉_z versus 〈Φ〉_z for the model B2 (B₀ = 1.0; see Table 1) at different time steps (t = 0, 10). The initial anti-correlation gives place to a tight correlation between the passive scalar and density. Right side of Figure 15 shows the evolution of 〈δΦ, δρ〉 (defined by Equation (8), where B must be substituted by Φ).

Zoom In Zoom Out Reset image size

Figure 16 shows the comparison of the evolution of 〈B_z〉_R=0.35L for models C1, C3, C4, C5, C6, and C7 of Table 2 with similar resistive models without turbulence of Table 4, with different values of Ohmic diffusivity η_Ohm. From these comparisons (and others where 〈B_z〉_R is measured for different values of R, like in Figure 10), we have estimated the values of η_turb shown in Table 2.

Zoom In Zoom Out Reset image size