MAGNETIC FIELD LINE RANDOM WALK IN ISOTROPIC TURBULENCE WITH ZERO MEAN FIELD

We consider the FLRW for homogeneous, isotropic turbulence with no mean magnetic field. With isotropy, we can assume the fluctuations in the x, y, and z directions to be statistically identical, and with no mean field, the diffusion coefficients in each direction are the same:

$$\begin{equation} D = D_x = D_y = D_z. \end{equation} \tag{ 13 }$$

When we write D or $\cal D$ with no functional dependence, we are referring to the asymptotic limit at large τ or s, respectively. The following derivation is similar to that of Ruffolo et al. (2004).

For zero mean field, we use B₀ = 0 in Equations (4). Then following Jokipii & Parker (1968), we express the displacement in, say, the x coordinate of a field line over τ by

$$\begin{equation} x = \int _0^{\tau }b_x[x(\tau ^{\prime }),y(\tau ^{\prime }),z(\tau ^{\prime })]d\tau ^{\prime }. \end{equation} \tag{ 14 }$$

The ensemble average variance is then given by

$$\begin{equation} \langle x^2\rangle = \int _0^{\tau } \int _0^{\tau } \langle b_x(x^{\prime },y^{\prime },z^{\prime })b_x(x,^{\prime \prime }y,^{\prime \prime }z^{\prime \prime })\rangle d\tau ^{\prime }d\tau,^{\prime \prime } \end{equation} \tag{ 15 }$$

where we introduce the notation x' for x(τ'), x'' for x(τ''), etc. Setting Δτ ≡ τ'' − τ', and with the assumption of homogeneity,

$$\begin{eqnarray} \langle x^2 \rangle &=& \int _0^{\tau }\negthinspace \negthinspace \int _{-\tau ^{\prime }}^{\tau - \tau ^{\prime }}\negthinspace \negthinspace {\cal L}_{xx}(\Delta \tau) d \Delta \tau d\tau ^{\prime }. \end{eqnarray} \tag{ 16 }$$

We use the symbol ${\cal L}_{xx}(\Delta \tau)$ to denote the ensemble average Lagrangian correlation of b_x at two locations along the same magnetic field line separated by Δτ.

The Lagrangian correlation function differs from a standard (Eulerian) correlation function in that the positions themselves depend on the field line trajectory. The Eulerian correlation function is the inverse Fourier transform of the power spectrum, which can often be specified based on physical considerations. It is possible to relate the Lagrangian correlation function ${\cal L}_{xx}$ to the Eulerian correlation function R_xx (and thence the power spectrum) using the following approximation:

$$\begin{equation} {\cal L}_{xx}(\Delta \tau) = \int R_{xx}({\bf \Delta x})P({\bf \Delta x}|\Delta \tau)d{\bf \Delta x}, \end{equation} \tag{ 17 }$$

where P(Δx|Δτ) is the distribution of the displacement $\bf \Delta x$ of field lines after Δτ. This approximation is known as Corrsin's independence hypothesis (Corrsin 1959; Salu & Montgomery 1977; see also McComb 1990). Computer simulations have been used to verify that theories based on Corrsin's hypothesis are successful for the FLRW (Gray et al. 1996; Ghilea et al. 2011) and field line separation (Ruffolo et al. 2004) in two-component turbulence and the FLRW in reduced magnetohydrodynamic (RMHD) turbulence (Snodin et al. 2013b), except that limiting quasi-2D cases can be problematic. The limitation of Corrsin's hypothesis is that it accounts for Lagrangian trajectories in terms of the two-point Eulerian correlation, so it does not account for coherent structures or for "memory" in the random walk, which occurs for nearly 2D turbulence because the FLRW can be temporarily trapped in "islands" with closed 2D trajectories (Ruffolo et al. 2003; Chuychai et al. 2007; Seripienlert et al. 2010).

We employ Corrsin's hypothesis and use P(Δx|Δτ) = P(Δx|Δτ)P(Δy|Δτ)P(Δz|Δτ) because of the statistical independence of Δx, Δy, and Δz, which is exact for our present case of isotropic turbulence and mirror symmetry, and also for some other symmetries. Then substituting Equation (17) into Equation (16) yields

$$\begin{eqnarray} \langle x^2\rangle &=& \int _0^{\tau }\negthinspace \negthinspace \int _{-\tau ^{\prime }}^{\tau - \tau ^{\prime }}\negthinspace \negthinspace \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } R_{xx}(\Delta x,\Delta y,\Delta z) \nonumber \\ && \!\!\!\!\times P(\Delta x|\Delta \tau)P(\Delta y|\Delta \tau)P(\Delta z|\Delta \tau)\, d\Delta x\, d\Delta y\, d\Delta z\, d \Delta \tau \, d\tau ^{\prime }. \nonumber\\ \end{eqnarray} \tag{ 18 }$$

There are different versions of the theory. Here we start by using diffusive decorrelation (DD; Taylor & McNamara 1971; Salu & Montgomery 1977; Matthaeus et al. 1995). The key assumptions of DD are that the magnetic field lines spread diffusively over the decorrelation scale of the random walk and their distributions are Gaussian. Thus the displacement distributions P of the DD model are Gaussian with variances

$$\begin{equation} \sigma _x^2 = \sigma _y^2 = \sigma _z^2 = 2D|\Delta \tau |, \end{equation} \tag{ 19 }$$

where D is the asymptotic diffusion coefficient in terms of τ (see Equations (6) and (13)).

Before delving deeper into the mathematical derivation, let us consider a simple conceptual derivation (adapted from Ruffolo et al. 2004):

$$\begin{equation} D = \frac{\langle x^2\rangle }{2\Delta \tau } \sim \left\langle \left(\frac{dx}{d\tau } \right)^2\,T \right\rangle = \langle b_x^2\,T\rangle, \end{equation} \tag{ 20 }$$

where T represents the τ scale over which b_x decorrelates and thus can be considered a "mean free τ." This FLRW is entirely "intrinsic," i.e., the decorrelation depends entirely on the FLRW itself. Suppose the decorrelation occurs over a spatial scale ℓ. Then for DD, we use

$$\begin{eqnarray} T &\sim & \frac{{\ell }^2}{2D} \nonumber \\ D &\sim & \frac{1}{2}\langle b_x^2\rangle \frac{{\ell }^2}{D} \nonumber \\ D &\sim & \frac{b}{\sqrt{6}}{\ell }, \end{eqnarray} \tag{ 21 }$$

where b² ≡ 〈|b|²〉. Note that in the above, T is expressed in terms of ℓ and D, which are constant for all field line trajectories in the ensemble. We obtain a diffusion coefficient in the form of Bohm diffusion as derived in Section 2 (see Equation (11)).

Similarly, for $\cal D$, the diffusion coefficient in terms of s,

$$\begin{equation} {\cal D}\sim \left\langle \left(\frac{dx}{ds} \right)^2\,S \right\rangle = \left\langle \frac{b_x^2}{|{\bf b}|^2}\,S\right\rangle, \end{equation} \tag{ 22 }$$

where S is a "mean free s." For DD, we use

$$\begin{eqnarray} S &\sim & \frac{{\ell }^2}{2{\cal D}} \nonumber \\ \cal D &\sim & \frac{1}{2} \left\langle \frac{b_x^2}{|{\bf b}|^2}\right\rangle \frac{{\ell }^2}{\cal D} \nonumber \\ {\cal D} &\sim & \frac{{\ell }}{\sqrt{6}}. \end{eqnarray} \tag{ 23 }$$

Then the ratio between the two diffusion coefficients is expected to be

$$\begin{equation} D/{\cal D}=b, \end{equation} \tag{ 24 }$$

the rms fluctuation, a prediction that is independent of the form of the power spectrum.

To continue the mathematical derivation of D, we make use of the power spectrum S_xx(k):

$$\begin{eqnarray} R_{xx}(\Delta x, \Delta y, \Delta z) &=& \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } \int _{-\infty }^{\infty } S_{xx}({\bf k}) \nonumber \\ && {} \times e^{ik_x\Delta x}e^{ik_y\Delta y}e^{ik_z\Delta z} dk_x dk_y dk_z. \end{eqnarray} \tag{ 25 }$$

Note that for Gaussian P,

$$\begin{eqnarray} \int _{-\infty }^{\infty } e^{ik_x\Delta x} P(\Delta x|\Delta \tau) d\Delta x &=& \int _{-\infty }^{\infty } \frac{e^{ik_x\Delta x}}{\sqrt{2\pi \sigma ^2}} \exp \left[-\frac{(\Delta x)^2}{2\sigma ^2}\right] d\Delta x \nonumber \\ &=& e^{-\frac{1}{2}\sigma ^2k_x^2} \nonumber \\ &=& e^{-D k_x^2|\Delta \tau |}, \end{eqnarray} \tag{ 26 }$$

which makes use of Equation (19), and with the analogous formulae for the Δy and Δz integrals, we obtain

$$\begin{equation} \langle x^2 \rangle = \int \int _0^{\tau } \int _{-\tau ^{\prime }}^{\tau - \tau ^{\prime }} S_{xx}({\bf k}) e^{-D k^2 |\Delta \tau |} d \Delta \tau d\tau ^{\prime }d{\bf k}. \end{equation} \tag{ 27 }$$

Note that S_xx(k) as used here follows a different Fourier transform convention than P_xx(k) as used in some of our previous work. The physical results do not depend on the convention used. Effectively, for the determination of the asymptotic diffusion coefficient at large τ, the limits of integration over Δτ in Equation (27) can be taken to be −∞ and ∞ (Jokipii & Parker 1968). In this diffusive regime, we have <x² > =2Dτ and

$$\begin{equation} D = \sqrt{\int \frac{S_{xx}(\bf k)}{k^2} d{\bf k}}. \end{equation} \tag{ 28 }$$

In terms of the modal energy spectral density S(k) = S_xx(k) + S_yy(k) + S_zz(k), which for isotropic turbulence is a function only of the wavevector magnitude k, we have

$$\begin{equation} D = \sqrt{\frac{1}{3}\int \frac{S(k)}{k^2} d{\bf k}}. \end{equation} \tag{ 29 }$$

In terms of the arclength, parameterized by ds = Bdτ, we can estimate the diffusion coefficient $\cal D$ using Equation (24):

$$\begin{equation} {\cal D} = \frac{D}{b} = \sqrt{\frac{1}{3}\frac{1}{b^2}\int \frac{S(k)}{k^2} d{\bf k}} = \sqrt{\frac{1}{3}\frac{\int S(k)/k^2 d{\bf k}}{\int S(k)d{\bf k}}}. \end{equation} \tag{ 30 }$$

These diffusion coefficients are quite similar to that obtained by Matthaeus et al. (1995) for a 2D axisymmetric component of transverse turbulence with a mean field B₀:

$$\begin{equation} D_{2D}\equiv \frac{\left< x^2 \right>}{2z} = \sqrt{\frac{1}{2} \frac{1}{B_0^2} \int \frac{S(k_\perp)}{k_\perp ^2} d{\bf k}_\perp }, \end{equation} \tag{ 31 }$$

where z serves to parameterize the FLRW, k_⊥ ≡ (k_x, k_y), and $k_\perp ^2=k_x^2+k_y^2$. The only differences in these expressions are due to the different dimensionality n (n = 3 for isotropic 3D turbulence, n = 2 for 2D turbulence, which appears in the factor 1/n inside the square root), and the magnetic field component associated with the independent variable (not used for isotropic turbulence in terms of τ; b for isotropic turbulence in terms of s; and B₀ for 2D turbulence in terms of z). The reason for the similar expressions is that both diffusion processes are intrinsic, involving random flights in which the magnetic fluctuation along the trajectory decorrelates stochastically due to the random motion itself (not due to B₀).

Matthaeus et al. (1995) pointed out that the integral in Equation (31) introduces a new length scale, the ultrascale $\tilde{\lambda }$. That equation can be expressed as (Ruffolo et al. 2004)

$$\begin{equation} D_{2D} = \frac{\tilde{\lambda }}{\sqrt{2}}\frac{b}{B_0} \end{equation} \tag{ 32 }$$

with

$$\begin{equation} \tilde{\lambda }\equiv \sqrt{\frac{\int S(k)/k^2 d{\bf k}_\perp }{\int S(k)d{\bf k}_\perp }}. \end{equation} \tag{ 33 }$$

These concepts can be generalized to 3D turbulence. For isotropic turbulence, the diffusion coefficients in Equations (28) and (30) can be expressed as

$$\begin{equation} D=\tilde{\lambda }\sqrt{\langle b_x^2\rangle } = \frac{\tilde{\lambda }}{\sqrt{3}}b \quad {\rm and}\quad {\cal D}=\frac{\tilde{\lambda }}{\sqrt{3}}, \end{equation} \tag{ 34 }$$

where

$$\begin{equation} \tilde{\lambda }\equiv \sqrt{\frac{\int S(k)/k^2 d{\bf k}}{\int S(k)d{\bf k}}}, \end{equation} \tag{ 35 }$$

i.e., $\tilde{\lambda }^2$ is the k⁻² moment of the power spectrum. This form of Bohm diffusion is similar to that anticipated from our heuristic derivation (Equation (21)).

In the random ballistic decorrelation (RBD) version of the theory, instead of the assumption of diffusive spreading as in Equation (19), we assume that the magnetic field lines ballistically spread in random directions over the decorrelation distance. Adapted from Ghilea et al. (2011),

$$\begin{eqnarray} &&\sigma _x^2 = \sigma _y^2 = \sigma _z^2 = {\left\langle {b_x^2} \right\rangle } \Delta {\tau ^2}. \end{eqnarray} \tag{ 36 }$$

In the conceptual derivation, T depends on the magnitude of the local field b as

$$\begin{eqnarray*} T &\sim & \frac{{\ell }}{|{\bf b}|}, \nonumber \end{eqnarray*}$$

so

$$\begin{eqnarray} D &\sim & \left< \left(\frac{dx}{d\tau } \right)^2 T \right> \sim \frac{1}{3}\left\langle \frac{|{\bf b}|^2}{|{\bf b}|}\right\rangle \,{\ell } = \frac{1}{3}\langle |{\bf b}|\rangle \,{\ell }. \end{eqnarray} \tag{ 37 }$$

Since S is determined for a ballistic trajectory over a distance ℓ, we simply use

$$\begin{eqnarray} S &\sim & {\ell } \nonumber \\ \cal D &\sim & \left< \left(\frac{dx}{d s} \right)^2 S \right> \sim \frac{{\ell }}{3}. \end{eqnarray} \tag{ 38 }$$

Then the ratio $D/{\cal D}$ for RBD is expected to be

$$\begin{equation} \frac{D}{\cal D} \sim \langle |{\bf b}|\rangle. \end{equation} \tag{ 39 }$$

The ratio $D/{\cal D}$ can help distinguish between DD and RBD, as will be discussed later in this section.

Let us turn to a full mathematical derivation of the asymptotic diffusion coefficient for RBD. With the RBD variance from Equation (36), Equation (26) becomes

$$\begin{eqnarray} \int \nolimits _{ - \infty }^\infty {{e^{ i{k_x}\Delta x}}P\left({\Delta x\left| {\Delta \tau } \right.} \right)} d\Delta x &=& \int \nolimits _{ - \infty }^\infty {\frac{{{e^{ i{k_x}\Delta x}}}}{{\sqrt{2\pi \left\langle {b_x^2} \right\rangle \Delta \tau {^2}} }}}\nonumber\\ && \times \exp {\left[ {\frac{{ - {{\left({\Delta x} \right)}^2}}}{{2\left\langle {b_x^2} \right\rangle \Delta \tau {^2}}}} \right]} d\Delta x\nonumber \\ &=&{e^{ - \frac{1}{2}\left\langle {b_x^2} \right\rangle k_x^2\Delta \tau {^2}}}. \end{eqnarray} \tag{ 40 }$$

Using the analogous results for the Δy and Δz integrals and $\left\langle {b_x^2} \right\rangle = \langle b_y^2\rangle = \langle b_z^2\rangle = b^2/3$ from isotropy, we obtain 〈x²〉 for large τ as

$$\begin{eqnarray} \left\langle { {x^2}} \right\rangle &=& \int {\left[ {\int _0^{ \tau } {\int _{ - \infty }^{\infty } {{S_{xx}}(\mathbf {k} }){e^{ - \frac{1}{6}{b^2}{k^2}\Delta \tau {^2}}}d\Delta \tau d\tau ^{\prime }} } \right]} d\mathbf {k} \nonumber \\ &=&\frac{\sqrt{6\pi } \tau }{b} \int {\frac{{ {S_{xx}}(\mathbf {k}) }}{{k}}}d\mathbf {k}. \end{eqnarray} \tag{ 41 }$$

Now we can obtain the random ballistic diffusion coefficient as

$$\begin{eqnarray} D = {D_x} &\equiv & \frac{{\left\langle { {x^2}} \right\rangle }}{{2 \tau }} \nonumber \\ &=& \frac{1}{b}\sqrt{\frac{\pi }{{6 }}}\int {\frac{{ {S}(k)}}{{k}}} {d}\mathbf {k} \nonumber \\ &=&b\sqrt{\frac{\pi }{6}} \frac{\int _0^{\infty } k S(k) dk}{\int _0^{\infty } k^2 S(k) dk}. \end{eqnarray} \tag{ 42 }$$

The integrals in Equation (42) are related to the total correlation length λ_c. To see this, first note that for an isotropic field, we can consider the total correlation function R(r) = 〈b(r₀) · b(r₀ + r)〉, where r₀ is an arbitrary position. The correlation length is then

$$\begin{equation} \lambda _c = \frac{\int _0^{\infty } R(r) dr}{R(0)}. \end{equation} \tag{ 43 }$$

If we define R(r) in terms of its Fourier transform

$$\begin{equation} R(r) = 4 \pi \int _{0}^{\infty } \frac{\sin (kr)}{kr} k^2 S(k) dk, \end{equation} \tag{ 44 }$$

then we get R(0) = 4π∫k²S(k)dk and $\int _0^{\infty } R(r) dr = 2 \pi ^2 \int _0^{\infty } k S(k) dk$, which when combined give

$$\begin{equation} \lambda _c = \frac{\pi }{2} \frac{\int _0^{\infty } k S(k) dk}{\int _0^{\infty } k^2 S(k) dk}. \end{equation} \tag{ 45 }$$

A similar result for isotropic 2D turbulence was presented by Matthaeus et al. (2007). Thus we obtain

$$\begin{eqnarray} D&=&\sqrt{\frac{2}{3 \pi }} \lambda _c b. \end{eqnarray} \tag{ 46 }$$

Guided by the heuristic result in Equation (37), we consider how to express Equation (46) in terms of 〈|b|〉. Assuming that the local distribution in b is Gaussian with width σ in each component, which is the case in our simulations, we have

$$\begin{eqnarray} \langle |{\bf b}|\rangle &=& \frac{\int _0^\infty b^3 e^{-b^2/(2\sigma ^2)}db}{\int _0^\infty b^2 e^{-b^2/(2\sigma ^2)}db} =\frac{2\sigma ^4}{(\sqrt{\pi /2})\sigma ^3} \nonumber \\ &=& (2\sqrt{2/\pi })\sigma. \end{eqnarray} \tag{ 47 }$$

Note that b² = 3σ², so $\sigma = b/\sqrt{3}$. Then

$$\begin{eqnarray} \langle |{\bf b}|\rangle &=& (2\sqrt{2/\pi })\frac{b}{\sqrt{3}}. \end{eqnarray} \tag{ 48 }$$

Thus we can rewrite Equation (46) in terms of 〈|b|〉 as

$$\begin{eqnarray} D&=&{\frac{ \lambda _c}{2}} \langle |{\bf b}|\rangle, \end{eqnarray} \tag{ 49 }$$

an expression of Bohm diffusion that is similar to the result of our heuristic derivation (Equation (37)).

Finally, we note that DD and RBD give different predictions for $D/{\cal D}$. For DD,

$$\begin{equation} \frac{D}{\cal D} \sim b,\quad {\rm (DD)} \end{equation} \tag{ 50 }$$

and for RBD,

$$\begin{equation} \frac{D}{\cal D} \sim \langle |{\bf b}|\rangle = 0.9213\, b\quad {\rm (RBD).} \end{equation} \tag{ 51 }$$

Thus in addition to comparing D versus τ between simulations and different versions of the theory, we can compare the ratio of asymptotic values $D/{\cal D}$ between the simulations and DD and RBD predictions.

In the previous two subsections, we determined the magnetic field line diffusion coefficients in the asymptotic limit τ → ∞. Here we consider the evolution of diffusion coefficients with τ by using an ordinary differential equation (ODE). To formulate the ODE, we start by considering the running diffusion coefficient, defined by

$$\begin{eqnarray} &&D(\tau) =\frac{1}{2} \frac{{d\left\langle { {x^2}} \right\rangle }}{{d\tau }}. \end{eqnarray} \tag{ 52 }$$

Substituting Equation (16) into Equation (52),

$$\begin{eqnarray} &&\frac{1}{2} \frac{{d\left\langle { {x^2}} \right\rangle }}{{d\tau }}= D(\tau) = \int \nolimits _0^\tau \mathcal {L}_{xx }\left(\Delta \tau \right) d \Delta \tau. \end{eqnarray} \tag{ 53 }$$

We define V(τ) = 〈x²〉, and the initial condition for Equation (53) is V(0) = 0. Then, differentiating this equation, we obtain

$$\begin{eqnarray} &&\frac{{{d^2}V(\tau)}}{{d{\tau ^2}}} = 2\frac{{dD(\tau)}}{{d\tau }} = 2 \mathcal {L}_{xx }\left(\tau \right), \end{eqnarray} \tag{ 54 }$$

with the initial condition D(0) = 0. Assuming Corrsin's hypothesis and a Gaussian displacement distribution with variance σ² along each coordinate, the ODE (Equation (54)) becomes

$$\begin{eqnarray} \frac{{{d^2}V(\tau)}}{{d{\tau ^2}}} = 2\frac{{dD(\tau)}}{{d\tau }} &=& 2\int S_{xx}(\mathbf {k})e^{- \frac{k^2\sigma ^2(\tau)}{2}} d\mathbf {k}. \end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray} &=& 2\int \frac{S(k)}{3} e^{ - \frac{k^2 \sigma ^2(\tau)}{2}} d\mathbf {k}. \end{eqnarray} \tag{ 56 }$$

In order to determine the FLRW evolution, one can calculate the running diffusion coefficient D(τ) by integrating Equation (56) from 0 to τ. The asymptotic diffusion coefficient D is obtained by integrating to ∞. For the DD model, one can start by setting σ²(τ) = 2Dτ in Equation (56), using the asymptotic diffusion coefficient D, where τ is assumed to be positive. Then we can determine D by integrating Equation (56) over τ from 0 to ∞ and solving the resulting implicit equation for D, as in Section 3.2. After determining D, one can obtain the evolution of the DD running diffusion coefficient by integrating Equation (56) from 0 to any τ. For the RBD model, we substitute $\sigma ^2(\tau)={\left\langle {b_x^2} \right\rangle }{\tau ^2}$ from Equation (36) into Equation (56). Then we can determine the RBD running diffusion coefficient and RBD asymptotic diffusion coefficient by integrating Equation (56). In addition to the DD and RBD models, in Equation (56), we can identify σ²(τ) = V(τ), which provides self-closure of the ODE. We refer to this as the "ODE model," which was suggested by Saffman (1963) and Taylor & McNamara (1971), and is sometimes simply referred to in the literature as the "Corrsin approximation." We obtain a result equivalent to those obtained by Lundgren & Pointin (1976), Wang et al. (1995), Shalchi & Kourakis (2007), and Snodin et al. (2013a):

$$\begin{eqnarray} \frac{{{d^2}V(\tau)}}{{d{\tau ^2}}} &=& 2\int \frac{{S(k)}}{3} e^{{ - \frac{{{k^2}{V(\tau)}}}{2}}} d\mathbf {k}. \end{eqnarray} \tag{ 57 }$$

Using the definition of the diffusion coefficient in Equation (52), Equation (57) can be rewritten as a system of two first-order ODEs:

$$\begin{eqnarray} \frac{{dV(\tau)}}{{d\tau }} &=& 2D(\tau) \end{eqnarray} \tag{ 58 }$$

$$\begin{eqnarray} \frac{{{d}D(\tau)}}{{d{\tau }}} &=& \int \frac{{S(k)}}{3} e^{{ - \frac{{{k^2}{V(\tau)}}}{2}}} d\mathbf {k}. \end{eqnarray} \tag{ 59 }$$

Solving Equations (58) and (59), subject to the initial conditions V(0) = 0 and D(0) = 0, the evolution of the ODE running diffusion coefficient is obtained. At low τ, we obtain V(τ)∝τ², which is consistent with ballistic field line trajectories. One can obtain the asymptotic diffusion coefficient via a first integral, following Taylor & McNamara (1971). From Equation (57),

$$\begin{eqnarray*} \frac{ { d }^{ 2 }V(\tau) }{ d{ \tau }^{ 2 } } &=&2\frac{ 4\pi }{ 3 } \int _{0}^{\infty } { { k }^{ 2 }S(k){ e^{ -\frac{ { { k^{ 2 } }V (\tau)} }{ 2 } } } } dk.\nonumber \end{eqnarray*}$$

Multiplying by dV(τ)/dτ, we obtain

$$\begin{eqnarray*} \frac{ { d } }{ d{ \tau } } \left[ \frac{ 1 }{ 2 } { \left(\frac{ { dV(\tau) } }{ d{ \tau } } \right) }^{ 2 } \right] &=&\frac{ { d } }{ d{ \tau } } \left[ -\frac{ 16\pi }{ 3 } \int _{0}^{\infty } { S(k){ e^{ -\frac{ { { k^{ 2 } }V (\tau) } }{ 2 } } } } dk \right] \nonumber \\ \frac{ 1 }{ 2 } { \left(\frac{ { dV(\tau) } }{ d{ \tau } } \right) }^{ 2 } = 2{ D }^{ 2 }\left(\tau \right) &=&C-\frac{ 16\pi }{ 3 } \int _{0}^{\infty } { S(k){ e^{ -\frac{ { { k^{ 2 } }V (\tau) } }{ 2 } } } } dk.\nonumber \end{eqnarray*}$$

Let us consider the difference between D²(τ = ∞) and D²(τ = 0):

$$\begin{eqnarray} { D}^{ 2 }\left(\infty \right) -{ D }^{ 2 }\left(0 \right) &=&\frac{ 8\pi }{ 3 } \int _{0}^{\infty } { S(k) } dk\nonumber \\ { D }^{ 2 }\left(\infty \right)& =&\frac{ 2 }{ 3 }{ b }^{ 2 } \frac{ \int _{0}^{\infty } { S(k) } dk }{ \int _{0}^{\infty } { { k }^{ 2 }S(k) } { d }k } \nonumber \\ { D }\left(\infty \right) &=&\sqrt{\frac{ 2 }{ 3 }} { b }\tilde{ \lambda }. \end{eqnarray} \tag{ 60 }$$

We note that this diffusion coefficient is a factor of $\sqrt{2}$ greater than the one obtained for DD. This factor of $\sqrt{2}$ was also obtained by Snodin et al. (2013a) and Ruffolo & Matthaeus (2013) for 2D or quasi-2D axisymmetric transverse turbulence.

To test the validity of our analytic models of magnetic field line diffusion for isotropic turbulence with zero mean field, we performed direct computer simulations to calculate the running diffusion coefficients without using any of the assumptions of the analytic models. The isotropic 3D fluctuating magnetic field is generated on a regular Fourier grid in wave number space (k space). The solenoidal property ∇ · b = 0 in real space implies k · b(k) in wave number space. Therefore, the component b(k) at a given k is given by

$$\begin{equation} {\mathbf {b}}({\mathbf {k}}) = \frac{1}{2}\sqrt{S(k)}\left[{\mathbf {b}}_1({\mathbf {k}})e^{{\rm i}\phi _1({\mathbf {k}})} + {\rm i} {\mathbf {b}}_2({\mathbf {k}})e^{{\rm i}\phi _2({\mathbf {k}})}\right], \end{equation} \tag{ 61 }$$

where b₁ and b₂ are unit vectors perpendicular to each other and also to k, which ensures that the constructed field is divergence free. Here ϕ₁ and ϕ₂ are independent random phases at every k for these two polarization directions and S(k) is the power spectrum, given by

$$\begin{eqnarray} &&S(k) = C\frac{{{k^2}{\lambda ^2}}}{{{{\left({1 + {k^2}{\lambda ^2}} \right)}^{{\Gamma /2+2}}}}}, \end{eqnarray} \tag{ 62 }$$

with a normalization constant C, used to control the magnetic field strength b, and where λ is the bendover scale. Note that for 3D isotropic fluctuations at a given magnitude k, the omnidirectional power spectrum is E(k) = 4πk²S(k), which at large k satisfies the scaling E(k)∝k^−Γ. For example, when we set Γ = 5/3 we obtain the Kolmogorov scaling E(k)∝k^−5/3. At low k, Equation (62) yields S(k)∝k² as required by homogeneity (Batchelor 1970). To yield a good resolution of the discrete power spectrum, we use a zero padding technique, setting the power to zero for all k above half of the Nyquist value. Taking an inverse fast Fourier transform of b(k) yields the real space magnetic field representation. To simulate the power spectrum covering the energy containing range and the inertial range of turbulence with zero padding, we generate the 3D isotropic field for b = 1 and λ = 1 on a periodic box of size L_x = L_y = L_z = 100 with N_x × N_y × N_z = 1024 × 1024 × 1024 grid points. In this work it is interesting to note that simulation results for diffusion coefficients are essentially unchanged (at least for Γ = 5/3) if we simply use ${\bf b}({\bf k}) = \sqrt{S(k)} {\bf b}_1({\bf k})$, instead of the expression in Equation (61).

The above random prescription allows us to generate as many different realizations of a magnetic field with a given spectrum as we wish. In each realization, we trace multiple magnetic field lines, where each line is obtained by solving the system of magnetic field line equations from a random initial location, using Equations (4) up to τ = 125 and Equations (2) up to s = 125. We solve the magnetic field line equations numerically using a fifth-order Runge–Kutta method with adaptive time stepping regulated by a fourth-order error estimate step. We traced 500,000 field lines, using 12,500 lines per field representation, to help ensure sufficient sampling. Collecting displacements for all pairs of points along each field line, the running diffusion coefficients were calculated.

The theories presented earlier predict $D\propto b\tilde{\lambda }$ or bλ_c, and ${\cal D}\propto \tilde{\lambda }$ or λ_c. When deriving $\tilde{\lambda }$ and λ_c from the above spectrum of Equation (62), they are found to be proportional to the bendover scale λ. Therefore, one expects that the diffusion coefficients will scale proportionally with λ. As a basic test of our code, we performed simulations with various values of λ, and confirmed that D∝bλ and ${\cal D}\propto \lambda$.

We explore the effect of the power spectrum on the FLRW by considering three values of Γ corresponding to different models of magnetohydrodynamic turbulence phenomenology (for a review, see Zhou et al. 2004). These are Iroshnikov–Kraichnan scaling with Γ = 3/2 (Iroshnikov 1963; Kraichnan 1965), Kolmogorov scaling with Γ = 5/3 (Kolmogorov 1941; Batchelor 1970), and "weak turbulence" scaling with Γ = 2 (Ng & Bhattacharjee 1997; Galtier et al. 2000). The value of Γ is used to specify the power spectrum of fluctuations according to Equation (62), both in the computer simulations and for evaluation of the analytic theory.

After tracing 500,000 field lines for b = 1 and λ = 1, the running diffusion coefficients D_i(τ) and ${\mathcal {D}}_i(s)$ were determined, as shown in Figure 2 for the example case of Γ = 5/3. The simulations are sufficiently accurate that the estimates of diffusion coefficients in each direction are nearly coincident, according to the symmetry of the problem, indicating that a sufficient number of field lines and realizations of turbulence have been simulated. At large values of τ or s, the diffusion coefficients approach asymptotic values. For Γ = 3/2, 5/3, and 2, Table 1 shows the asymptotic diffusion coefficients D and $\mathcal {D}$, which were calculated by averaging D_i(τ) and $\mathcal {D}_i(s)$ for i = x, y, and z, weighting according to the number of available pairs of points, over the interval 60 < τ < 80 or 60 < s < 80. The error of the mean, indicated by parentheses, was determined based on the spread among the values along the three directions. Table 1 also shows the diffusion coefficients calculated for the three versions of the theory. For a more precise comparison with computer simulation results, we evaluated the integrals in Equation (33) over a finite k  range, from the minimum to the maximum k values in the Fourier grid used in the simulations. Recall that all three versions of the theory predict Bohm diffusion (as generalized to the case of zero mean field in Section 2) as is usual for theories using Corrsin hypothesis at high Kubo number R. The computer simulation results confirm that for this case of R = ∞, the asymptotic diffusion coefficients are consistent with Bohm diffusion (the diffusion coefficients are finite). The simulation results are also numerically close to the results of all three versions of the Corrsin-based theory. We can see that the ratios between all asymptotic diffusion coefficients are very similar when changing the power spectral index in Equation (62) to change the omnidirectional power spectrum at high k from k^−5/3 to k^−3/2 or k⁻², indicating that our results are not sensitive to such details of the power spectrum.

Zoom In Zoom Out Reset image size

Table 1. Asymptotic Field Line Diffusion Coefficients and Their Ratios (b = 1, λ = 1)

Quantity	Parameters	Evaluation	Value
			Γ = 3/2	Γ = 5/3	Γ = 2
D	τ	Theory:DD	0.2881	0.3094	0.3522
		Theory:RBD	0.2680	0.2916	0.3405
		Theory:ODE	0.4075	0.4376	0.4981
		Simulation	0.3357(9)	0.3600(5)	0.4069(7)
$\mathcal {D}$	s	Simulation	0.3759(1)	0.3998(5)	0.4517(9)
D/$\mathcal {D}$	τ and s	Theory:DD	1.0000	1.0000	1.0000
		Theory:RBD	0.9213	0.9213	0.9213
		Simulation	0.8929(24)	0.9005(17)	0.9008(23)

Download table as:  ASCIITypeset image

Let us now consider whether the simulation results favor a particular version of the theory. Among the DD, RBD, and ODE versions, the DD result is closest to the simulation result (DD is 14% lower). However, the evolution of running diffusion coefficients as a function of τ as shown in Figure 3 indicates that RBD is more accurate than DD in comparison with computer simulations at small τ and thus better tracks the transition from free streaming to asymptotic diffusion. Moreover, the asymptotic diffusion coefficient predicted by RBD is only slightly lower than the DD result and is 18% lower than the simulation result. In addition to comparing theory and simulation results for D(τ), we also compare the ratio of asymptotic values $D/\mathcal {D}$ between simulations and DD and RBD models. It came as a surprise that tracing the field line diffusion in τ and s at the same time can shed some light on the physical processes. After all, τ and s provide two parameterizations of the same field line trajectories. However, in the RBD calculation of diffusion coefficients, the mean free value T weights the b ensemble differently than the mean free value S, while in the DD calculation the weights are independent of the local field, so the comparison between D and $\mathcal {D}$ helps indicate the dominant FLRW mechanism. The $D/\mathcal {D}$ ratios from simulations are close to 0.90, and are much closer to the RBD prediction of 0.9213 than to the DD prediction of 1. This may indicate that the RBD version best captures the physics of the FLRW in addition to being physically simplest and easiest to implement.

Zoom In Zoom Out Reset image size