Anomalous diffusion of field lines and charged particles in Arnold-Beltrami-Childress force-free magnetic fields

The possibility that the cosmic magnetic fields are force-free was initially pointed out by Lüst and Schlüter.¹ In plasma regions of high conductivity, large current flows give rise to Lorentz forces that cannot be balanced by the observed pressure gradients, or gravitational forces, or inertial forces.^2,3 Consequently, the Lorentz force has to vanish, which constraints the currents to flow parallel to the magnetic field. In these plasmas, the kinetic energy density of the ionized gas is much less than the energy density in the magnetic fields.

In a plasma, the force-free magnetic field (or magnetic induction) obeys

where B = μ₀H is the magnetic field, H is the magnetic intensity, μ₀ is the free-space permeability, and λ, in general, is an arbitrary function of position. Since ∇.B = 0, it follows that B$⋅$∇λ = 0 for force-free fields.

An equation similar to (1) arises for steady state, incompressible Euler fluid flows

where u is the fluid velocity, $u=|u|$⁠, p is the pressure, and ρ is the (constant) fluid density. Then,

where $H(r)=(p/ρ)+(u2/2)$ is the Bernoulli function. Consequently, the flow lines of both u and the vorticity ∇ × u lie in the surfaces of constant $H(r)$⁠. This requires that such surfaces exist, and that the function $H(r)$ itself is not a constant independent of r. In two seminal papers,^4,5 Arnold proved that the three-dimensional, incompressible, steady state Euler flows are integrable except when the vorticity ∇ × u is parallel to u, i.e.,

The above equations represent Beltrami flows.⁶ While λ can, in general, be a function of r, Arnold^5,7 noted that the Beltrami flows are non-integrable only when λ is a constant. It follows that the force-free magnetic field equation (1), being identical to the Beltrami flow (4), will have non-integrable magnetic field lines when λ is a constant independent of space.

Based on an earlier study by Chandrasekhar and Woltjer³ on force-free magnetic fields in cosmic plasmas, Woltjer⁸ proved that the state with constant λ is the lowest magnetic energy state in a closed system.

There are several different astrophysical plasma environments where force-free magnetic fields are believed to be present. These include the solar photosphere⁹ and the solar corona.^10,11 The force-free plasma state has also been studied within the context of laboratory plasmas.^12,13

In this paper, we study the nonlinear spatial evolution of the force-free magnetic field lines and the spatial transport of charged particles moving in these fields. The form of the force-free magnetic field chosen for these studies satisfies (1) for constant λ. It is identical to the velocity field which was first proposed by Arnold⁴ in the context of Euler flows. This particular vector field is referred to as the Arnold-Beltrami-Childress (ABC) flow,¹⁴ or, in our case, the ABC force-free magnetic field.¹⁵ The ABC flow is helical and the properties of the Lagrangian structure of these flows have been extensively studied in Refs. 14 and 16. Such flows in conducting fluids can lead to the kinematic dynamo action responsible for generating astrophysical magnetic fields.^17–20 The most common examples being galactic magnetic fields, solar magnetic fields, and planetary magnetic fields. Observations of solar magnetic flares point to the helical nature of the magnetic field ejected from the Sun and carried away by the solar wind.²¹ 

The force-free magnetic fields also affect the motion of charged particles. The transport of cosmic rays across the magnetic field is believed to be responsible for regulating the high-energy particle density in galactic magnetic fields.^22,23 The propagation of solar cosmic rays in the interplanetary space cannot be explained unless there is transport of energetic particles across the interplanetary magnetic field.²⁴ Since these plasmas are essentially collisionless, the existence of chaotic magnetic fields is postulated as a means for explaining the diffusion of charged particles.^22,25

In view of the omnipresence of helical force-free magnetic fields in plasmas, a study of the properties of field lines, and of spatial diffusion of charged particles in these fields, is of fundamental importance.

The magnetic field line equations can be formulated as a Hamiltonian system.²⁶ In three spatial dimensions, the field line equations are non-integrable.²⁷ For such Hamiltonian systems, the spatial structure of the field lines is a mix of regular and chaotic trajectories.²⁸ The phase space structure of ABC flows is indeed such a blend of flow lines.¹⁴ For the ABC magnetic field, we quantify the spatial divergence of nearby field lines in terms of Lyapunov exponents (LEs).²⁸ Only the largest LE is analyzed as it represents the dominant growth rate, and the orientations of the other LEs tend toward the direction of the largest LE as the trajectory evolves over longer distances.

The nonlinear evolution of the ABC field lines as a function of s, the distance along the field line, is analyzed by evaluating the finite-distance LEs^29,30 and the asymptotic-distance LEs.^31,32 The analogy here is with finite-time and time-asymptotic LEs for systems in which the evolution parameter is time. The finite-distance LE gives a measure of the local nonlinearity, while the asymptotic-distance LE is a measure of the intensity of chaos. For the magnetic field lines, the contours of constant finite-distance LEs do not line up with the contours of constant asymptotic-distance LEs. The former are aligned more closely with the contours of constant magnetic field amplitude, while the latter are aligned more closely with the Poincaré surface-of-section. In regions where the field lines are chaotic, the spatial diffusion of field lines is determined by considering an ensemble of initial points that are located within a close vicinity of a reference point. The field lines originating from all these points are followed along s, and the ensemble-averaged variance with respect to the reference field line is evaluated. There are two stages to the spatial evolution of the variance. The first stage corresponds to mixing of the field lines in which the ensemble-averaged mean position of the field lines is anisotropic. In the second stage, the ensemble-averaged mean position evolves isotropically, and the variance follows a power law in s. The spatial diffusion of the field lines, which follows directly from the variance in the second stage, scales like s^α with α > 1; i.e., the field lines are super-diffusive in space.³³ It is in this region of chaotic, and diffusive, field lines that the charged particles are also chaotic and diffusive. And, not surprisingly, the diffusion of particles is anomalous; in fact, the particles experience superdiffusion.

The analysis of the particle motion follows that of the field lines. The nonlinear evolution of particles is determined by evaluating the finite-time and time-asymptotic LEs. Again, we concentrate only on the largest Lyapunov exponent. In contrast to the nonlinear evolution of the field lines, for the particles, the finite-time LEs not only depend on the initial spatial location of the particles but also on their energy and their initial pitch-angle. The pitch-angle is defined as the angle between the velocity vector and the magnetic field vector at the spatial location of the particle. As in the case of field lines, the contours of finite-time LEs are aligned more closely with contours of constant magnetic field amplitude when the pitch-angle is zero; i.e., the initial velocity vector is aligned with the magnetic field line. This is not the case when the initial pitch-angle is different from zero, and the particles drift due to the spatial inhomogeneity of the magnetic field.³⁴ However, importantly, the time-asymptotic LEs are independent of the pitch-angle, but depend on the initial spatial location of the particles. The largest LEs are confined to the same regions where the field lines are chaotic. For an ensemble of particles initially located in the chaotic phase space, the ensemble-averaged mean square displacement is evaluated as a function of the time lag. There are three stages in the evolution of the mean square displacement. The first stage is the ballistic stage where the particles are, essentially, free streaming. During the second stage, the ensemble-averaged velocity autocorrelation function decays to near zero. In the third stage, the evolution of the particles, as a function of the time lag, is super-diffusive.³³ 

Diffusion of magnetic field lines has been considered in various studies of highly turbulent plasmas.^35,36 In these studies, as the distance along the magnetic field lines increases, the separation between field lines follows ordinary diffusion and scales linearly with s. On the contrary, in this paper, the diffusion of field lines for the completely deterministic ABC force-free magnetic field is studied without imposing any ad hoc stochastic fields. For s larger than the correlation length, the separation between field lines is found to be super-diffusive.

The presumption of stochastic, or chaotic, magnetic fields in the previous studies has been necessary in order to explain the observed transport of particles in various astrophysical plasmas.^22,23,25 The particle diffusion is ordinary in the sense that the mean square displacement scales linearly as a function of the time lag.²³ In this paper, it is shown that there is an inherent spatial diffusion or transport of charged particles in the three-dimensional, deterministic, ABC magnetic fields in regions where the field lines are chaotic. However, the diffusion is not ordinary—the particles, just like the field lines, undergo superdiffusion.

The rest of the paper is divided into two main parts. The first part is on the properties of the Lyapunov exponents for the ABC magnetic field lines. It also includes a description of the spatial diffusion of field lines. The second part is on the dynamics of charged particles in the ABC fields. This part includes sections on the finite-time and time-asymptotic Lyapunov exponents, the velocity autocorrelation function, the mean square displacement, and the spatial diffusion of particles in the chaotic ABC magnetic fields.

The general solution of the force-free equation (1) for constant λ has been discussed in Ref. 37. Of the large class and variety of solutions of (1), we will choose a particular form of the three-dimensional, force-free, helical magnetic field—the ABC field—that was initially suggested by Arnold⁴ and shown to have chaotic field lines.^14,38

The ABC field is expressed in the Cartesian coordinate system and takes on the form

where A, B, and C are arbitrary constants, and the unit vectors $x̂, ŷ, ẑ$⁠, are along the corresponding axes. The ABC field is spatially periodic in the three principal directions.

The magnetic field line trajectories are obtained from

where $r=xx̂+yŷ+zẑ, B=Bxx̂+Byŷ+Bzẑ, |B|$ is the magnitude of B, and s is the path length along the magnetic field line.

We normalize distances by the constant λ and B by a constant magnetic field B₀. λ⁻¹ basically determines the spatial scale length for the variation of the magnetic field. The equations are expressed in a normalized form by replacing (λx, λy, λz) and (A/B₀, B/B₀, C/B₀) by the dimensionless variables (x, y, z) and (A, B, C), respectively. Figure 1 shows the z = 0 Poincaré surface-of-section for the magnetic field lines in the x–y plane. The parameters A, B, and C are the same as used in Ref. 14, and this figure is similar to Fig. 5 in Ref. 14. The intermixing of chaotic and regular field lines is quite evident. The sparsely populated regions in Fig. 1 are along the contour $Bz≡C sin y+B cos x=0$⁠. From (5), the magnetic field lines along this contour are tangential to the x–y plane. Since the field is force-free, the magnetic field vorticity is also tangential to the x–y plane along this contour.

Superimposed on the z = 0 Poincaré surface-of-section, in Fig. 1, are contours of constant magnitude of the magnetic field $|B|$⁠. While the field lines can wander chaotically in space, the contours of constant $|B|$ are well-defined and fixed. These contours can be useful in understanding some of the kinematics of charged particles in the ABC magnetic fields. The structure of field lines has been analyzed previously, for example, in Refs. 14 and 16. In the following subsections, we study the spatial evolution of magnetic field lines by evaluating the finite-distance Lyapunov exponents^29,30 and the asymptotic-distance Lyapunov exponents.^31,32 These LEs are useful for identifying, and characterizing, regular and chaotic behavior in different regions of space. The nonlinear spreading of magnetic fields lines, which is a measure of their diffusion, is determined from the spatial evolution of a bundle of field lines that are initially close to a reference field line.

In the equations for the field lines (6), the dimensionless path length s is taken to be the spatial “evolution” parameter. The three coupled differential equations can be linearized around any point (x, y, z), to obtain the tangent space equations in terms of the Jacobian matrix. The solutions of the tangent space equations determine the LEs, which give a measure of the rate of divergence of magnetic field lines which are initially close to each other. As mentioned in the introduction, only the largest LE will be considered as it dominates the rate of divergence. There are two sets of LEs that are calculated: one set is evaluated after a short distance of evolution, while the other set is for a much longer distance of evolution. The elements of the former set are referred to as the finite-distance LEs^29,30 or finite-time LEs when the evolution parameter is time. The elements of the latter set are the asymptotic-distance LEs,^31,32 or correspondingly, time-asymptotic LEs.

Figure 2 is a contour plot of the largest finite-distance LE for a 300 × 300 ensemble of initial conditions that are uniformly distributed in x–y space with z = 0. The contours are superimposed on top of the z = 0 Poincaré surface-of-section shown in Fig. 1. For this plot, all the initial conditions are spatially evolved over a distance Δs = π/8. If we double Δs, the contours essentially retain their shape. The contours of LEs in Fig. 2 are a measure of the nonlinearity of the system. Even in the region where the field lines are regular, the evolution of the field lines is still nonlinear, and nearby field lines diverge, albeit at a slower rate. The rate of divergence, in certain regions of space, is the same whether the field lines are regular or chaotic. However, the larger values of the finite-distance LE occur in the region where the magnetic field lines are chaotic.

Figure 3 is a contour plot of the largest asymptotic-distance LE as a function of the initial conditions which are uniformly distributed in the x–y plane. The details are indicated in the figure caption. The contours are jagged as we have used fewer initial conditions compared to the number used for the evaluation of the finite-distance LEs. It is easy to correlate this plot with the Poincaré surface-of-section in Fig. 1. The magnitude of the asymptotic-distance LE is large, by several orders of magnitude, in regions where the magnetic field lines are chaotic, compared to the regions where the field lines are regular. It should be noted that the maximum of the asymptotic-distance LE is much smaller than the maximum of the finite-distance LE.

The asymptotic-distance LE shows that the largest rates of separation, between neighboring lines, are in the same region of space where the field lines are chaotic. The LEs are determined from a linearization of the full field line equation (6), and are, consequently, representative of the evolution of the field lines over a short range in s. The LEs provide a measure of the rate of decorrelation, or loss of memory, between neighboring field lines. Beyond this short range, the full set of Eq. (6) has to be solved in order to study the spatial transport, or diffusion, of field lines.

Consider the spatial evolution of an ensemble of magnetic field lines that are, initially, located in close proximity of a reference field line; all the initial positions being in the chaotic phase space shown in Fig. 1. Let r_i(s) = (x_i(s), y_i(s), z_i(s)) define the trajectory of the i-th field line; i = 0, 1, 2,…N, where N + 1 is the total number of field lines that are tracked. The initial conditions are prescribed at s = 0. If i = 0 is the reference field line, then the initial starting positions of the other N field lines are chosen to be such that $|rj−r0|≪1$ for j = 1, 2,…N.

In the numerical results that follow, the initial location of the reference field line is marked by center of the cross in Fig. 2. Setting N = 69, the starting positions of the other 69 field lines are within a region which is smaller than the size of the cross. The cross is located entirely in the chaotic phase space. Figure 4 shows the ensemble-averaged position, $⟨r(s)⟩B=(⟨x(s)⟩B,⟨y(s)⟩B,⟨z(s)⟩B)$⁠, as a function of s for the 70 field lines. $⟨r(s)⟩B$ is defined as

Initially, for $s/2π≤10, ⟨x(s)⟩B, ⟨y(s)⟩B$⁠, and $⟨z(s)⟩B$ are quite different from each other. But, as the distance along the field lines increases, the ensemble-averaged coordinates track each other; the evolution of the field lines is essentially isotropic. This is not surprising, as the ABC magnetic field (5) does not have any preferential direction.

In order to get a sense of the spatial diffusion of magnetic field lines, the following two statistical measures are evaluated as a function of s:

and

where

is the distance between the i-th field line and the reference field line, after a distance s along each field line. Here, (8) is an expression for the mean, or ensemble average, of the distance between all field lines and the reference field line. Equation (9) gives the variance or the square of the standard deviation, with respect to the reference field line, and is a measure of the spatial diffusion of the magnetic field lines.

Figure 5 illustrates μ_B(s) as a function of s for the field lines that are part of Fig. 4. Clearly, the mean increases as a function of distance along the field lines. For s/2π  ≤ 10, μ_B(s) changes rapidly; it indicates a mixing of the field lines during the early stages of evolution in s. This, as shown in Fig. 4, is directly related to the initial anisotropy in the spatial evolution of the field lines. For the range 10 ≤ s/2π≤ 200 000, a least squares fit, indicated by the dashed line in Fig. 5, yields μ_B(s) ≈ 1.09s^0.77. The least squares fit is obtained by a linear fit to the data relating $log (μB)$ to $log (s)$⁠. The coefficient of determination for the least squares fit is R² ≈ 0.98 − R² is a statistical measure of the regression line being an approximation of the actual data. From Fig. 4, in this range of s, the field lines evolve isotropically in space.

Figure 6 is a plot of $σB2(s)$⁠, as function of s, for the field lines. This gives the spatial “rate” at which the field lines spread away from the reference field line. The least squares fit to the data in the range 10 ≤ s/2π ≤ 200 000, indicated by the dashed line in the figure, yields $σB2(s)≈0.18s1.425$⁠. The coefficient of determination is R² ≈ 0.9. Since $σB2(s)$ is a measure of the diffusion coefficient, the least squares fit reveals that the spatial diffusion of the magnetic field lines is anomalous; in fact, the field lines are super-diffusive.³³ 

The separation of field lines, in magnetic fields in which randomness is prescribed ab initio, has been discussed before.³⁵ The separation of field lines is considered to be an important step in explaining the transport of charged particles in astrophysical plasmas.²² For the completely deterministic, force-free, ABC magnetic field, the field line separation is super-diffusive in regions where the field lines are chaotic. Consequently, one has to be careful in formulating the diffusion equation for the field lines.³³ 

We are particularly interested in the motion of charged particles which are initially located in the region of space where the field lines are chaotic. Passive scalars in fluid flows follow the flow lines given by (6). The motion of charged particles in a magnetic field is quite different. Instead of the field line equation (6), the particle motion is described by the Lorentz equations

where v is the velocity, at time t, of the particle at the position r, q is the charge of the particle, and m is its mass. In terms of the dimensionless parameters discussed above, the second equation in (11) takes on the form

where τ = Ωt, Ω = qB₀/m, and $v=|v|$ is normalized by Ω/λ.

In a spatially homogeneous magnetic field, the charged particles gyrate with the center of gyration, the gyrocenter, following a magnetic field line. This would be equivalent to the motion of passive scalars in homogeneous flows. In spatially inhomogeneous magnetic fields, particle orbits drift off the field lines due to gradient and curvature drifts.³⁴ The particles are also affected by the change in the magnetic field strength along the field line; they can be reflected as they move toward higher field strengths. These effects bring about significant differences between the flow of passive scalars and the motion of charged particles in force-free fields (5).

The normalized energy of a particle, $E=v2=(v⊥2+v∥2)$⁠, moving in a spatially dependent magnetic field is conserved. Here, v_⊥ and v_∥ are, respectively, components of the particle velocity perpendicular and parallel to the direction of the magnetic field, at the location of the particle. Corresponding to the normalizations we have used, $E$ is a measure of the ratio of the Larmor radius of the particle to the scale length of the magnetic field inhomogeneity.

A previous study has shown that the dynamics of charged particles in spatially chaotic field lines depends not only on the energy of the particle but also on the direction of the initial velocity vector with respect to the local direction of the field line.³⁹ We define a pitch-angle parameter, $ξ=v∥/E$⁠. If ξ = 1, the velocity of the particle is along the local magnetic field line. For ξ = 0, the velocity is perpendicular to the local magnetic field line.

In the following discussions, it is assumed that the energy of the particles is small, so that the ratio of the Larmor radius, for ξ = 0, to the magnetic field scale length is small. Then, it is easier to identify the field line associated with the particle. For larger energies, the particle will sample a variety of inhomogeneous magnetic field lines during a single cyclotron orbit, thereby, adding to the complexity of the dynamics. Also, a particle started in the region where the field lines are chaotic could sample the region where the field lines are regular.

For the range of dimensionless parameters which will be used in the ensuing simulations, it is useful to put into context the corresponding physically relevant quantities. From the approximate parameter range for astrophysical plasmas in interstellar space, tabulated in Ref. 40, assume that the electron density is n_e = 10³m⁻³, the temperature of the plasma is T_e = 10²K, and the magnetic field strength is B₀ = 0.1 nT. The ratio of the kinetic pressure to the magnetic pressure is approximately 7 × 10⁻⁴, if the primary constituents of the plasma are electrons and protons. For the electrons, Ω_e ≈ 17.6 rad s⁻¹, and their thermal speed is v_Te ≈ 3.9 × 10⁴m s⁻¹. For thermal electrons, $E=ρeλ$⁠, where ρ_e is the electron Larmor radius, and $E$ is the ratio of the Larmor radius to the scale length of the magnetic field. If $E=10−2$⁠, then λ ≈ 4.5 × 10⁻⁶ m⁻¹.

Another example is the cool solar corona, for which the parameters are tabulated in Ref. 41. The typical parameters are n_e = 10¹⁵m⁻³, T_e = 10⁶K, and B₀ = 10³T. If the main constituents are assumed to be electrons and protons, the ratio of the kinetic pressure to the magnetic pressure is, approximately, 7 × 10⁻². For electrons, Ω_e ≈ 1.76 × 10⁸rad s⁻¹ and v_Te ≈ 3.9 × 10⁶m s⁻¹ so that, for $E=10−2$⁠, we obtain λ ≈ 0.45 m⁻¹.

The tangent-space equations, necessary to study the behavior of Lyapunov exponents, for particle trajectories can be easily constructed by linearizing Eqs. (11). The resulting equations are given in terms of the Jacobian matrix. For the six-dimensional phase space, there are six distinct LEs. As stated earlier, only the behavior of the largest of these Lyapunov exponents will be studied. In order to illustrate the differences in the particle dynamics, the finite-time LEs for two different values of ξ, 1 and 0, are studied. For $E=10−2$⁠, Figs. 7 and 8 are, respectively, the contour plots for the two values of ξ in the x–y plane. The particles are, initially, distributed in the x–y plane with z = 0, and the finite-time LEs are determined after a time lapse of δt = 3π/20. Doubling δt has no effect on the contour plots, so that the results are well within the range of validity of the tangent-space equations. From the two figures, it can be observed that, while the range for the largest LE is essentially the same in the two figures, the contours of constant LEs occupy different regions of configuration space.

In Fig. 7, the maximum range between the displayed Lyapunov contours is, approximately, 0.004, while in Fig. 8, it is 0.007. Even though these ranges are small, they are still meaningful. The technique used for solving the magnetic field line equation (6), and particle orbit equations (11) is highly accurate. The symplectic algorithm that is implemented in the numerical codes has been discussed in Ref. 39.

Figure 9 shows some of the contours of the largest finite-time LE, superimposed on a selection of contours of $|B|$⁠. The contours of constant LEs almost follow surfaces of constant $|B|$⁠. There is some noticeable drift off the $|B|$ surfaces, but the large LEs are confined to within $|B|>1.8$ and $|B|<1$⁠. Upon comparing this figure with the surface-of-section plot in Fig. 1, it is evident that the largest LEs are not wholly in the region where the field lines are chaotic.

Figure 10 amplifies the distinction between ξ = 0 (dashed lines) and ξ = 1 (solid lines), where the same amplitude level contours are displayed on top of each other for the two values of ξ. There is little overlap between the two sets, and, for ξ = 0, a larger region of space is covered by the larger LEs. Figure 11 shows contours of the largest finite-time LEs when, initially, $ξ=1/2$⁠, i.e., initially, v_⊥ = v_∥. The overall structure is almost similar to that for ξ = 0.

The range and magnitude of the largest, finite-time LEs increases with the energy of the particles. To make a reasonably fair comparison for different energies, the product $Eδt$ is kept fixed at the same value as for the case of $E=10−2$⁠. Linearly, the speed of the particles is proportional to $E$⁠. By keeping $Eδt$ constant, the particles move approximately the same distance away from their initial positions. The corresponding results, displayed in Table I, show an increase in the range, as well as the magnitude, of the finite-time LEs.

While the finite-time LEs, for a fixed energy, display interesting features as a function of ξ, the time-asymptotic LEs are essentially independent of the initial value of ξ. This is primarily because, in the long time limit, particles get to sample the entire range of ξ as they move through the chaotic region. In the long time limit, the particles sample the entire region of chaotic phase space. Figure 12 is a contour plot of the largest time-asymptotic LE as a function of initial conditions. The particles are, initially, uniformly spread in x and y, with z = 0. The contours are jagged as the number of initial conditions is fewer than for finite-time LEs. Comparing with the surface-of-section in Fig. 1, it is clear that the large values of the time-asymptotic LE for the particles are in regions where the field lines are chaotic. In regions where the field lines are regular, the time-asymptotic LEs for the particles are small. It should be noted that, the contours in Fig. 12 are for the logarithm of the time-asymptotic LEs. Thus, the variations range over 5 orders of magnitude. In regions where the field lines are regular, the particle time-asymptotic LEs are, at least, two orders of magnitude smaller than those in the regions where the field lines are chaotic. The magnitudes of time-asymptotic LEs are considerably smaller in comparison to the finite-time LEs.

The time-asymptotic LE in Fig. 12 shows that the divergence of trajectories of two initially nearby particles is most prominent in the region where the field lines are chaotic. The nonlinear long time behavior of particles in this region will be a diffusive process, which occurs after particle orbits are no longer correlated. The LEs give the rate at which nearby trajectories diverge, but do not yield the decorrelation time for these trajectories. The correlation time is determined from the velocity autocorrelation function.⁴² The diffusive process manifests itself only after correlations have decayed away. The study of the velocity autocorrelation function, and the diffusion process, is carried out numerically for a set of particles which are in the region where the field lines are chaotic. The trajectories of an ensemble of particles, initially located at the point marked by a cross in Fig. 2, which have the same energy but different values of the pitch-angle ξ, are following for long times. The cross in Fig. 2 also indicated the initial location of the field lines for studying the spatial diffusion of field lines.

In time-independent magnetic fields, the energy of a particle is conserved, so that $|vx|,|vy|,|vz|≤E$⁠. For long times, each velocity component evolves chaotically in time within the given range. We will focus our attention on the x-component of the velocity vector for an analysis of the autocorrelation function. The other components give similar results for the autocorrelation function.

Consider an ensemble of N_p particles. For each particle, the Lorentz equations (11) are solved to construct, at regular time intervals δt, an array for the x-component of the velocity vector along the particle trajectory. For a particular particle p, let the array be of the form ${ v1(p),v2(p),v3(p),…vN(p) }$⁠, where (N − 1) δt is the total time for which the orbit is evaluated. The ensemble-averaged, mean of the x-component of the velocity is

where ⟨…⟩_p indicates ensemble averaging. For a lag time Δt = kδt, the ensemble-averaged autocorrelation function is⁴³ 

where k = 0, 1, 2,… K with K < N. The outer sum is the ensemble average over N_p particles, while the inner sum is the discrete autocorrelation function with a lag time represented by k. Thus, ⟨C_k⟩_p is the autocorrelation function for the k^th lag, with the implicit assumption that the corresponding time is Δt.

Figure 13 is a plot of $⟨Ck⟩p/⟨C0⟩p$ as a function of the lag time Δt/2π. The energy of the particles, E = 10⁻⁴, is taken to be small enough so that, for ξ = 0, the Larmor radius is small enough that the particles do not cross into the spatial region of regular field lines. The correlation function for, either the y or z component of the velocity vector, is similar to that for the x component. The correlation function decays to zero for Δt ≈ 120, so that the decorrelation time is ΔT_c ≈ 240π.

The ensemble-averaged, velocity autocorrelation function changes with energy. Figure 14 shows the correlation function when the energy of the particles is E = 10⁻². The decay of ⟨C_k⟩_p is much more rapid than in Fig. 13, but there persists a tail which indicates the presence of long-time correlations.

For a particle p, an array of position vectors ${ r1(p),r2(p),r3(p),…rN(p) }$ is evaluated from the trajectory of the particle. Using the same notation as in the previous subsection, the ensemble-averaged, mean square displacement⁴² is expressed as

where k = 0, 1, 2,… K, with K < N. The outer sum, corresponding to $⟨…⟩p$⁠, represents ensemble-averaging, while the inner sum, corresponding to $⟨…⟩$⁠, is the mean square displacement obtained by averaging over all pairs of points that are separated in time by Δt = kδt.⁴⁴ 

The mean square displacement is a measure of the spatial transport of particles.³³ For example, if $⟨⟨Δr2⟩⟩p∼(Δt)2$⁠, then the spatial transport is normal, and the particles have a uniform rectilinear, or ballistic, motion. For classical random walk diffusion, the mean square displacement scales like $⟨⟨Δr2⟩⟩p∼Δt$⁠. Figure 15 is a log-log plot of the ensemble-averaged, mean square displacement as a function of the lag time, Δt/2π, for the same particles and conditions used to evaluate the autocorrelation function (14). A closer analysis of this figure yields three regions, distinguished by different scalings of $⟨⟨Δr2⟩⟩p$ with respect to Δt. The first region, shown in Fig. 16, spans the time lags Δt/2π ∈ [1, 10]. The dotted line is a least squares fit to the mean square displacement, and gives the scaling $⟨⟨Δr2⟩⟩p≈5×10−5(Δt)1.93$⁠. The exponent of Δt indicates that, over the given range of time lag, the particles are essentially undergoing ballistic motion. In the second region, shown in Fig. 17, covers the range Δt/2π ∈ [11, 200]. The least squares fit yields a scaling $⟨⟨Δr2⟩⟩p≈1.92×10−4(Δt)1.67$⁠. From the autocorrelation function in Fig. 13, this is the time range over which the correlations decay away. The upper limit of Δt in this range is, somewhat, flexible. The third, and final, region spans the range Δt/2π ∈ [201, 10⁵] in which the least squares fit, indicated by the dashed line in Fig. 15, yields $⟨⟨Δr2⟩⟩p≈2.44×10−4(Δt)1.575$⁠. Over this range of time lags, the particles diffuse in space; in fact, they are super-diffusive, as the exponent of Δt is greater than 1.³³ Since the super-diffusive region spans nearly three decades in time lags, it does not represent transient phenomenon.

At higher energies, the separation of the mean square displacement into three regions becomes less clear. Figure 17 shows the mean square displacement as a function of the time lag for $E=0.1$⁠. If the same analysis as for the results in Fig. 15 is followed, then $⟨⟨Δr2⟩⟩p≈5.93×10−3(Δt)1.75$ in the range Δt/2π ∈ [1, 10]. The autocorrelation function, shown in Fig. 14, decays to below the 1% level for Δt/2π ≳ 4900. In the range Δt/2π ∈ [11, 5000], the least squares fit gives $⟨⟨Δr2⟩⟩p≈4.62×10−3(Δt)1.72$⁠. It is noteworthy that, unlike the case for $E=0.01$⁠, the scaling in the two sub-domains Δt/2π ∈ [1, 10] and Δt/2π ∈ [11, 5000] is very similar. In the diffusive range, Δt/2π ∈ [5001, 100000], the least squares fit to the mean square displacement gives $⟨⟨Δr2⟩⟩p≈0.01(Δt)1.65$ (Fig. 18). The spatial transport of particles is super-diffusive. Consequently, for the higher energy, the ensemble-averaged, mean square displacement shows less of a variation with respect to time lag over the entire range of time lags, Δt/2π ∈ [1, 100000].

The ratio of the mean square displacement, evaluated in the diffusive regime for the two energies, is $⟨⟨Δr2⟩⟩p(E=0.1)/⟨⟨Δr2⟩⟩p(E=0.01)≈82$ for Δt = 10⁴. Hence, the spatial superdiffusion of particles is weakly dependent on the time lag for large values of the time lag.

The Arnold-Beltrami-Childress force-free magnetic field provides an interesting paradigm for studying the spatial diffusion of magnetic field lines, and the spatial transport of charged particles. The ABC field is also a paradigm, in astrophysical plasmas, for the magnetic field in regions of high plasma conductivity, where the Lorentz forces cannot be balanced by pressure gradients. The equations for the evolution of the ABC magnetic field lines, and for the particle trajectories in these fields, are Hamiltonian systems. Consequently, the corresponding phase space is a mix of regular and chaotic trajectories.

The nonlinear spatial evolution of field lines is determined from the properties of the Lyapunov exponents. Only the largest LEs are followed since they represent the dominant growth rate. Over short distances along the field lines, the finite-distance LEs are large in the region where the field lines are nearly tangential to the x–y plane. The asymptotic-distance LEs essentially follow the z = 0 Poincare surface-of-section for the magnetic field lines. The largest asymptotic-distance LEs are in regions where the magnetic field lines are chaotic. In these regions, the variance of an ensemble of field lines, with respect to a reference field line, scales like s^1.425, where s is the path length along the field lines. Hence, the spatial spreading of chaotic field lines is super-diffusive.

There have been prior studies on chaotic magnetic field lines that pertain to laboratory fusion plasmas.^45,46 While the magnetic field in fusion plasmas is not force-free, it is interesting to compare the diffusion characteristics of chaotic fields. In the Rechester-Rosenbluth paper,⁴⁵ the diffusion of chaotic field lines is ordinary—the mean square displacement is a linear function of the distance along the toroidal direction. This is in stark contrast to the superdiffusion of magnetic field lines associated with the ABC fields. The transport of heat due to the chaotic magnetic fields in Ref. 45 is due to the classical random walk diffusion;³³ the transport of particles in the chaotic ABC fields is super-diffusive. In Ref. 46, the field lines are initially super-diffusive and, subsequently, after a few iterations, sub-diffusive. The sub-diffusive behavior is attributed to the existence of island chains in the Poincaré surface-of-section for the field lines. The results in our paper differ in two salient ways. First, the initial evolution of the ABC field lines is correlated and anisotropic; there is mixing of the field lines but no diffusion. Second, after the field lines have decorrelated and are evolving isotropically in space, the field lines are in the super-diffusive regime. The existence of island chains in the Poincaré surface-of-section in Fig. 1, does not impede the spatial super-diffusive behavior of the field lines—a consequence of the three-dimensional structure of the ABC magnetic fields.

From Fig. 6, and the discussion associated with this figure, it is meaningful to evaluate the diffusive aspect of field lines only after the decay of the correlation function. As shown in Fig. 4, the decay of the correlation function is related to the isotropic evolution of the field lines, which starts around s/2π ≈ 20. Consequently, the magnitude of the mean free path for the field lines, which is the ensemble-averaged correlation length, is s_mfp ≈ 40π.

For charged particles in ABC fields, the finite-time LEs depend on the initial pitch-angle of the particles. For particles started off with their velocities along the field lines, the finite-time LEs are large in regions where the magnitude of the magnetic field is large. For particles with their initial velocities perpendicular to the field lines, the large finite-time LEs are more spread out in space. These particles drift off the magnetic surfaces in the inhomogeneous magnetic field. The spatial distribution of the largest time-asymptotic LEs is, essentially, independent of the pitch-angle. The exponents are large in regions where the magnetic field lines are chaotic. In this region, the particle motion is also chaotic.

The decay of the ensemble-averaged, velocity autocorrelation function, for particles in the chaotic phase space, depends on the energy of the particles. For higher energy, the correlations persist for long time lags. The ensemble averaged, mean square displacement can be divided into essentially three different stages of evolution with respect to the time lag. This is particularly true for low energies where, in the first stage, the mean square displacement is proportional to (Δt)². This corresponds to ballistic motion or normal spatial transport.³³ The second stage is the mixing stage representing the time period over which the velocity autocorrelation function decays to near zero. The third and final stage is where the particles experience diffusion–anomalous superdiffusion. The mean square displacement is not a linear function of time. The distinction among the three regions is less pronounced as the energy of the particles is increased. However, the distribution of particles continues to evolve super-diffusively.

The magnitude of the mean free path for the particles, just like for the magnetic field lines, is the ensemble-averaged distance that the particles travel before the autocorrelation function decays away. From Fig. 13, the decorrelation time is ΔT_c ≈ 240π for particles with speed $v=E=10−2$⁠. Then, it follows that the normalized, ensemble-averaged, magnitude of the mean free path is r_mfp ≈ 2.4π. It is interesting to compare this with the magnitude of the mean free path for the field lines s_mfp ≈ 40π. The particles get decorrelated over a distance which is much shorter than the decorrelation length for the field lines.

For the case of particles with higher energy, as shown in Fig. 14, v = 0.1 and ΔT_c ≈ 9800π. Accordingly, r_mfp ≈ 980π and the particles orbits get decorrelated after a distance which is much longer than the decorrelation length for the field lines. It can be safely concluded that, the mean free path of particles is not related to the mean free path of the magnetic field lines.

The separation of time scales between the evaluation of the Lyapunov exponents and the evaluation of the diffusion scaling is distinct. The evaluation of the Lyapunov exponents, whether finite time (finite distance) or asymptotic, occurs over a limited range in space or time. The limitation is imposed by the linearization of the trajectory that varies nonlinearly in space and/or time; the linearized equations leading to either the contraction or expansion of a set of initial conditions. The diffusive process is a consequence of the full, nonlinear equations and occurs after the correlation function has decayed away. The persistence of correlations over extended spatial and temporal scales separates the diffusive process from the local, short scale, dynamics that is the basis for calculating the Lyapunov exponents. Subsequently, it is not possible to determine a connection between the Lyapunov exponents and the anomalous diffusion of either the magnetic field lines or the particles.

The spatial superdiffusion of the ABC magnetic field lines and of the particles is likely linked. But, the rate of superdiffusion of particles depends on their energies while that of the magnetic field lines is fixed. Regardless, the fact that both the field lines and the particles do not follow regular diffusion in space, it is imperative that the conventional approach to understanding the transport of particles in plasmas has to be modified. The conventional approach follows the classical theory of random walks in which the mean square displacement is proportional to the evolution parameter; in our case, s for magnetic field lines and t for the particles. It does not take into account the persistence of correlations over long distances or for long times. In order to understand the transport of particles, the results in this paper reinforce the need to go beyond the quasilinear model, which is based on the simple random walk processes represented by Brownian motion.^33,47,48

The work was supported by U.S. DOE Grant No. DE-FG02-91ER-54109. B.D. was supported by NSF Grant No. AGS-1062050 and through the Individual Investigator Distinguished Research award of the University of Alabama in Huntsville. B.D. is thankful to Gary Zank for encouragement and support. D.M. was supported in part by the European Research Council under the AstroDyn Research Project and the Swedish Research Council.