PARTICLE ACCELERATION VIA RECONNECTION PROCESSES IN THE SUPERSONIC SOLAR WIND

Magnetic reconnection has been widely invoked to explain the energization of ions and electrons throughout the heliosphere, from solar flares and Earth's magnetosphere (e.g., Drake et al. 2006a; Oka et al. 2010; Birn et al. 2012; Cargill et al. 2006, 2012 and references therein) to even the heliopause (Lazarian & Opher 2009). An interesting series of studies has also considered particle energization in astrophysical flows, identifying related first- and second-order Fermi acceleration-like mechanisms due to magnetic reconnection (de Gouveia dal Pino & Lazarian 2005; Kowal et al. 2011, 2012). Despite the wide application of some form or another of reconnection mechanisms to a variety of heliospheric and astrophysical environments, the precise physical mechanism by which magnetic reconnection energizes charged particles is still not fully understood. Typically, particle acceleration within reconnecting current sheets has been invoked because this can generate direct current (DC) electric fields that accelerate (or decelerate) charged particles. The idea that a statistical acceleration mechanism related to magnetic reconnection was responsible for accelerating particles to high energies was first advanced quantitatively by Matthaeus et al. (1984) and Ambrosiano et al. (1988). They investigated the effect of turbulence on particle acceleration in an MHD field by computing test particle trajectories in turbulent MHD reconnecting fields, examining the acceleration mechanism in some detail. Matthaeus et al. (1984) and Ambrosiano et al. (1988) found that turbulence influences the acceleration in two ways. It enhances the reconnection electric field while producing a stochastic electric field that gives rise to momentum diffusion; and it produces magnetic "bubbles" and other irregularities that can temporarily trap test particles in the strong reconnection electric field for times comparable to the magnetofluid characteristic time. A power-law distribution for the energetic particle distribution was found from their test particle simulations.

Here we consider the acceleration of particles assuming an emerging paradigm for the dissipation of magnetic turbulence via localized small-scale reconnection processes. We examine the physics of magnetic island merging and include this statistically in deriving a transport equation for particles experiencing energization in a "sea" of dynamically interacting magnetic islands embedded in a super-Alfvénic flow such as the supersonic solar wind. We solve the transport equation for an incompressible background flow and show that the accelerated particle distribution is a power law in particle speed c. The spectral index of the power-law distribution is found to depend on the Alfvén Mach number, the characteristic diffusion time of charged particles, and the characteristic time for magnetic island contraction.

Below, we discuss reconnection in a turbulent magnetofluid with magnetic island merging. Two forms of the transport equation for particles experiencing pitch-angle scattering and interacting with reconnecting magnetic islands are derived; the first, a gyrophase-averaged transport equation appropriate to nonisotropic particle distributions and a second that assumes near isotropy of the particles, which yields a simpler transport equation. The isotropic 1D form of the transport equation is solved in the steady-state for an incompressible flow.

2.1. Simulations of Multiple Magnetic Islands, Reconnection, and Turbulence

2.2. Observations of Magnetic Islands and Merging

2.3. Summary of the Simulations and Observations

2.4. Magnetic Island Merging and Particle Energization

Although we have emphasized above that reconnection is intrinsic to MHD turbulence and a natural consequence of the nonlinear cascade process, considerable insight into the interaction and formation of magnetic islands is gained from more typical symmetric, isolated models of reconnection (both kinetic and MHD; e.g., Sato et al. 1982; Ambrosiano et al. 1988; Hoshino et al. 2001; Drake et al. 2005, 2006a; Pritchett 2006, 2008; Oka et al. 2010; Huang & Bhattacharjee 2013). Of particular interest is multimagnetic island coalescence (Finn & Kaw 1977; Pritchett & Wu 1979; Tajima et al. 1987; Pritchett 2007, 2008; Wan et al. 2008; Oka et al. 2010), especially in view of the observations presented above that showed active merging of magnetic islands occurring in the solar wind.

In the context of isolated reconnection simulations, a tearing-mode instability typically generates localized currents and multiple x lines. The tearing-mode instability is not necessary for the formation of multiple islands separated by x points because these arise naturally in freely decaying turbulence—see, e.g., Figure 4 of Servidio et al. (2010). In either case, neighboring islands can be attracted to one another by the Lorentz force. As illustrated in Figure 1, the merging of the two islands introduces a reconnection event, marked by the heavy X.

Zoom In Zoom Out Reset image size

In the context of isolated reconnection models, the reconnection due to the merging of magnetic islands is termed "antireconnection" (Pritchett 2008), and the electric field induced by the antireconnection event is opposite in direction to the electric field of the primary reconnection event. Such terminology is of course inappropriate to reconnection associated with island merging in the context of turbulence because there is no "primary reconnection event." Nonetheless, the physics of the island merger is precisely the same for both MHD turbulence simulations and simulations of isolated reconnection, and we will sometimes use the "antireconnection" electric field terminology where its use is unavoidable. As illustrated in Figure 1, the two magnetic islands eventually merge to form a single large island. Illustrated in Figure 2 is a cartoon, drawn from the simulations of Servidio et al. (2009, 2010, 2011), that shows the in-plane magnetic field. A "sea" of different sized magnetic islands is clearly apparent, surrounded by in-plane magnetic field lines. A magnetic field line component orthogonal to the plane depicted in Figure 2, the guide field, is also present. Particles propagate along magnetic field lines, sometimes experiencing trapping in magnetic islands and other times propagating along the guide magnetic field component. In this way, charged particles experience energization by being trapped in and escaping from multiple magnetic islands.

Zoom In Zoom Out Reset image size

Test particle simulations in a dynamical "sea" of magnetic islands show that particles experience energization, with spectra that frequently resemble a power law. Three processes essentially govern the gain in energy that a particle experiences. The first, discussed by Drake et al. (2006a), is based on the idea that magnetic islands formed during magnetic reconnection will contract. Particles trapped within the collapsing island experience repeated reflections at the ends of converging mirrors, gaining energy in a first-order Fermi manner (Drake et al. 2006a). PIC simulations presented by Oka et al. (2010) showed that magnetic islands could "oscillate" or " bounce," i.e., undergo both contraction and expansion. Particles trapped in an oscillating island may experience a second-order Fermi energization. A second process for particle energization is also related to the shortening of the magnetic island field line length but is due now to the merging of two adjacent islands. The merged magnetic field lines (Figure 1) shorten as the merging progresses. The contraction of the magnetic field leads to an increase in the particle velocity component parallel to the local magnetic field and a decrease in the perpendicular energy. By sampling multiple merging magnetic islands, particles can be energized via a second-order Fermi process (Drake et al. 2013). A third energization mechanism is discussed by Oka et al. (2010). In this case, particles trapped in a merging magnetic island (Figure 1) can experience multiple interactions with the reconnection electric field generated by the merging of two islands. This is quite distinct from the usual energization of particles by an electric field generated by the reconnection of two long antiparallel magnetic field lines. Unlike the usual particle energization process via reconnection, particles trapped in a merging magnetic island experience an extended period of interaction with the reconnection electric field, and thus significant energy gain is possible (Oka et al. 2010; Pritchett 2008; Tanaka et al. 2010). The test particle simulations of Oka et al. (2010) suggest that this third energization mechanism may be the dominant energization process for particles in reconnection layers.

In this paper, we adopt the emerging paradigm that turbulence in the supersonic solar wind is dissipated via quasi-2D coherent structures or magnetic islands that are created through dynamical processes. The islands interact nonlinearly, merge, contract, stretch, and attract and repel one another. We will consider the three processes above that have been identified as particle energization mechanisms related to the merging of magnetic islands, i.e., magnetic island contraction, magnetic field shortening due to merging, and the antireconnection electric field.

We therefore hypothesize that medium- and small-scale magnetic islands exist throughout the solar wind, probably as the end-point of turbulent processes, and that they experience multiple mergers. In the process of merging, particles are trapped within the dynamically evolving islands. Charged particles gain energy via magnetic island contraction, merging, or repeated interactions with the electric field induced by the merging process itself. Some particles gain energy, and others lose energy. We will show that the accelerated particle distribution is essentially a power law with an index determined by the Alfvén Mach number and the microphysical timescales associated with magnetic island contraction and charged particle diffusion.

A charged particle propagating through a sea of dynamically interacting magnetic islands (Figure 2) surrounded by magnetic field lines will experience a total force

$$\begin{equation} {\bf F} = {\bf F}_L + {\bf F}_c + {\bf F}_m + {\bf F}_E, \end{equation} \tag{ 1 }$$

where F_L denotes the usual Lorentz force, F_c particle acceleration due to island contraction (Drake et al. 2006a), F_m particle energization due to island merging (magnetic field line shortening) (Drake et al. 2013), and F_E is the particle acceleration due to energization by the reconnection electric fields induced by magnetic island coalescence or other events. Here F_L = q(E + v × B), and q and m denote particle charge and mass, respectively, v is the (nonrelativistic) particle velocity, E is the large-scale or mean electric field, and B the large-scale or mean magnetic field vector.

Consider first the acceleration term due to magnetic island contraction F_c = m.dv/dt|_c. Drake et al. (2006a) show the trajectory of a charged particle trapped in a contracting island (their Figure 2). Particularly revealing are Figures 2(b) and (c) of Drake et al. (2006a), which show that the trapped particle gains energy in a monotonic manner, with increases occurring at the converging island ends. These figures illustrate very clearly that the process of energy gain is a standard first-order Fermi process. The acceleration term is most simply understood on the basis of the conservation of the first and second adiabatic invariants. Because parallel action is constant, i.e., ∫v_∥dℓ' = const., where ℓ' is the field line length, we have ∫v_∥dℓ' ∼ v_∥ℓ ∼ const. on assuming a characteristic length scale ℓ for an island. This therefore yields

$$\begin{eqnarray*} &&\frac{d v_{\parallel } }{{dt}} = -\frac{1}{\ell } \frac{d \ell }{{dt}} v_{\parallel } \equiv \eta _{c \parallel } v_{\parallel }. \end{eqnarray*}$$

For field line shortening, dℓ/dt < 0 and thus η_c∥ > 0, implying that the parallel velocity increases in time during magnetic island contraction. Conservation of the first adiabatic invariant $v_{\perp }^2 /B = {\rm const.}$ shows that

$$\begin{eqnarray*} &&\frac{d v_{\perp }^2 }{{dt}} = \frac{1}{B} \frac{{dB}}{{{dt}}} v_{\perp }^2 \equiv \eta _{c \perp } v_{\perp }^2. \end{eqnarray*}$$

Because the island is contracting, dB/dt > 0, and so the perpendicular energy increases in time as the magnetic island contracts. Recall that for a 2D magnetic island, the magnetic field decreases with radial distance as B = Ψ(r)/r, where Ψ(r) is the magnetic flux,⁴ so that for a characteristic island size, Ψ ∼ Bℓ. The contracting island flux is constant so that

$$\begin{eqnarray*} &&\frac{{dB}}{{{dt}}} = -\frac{\Psi }{\ell ^2} \frac{d \ell }{d t} \Rightarrow \frac{1}{B} \frac{{dB}}{{dt}} = - \frac{1}{\ell } \frac{d \ell }{{dt}}, \end{eqnarray*}$$

from which we infer that

$$\begin{eqnarray*} &&\eta _{c \parallel } = \eta _{c \perp } \equiv \eta _c. \end{eqnarray*}$$

For a contracting island, we therefore obtain the acceleration terms (Drake et al. 2006a, 2013)

$$\begin{equation} \left. \frac{d v_{\parallel } }{{dt}} \right|_c = \eta _c v_{\parallel } \qquad {\rm and} \qquad \left. \frac{d v_{\perp } }{{dt}} \right|_c = \frac{1}{2} \eta _c v_{\perp }. \end{equation} \tag{ 2 }$$

The island contracts at a speed V = dℓ/dt. Dimensionally, the rate of energy gain due to conservation of longitudinal action is then simply η_c = −αV/ℓ, where the coefficient α is taken to be the ratio of the reconnecting magnetic island magnetic field energy density to the total magnetic energy density, i.e., b²/B² (Drake et al. 2006a, 2013; Bian & Kontar 2013). A characteristic speed of contraction is V ∼ V_A, the local Alfvén speed. Of course, if the magnetic island is expanding, a particle will lose energy at the corresponding rate. If contraction is favored, a particle will experience a gain in energy. Bian & Kontar (2013) use an identical argument, identifying η_c as proportional to the mean compression associated with an average over the ensemble of magnetic islands in the region of interest.

Consider now magnetic island merging (Figure 1) while neglecting the induced antireconnection electric field. Here we follow the analysis of Fermo et al. (2012) and Drake et al. (2013). On crudely approximating the merged island as a circle, we assume that the merging of two islands of initial areas A₁ and A₂ yields a single island of area A = A₁ + A₂ = const. Thus, the merged island scales according to $r \sim \sqrt{A / \pi }$. Again, we assume that the first and second adiabatic invariants hold, so that for an island of characteristic magnetic field line length ℓ,

$$\begin{eqnarray*} &&\frac{d v_{\parallel } }{{dt}} = -\frac{1}{\ell } \frac{d \ell }{{dt}} v_{\parallel } \equiv \eta _{m \parallel } v_{\parallel }. \end{eqnarray*}$$

Because the island area is assumed constant, $2\pi r = 2 \pi (r_1^2 + r_2^2)^{1/2} < 2\pi ((r_1 + r_2)^2)^{1/2}$, implies that the characteristic length ℓ decreases with merging, or dℓ/dt < 0. The second assumption is that magnetic flux is constant (because reconnection simply changes magnetic connectivity, and the flux will be the larger of the fluxes in the two initial islands). Thus, combining the constant area of the merged island A ∼ 2πrℓ with the constant flux Ψ ≃ B2πr yields

$$\begin{eqnarray*} &&\frac{B}{\ell } \sim {\rm const.} \qquad {\rm or} \qquad \frac{1}{B} \frac{{dB}}{{dt}} = \frac{1}{\ell } \frac{d \ell }{{dt}} < 0. \end{eqnarray*}$$

Consequently, during island merging, the perpendicular particle energy changes according to

$$\begin{eqnarray*} &&\frac{d v_{\perp }^2 }{{dt}} = \frac{1}{B} \frac{{dB}}{{dt}} v_{\perp }^2 = \frac{1}{\ell } \frac{d \ell }{{dt}} v_{\perp }^2 \equiv \eta _{m \perp } v_{\perp }^2. \end{eqnarray*}$$

Thus, η_m∥ = −η_m⊥ ≡ η_m, and the particle acceleration terms for merging islands (neglecting the induced antireconnection electric field) become (Drake et al. 2013)

$$\begin{equation} \left. \frac{d v_{\parallel } }{{dt}} \right|_m = \eta _m v_{\parallel }, \qquad {\rm and} \qquad \left. \frac{d v_{\perp } }{{dt}} \right|_m = -\frac{1}{2} \eta _m v_{\perp }. \end{equation} \tag{ 3 }$$

In this case, a particle gains in parallel energy and loses in perpendicular energy.

Denote the fluctuating reconnection electric field by δE' and regard this as a perturbation on the electric field E measured in the stationary frame. We can therefore write

$$\begin{equation} \left. {\bf F}_E = m \frac{d {\bf v} }{{dt}} \right|_E = q \delta {\bf E}^{\prime }. \end{equation} \tag{ 4 }$$

Note that the characteristic timescale for the existence of the fluctuating reconnection electric field for a characteristic island size ℓ and Alfvén speed V_A is τ_E ∼ ℓ/V_A, yielding $\langle \delta {\bf E}^{\prime } \rangle \sim \delta E_0 e^{-t / \tau _E}$ for a typical individual event.

We begin with the 6D Liouville equation and assume that fluctuating magnetic fields define a scattering term δf/δt)_s and then rewrite the equation in the Vlasov form. The acceleration term is separated into two parts, one corresponding to the usual Lorentz electromagnetic force term F_L associated with the large-scale fields and the other to energization by coalescing magnetic islands, as discussed above. In principle, all electric fields generated by reconnection processes are included but in dissipative turbulence; this amounts to electric fields induced by magnetic islands. We restrict our attention to nonrelativistic particles for the present. We therefore have

$$\begin{eqnarray} \frac{\partial f}{\partial t} &+& v_i \frac{\partial f}{\partial x_i} + \frac{{\bf F}_L }{m} \cdot \nabla _v f\nonumber\\ & =& \left. \frac{\delta f}{\delta t} \right)_s - \frac{{\bf F}_c }{m} \cdot \nabla _v f - \frac{{\bf F}_m }{m} \cdot \nabla _v f\nonumber\\ && - \frac{{\bf F}_E }{m} \cdot \nabla _v f - \nabla _v \cdot\left(\frac{{\bf F}_L }{m} + \frac{{\bf F}_c }{m} + \frac{{\bf F}_m }{m} + \frac{{\bf F}_E }{m} \right)\nonumber\\ f& =& \left. \frac{\delta f}{\delta t} \right)_s - \left. \frac{d {\bf v} }{{dt}} \right|_c \cdot \nabla _v f - \left. \frac{d {\bf v} }{{dt}} \right|_m \cdot \nabla _v f\nonumber\\ && - \left. \frac{d {\bf v} }{{dt}} \right|_E \cdot \nabla _v f - \nabla _v \cdot \left(\frac{{\bf F}_c }{m} \right) f, \end{eqnarray} \tag{ 5 }$$

where ∇_v is the gradient operator with respect to the particle velocity. Ordinarily, the velocity divergence terms on the right-hand side of Equation (5) are absent in the Vlasov equation because (typically) the force terms are independent of the particle velocity (except the Lorentz term v × B for which ∇_v · (v × B) = 0). However, both F_c/m = η_c(v_∥, v_⊥/2) and F_c/m = η_c(v_∥, −v_⊥/2) are functions of particle velocity, for which ∇_v · F_c/m = 2η_c and ∇_v · F_m/m = η_m − η_m = 0.

Our purpose is to derive an equation that describes the transport of charged particles in a region of characteristic size L that is much larger than the characteristic size of magnetic islands, i.e., L ≫ ℓ, with a timescale characterized by the large-scale super-Alfvénic flow velocity U within which the reconnecting region is embedded. Consequently, the characteristic timescale for the system is T ∼ L/U. If we normalize Equation (5) with respect to the transport time and length scales and suppose that the characteristic pitch-angle scattering time is ∼τ_s and the characteristic fluctuating merging magnetic island timescale is τ_{c, m, E} ∼ ℓ/V_A, then it is straightforward to see that the scattering term is proportional to T/τ_s, and the merging terms are proportional to T/τ_{c, m, E} ∼ (V_A/U)(L/ℓ). Because T ≪ τ_s and L ≪ ℓ, both terms on the right-hand side of Equation (5) may justifiably be regarded as O(1/ε), where ε(≪ 1) is the measure of the ratio of the scattering, contracting, merging, or fluctuating reconnection electric field timescale to the transport timescale. The right-hand side can therefore be separated from the slowly varying left-hand side.

In the absence of an applied electric field, E = −U × B is the motional electric field associated with the large-scale flowing plasma. To eliminate the electric field and move into the frame that ensures that particle scattering does not result in momentum or energy changes, we transform to the flow frame and rewrite v = c + U or c = v − U. We further assume that the flow is super-Alfvénic, specifically, that U ≫ V_A, V_A the Alfvén speed. In principle, we could alternatively assume that the cross-helicity of the background turbulent fluctuations responsible for scattering the particles is zero. In either case, we can avoid boosting to the background flow plus Alfvén speed. Consequently, the transport equation below is appropriate only to the super-Alfvénic solar wind and not to the sub-Alfvénic regions such as low in the solar corona. For simplicity, we assume that flow accelerations are absent in the large-scale flow field. The transformation of Equation (5) to the plasma flow frame then yields immediately

$$\begin{eqnarray} \frac{\partial f}{\partial t} &+& (U_i + c_i) \frac{\partial f}{\partial x_i} + \frac{q}{m} ({\bf c} \times {\bf B})_i \frac{\partial f}{\partial c_i} = \left. \frac{\delta f}{\delta t} \right)_s\nonumber\\ & -& \left. \frac{d {\bf c} }{{dt}} \right|_c \cdot \nabla _c f - \left. \frac{d {\bf c} }{{dt}} \right|_m \cdot \nabla _c f - \frac{q}{m} \delta {\bf E} \cdot \nabla _c f\nonumber\\ &-& \nabla _v \cdot \left(\left. \frac{d {\bf c} }{{dt}} \right|_c \right) f, \end{eqnarray} \tag{ 6 }$$

where δE is the induced electric field in the moving frame. Because we are not interested in the gyromotion of the particles, we gyrophase average Equation (6), making the assumption that f(x, t, c, μ, ϕ) ≃ f(x, t, c, μ), where μ = cos θ is the cosine of the particle pitch angle.

Let b ≡ B/|B| and introduce the Cartesian basis set ${\bf e} = (\hat{\bf e}_1, \hat{\bf e}_2, \hat{\bf e}_3 = {\bf b})$. The magnetic field, no matter how complex, defines the z or $\hat{\bf e}_3$ coordinate through the directional vector b. δE₃ is the corresponding electric field fluctuation in this direction. The gradient operator is given by

$$\begin{eqnarray*} &&\frac{\partial }{\partial {\bf c} } = \hat{\bf e}_c \frac{\partial }{\partial c} + \hat{\bf e}_{\theta } \frac{1}{c} \frac{\partial }{\partial \theta } + \hat{\bf e}_{\phi } \frac{1}{c \sin \theta } \frac{\partial }{\partial \phi }, \end{eqnarray*}$$

where

$$\begin{eqnarray*} {\bf c} &=& c (\sin \theta \cos \phi \hat{\bf e}_1 + \sin \theta \sin \phi \hat{\bf e}_2 + \cos \theta {\bf b}); \\ {\bf c}_{\parallel } &=& c \cos \theta {\bf b} = c\mu {\bf b} = c_{\parallel } {\bf b}; \\ {\bf c}_{\perp } &=& c_{\perp } \left(\cos \phi \hat{\bf e}_1 + \sin \phi \hat{\bf e}_2 \right); \\ \hat{\bf e}_c &=& (\sin \theta \cos \phi \hat{\bf e}_1 + \sin \theta \sin \phi \hat{\bf e}_2 + \cos \theta {\bf b}); \\ \hat{\bf e}_{\theta } &=& (\cos \theta \cos \phi \hat{\bf e}_1 + \cos \theta \sin \phi \hat{\bf e}_2 - \sin \theta {\bf b}); \\ \hat{\bf e}_{\phi } &=& -\sin \phi \hat{\bf e}_1 + \cos \phi \hat{\bf e}_2. \end{eqnarray*}$$

The expressions for energy gain in a contracting magnetic island, Equation (2), can then be written as

$$\begin{eqnarray*} \left\langle \left. \frac{d {\bf c}_{\parallel } }{{dt}} \right|_c \cdot \nabla _c f \right\rangle _{\phi } &=& \eta _c c \mu \left(\mu \frac{\partial f}{\partial c} + \frac{ 1 - \mu ^2}{c} \frac{\partial f}{\partial \mu } \right); \\ \left\langle \left. \frac{d {\bf c}_{\perp } }{{dt}} \right|_c \cdot \nabla _c f \right\rangle _{\phi } &=& -\frac{1}{2} \eta _c c(1 - \mu ^2) \left(\frac{\partial f}{\partial c} - \mu \frac{1}{c} \frac{\partial f}{\partial \mu } \right); \\ \left\langle \nabla _c \cdot \left(\eta _c {\bf c}_{\parallel } \right) \right\rangle _{\phi }& =& \left\langle \nabla _c \cdot \left(\eta _c {\bf c}_{\perp } /2 \right) \right\rangle _{\phi } = \eta _c; \\ \left\langle \nabla _c \cdot \left(\eta _m {\bf c}_{\parallel } \right) \right\rangle _{\phi } &=& \eta _m = - \left\langle \nabla _c \cdot \left(\eta _m {\bf c}_{\perp } /2 \right) \right\rangle _{\phi }, \end{eqnarray*}$$

where 〈⋅⋅⋅〉_ϕ denotes an average over gyrophase. We therefore obtain

$$\begin{equation} \left. \left\langle \frac{d {\bf c} }{{dt}} \right|_c \cdot \nabla _c f \right\rangle _{\phi } = \eta _c c \frac{1 + \mu ^2}{2} \frac{\partial f}{\partial c} + \frac{1}{2} \eta _c (1 - \mu ^2) \mu \frac{\partial f}{\partial \mu } + 2 \eta _c f. \end{equation} \tag{ 7 }$$

Similarly, the merging magnetic island energization terms can be expressed as

$$\begin{equation} \left. \left\langle \frac{d {\bf c} }{{dt}} \right|_m \cdot \nabla _c f \right\rangle _{\phi } = \eta _m c \frac{3 \mu ^2 - 1}{2} \frac{\partial f}{\partial c} + \frac{3}{2} \eta _m (1 - \mu ^2) \mu \frac{\partial f}{\partial \mu }. \end{equation} \tag{ 8 }$$

Finally, the fluctuating electric field δE term, expressed in terms of the Cartesian basis $\delta {\bf E} = \delta E_1 \hat{e}_1 + \delta E_2 \hat{e}_2 + \delta E_3 \hat{e}_3$ becomes, after gyrophase averaging,

$$\begin{equation} \left\langle \delta {\bf E} \cdot \nabla _c f \right\rangle _{\phi } = \delta E_3 \left(\mu \frac{\partial f}{\partial c} + (1 - \mu ^2) \frac{1}{c} \frac{\partial f}{\partial \mu } \right). \end{equation} \tag{ 9 }$$

Gyrophase averaging will obviously eliminate particle energization by the fluctuating electric field components orthogonal to the b direction because the gain and loss of energy as the particle gyrates will gyrophase average to zero. Consequently, only δE₃ defined by the magnetic field direction can contribute to the energy gain or loss of the particle.

The gyrophase averaging of the left-hand side of Equation (6) proceeds in the usual way (Skilling 1975; Isenberg 1997; le Roux et al. 2007; Webb et al. 2009; see Zank 2013 for details). Including the magnetic island energization terms then yields

$$\begin{eqnarray} &&\frac{\partial f}{\partial t} + \left(U_i + c \mu b_i \right) \frac{\partial f}{\partial x_i} + \frac{1 - \mu ^2}{2}\nonumber\\\ms &&\quad \times \left[ c \nabla \cdot {\bf b} + \mu \nabla \cdot {\bf U} 3 b_i b_j \frac{\partial U_j}{\partial x_i} - \frac{2 b_i}{c} \frac{D U_i}{{dt}} \right] \frac{\partial f}{\partial \mu }\nonumber\\\ms &&\quad + \left[ \frac{1 - 3 \mu ^2 }{2} b_i b_j \frac{\partial U_j}{\partial x_i} - \frac{1 - \mu ^2}{2} \nabla \cdot {\bf U} - \frac{\mu b_i}{c} \frac{D U_i}{{dt}} \right]\nonumber\\\ms &&\quad \times c \frac{\partial f}{\partial c} = \frac{\partial }{\partial \mu } \left(\nu _s (1 - \mu ^2) \frac{\partial f}{\partial \mu } \right) - \frac{q}{m} \delta E_3\nonumber\\ &&\quad\times \left(\mu \frac{\partial f}{\partial c} + \frac{1 - \mu ^2}{c} \frac{\partial f}{\partial \mu } \right) - \left(\eta _c \frac{1 + \mu ^2}{2} + \eta _m \frac{3\mu ^2 - 1}{2} \right) c \frac{\partial f}{\partial c}\nonumber\\ &&\quad - 2 \eta _c f - \left(\frac{1}{2} \eta _c + \frac{3}{2} \eta _m \right) \mu (1 - \mu ^2) \frac{\partial f}{\partial \mu }, \end{eqnarray} \tag{ 10 }$$

where we have assumed the simplest possible diffusion form for the pitch-angle scattering term

$$\begin{eqnarray*} &&\left. \frac{\delta f}{\delta t} \right)_s = \frac{\partial }{\partial \mu } \left(\nu _s (1 - \mu ^2) \frac{\partial f}{\partial \mu } \right), \end{eqnarray*}$$

ν_s = 1/τ_s, and τ_s ≃ Ω is a characteristic pitch-angle scattering time, Ω = eB/m, the particle gyrofrequency. Here, D/Dt ≡ ∂/∂t + U_j∂/∂x_j is the convective derivative. The gyrophase-averaged Equation (10) describes the transport of energetic particles experiencing pitch-angle scattering in a turbulent plasma embedded with merging and contracting magnetic islands. The distribution function f(x, v, μ, t) is gyrotropic but not isotropic. Equation (10) therefore allows for the description of arbitrarily anisotropic distributions, including beams, etc. Such a description will undoubtedly be important in considering the possible development of firehose or Weibel instabilities, as suggested by Drake et al. (2006b, 2010); Schoeffler et al. (2011), in the reconnection regions due to the energization of particles. Several points about the gyrophase-averaged Equation (10) are noteworthy. Two particle energization terms are present (i.e., ∝c∂f/∂c), one due to the divergence of the large-scale flow velocity, ∇ · U, which is the familiar term responsible for the energization of particles at shocks, and the other due to magnetic island contraction (η_c) and magnetic field line shortening due to island merging (η_m). Notice the close correspondence between the ∇ · U and the contracting island terms. Both enter as first-order terms proportional to (1 − μ²)/2, the first describing the convergence of the bulk flow on the timescale L/U and the latter the convergence/contraction of an island on the timescale ℓ/V_A. By contrast, the energization term associated with magnetic island merging has the coefficient

$$\begin{eqnarray*} &&\eta _m \frac{3 \mu ^2 - 1}{2}, \end{eqnarray*}$$

which is proportional to the second-order Legendre polynomial. This indicates immediately that energization via island merging (neglecting for the present the role of the antireconnection electric field) will yield only second-order Fermi energization, in agreement with the result of Drake et al. (2013). Physically, it is clear that assuming conservation of adiabatic moment within the contracting magnetic island will simply redistribute energy between the parallel and perpendicular components. Bian & Kontar (2013) arrive at a similar conclusion although via a quite different approach.

Numerous simulations of test particle energization in reconnecting current sheets yield particle distributions that are essentially isotropic. We invoke this assumption and consider the Legendre polynomial expansion of the gyrophase-averaged Equation (10). This allows us to derive a simpler transport equation analogous to the energetic particle transport equation used in, e.g., cosmic ray investigations. Accordingly, we expand the gyrophase-averaged particle distribution function f as

$$\begin{eqnarray*} \hspace{-2pc}f({\bf x}, t, c, \mu) &=& \sum _{n=0}^{\infty } \frac{1}{2} (2n + 1) P_n (\mu) f_n ({\bf x}, t, c), \\ {\rm where} \quad f_n ({\bf x}, t, c) &=& \int _{-1}^1 f P_n (\mu) d\mu. \end{eqnarray*}$$

Following the derivation in Zank (2013) and neglecting flow shear and acceleration terms yields an infinite set of partial differential equations in the coefficients f_n of the Legendre polynomials. The expression is lengthy (Equation (A1)), and we relegate a sketch of the derivation and the fully expanded equation to the Appendix. Applying the Legendre polynomial expansion to the electric field term on the right-hand side of Equation (10) yields

$$\begin{eqnarray*} &&-\frac{q}{m} \frac{\delta E_3}{2m + 1} \left[ (m + 1) \frac{\partial f_{m+1} }{\partial c} + m \frac{\partial f_{m-1} }{\partial c} \right]\\ &&\quad - \frac{q}{m} \frac{\delta E_3}{c} \left[ \frac{(m + 1)(m + 2)}{2m + 1} f_{m+1} - \frac{m(m - 1)}{2m + 1} f_{m - 1} \right]. \end{eqnarray*}$$

The third term on the right-hand side of Equation (10) maps to

$$\begin{eqnarray*} &&-\frac{\eta _c }{2} \left(c\frac{\partial f_m}{\partial c} + \frac{(m - 1) m}{(2m + 1)(2m - 1)} c\frac{\partial f_{m+2} }{\partial c}\right.\\ &&\quad + \left(\!\frac{(m + 1)^2 }{2m + 3} + \frac{m^2}{2m - 1} \!\right) c\frac{\partial f_m}{\partial c} \left. {} + \frac{(m + 2)(m + 1) }{(2m + 3)(2m + 1) } c\frac{\partial f_{m+2}}{\partial c} \right)\\ &&\quad + \frac{3 \eta _m}{2} \left(\frac{1}{3} c \frac{\partial f_m }{\partial c} - \frac{(m - 1)m }{(2m + 3)(2m + 1) } c \frac{\partial f_{m-2} }{\partial c} \right. \\ &&\quad - \left(\frac{(m + 1)^2 }{2m + 3} + \frac{m^2}{2m - 1} \right) \frac{1}{2m + 1} c \frac{\partial f_m }{\partial c}\\ &&\quad -\left. \frac{(m + 2)(m + 1) }{(2m + 3)(2m + 1) } c \frac{\partial f_{m+2} }{\partial c} \right), \end{eqnarray*}$$

and the final two terms map to

$$\begin{eqnarray*} &&-2 \eta _c f_m - \left(\frac{1}{2} \eta _c + \frac{3}{2} \eta _m \right) \left(\frac{(m + 1)(m + 2)(m + 3) }{(2m + 1)(2m + 3) } f_{m+2}\right. \\ &&\quad + \frac{m (m + 1)}{(2m - 1)(2m + 3)} f_m \left. {} - \frac{m(m - 1)(m - 2) }{(2m + 1)(2m - 1) } f_{m-2} \right). \end{eqnarray*}$$

The full set of infinite partial differential equations corresponding to the Legendre expansion of the gyrophase-averaged Equation (10) is given by Equation (A1) in the Appendix, including the above terms.

To include the possible effects of second-order Fermi energization in a reduced model describing a nearly isotropic distribution function, we need to consider the f₂ approximation (i.e., consider m = 0, 1, 2 or f₀, f₁, f₂ ≠ 0 and f_n = 0 ∀ n ⩾ 3). For simplicity, we henceforth neglect the shear terms ∂U_j/∂x_i and convective derivatives DU_i/Dt—these are easily if tediously incorporated.

The m = 0 expansion yields

$$\begin{eqnarray} &&\frac{\partial f_0}{\partial t} + U_i \frac{\partial f_0}{\partial x_i} - \frac{c}{3} \frac{\partial U_i}{\partial x_i} \frac{\partial f_0}{\partial c} + \eta _c \frac{2}{3} c \frac{\partial f_0}{\partial c} - 2 \eta _c f_0\nonumber\\ &&\quad = -\frac{\partial }{\partial x_i} \left(c b_i f_1 \right) - \frac{q}{m} \frac{\delta E_3}{c^2} \frac{\partial }{\partial c} \left(c^2 f_1 \right)\nonumber\\ &&\qquad - \left[ \frac{1}{3} \frac{\partial U_i}{\partial x_i} + \frac{\eta _c}{3} + \eta _m \right] \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^3 f_2 \right), \end{eqnarray} \tag{ 11 }$$

which shows explicitly the important result that island contraction yields a first-order energization term for the isotropic or zeroth-order distribution function f₀ proportional to (c/3)∂f₀/∂c. The first of the energization terms is familiar from the transport of cosmic rays and describes energy gain or loss due to a compressive or divergent large-scale flow field (e.g., diffusive shock acceleration). The second energization term on the left results from contracting magnetic islands. The right-hand side of Equation (11) contains the higher-order moments of the expansion, f₁ and f₂. To compute these terms, subject to the closure assumption that f_m = 0 ∀ m ⩾ 3, we need to consider the next-order expansions. On considering the m = 1 terms, we find

$$\begin{eqnarray*} &&\frac{\partial f_1}{\partial t} + U_i \frac{\partial f_1}{\partial x_i} + \frac{2}{3} cb_i \frac{\partial f_2}{\partial x_i} + \frac{c}{3} b_i \frac{\partial f_0}{\partial x_i} - \frac{c}{5} \frac{\partial U_i}{\partial x_i} \frac{\partial f_1}{\partial c} + c \frac{\partial b_i}{\partial x_i} f_2\\ && \quad + \frac{1}{5} \frac{\partial U_i}{\partial x_i} f_1 = -\nu _s f_1 - \frac{q}{m} \frac{\delta E_3}{3} \left[ 2 \frac{\partial f_2}{\partial c} + \frac{\partial f_0}{\partial c} \right] - 2 \frac{q}{m} \frac{\delta E_3}{c} f_2 \\ &&\quad - \frac{4}{5} c \eta _c \frac{\partial f_1}{\partial c} - \eta _m \frac{2 c}{5} \frac{\partial f_1}{\partial c} - \frac{1}{5} \left(\eta _c + 3 \eta _m \right) f_1 - 2 \eta _c f_1. \end{eqnarray*}$$

By assuming that f₁ ≪ f₀ and f₂ ≪ f₀ and that the scattering frequency ν_s = (τ_s)⁻¹ is large (fast scattering), we obtain the approximate result

$$\begin{equation} f_1 \simeq -\frac{cb_i \tau _s}{3} \frac{\partial f_0}{\partial x_i} - \frac{q}{m} \frac{\delta E_3}{3} \tau _s \frac{\partial f_0}{\partial c}. \end{equation} \tag{ 12 }$$

Thus, the first-order correction may be approximated by (1) the diffusion description (as is well known, and the spatial diffusion coefficient is proportional to the particle scattering frequency ν_s and independent of the reconnection terms) and (2) the induced antireconnection electric field energization term.

To obtain the solution for f₂, we consider the m = 2 order expansion. The appropriate expansion yields the equation

$$\begin{eqnarray*} &&\frac{\partial f_2}{\partial t} + U_i \frac{\partial f_2}{\partial x_i} - \frac{5c}{21} \frac{\partial U}{\partial x_i} \frac{\partial f_2}{\partial c} + \frac{c}{15} \frac{\partial U}{\partial x_i} \frac{\partial f_0}{\partial c} + \frac{2}{5} cb_i \frac{\partial f_1}{\partial x_i} - \frac{c}{5} \frac{\partial b_i}{\partial x_i} f_1\\ &&\quad + \frac{1}{7} \frac{\partial U_i}{\partial x_i} f_2 = -3\nu _s f_2 - \frac{2}{5} \frac{q}{m} \delta E_3 \frac{\partial f_1}{\partial c} + \frac{q}{m} \frac{2}{5} \frac{\delta E_3}{c} f_1\\ &&\quad - \eta _c \frac{c}{15} \frac{\partial f_0}{\partial c} - \eta _c \frac{16 c}{21} \frac{\partial f_2}{\partial c} - \eta _m \frac{c}{5} \frac{\partial f_0}{\partial c} \\ &&\quad {} - \eta _m \frac{2c}{7} \frac{\partial f_2}{\partial c} - \frac{1}{7} \left[ \eta _c + 3 \eta _m \right] f_2 - 2 \eta _c f_2. \end{eqnarray*}$$

Again making the assumption of rapid scattering and the subsequent smallness of the higher-order corrections, we find that the approximate solution for f₂ is

$$\begin{equation} f_2 \simeq - \frac{1}{15} \frac{\partial U_i}{\partial x_i} \frac{c \tau _s}{3} \frac{\partial f_0}{\partial c} - \eta _c \frac{1}{15} \frac{c \tau _s}{3} \frac{\partial f_0}{\partial c} - \eta _m \frac{1}{5} \frac{c \tau _s }{3} \frac{\partial f_0}{\partial c}. \end{equation} \tag{ 13 }$$

The first term introduces diffusion in energy due to particle pitch-angle scattering, and the second and third terms are energy diffusion terms due to contracting and merging magnetic islands, respectively. Notice also that the reconnection energization terms also include the effect of particle scattering. Scattering produces two interesting effects for second-order Fermi energization. First, one can obtain second-order Fermi heating from a single island because particle pitch-angle scattering can sometimes scatter particle pitch angles close to 90°, in which case the particle will be lost to the island, and sometimes particle pitch angles will be scattered close to 0°, leading to particle trapping and therefore heating. Multiisland interaction is not therefore a requirement for second-order Fermi heating to occur. Second, in a multiisland context, particle scattering can ensure particle escape and entrance of particles into multiple contracting or merging magnetic islands and therefore facilitates a second-order Fermi energization. Notice also that the antireconnection electric field does not enter at this order.

The approximate solutions for f₁ and f₂ can be used to determine a closed transport equation in f₀. Let us introduce the notation

$$\begin{eqnarray*} &&\kappa _s \equiv \frac{c^2 \tau _s}{3}; \qquad K_{ij} \equiv b_i \kappa _s b_j, \end{eqnarray*}$$

which will describe the spatial diffusion coefficient in the final transport equation. The most general form of the transport equation for a distribution that is close to isotropic is then given by

$$\begin{eqnarray} &&\frac{\partial f_0}{\partial t} + U_i \frac{\partial f_0}{\partial x_i} - \frac{c}{3} \frac{\partial U_i}{\partial x_i} \frac{\partial f_0}{\partial c} + 2\eta _c \frac{c}{3} \frac{\partial f_0}{\partial c} + 2 \eta _c f_0\nonumber\\ &&\quad = \frac{\partial }{\partial x_i} \left(K_{ij} \frac{\partial f_0}{\partial x_j} \right) + \frac{1}{3} \frac{q}{m} \delta E_3 \frac{\partial }{\partial x_i} \left(c b_i \tau _s \frac{\partial f_0}{\partial c} \right) \nonumber \\ &&\qquad + \frac{q}{m} \frac{\delta E_3}{c^2} \frac{\partial }{\partial c} \left(\frac{c^3}{3} b_i \tau _s \frac{\partial f_0}{\partial x_i} \right) + \left(\frac{q}{m} \right)^2 \frac{\delta E_3^2}{c^2} \frac{\partial }{\partial c} \left(\frac{c^3}{3} \tau _s \frac{\partial f_0}{\partial c} \right) \nonumber \\ &&\qquad + \frac{1}{5} \left(\frac{1}{3} \frac{\partial U_i}{\partial x_i} + \frac{\eta _c}{3} + \eta _m \right)^2 \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^2 \kappa _s \frac{\partial f_0}{\partial c} \right). \end{eqnarray} \tag{ 14 }$$

We can rewrite Equation (14) slightly by introducing a scattered antireconnection electric-field-induced velocity

$$\begin{equation} V_E \equiv \frac{1}{3} \frac{q}{m} \delta E_3 \tau _s, \end{equation} \tag{ 15 }$$

which has a sign dependence because of the charge. We use V_Ei = −|V_E|b_i ≡ |V_E|_i and assume that (q/m)δE₃ is independent of x to obtain

$$\begin{eqnarray} &&\frac{\partial f_0}{\partial t} + \left(U_i + 3 |V_E |_i \right) \frac{\partial f_0}{\partial x_i} - \frac{c}{3} \frac{\partial U_i}{\partial x_i} \frac{\partial f_0}{\partial c} + \frac{1}{c^2} \frac{\partial }{\partial c} \left(\frac{c^3}{3} 2 \eta _c f_0 \right) \nonumber\\ &&\quad = \frac{\partial }{\partial x_i} \left(K_{ij} \frac{\partial f_0}{\partial x_j} \right) - c \frac{\partial }{\partial x_i} \left(|V_E |_i \frac{\partial f_0}{\partial c} \right) - |V_E |_i c \frac{\partial ^2 f_0}{\partial x_i \partial c}\nonumber\\ &&\qquad + 3\frac{V_E^2}{\tau _s} \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^2 \frac{\partial f_0}{\partial c} \right) + \frac{1}{5} \left(\frac{1}{3} \frac{\partial U_i}{\partial x_i} + \frac{\eta _c}{3} + \eta _m \right)^2\nonumber\\ &&\qquad\times \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^2 \kappa _s \frac{\partial f_0}{\partial c} \right). \end{eqnarray} \tag{ 16 }$$

In many respects, the transport Equation (16) is very similar to the standard Parker–Axford–Gleeson transport equation, elaborated further below. Equation (16) illustrates that particles propagating within a region containing dynamically merging magnetic islands and antireconnection electric fields can experience first-order energization or cooling. Related transport equations that include the effects of fluctuating electric fields have been derived previously by le Roux et al. (2002); Kichatinov (1983); Fedorov et al. (1992); Dorman et al. (1987), although not related to the concept of merging magnetic islands.

Equation (14) can be rewritten in a particularly revealing and useful phase space conservation form,

$$\begin{equation} \frac{\partial f_0}{\partial t} + \frac{\partial S_i }{\partial x_i} + \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^2 J_p \right) + \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^2 D_{cc} \frac{\partial f_0}{\partial c} \right) = 0, \end{equation} \tag{ 17 }$$

where

$$\begin{eqnarray*} &&{\bf S} \equiv - {\bf K} \cdot \nabla f_0 - \frac{c}{3} \frac{\partial f_0}{\partial c} \left({\bf U} - 3 |{\bf V}_E |\right), \end{eqnarray*}$$

is the energetic particle streaming in space,

$$\begin{eqnarray*} &&J_p \equiv \frac{c}{3} \left({\bf U} + |{\bf V}_E |\right) \cdot \nabla f_0 + \frac{c}{3} 2 \eta _c f_0, \end{eqnarray*}$$

is the streaming in momentum space, and the diffusion coefficient in velocity space is given by

$$\begin{eqnarray*} &&D_{cc} = \kappa _s \left[ 3 |V_E|^2 + \frac{1}{5} \left(\frac{1}{3} \frac{\partial U_i}{\partial x_i} + \frac{\eta _c}{3} + \eta _m \right)^2 \right]. \end{eqnarray*}$$

By estimating δE₃ ∼ V_AB and introducing a characteristic particle speed c₀, we can order the terms in Equation (16) using $\bar{c} \equiv c/c_0$. The second and third terms on the right-hand-side scale as $\Omega V_A \tau _s \bar{c} \partial ^2 f_0 /\partial x_i\partial c$ and $\Omega V_A \tau _s \bar{c}^2 \partial /\partial \bar{c} \left(\bar{c}^3 b_i \partial f_0/\partial x_i \right)$, whereas the fourth term scales as $\Omega ^2 \left(V_A^2 / c_0^2 \right) \tau _s \bar{c}^2 \partial /\partial \bar{c} \left(\bar{c}^2 \partial f_0 /\partial \bar{c} \right)$. Thus, the first-order correct transport equation for $V_A^2 / c_0^2 \ll 1$ is

$$\begin{eqnarray} &&\frac{\partial f_0}{\partial t} + \left(U_i + 3 |V_E | b_i \right) \frac{\partial f_0}{\partial x_i} - \frac{c}{3} \frac{\partial U_i}{\partial x_i} \frac{\partial f_0}{\partial c} + \frac{1}{c^2} \frac{\partial }{\partial c} \left(\frac{c^3}{3} 2 \eta _c f_0 \right)\nonumber\\ &&\quad = \frac{\partial }{\partial x_i} \left(K_{ij} \frac{\partial f_0}{\partial x_j} \right) - c \frac{\partial }{\partial x_i} \left(|V_E |_i \frac{\partial f_0}{\partial c} \right) - |V_E |_i c \frac{\partial ^2 f_0}{\partial x_i \partial c}.\nonumber\\ && \end{eqnarray} \tag{ 18 }$$

Evidently, only island contraction and the induced antireconnection electric field associated with island merging are important for particle energization at the leading order. The transport equation that describes energization by contracting islands alone (i.e., η_c ≠ 0, η_m = 0, δE₃ = 0) is given by

$$\begin{eqnarray} &&\frac{\partial f_0}{\partial t} + U_i \frac{\partial f_0}{\partial x_i} - \frac{c}{3} \frac{\partial U_i}{\partial x_i} \frac{\partial f_0}{\partial c} + \frac{1}{c^2} \frac{\partial }{\partial c} \left(\frac{c^3}{3} 2 \eta _c f_0 \right) \nonumber\\ &&\quad = \frac{\partial }{\partial x_i} \left(K_{ij} \frac{\partial f_0}{\partial x_j} \right) + \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^2 \frac{\kappa _s }{45} \left(\frac{\partial U_i}{\partial x_i} + \eta _c \right)^2 \frac{\partial f_0}{\partial c} \right),\nonumber\\ && \end{eqnarray} \tag{ 19 }$$

where the second term on the right-hand side of Equation (19) describes second-order Fermi energization associated with both the flow divergence and magnetic island contraction. Particle diffusion, expressed through the spatial diffusion coefficient κ_s enables particles to experience on average more compressive/contraction energy gains than losses. By contrast, in the presence of merging alone (i.e., neglecting both magnetic island contractions η_c = 0 and the antireconnection electric field δE₃ = 0), the transport equation reduces to

$$\begin{eqnarray} &&\frac{\partial f_0}{\partial t} + U_i \frac{\partial f_0}{\partial x_i} - \frac{c}{3} \frac{\partial U_i}{\partial x_i} \frac{\partial f_0}{\partial c} = \frac{\partial }{\partial x_i} \left(K_{ij} \frac{\partial f_0}{\partial x_j} \right)\nonumber\\ &&\quad + \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^2 \frac{\kappa _s}{45} \left(\frac{\partial U_i}{\partial x_i} + \eta _c \right)^2 \frac{\partial f_0}{\partial c} \right), \end{eqnarray} \tag{ 20 }$$

showing that the reconnection through merging contributes to particle energization only through a second-order Fermi process.

Finally, retaining only the antireconnection electric field and setting η_c = η_m = 0 yields

$$\begin{eqnarray} &&\frac{\partial f_0}{\partial t} + \left(U_i + 3 |V_E | b_i \right) \frac{\partial f_0}{\partial x_i} - \frac{\partial U_i}{\partial x_i} \frac{c}{3} \frac{\partial f_0}{\partial c} = \frac{\partial }{\partial x_i} \left(K_{ij} \frac{\partial f_0}{\partial x_j} \right) \nonumber \\ &&\quad - c \frac{\partial }{\partial x_i} \left(|V_E |_i \frac{\partial f_0}{\partial c} \right) - |V_E|_i c \frac{\partial f_0^2}{\partial x_i \partial c} + \frac{1}{c^2} \frac{\partial }{\partial c} \left(c^2 D^E_{cc} \frac{\partial f_0}{\partial c} \right),\nonumber\\ \end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray*} &&D^E_{cc} \equiv \kappa _s \left(\frac{1}{45} \left(\frac{\partial U_i}{\partial x_i} \right)^2 + 3 |V_E|^2 \right), \end{eqnarray*}$$

is the second-order velocity diffusion coefficient.

Were we to neglect the second-order Fermi acceleration terms on the right-hand side of Equation (16), the transport equation would be very similar to the standard Gleeson–Axford–Parker transport equation used in the study of cosmic ray transport (e.g., Zank 2013) except for the reconnection terms. Equation (16) illustrates that particles propagating within a region, such as a reconnecting current sheet that produces numerous merging and contracting magnetic islands and antireconnection electric fields, can experience first-order and second-order energization. The second-order energization of particles relies on the presence of scattering (the second Fermi acceleration term is the product of both reconnection and spatial diffusion), which allows the particle to scatter out of an energization region before losing energy.

Finally, we mention that we have also derived the corresponding focused transport equation for relativistic charged particles and its isotropic, advection–diffusion counterpart. The equations are structurally identical except for a relativistic factor in the time derivative and the use of the momentum rather than the velocity variable.

Let us consider perhaps the simplest of problems: an incompressible super-Alfvénic flow, ∇ · U = 0, with a steady injection of particles of fixed initial speed c₀ at the origin, and consider a 1D problem. We further assume for analytic convenience that the spatial diffusion coefficient is constant, and we neglect the second-order energy diffusion terms. We will consider several examples to illustrate the effect of the different energization mechanisms.

4.1. Particle Acceleration by the Antireconnection Electric Field Alone

After rearranging (21) slightly and discarding the second-order velocity diffusion term, introducing ξ = ln (c/c₀) yields the governing equation as (neglecting the "0" subscript on f₀)

$$\begin{equation} \frac{\partial ^2 f}{\partial x^2} - \frac{2 |V_E | }{\kappa } \frac{\partial ^2 f}{\partial x \partial \xi } - \frac{U + 3 |V_E | }{\kappa } \frac{\partial f}{\partial x} = - \frac{n_0 \delta (c - c_0) }{4\pi \kappa c^2} c_0 \delta (x), \end{equation} \tag{ 22 }$$

where all the coefficients are constant, including the particle number density n₀. The principal part of the differential operator in Equation (22) shows that Equation (22) is hyperbolic with characteristic curves ξ⁺ = −2(|V_E|/κ)x + const. and ξ⁻ = const. To solve (22), we introduce a Laplace transform in the variable ξ, $\bar{f} (x,s) = {\cal L} [f(x,\xi)]$, which reduces (22) to the second-order nonhomogeneous ordinary differential equation,

$$\begin{eqnarray} &&\frac{d^2 \bar{f} }{dx^2} - \frac{g(s) }{\kappa } \frac{d \!\bar{f} }{dx} = -\frac{Q}{\kappa } \delta (x); \\ &&g(s) = U + 3 |V_E | + 2 |V_E | s = 2 |V_E | (s - s_0); \nonumber \\ &&s_0 = -(3 + U/|V_E |)/2; \qquad {} \qquad Q = \frac{n_0}{4\pi c_0^2}. \nonumber \end{eqnarray} \tag{ 23 }$$

The complementary solution to Equation (23) is

$$\begin{eqnarray*} \bar{f}\! (x,s) = \left\lbrace \begin{array}{ll}A, & x > 0 \\ A e^{g(s) s / \kappa }, & x < 0 \end{array} \right.\!\!, \end{eqnarray*}$$

after demanding $\bar{f} (x,s)$ be bounded as x → +∞, assuming that $\bar{f} (x,s) \rightarrow 0$ as x → −∞, and assuming continuity of $\bar{f} (0,s)$ i.e., $\bar{f} (-\varepsilon,s) = \bar{f} (\varepsilon,s)$ as ε → 0. The constant A can be obtained by integrating Equation (23) over the interval [ − ε, ε] and taking the limit ε → 0, which yields the jump condition

$$\begin{eqnarray*} &&\left[ \frac{d \bar{f} }{dx} \right] = -\frac{Q}{\kappa }, \end{eqnarray*}$$

where [F] ≡ F_{0 −} − F_{0 +}. Evaluating the derivatives across x = 0 yields A = Q/g(s), or

$$\begin{equation} \bar{f} (x,s) = \left\lbrace \begin{array}{ll}\displaystyle \frac{Q}{g(s)}, & x > 0 \\ \displaystyle \frac{Q}{g(s)} e^{g(s) s / \kappa }, & x < 0 \end{array} \right.. \end{equation} \tag{ 24 }$$

The solution ahead (x < 0) of the particle injection point decays exponentially from the source with a scale length determined essentially by the spatial diffusion coefficient κ and the bulk velocity U and is constant behind (x > 0) the source.

The inversion of $\bar{f} (x,s)$ for x > 0 is straightforwardly obtained by the method of residues for the simple pole at s = s₀. For x < 0, some care needs to be exercised in choosing the contour in the complex s plane because one needs to ensure the convergence of the integral

$$\begin{eqnarray*} &&\frac{1}{2 \pi i} \int _{a - i \infty }^{a + i \infty } \frac{\exp \left[ \left(\xi + 2 |V_E | x /\kappa \right) s \right] }{s - s_0} ds. \end{eqnarray*}$$

Only if ξ + 2|V_E|x/κ > 0 can an enclosed contour containing s₀ be chosen, whereas if ξ + 2|V_E|x/κ < 0, the closed contour cannot include s₀. When the latter condition holds, f(x, c) = 0. The solution for x < 0 is then obtained, yielding the full solution for c > c₀ as

$$\begin{eqnarray} &&f (x,c) \nonumber\\ &&=\! \left\lbrace\!\! \begin{array}{ll} \displaystyle \frac{n_0}{8\pi c_0^2 |V_E |} \left(\!\frac{c}{c_0}\! \right)^{-(3 + U/|V_E)/2} H(c - c_0), &\! x > 0 \\ \displaystyle \frac{n_0}{8\pi c_0^2 |V_E |} \left(\!\frac{c}{c_0}\! \right)^{-(3 + U/|V_E)/2} H(\xi + 2 |V_E | x /\kappa), & \!x < 0 \end{array} \right.\!\!\!\!,\quad \end{eqnarray} \tag{ 25 }$$

where H(x) is the Heaviside step function. The solution for x < 0 essentially reflects "causality" because recall that ξ + 2|V_E|x/κ is one of the characteristic variables of the hyperbolic Equation (22), and the step function implies ln (c/c₀) + 2|V_E|x/κ > 0. The solution for c > c₀ is a simple power law in particle speed. By integrating Equation (25) over c, one can obtain the spatial distribution of the number density n for the energized particles, thus

$$\begin{equation} n(x) = \left\lbrace\!\! \begin{array}{ll} \displaystyle \frac{n_0}{8\pi c_0 \left(U + |V_E | \right) }, & x > 0 \\ \displaystyle \frac{n_0}{8\pi c_0 \left(U + |V_E | \right) } \exp \left[ \frac{ \left(U + |V_E | \right) x }{\kappa } \right] & x < 0 \end{array} \right.\!\!. \end{equation} \tag{ 26 }$$

The number density is constant downstream of the particle source and decays exponentially upstream of the source. The latter is the result of particles having to diffuse upstream against the flow.

The solution for c < c₀ can be derived directly from the analysis above by noting that defining |ξ| = |ln (c/c₀)| yields the same Equation (22) if we use ξ↦|ξ|. The corresponding solution for c < c₀ then follows from Equation (25) if we substitute |ln (c/c₀)| for ξ. Thus, for particle speeds c < c₀, we obtain

$$\begin{eqnarray} &&f (x,c) \nonumber\\ &&=\! \left\lbrace\!\!\! \begin{array}{ll} \displaystyle \frac{n_0}{8\pi c_0^2 |V_E |} \left(\!\frac{c}{c_0}\! \right)^{(3 + U/|V_E)/2} H(c_0 - c), &\! x > 0 \\ \displaystyle \frac{n_0}{8\pi c_0^2 |V_E |} \left(\!\frac{c}{c_0} \!\right)^{(3 + U/|V_E)/2} H(|\xi | + 2 |V_E | x /\kappa), &\! x < 0 \end{array} \right.\!\!\!\!\!, \nonumber\\ \end{eqnarray} \tag{ 27 }$$

which when combined with (25) yields a (logarithmically) triangular particle distribution bounded by two symmetric power laws. The spatial dependence of the number density for cooled particles is similar to that of Equation (26). A sketch of the solution is illustrated in Figure 3. The step function in Equation (27) shows that 2|V_E|x/κ > ln (c/c₀), illustrating that some particles experience deceleration (cooling) through interacting with the reconnection electric field. This is not surprising because particles can as easily lose energy as gain energy in the presence of a direct electric field for an isotropic distribution of charged particles. This effect is contained in the transport equation that we derived for the case of a nearly isotropic distribution. The power-law solutions (25) and (27) yield the complete solution to the stationary Equation (22), and the normalization constant shows that half of the injected particles are accelerated, and half are decelerated. Thus, in terms of net acceleration, there is none because half the particles are cooled, and the other half are energized symmetrically for an isotropic distribution; there is, however, an accelerated component that has energies greater than c₀, and this component admits a power-law distribution for accelerated ions or electrons. Any approach to accelerating particles using a direct electric field, whether simulations, as in Oka et al. (2010), or theoretically, as in this work, must always result in a distribution of heated and cooled particles.

Zoom In Zoom Out Reset image size

The power-law index can be estimated from

$$\begin{eqnarray*} &&|V_E | = \frac{ |q| }{m} |\delta E_3| \tau _s \sim \frac{ |q| }{m} \frac{V_A B}{\Omega } = V_A, \end{eqnarray*}$$

so that the index (3 + U/|V_E|)/2 ∼ (3 + M_A)/2, where M_A = U/V_A is the Alfvén Mach number. For Alfvén Mach numbers M_A < 5, the spectral index is hard, being less than 4.

4.2. Particle Acceleration by Magnetic Island Contraction Alone

The original discussion of particle acceleration by magnetic island reconnection (Drake et al. 2006a) focused on the effects of magnetic island contraction exclusively. The subsequent numerical simulations by Oka et al. (2010) separated the contributions of the various energization mechanisms associated with particle acceleration. To isolate particle acceleration by magnetic island contraction, we consider the solution of the transport Equation (19), i.e., only η_c ≠ 0, in the limit of an incompressible flow. As with the previous example, we introduce ξ = ln (c/c₀) to obtain (neglecting the 0 subscript on f₀)

$$\begin{equation} \frac{\partial ^2 f}{\partial x^2} - \frac{U}{\kappa } \frac{\partial f}{\partial x} - \frac{2\eta _c}{3 \kappa } \frac{\partial f}{\partial \xi } - \frac{2\eta _c}{\kappa } f = - \frac{Q}{\kappa } \delta (c - c_0) \delta (x), \end{equation} \tag{ 28 }$$

where Q = n₀/4πc². Unlike Equation (22), Equation (28) is strictly parabolic. The Laplace transform in ξ, $\bar{f} (x,s) = {\cal L}[f(x,s)]$, yields a second-order linear nonhomogeneous ordinary differential equation,

$$\begin{eqnarray} \frac{d^2 \bar{f} }{dx^2} - \frac{U}{\kappa } \frac{d\! \bar{f} }{{dx}} - g(s) \bar{f} &=& -\frac{Q}{\kappa } \delta (x); \\ g(s) &=& \frac{2 \eta _c}{3 \kappa } (s + 3); \nonumber \\ Q &=& \frac{n_0}{4\pi c_0^2}. \nonumber \end{eqnarray} \tag{ 29 }$$

Equation (29) is solved straightforwardly using the method of variation of parameters. The complementary solutions of Equation (29) are

$$\begin{eqnarray} y_1 (x) &=& e^{\lambda _1 x}, \qquad y_2 (x) = e^{\lambda _2 x}; \\ \lambda _{1,2}& \equiv & \frac{ U/\kappa \pm \sqrt{\Delta } }{2}, \quad \sqrt{\Delta } = \left[ \left(\frac{U}{\kappa } \right)^2 + \frac{8 \eta _c}{3\kappa } (s + 3) \right]^{1/2}\nonumber\\ & =& \lambda _1 - \lambda _2. \nonumber \end{eqnarray} \tag{ 30 }$$

On noting that y₁ → 0 and y₂ → ∞ as x → −∞, and y₁ → ∞ and y₂ → 0 as x → ∞, the solution of Equation (29) is given by

$$\begin{eqnarray*} &&\bar{f} (x,s) = \frac{n_0}{4\pi c_0^2 \kappa } \frac{1}{\lambda _1 - \lambda _2} \left[ e^{\lambda _1 x} \left(1 - H(x) \right) + e^{\lambda _2 x} H(x) \right], \end{eqnarray*}$$

or, equivalently,

$$\begin{equation} \bar{f}(x,s) = \frac{n}{8\pi c_0^2 } \sqrt{\frac{3}{2\kappa \eta _c} } \frac{1}{\sqrt{s - s_1} } \left\lbrace \begin{array}{ll}e^{\lambda _1 x}, & x > 0 \\ e^{\lambda _2 x}, & x < 0 \end{array} \right.\!\!\!, \end{equation} \tag{ 31 }$$

after expressing $\kappa (\lambda _1 - \lambda _2) = \kappa \sqrt{\Delta } = \sqrt{8\kappa \eta _c /3} \sqrt{s - s_1}$, where

$$\begin{eqnarray*} &&s_1 = -3 \left(1 + \frac{U^2}{8\kappa \eta _c} \right). \end{eqnarray*}$$

On rewriting (31) slightly, we can use a tabulated Laplace integral transform (Erdélyi et al. 1954, p. 246. formula (6)) to obtain

$$\begin{eqnarray} f(x,c/c_0) &=& \frac{n_0}{8\pi c_0^2} \sqrt{\frac{3}{2\kappa \eta _c} } \left(\frac{c}{c_0} \right)^{-3 (U^2/(8\kappa \eta _c) + 1)} \nonumber \\ &&\times \frac{1}{\sqrt{\pi \ln (c/c_0)} } \exp \left[ \left(\frac{U}{2\kappa } - \frac{\eta _c}{6\kappa \ln (c/c_0) } x \right) x \right]\nonumber\\ &&\times H \left(|\ln (c/c_0) | \right]). \end{eqnarray} \tag{ 32 }$$

It is instructive to express the solution (32) in terms of the diffusive timescale τ_diff ≡ κ/U², the magnetic island contraction timescale $\tau _c \equiv \eta _c^{-1}$, and the diffusion length scale L_diff ≡ κ/U. The distribution function (32) for the accelerated particles then becomes

$$\begin{eqnarray} f(x,c/c_0) &=& \frac{n_0}{8\pi c_0^2 U} \sqrt{\frac{3\tau _c}{2\tau _{{\rm diff}} } } \left(\frac{c}{c_0} \right)^{-3 (\tau _c/(8\tau _{{\rm diff}}) + 1)} \frac{1}{\sqrt{\pi \ln (c/c_0)} }\nonumber\\ && \times \exp \left[ \left(\frac{1}{2} - \frac{1}{6 \ln (c/c_0) } \frac{\tau _{{\rm diff}} }{\tau _c} \frac{x}{L_{{\rm diff}} } \right) \frac{x}{L_{{\rm diff}} } \right]\nonumber\\ &&\times H \left(|\ln (c/c_0) | \right]). \end{eqnarray} \tag{ 33 }$$

Evidently, solution (33) is essentially a power law for c > c₀, which is consistent with the numerous test particle simulations cited above that yield isotropic power-law distribution functions for the energized particles. However, it is apparent that the underlying microphysics of the energization process in contracting magnetic islands and the physics of charged particle pitch-angle scattering plays a critical role in determining the actual power-law exponent. If the characteristic magnetic island contraction time τ_c ≪ 8τ_diff, the power-law spectral index is very hard with f ∼ c⁻³. Plotted in Figure 4 are several solutions (33), of which the hardest spectrum corresponds to assuming a ratio of τ_c/τ_diff = 0.5 at x/L_diff = 1. For c > c₀, the distribution is essentially a power law with slope ∼−3. Also illustrated in Figure 4 are three examples of the solution to Equation (33) for different choices of the ratio τ_c/τ_diff = 5, 7, 9. These slopes are steeper and correspond to the power-law index of Equation (33).

Zoom In Zoom Out Reset image size

Precise estimates of both τ_c and τ_diff are beyond the scope of this paper, but order-of-magnitude estimates are possible. We may take τ_c = (η_c)⁻¹ = ℓ/V_A (Drake et al. 2006a; Bian & Kontar 2013; neglecting a possible factor of b²/B²). Considering the relation τ_c ⩽ 8τ_diff or ℓ/V_A ⩽ 8c²τ_s/(3U²) and rewriting τ_s = λ_mfp/c, where λ_mfp is the particle scattering mean free path, yields

$$\begin{equation} \frac{\ell }{V_A} \le \frac{8c \lambda _{mfp} }{3U^2} \quad {\rm or} \quad \frac{3M_A}{8} \frac{U}{c} \le \frac{\lambda _{mfp} }{\ell }. \end{equation} \tag{ 34 }$$

The inequality (34) delineates the conditions under which a relatively hard power-law index (α = −3 to say −5) can be expected. The inequality (34) is easily met for typical quiet-time solar wind conditions because the left-hand side is of O(1). Of course, for particle energization to occur relatively efficiently via the magnetic island contraction mechanism, particle scattering within the island should be relatively weak, and so the condition ℓ ≪ λ_mfp is implicitly assumed. We note that nonlinear particle feedback effects might well damp the magnetic islands' ability to accelerate particles, leading to a steeper spectrum than suggested here.

Acceleration by the island contraction mechanism can therefore be expected to produce rather hard power-law spectra, as illustrated in Figure 4. However, unlike the induced or antireconnection electric field case above or diffusive shock acceleration, the properties of the power-law solution are not determined from macroscopic flow boundary conditions but rather by the detailed microphysics of charged particle scattering and magnetic island contraction.

4.3. Complete First-order Correct Solution

Having identified the individual contributions of the first-order antireconnection electric field and magnetic island contraction to the energization of particles in regions of reconnecting current sheets, we consider briefly solutions to the full first-order correct steady-state transport Equation (18) for an incompressible flow. The steady-state 1D equation may be expressed as (neglecting the 0 subscript again)

$$\begin{eqnarray} &&\frac{\partial ^2 f}{\partial x^2} - 2 \frac{|V_E|}{\kappa } \frac{\partial ^2 f}{\partial x \partial \xi } - \frac{U + 3 |V_E|}{\kappa } \frac{\partial f}{\partial x} - \frac{2\eta _c}{3\kappa } \frac{\partial f}{\partial \xi }- \frac{2\eta _c}{\kappa } f\nonumber\\ &&\quad = -\frac{Q}{\kappa } \delta (c - c_0) \delta (x), \end{eqnarray} \tag{ 35 }$$

where the various terms have been defined above. Like the anti-reconnection electric field case, Equation (35) is a hyperbolic equation with the same characteristic curves. The discussion of Section 4.1 therefore carries over. We follow the analysis of Section 4.2, Laplace transforming (35) to obtain the second-order, nonhomogeneous linear ordinary differential equation

$$\begin{eqnarray} &&\frac{d^2 \bar{f} }{dx^2} - \frac{g(s) }{\kappa } \frac{d \bar{f} }{dx} - \frac{h(s) }{\kappa } \bar{f} = -\frac{Q}{\kappa } \delta (x); \\ &&g(s) \equiv U + 3 |V_E| + 2 |V_E|s; \nonumber \\ &&h(s) \equiv \frac{2\eta _c}{3} (s + 3), \nonumber \end{eqnarray} \tag{ 36 }$$

to obtain

$$\begin{eqnarray} \bar{f} (x,s) &=& \frac{n}{4\pi c_0^2} \frac{1}{\kappa \sqrt{\Delta } } \left\lbrace \begin{array}{ll}e^{\lambda _1 x}, & x > 0 \\ e^{\lambda _2 x}, & x < 0 \end{array} \right.\!\!\!, \\ \lambda _{1,2} &\equiv & \frac{g(s) /\kappa \pm \sqrt{\Delta } }{2}; \nonumber \\ \kappa \sqrt{\Delta } &\equiv & 2|V_E| \left[\! (s \!-\! s_1)^2 \!-\! \frac{\eta _c \kappa }{3 |V_E|^3 } \left(\!U - 3|V_E| + \frac{\eta _c \kappa }{3 |V_E|}\! \right)\! \right]^{1/2}; \nonumber \\ s_1 &\equiv & -\frac{1}{2} \left(\frac{U + 3|V_E|}{|V_E| } + \frac{2\eta _c \kappa }{3 |V_E|^2 } \right). \nonumber \end{eqnarray} \tag{ 37 }$$

By using formula (36), p. 249 of Erdélyi et al. (1954), the inverse Laplace transform of Equation (37) yields the solution

$$\begin{eqnarray} f(x,c/c_0) &=& \frac{n_0}{8\pi c_0^2 |V_E| } \left(\frac{c}{c_0} \right)^{-(3 + U/|V_E| + 2\eta _c \kappa /(3|V_E|^2)/2}\nonumber\\ &&\times \exp \left[ -\frac{2\eta _c }{3|V_E| } x \right] I_0 \left(\Phi \right) H\left[\ln (c/c_0) \right]\nonumber\\ &&\times H \left[ \ln (c/c_0) + \frac{2 |V_E| }{\kappa x} \right]; \\ \Phi &\equiv& \sqrt{ \frac{\eta _c \kappa }{3 |V_E|^3 } \left(U - 3|V_E| + \frac{\eta _c \kappa }{3|V_E| } \right) } \nonumber\\ &&\times\left(\left(\ln (c/c_0) \right)^2 + \frac{4 |V_E| }{\kappa } \ln (c/c_0) x \right)^{1/2}, \nonumber \end{eqnarray} \tag{ 38 }$$

where I₀ is the modified Bessel function of order 0. Note again the presence of the characteristic curves ξ⁺ = −2|V_E|/κ + const. and ξ⁻ = const. We do not repeat the discussion about accelerated and cooled particles. Equation (38) describes the particles that have acquired speeds c > c₀, and the power-law-like character of the solution is evident. As before, the nature of the solution (38) is most clearly revealed when expressed in terms of characteristic scales and the Alfvén Mach number. Hence

$$\begin{eqnarray} f(x,c/c_0)& =& \frac{n_0}{8\pi c_0^2 V_A } \left(\frac{c}{c_0} \right)^{-(3 + M_A + 2\tau _{{\rm diff}} /(3\tau _c) M_A^2)/2}\nonumber\\ &&\times \exp \left[ -\frac{2\tau _{{\rm diff}} }{3\tau _c} M_A \frac{x}{L_{{\rm diff}} } \right] I_0 \left(\Phi \right) H \left[ \ln (c/c_0) \right]\nonumber \\ &&\times H \left[ \ln (c/c_0) + \frac{2}{M_A} \frac{x}{L_{{\rm diff}} } \right]; \\ \Phi &=& \sqrt{ \frac{\tau _{{\rm diff}} }{3\tau _c} M_A^2 \left(M_A - 3 + \frac{\tau _{{\rm diff}} }{3\tau _c} M_A^2 \right) }\nonumber\\ &&\times \left((\ln (c/c_0)^2 + \frac{4}{M_A} \ln (c/c_0) \frac{x}{L_{{\rm diff}} } \right)^{1/2}. \nonumber \end{eqnarray} \tag{ 39 }$$

Solutions (39), i.e., for c > c₀, are plotted in Figure 5. Figure 5 (top) assumes M_A = 7, and Figure 5 (bottom) uses M_A = 11, together with the four values of τ_c/τ_diff used in Figure 4. The solutions for both values of M_A exhibit power laws for virtually all values of c > c₀, and the slope depends on the choice of M_A and τ_c/τ_diff. However, it is interesting that the power-law exponent α in Equation (39) for (c/c₀)^−α does not accurately represent the slope of the solutions plotted in Figure 5. The parameters used in the plot would give an exponent α much steeper than plotted. Instead, the modified Bessel function has a positive slope and contributes to the (c/c₀) distribution by flattening the overall slope. The smaller Mach number plots (Figure 5, top) have harder spectra, with power-law exponents ranging from ∼ − 3.2, −3.53, −3.63, to ∼ − 3.7, compared to the higher Mach number plots (Figure 5 bottom), where the exponents now range from ∼ − 3.2, −3.8, and −4, to ∼ − 4.1 for τ_c/τ_diff = 0.5, 5, 7, and 9, respectively. The argument of the modified Bessel's function (39) is large, allowing it to be approximated as

$$\begin{eqnarray*} &&I_0 (x) \simeq \frac{e^x}{\sqrt{2\pi x} }. \end{eqnarray*}$$

On neglecting the factor of 2π and assuming that aln ²(c/c₀) ≫ abln (c/c₀), where $a \equiv (M_A^2 /3 \bar{\tau }) (M_A - 3 + M_A^2 /3 \bar{\tau })$, $\bar{\tau } \equiv \tau _c /\tau _{{\rm diff}}$, b ≡ 4(x/L_diff)/M_A, yields

$$\begin{eqnarray} f &\sim& \left(\frac{c}{c_0} \right)^{-(3 + M_A + 2\tau _{{\rm diff}} /(3\tau _c) M_A^2)/2} \nonumber\\ &&\times\left(\frac{c}{c_0} \right)^{\sqrt{ M_A^2 /(3\bar{\tau }) (M_A - 3 + \tau _{{\rm diff}} /(3\tau _c) M_A^2) } - 0.14}. \end{eqnarray} \tag{ 40 }$$

Here the factor 0.14 results from approximating the term $1/\sqrt{\ln (c/c_0)}$ in the range c/c₀ ∈ [10, ..., 100]. It is easily verified that the expression (40) yields virtually identical values for the power-law slopes derived from the plots of Figure 5.

Zoom In Zoom Out Reset image size

The solutions given in Sections 4.1–4.3 are appropriate only locally in the solar wind because we neglected the divergence term in the transport Equation (18). However, we can estimate quite easily the relative importance of the divergence term compared to the acceleration timescales. The 1D time-dependent transport Equation (18) is simply

$$\begin{eqnarray*} &&\frac{\partial f_0}{\partial t} + \left(U + 3|V_E| \right) \frac{\partial f_0}{\partial x} - \frac{c}{3} \frac{\partial U}{\partial x} \frac{\partial f_0}{\partial c} + \frac{2\eta _c}{3} c \frac{\partial f_0}{\partial c} + 2\eta _c f_0\\ &&\quad = \frac{\partial }{\partial x} \left(K \frac{\partial f}{\partial x} \right) - 2 |V_E | c \frac{\partial ^2 f}{\partial x \partial c}. \end{eqnarray*}$$

If we neglect the spatial diffusion term and introduce a characteristic scale length ℓ_E so that the last term may be approximated as −2|V_E|/ℓ_Ec∂f/∂c, the transport equation may be approximated as a first-order quasi-linear equation. One of the characteristics is then given by

$$\begin{equation} \frac{1}{c} \frac{d c}{{dt}} = -\frac{1}{3} \frac{\partial U}{\partial x} + 2\frac{\eta _c}{3} + 2 \frac{|V_E |}{\ell _E}, \end{equation} \tag{ 41 }$$

which is the inverse of the particle acceleration timescale. The last two terms are the inverse acceleration timescales due to magnetic island contraction and the antireconnection-induced electric field, respectively. To understand the physical meaning of ℓ_E, we express the length scale in terms of a timescale that is related to the duration of the magnetic island merging event (because this is the time available to a trapped particle to experience repeated encounters with δE₃); thus ℓ_E ∼ V_Aτ_m, where τ_m is the time taken for the merging of two adjacent magnetic islands. If ℓ is the characteristic size of a magnetic island, then τ_m ∼ ℓ/V_A. Thus, because η_c ∼ V_A/ℓ and |V_E| ∼ V_A, the acceleration rate due to magnetic island merging (both island contraction and induced electric field) is $\tau _{{\rm acc}}^{-1} \simeq 2 V_A/(3 \ell) + 2V_A /\ell = 8V_A /(3\ell)$. We can now compare the size of the particle acceleration rate to that of the expansion or cooling term ∂U/∂x ≃ U/R, where R is the solar wind spatial scale ∼1 AU. If we assume that ℓ is less than the correlation length scale at 1 AU (i.e., ∼0.01 AU), then say $\tau _{{\rm acc}}^{-1} \sim 3 (50\, {\rm km\, s^{-1}}) / (10^{-3} \,{\rm AU})$ and U/R ∼ (400 km s⁻¹)/(1 AU) implies that the acceleration rate is at least 100 times greater than adiabatic cooling. This justifies our neglecting the divergence term to compute the accelerated particle spectrum.

4.4. The Alfvén Mach Number Throughout the Heliosphere

In view of the central role played by the Alfvén Mach number in the solutions described in Sections 4.1–4.3, in this subsection we present observations describing the variation of the Alfvén Mach number as a function of heliocentric distance and solar activity. In Figure 6, we plot the Alfvén Mach number as a function of heliocentric distance, together with the corresponding value of the spectral index (3 + M_A)/2 derived from the induced electric field acceleration example. Specifically, we used hourly data from Helios 2, Pioneer Venus, Orbiter, IMP8, and Voyager 2 for 1976–1979 as described by Khabarova & Obridko (2012), together with Voyager 1 and Helios 1 data sets for the corresponding period. The values were derived as described in http://omniweb.sci.gsfc.nasa.gov/ftpbrowser/bow_derivation.html, taking into account the multispecies content of the solar wind and averaged by distance and all spacecraft. From Figure 6, the Alfvén Mach number is approximately constant beyond ∼1.5 AU with a value of 10–11, yielding a spectral index of ∼6.5–7 for this region. By contrast, the Alfvén Mach number increases from a value of about 5 at 0.3 AU to ∼9 at 1 AU. The corresponding spectral index for the antireconnection case changes from 4 to just less than 6. Of course, the variance, as illustrated by the error bars of Figure 6, can be large. Figure 7 shows that the Alfvén Mach number at 1 AU in the ecliptic plane can vary significantly with solar cycle. During solar maximum, values of M_A ∼ 7–8 are common, after which M_A increases as the solar minimum approaches (see also Svalgaard & Cliver 2010). Current sheet crossings are associated with plasma characterized by the highest possible Alfvén Mach number because of the abruptly decreasing strength of the interplanetary magnetic field and the enhanced solar wind density. Therefore, even at 1 AU, the corresponding spectral index should reach values typical for the outer heliosphere (and much more) in the vicinity of the heliospheric current sheet.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is interesting that if we consider the antireconnection-induced electric field-alone solution (Section 4.1), the roughly constant value of M_A in the outer heliosphere between ∼1.5 and 6 AU suggests that f ∼ c^−6.5 or c⁻⁷ should hold throughout this region during quiet times. It is also intriguing that from about 0.5–1 AU, the value of the Alfvén Mach number yields a spectral index in the range of ∼4–5, or f ∼ c⁻⁴–c⁻⁵ (Figure 6), which is similar to quiet-time observations of suprathermal (1–50 keV) particle spectra with c⁻⁵ reported by Fisk & Gloeckler (2006, 2012) and Mewaldt et al. (2001) for the supersonic solar wind. Similarly, during solar maximum periods at ∼1 AU, M_A ∼ 7–8, for which f ∼ c⁻⁵–c⁻⁶. However, these spectra may be masked in part by particle acceleration due to diffusive shock acceleration during solar maximum. Nonetheless, particle acceleration by reconnection of the kind described here may be useful in providing a seed particle population during this time.

Multiple reconnection has been identified as a primary dissipation process in an evolving turbulent magnetized plasma. Dissipation proceeds through the dynamical interaction of quasi-2D magnetic islands that form a "sea" of coherent structures within which particle energization proceeds. This perspective is supported by both numerical simulations (kinetic and MHD) and observations in the supersonic solar wind. Merging magnetic islands produce a reconnection-induced electric field within a larger closed structure, magnetic islands experience contraction, and magnetic field lines decrease in length during the merging process. Particles trapped in the vicinity of merging islands (Figure 1) repeatedly experience the induced electric field, leading some to gain energy. Particles trapped within contracting magnetic islands gain energy via a first-order Fermi mechanism. Oka et al. (2010) find that the electric field created by magnetic merging is the most effective of the particle mechanisms associated with the merging of magnetic islands. This in part may be due to the possibility that magnetic islands can oscillate rather than uniformly contract. An interesting result that is supportive (but not uniquely definitive) of the idea that magnetic islands might energize particles was obtained by Tessein et al. (2013). Tessein et al. used ACE magnetic field data to identify coherent structures in the solar wind using an approach related to that developed by, e.g., Greco et al. (2009a) and Osman et al. (2011). Using ACE energetic particle data and employing a PVI (partial variance of increments) analysis, Tessein et al. found a strong association of energetic particles in the range 0.047–4.75 MeV with solar wind discontinuities. One possible explanation (although certainly not the only one) is that the particles seen at PVI events were energized locally, possibly by reconnection events.

Khabarova et al. (2014) have presented observations that show the presence of magnetic islands in the vicinity of the heliospheric current sheet. They find that (1) magnetic islands in the solar wind possess a range of spatial scales and (2) magnetic islands experience dynamical merging in the solar wind. This work extends a variety of solar wind observations that support the physical picture described here, viz., one of a turbulent solar wind plasma in which quasi-2D magnetic islands with a range of spatial scales experience multiple reconnection events.

We have derived a transport equation for a gyrotropic distribution of particles experiencing pitch-angle scattering and energization via electric fields generated in a dissipative multireconnection super-Alfvénic plasma. Simulations of particle acceleration in a plasma with multiple reconnection events frequently exhibit isotropic energetic particle distributions. We therefore simplify the gyrophase-averaged or focused transport equation by assuming an isotropic particle distribution. This yields an advection–diffusion transport equation that resembles the well known cosmic ray transport equation except for energization terms due to stochastically distributed reconnection electric fields, contracting magnetic islands, and magnetic field line shortening associated with magnetic island merging.

The total acceleration timescale is given by a balance of the adiabatic term associated with the divergence of the large-scale flow and the magnetic island merging acceleration timescale, i.e.,

$$\begin{eqnarray*} &&\tau _{{\rm acc}}^{{\rm total}} \simeq -\frac{1}{3} \nabla \cdot {\bf U} + \frac{8 V_A}{3 \ell }, \end{eqnarray*}$$

where ℓ is a characteristic magnetic island size. At 1 AU, the divergence term can be neglected locally because it is ∼100 times smaller. The acceleration timescale V_A/ℓ reflects the time available for trapped particles in merging magnetic islands to either sample repeatedly the reconnection electric field induced by the merging or be trapped in a contracting magnetic island. This is a vital distinction between the time available for charged particles to sample the reconnection electric field induced by the typical reconnection model that is generated by antiparallel magnetic fields in the absence of islands. In the latter case, the acceleration timescale will be determined by Alfvén speed and the reconnection diffusion length scale (distinct from the diffusion scale discussed here in the context of particle transport). Consequently, particle acceleration due to particle trapping in contracting magnetic islands and acceleration in the reconnection electric field induced by island merging is a far more efficient acceleration mechanism than particle acceleration by reconnection electric fields in the absence of magnetic islands.

By considering a simple incompressible background flow, we solved the steady-state isotropic form of the 1D transport equation with a fixed source and injection speed c₀. We assumed that the spatial diffusion coefficient was constant and neglected the second-order diffusion terms in energy. We showed that the adiabatic expansion term was some 100 times slower than the particle energization terms associated with reconnecting magnetic islands and could therefore be neglected locally. We considered in some detail the individual contributions to particle energization by the induced electric field and contracting magnetic islands separately. We then solved the first-order correct transport equation with both the electric field and island contraction terms present to obtain the general solution.

The solution for the transport equation in which only the induced or antireconnection electric field energization term is retained yields pure power-law distributions for the accelerated particles, i.e., c > c₀, where c₀ is the injected initial particle speed. Because the particles experience stochastic encounters with the induced electric field, half the particles are accelerated, and half are cooled. Mathematically, this is expressed through the hyperbolic character of the form of the transport equation in phase space. The half of the injected particles that are accelerated have a power-law distribution with index α = −(3 + U/|V_E|)/2, where |V_E| is a scattered electric-field-induced velocity (15). A plausible estimate shows that |V_E| = V_A, yielding a power-law exponent α = −(3 + M_A)/2, where M_A is the Alfvén Mach number for particle speeds c > c₀. The half of the injected particles that lose energy have a power-law spectrum with index α = (3 + M_A)/2 for particle speeds c < c₀. The full distribution, when plotted logarithmically, therefore has a triangular shape with the apex at c = c₀.

The transport equation is strictly parabolic when the energization terms are restricted to magnetic island contraction alone. The accelerated form of the accelerated distribution function is of a power-law-like nature, being a power law for particle speeds c > c₀. The power-law exponent α = −3(1 + τ_c/(8τ_diff)) is a function of the magnetic island contraction timescale τ_c and the particle diffusion timescale τ_diff = κ/U², where U is the large-scale background flow speed. Both of these timescales reflect the microphysics of the island contraction and particle scattering, and neither is particularly well known. If τ_c < τ_diff, then the accelerated particle distribution is very hard, tending to a power-law index α = −3. An inequality was derived, Equation (34), identifying the conditions under which a relatively hard power-law solution can be expected. Under typical quiet-time solar wind conditions, it is only necessary that the particle scattering mean-free path λ_mfp be larger than the characteristic scale size ℓ of the contracting magnetic island.

Finally, the full first-order correct transport equation was solved when both the induced electric field and the magnetic island contraction energization terms were included. As before, power-law-like solutions were obtained, with the index depending on both the Alfvén Mach number and the ratio τ_c/τ_diff. The value of the approximate power-law exponent is easily estimated for this general case because the exponent is determined primarily by the large argument approximation of the modified Bessel's function component of the solution. We derive a solution for particle speeds c greater than a few c₀ that shows the accelerated particle distribution is a power-law-like solution depending on only the Alfvén Mach number and the ratio of the magnetic island contraction scale to the particle diffusion timescale. The combination of induced electric field and magnetic island contraction yields hard spectra for typical parameters. The numerical solutions confirm that the spectra are hard, the flatness being determined by the size of the ratio τ_c/τ_diff for a given Alfvén Mach number. The larger M_A implies more efficient escape from acceleration regions and therefore leads to steeper power-law-like solutions for the general case. The governing differential equation was hyperbolic, and therefore many of the characteristics of the electric-field-alone solution were similarly present for this more general case.

Given the importance of the Alfvén Mach number in determining the spectral index or the conditions under which energization via the contracting island mechanism would generate a power law, we presented observations showing the averaged radial dependence of M_A for the period 1976–1979 from 0.3–∼6 AU. From 0.3–∼1 AU, M_A ranges from ∼5–9, and >1.5 AU, M_A ∼ 10–11. In the case of the induced electric-field-alone solution of the transport equation, this yields a predicted power-law spectral index that ranges from −4 to −6 in the inner heliosphere to −6.5 to −7 in the outer heliosphere, varying with solar cycle. At least in the inner heliosphere, this is intriguingly close to the c⁻⁵ quiet-time observations of suprathermal (1–50 keV) particle spectra in the supersonic solar wind observed by Fisk & Gloeckler (2006, 2012) and Mewaldt et al. (2001) at 1 AU. These values of M_A illustrate that the inequality (34) will be relatively easily met for both the inner and outer heliosphere, so that particle acceleration by island contraction in a reconnecting flow should produce relatively hard spectra. As we showed, the combination of the induced antireconnection electric field and island contraction conspire to generate hard power-law-like solutions of the nearly isotropic form of the transport equation.

Because dissipation via multiple reconnection in a turbulent plasma should also occur in subsonic flows, such as in the lower corona, coronal loops, and the inner heliosheath, it is possible that the mechanism presented here might apply to particle acceleration in these regions. However, the transport equation will need to be rederived to take into account the sub-Alfvénic flow. Finally, we speculate that particle energization through a reconnection mechanism might be applicable to collisionless shock waves. Shocks typically generate significant levels of vortical turbulence downstream, which may therefore be a site for the acceleration of electrons or protons. This may be important for the injection of particles into the diffusive shock acceleration mechanism, or perhaps even partly responsible for the acceleration of particles at collisionless shock waves.

We acknowledge the partial support of NASA grants NNX08AJ33G, Subaward 37102-2, NNX09AG70G, NNX09AG63G, NNX09AJ79G, NNG05EC85C, Subcontract A991132BT, NNX09AP74A, NNX10AE46G, NNX09AW45G, and NNX14AF43G, and NSF grant ATM-0904007. This work was also supported in part by a UAH IIDR award to G.P.Z. O.K. was supported by RFBR grant 14-02-00769. G.P.Z. thanks J. F. Drake, M. Swisdak, W. H. Matthaeus, R. Bruno, and L. Trenchi for useful discussions. G.P.Z. thanks R. Bruno for his kind hospitality while visiting the INAF.

A Legendre polynomial expansion of the gyrophase-averaged Equation (6) can be used to derive a simpler transport equation analogous to that used to describe the transport of cosmic rays. We expand the gyrophase-averaged particle distribution function f as

$$\begin{eqnarray*} f({\bf x}, t, c, \mu) &=& \sum _{n=0}^{\infty } \frac{1}{2} (2n + 1) P_n (\mu) f_n ({\bf x}, t, c), \\ {\rm where} \quad f_n ({\bf x}, t, c) &=& \int _{-1}^1 f P_n (\mu)\, d\mu. \end{eqnarray*}$$

The orthogonality condition is given by

$$\begin{eqnarray*} \int _{-1}^1 P_m (\mu) P_n (\mu) \,d\mu = \left\lbrace \begin{array}{ll}0 & m \ne n \\ \displaystyle \frac{2}{2n + 1} & m = n \end{array} \right.\!\!\!, \end{eqnarray*}$$

and some useful recurrence relations are

$$\begin{eqnarray*} &&(n + 1) P_{n+1} (\mu) = (2n + 1) \mu P_n (\mu) - n P_{n-1} (\mu); \\ &&(1 - \mu ^2) \frac{d}{d \mu } P_n (\mu) = n P_{n-1} (\mu) - n \mu P_n (\mu)\\ &&\quad = \frac{n(n + 1)}{2n + 1} \left[ P_{n-1} (\mu) - P_{n+1} (\mu) \right]; \\ &&\frac{d}{d \mu } P_{n+1} (\mu) - \mu \frac{d}{d \mu } P_n (\mu) = (n + 1) P_n (\mu); \\ &&\mu \frac{d}{d \mu } P_n (\mu) - \frac{d}{d \mu } P_{n-1} (\mu) = n P_n (\mu); \\ &&\frac{d}{d \mu } \left[ P_{n+1} (\mu) - P_{n-1} (\mu) \right] = (2n + 1) P_n (\mu). \end{eqnarray*}$$

We need

$$\begin{eqnarray*} \mu ^2 P_n &=& \frac{(n+1)(n+2)}{(2n+1)(2n+3)} P_{n+2} + \left[ \frac{(n+1)^2}{2n+3} + \frac{n^2}{2n-1} \right] \frac{P_{n}}{2n+1} \\ &{}& {} + \frac{n(n-1)}{(2n-1)(2n+1)} P_{n-2}, \end{eqnarray*}$$

to evaluate

$$\begin{eqnarray*} &&\int _{-1}^1 (3\mu ^2 - 1) P_m P_n d \mu - 3 \int _{-1}^1 \mu ^2 P_m P_n d \mu - \int _{-1}^1 P_m P_n d\mu, \end{eqnarray*}$$

and

$$\begin{eqnarray*} \mu P_{n-1} - \mu P_{n+1} &=& \left(\frac{n}{2n - 1} - \frac{n-1}{2n - 1} \right) P_n \\ &&+ \frac{n - 1}{2n - 1} P_{n-2} - \frac{n + 2}{2n + 3} P_{n+2}. \end{eqnarray*}$$

We systematically project and expand each of the terms in Equation (6) from left to right using the Legendre polynomial P_m(μ) and the expansion for f. The left-hand side of Equation (6) and the diffusion in μ term and the additional reconnection terms yield the complete transformed focused transport Equation (6), now, however, expressed as an infinite set of partial differential equations in the coefficients f_n of the Legendre polynomials,

$$\begin{eqnarray} &&\frac{\partial f_m}{\partial t} + U_i \frac{\partial f_m}{\partial x_i} + \frac{c b_i}{2m+1} \left[ (m+1) \frac{\partial f_{m+1}}{\partial x_i} + m \frac{\partial f_{m-1}}{\partial x_i} \right] \nonumber \\ &&\quad {}+ c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{1}{2} \frac{\partial f_m}{\partial c} - c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{3}{2} \frac{(m-1)m}{(2m+1)(2m-1)} \frac{\partial f_{m-2}}{\partial c} \nonumber \\ &&\quad {} - c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{3}{2} \left(\!\frac{(m+1)^2}{2m+3} + \frac{m^2}{2m-1} \!\right) \frac{1}{2m+1} \frac{\partial f_{m}}{\partial c} \nonumber \\ &&\quad - c b_i b_j \frac{\partial U_j}{\partial x_i} \frac{3}{2} \frac{ (m+2)(m+1) }{(2m+3)(2m+1)}\frac{\partial f_{m+2}}{\partial c} - \frac{c}{2} \frac{\partial U_i}{\partial x_i} \frac{\partial f_m}{\partial c} \nonumber \\ &&\quad + \frac{c}{2} \frac{\partial U_i}{\partial x_i} \frac{(m-1)m}{(2m+1)(2m-1)} \frac{\partial f_{m-2}}{\partial c} + \frac{c}{2} \frac{\partial U_i}{\partial x_i}\left(\!\frac{(m+1)^2}{2m+3}\right.\nonumber\\ &&\quad\left. + \frac{m^2}{2m\!-\!1} \!\right) \frac{1}{2m+1} \frac{\partial f_{m}}{\partial c} + \frac{c}{2} \frac{\partial U_i}{\partial x_i} \frac{ (m+2)(m+1) }{(2m+3)(2m+1)} \frac{\partial f_{m+2}}{\partial c} \nonumber \\ &&\quad - \frac{DU_i}{{dt}} \frac{b_i}{2m+1} \left[ (m+1) \frac{\partial f_{m+1}}{\partial c} + m \frac{\partial f_{m-1}}{\partial c} \right] + \frac{c}{2} \frac{\partial b_i}{\partial x_i} \nonumber\\ &&\quad\times\left[ \frac{(m+1)(m+2)}{2m+1} f_{m+1} - \frac{m(m-1)}{2m+1} f_{m-1} \right] + \frac{1}{2} \frac{\partial U_i}{\partial x_i} \nonumber \\ &&\quad \times\frac{(m\!+\!1)(m\!+\!2)(m\!+\!3)}{(2m\!+\!1)(2m\!+\!3)} f_{m\!+\!2} \!+\! \frac{1}{2} \frac{\partial U_i}{\partial x_i} \frac{m(m+1)}{(2m\!-\!1)(2m+3)} f_m \nonumber \\ &&\quad - \frac{1}{2} \frac{\partial U_i}{\partial x_i} \frac{m(m-1)(m-2)}{(2m+1)(2m-1)} f_{m-2} {}- \frac{3}{2} b_i b_j \frac{\partial U_j}{\partial x_i} \nonumber \\ &&\quad \times \frac{(m+1)(m+2)(m+3)}{(2m+1)(2m+3)} f_{m+2} \!-\! \frac{3}{2} b_i b_j \frac{\partial U_j}{\partial x_i}\frac{m(m+1)}{(2m\!-\!1)(2m+3)} f_m\nonumber\\ &&\quad + \frac{3}{2} b_i b_j \frac{\partial U_j}{\partial x_i} \frac{m(m-1)(m-2)}{(2m+1)(2m-1)} f_{m-2} - \frac{b_i}{c} \frac{DU_i}{{dt}}\nonumber\\ &&\quad\times \left[\! \frac{(m+1)(m+2)}{2m+1} f_{m+1} \!-\! \frac{m(m-1)}{2m+1} f_{m-1}\! \right] \!=\! - \nu _s \frac{ m (m+1) }{2} f_m \nonumber \\ &&\quad - 2 \eta _c f_m{} -\frac{q}{m} \frac{\delta E_3}{2m + 1} \left[ (m + 1) \frac{\partial f_{m+1} }{\partial c} + m \frac{\partial f_{m-1} }{\partial c} \right] \nonumber \\ &&\quad {} - \frac{q}{m} \frac{\delta E_3}{c} \left[ \frac{(m + 1)(m + 2)}{2m + 1} f_{m+1} - \frac{m(m - 1)}{2m + 1} f_{m - 1} \right] \nonumber \\ &&\quad {} -\frac{\eta _c }{2} \left(c\frac{\partial f_m}{\partial c} + \frac{(m - 1) m}{(2m + 1)(2m - 1)} c\frac{\partial f_{m+2} }{\partial c}\right.\nonumber\\ &&\quad + \left(\!\frac{(m + 1)^2 }{2m + 3} + \frac{m^2}{2m - 1} \!\right) c\frac{\partial f_m}{\partial c} \left. {} + \frac{(m + 2)(m + 1) }{(2m + 3)(2m + 1) } c\frac{\partial f_{m+2}}{\partial c} \!\right)\nonumber\\ &&\quad + \frac{3 \eta _m}{2} \left(\frac{1}{3} c \frac{\partial f_m }{\partial c} - \frac{(m - 1)m }{(2m + 3)(2m + 1) } c \frac{\partial f_{m-2} }{\partial c} \right. \nonumber \\ &&\quad \ {} - \left(\!\frac{(m + 1)^2 }{2m + 3} + \frac{m^2}{2m - 1} \!\right) \frac{1}{2m + 1} c \frac{\partial f_m }{\partial c} \!-\! \frac{(m + 2)(m + 1) }{(2m + 3)(2m + 1) } c\nonumber\\ &&\quad\left.\times \frac{\partial f_{m+2} }{\partial c} \right) {} -\left(\frac{1}{2} \eta _c + \frac{3}{2} \eta _m \right) \left(\frac{(m + 1)(m + 2)(m + 3) }{(2m + 1)(2m + 3) } f_{m+2}\right.\nonumber\\ &&\quad + \frac{m (m + 1)}{(2m - 1)(2m + 3)} f_m\left. {} - \frac{m(m - 1)(m - 2) }{(2m + 1)(2m - 1) } f_{m-2} \right). \end{eqnarray} \tag{ A1 }$$

The infinite set of partial differential Equations (A1) is equivalent to the focused transport Equation (10) and therefore as challenging to solve. At each order of the expansion, i.e., the partial differential equation for a Legendre coefficient of particular order, it is clearly seen that the equation possesses coefficients of a higher order. This is another expression of the closure problem. Closure is typically affected by simply truncating the Legendre polynomial expansion at a finite number of coefficients. In the body of the text, we consider the m = 1 and m = 2 approximations.