ELECTRON ACCELERATION IN CONTRACTING MAGNETIC ISLANDS DURING SOLAR FLARES

As described by Drake et al. (2006a, 2010), the mechanism underlying the electron acceleration can be understood in terms of conservation of the adiabatic invariants of particle motion (Northrop 1963) within an evolving plasmoid.

The first invariant is associated with the particle's gyration about magnetic-field lines. It is referred to as the magnetic moment of the particle, μ:

$$\begin{eqnarray}&&\mu =\displaystyle \frac{q\omega {\rho }^{2}}{2c}=\displaystyle \frac{{W}_{\mathrm{gyr}}}{B}.\end{eqnarray} \tag{ 1 }$$

Here ω is the gyrofrequency, ρ is the particle's gyroradius, q is its electric charge, W_gyr is its kinetic energy of gyration, B is the magnitude of the magnetic field, and c is the speed of light. As the particle moves within the plasmoid, the value of the local magnetic field may increase or decrease. Conservation of the magnetic moment implies that the gyration energy W_gyr simultaneously increases or decreases, respectively. In a static magnetic field, the total particle energy W is preserved, so the energy of parallel motion ${W}_{\parallel }$ decreases or increases to compensate for the change in W_gyr. If ${W}_{\parallel }$ decreases, this leads to the phenomenon of magnetic mirroring. In a time-dependent magnetic field, on the other hand, an increase in the field strength leads to an increase in the gyration energy; however, the parallel energy does not decrease. The increase in total energy in this case is due to the associated induced electric field. This phenomenon is referred to as betatron acceleration.

The second invariant is associated with the particle's motion along field lines in geometries that support closed orbits, either mirroring or transiting. It is referred to as the parallel action, J:

$$\begin{eqnarray}&&J=\oint {{mv}}_{\parallel }{ds},\end{eqnarray} \tag{ 2 }$$

where the integral of the parallel momentum ${{mv}}_{\parallel }$ is taken along the field line coordinate s over one mirroring period or one loop of a transiting spiral, respectively. In the case of the magnetic configuration considered in the present paper, the plasmoid contracts after formation, thus reducing the length of the particle's path. Under this condition, conservation of the action invariant implies an increase in the velocity magnitude and kinetic energy of parallel motion, ${v}_{\parallel }$ and ${W}_{\parallel }$, respectively. As for the first invariant, this momentum and energy transfer to the particle occurs through the intermediary of the electric field and generalizes Fermi's mechanism (Fermi 1954) for accelerating particles via reflection off moving magnetic mirrors.

These simple considerations of conservation of adiabatic invariants within a plasmoid show that the magnetic-field evolution can strongly affect the motion of individual particles. If the local plasma conditions within a plasmoid experience a significant change during its evolution, the energies of the particles trapped within the plasmoid will also change significantly. The primary particle species affected by this mechanism is the electron species (Guidoni et al. 2016). For that reason they are the focus of this work.

In order to conserve the second invariant, the electrons need to close their orbits (mirroring or transiting) much more rapidly than the ambient magnetic field changes. This condition is usually met for energetic electrons, since their thermal speed is much larger than the bulk speed of the electron-ion plasma.

The effect of the acceleration per individual electron within a single contracting island is expected to be modest (Drake et al. 2006a). For example, Guidoni et al. (2016) found an increase in energy below a factor of five. This result is extended, and generally confirmed, by the analysis presented in this paper.

During a solar flare, magnetic free energy is released via reconnection and energy is transferred from magnetic field to plasma. As mentioned in Section 2.1, the evolution of a single plasmoid does not produce a big impact on the electron population, and therefore cannot ultimately explain the observed energetic radiation.

The key idea is that reconnection occurs at many locations during a solar flare (Sheeley et al. 2004; Drake et al. 2006b; Fermo et al. 2010; Karpen et al. 2012), yielding a large number of plasmoids. Apart from simple contraction, these plasmoids can coalesce together as considered by Drake et al. (2013) and Zank et al. (2014). Interaction and merger of many plasmoids might eventually lead to a large energy gain by the electron population. In this work, however, we only consider the effect of acceleration in a single magnetic island.

The background plasma conditions have been simulated with the Adaptively Refined Magnetohydrodynamics (MHD) Solver (ARMS) code (DeVore 1991). A detailed description of this simulation can be found in Guidoni et al. (2016). For the sake of simplicity, a 2.5D setting is adopted, i.e., the system is assumed to have axial symmetry. Ideal MHD is assumed. The flare is triggered according to the principles discussed in Karpen et al. (2012) and earlier papers on the breakout mechanism for solar eruptions (e.g., Antiochos et al. 1999). The modeled current sheet, as expected, is populated by several X-points separated by plasmoids, which are usually referred to as magnetic islands under 2.5D conditions. Because the current sheet forms in the initially low-lying, sheared magnetic field that powers the eruptive flare, the magnetic field in the sheet has a strong component in the translationally invariant direction (out of the plane in Figure 1).

Zoom In Zoom Out Reset image size

A number of islands were well resolved on the adaptive mesh. Their shape, evolution, and life cycle were carefully tracked. The simulation revealed the presence of one persistent X-point with smaller ones generated near it, some moving Sunward until reaching and merging with the flare arcade, others moving outward and merging with the erupting coronal mass ejection.

Those islands that were present in the system for a significantly long time and had clear features were chosen for investigation to test the acceleration mechanism of Drake et al. (2006a). Since ARMS does not include a particle-tracking capability, only theoretical estimates could be made. We found that individual electrons should gain up to a factor of five in energy (Guidoni et al. 2016). The simulated MHD data for Island 2 of Guidoni et al. (2016) was exported and used as input for the particle simulation described below.

The electron population's evolution has been simulated with the Adaptive Mesh Particle Simulator (AMPS) code developed at the University of Michigan and successfully used by Tenishev et al. (2008, 2010, 2013).

Rather than directly integrating Newton's equations of motions for the electrons in the electromagnetic field, ${\boldsymbol{E}}$ and ${\boldsymbol{B}}$, which would be prohibitively expensive (see below), we elected to integrate the equations for guiding-center motion (Northrop 1963). In this model, the dominant motions are (1) the displacement of the guiding center parallel to the magnetic field with velocity ${v}_{\parallel }$, (2) the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ drift perpendicular to the magnetic field, with velocity ${{\boldsymbol{U}}}_{\perp }$, and (3) the fast gyration also perpendicular to the magnetic field, with velocity ${{\boldsymbol{v}}}_{\mathrm{gyr}}$. The electron kinetic energy W is then given by

$$\begin{eqnarray}W & = & {W}_{\parallel }+{W}_{\perp }+{W}_{\mathrm{gyr}}\\ & = & \displaystyle \frac{{{mv}}_{\parallel }^{2}}{2}+\displaystyle \frac{m{{\boldsymbol{U}}}_{\perp }^{2}}{2}+\displaystyle \frac{m{{\boldsymbol{v}}}_{\mathrm{gyr}}^{2}}{2}\\ & = & \displaystyle \frac{{{mv}}_{\parallel }^{2}}{2}+\displaystyle \frac{m{{\boldsymbol{U}}}_{\perp }^{2}}{2}+\mu B,\end{eqnarray} \tag{ 3 }$$

where m is the electron's mass and μ is its magnetic moment introduced in Equation (1). Additional terms, such as the gradient-B and curvature drifts, contribute small corrections to the perpendicular velocity and kinetic energy of the guiding center (see Northrop 1963, and the estimates below). This implies that the guiding centers perfectly follow the magnetic-field lines to first order, since the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ velocity is the cross-field advection rate of the field lines in the ideal-MHD simulation.

The total rate of change of ${v}_{\parallel }$, following the electron's orbit, is

$$\begin{eqnarray}\displaystyle \frac{{{dv}}_{\parallel }}{{dt}} & = & -\displaystyle \frac{\mu }{m}\displaystyle \frac{\partial B}{\partial s}+{{\boldsymbol{U}}}_{\perp }\cdot \displaystyle \frac{d{\boldsymbol{b}}}{{dt}}\\ & = & -\displaystyle \frac{\mu }{m}\displaystyle \frac{\partial B}{\partial s}+{{\boldsymbol{U}}}_{\perp }\cdot \left[\displaystyle \frac{\partial {\boldsymbol{b}}}{\partial t}+{v}_{\parallel }\displaystyle \frac{\partial {\boldsymbol{b}}}{\partial s}+({{\boldsymbol{U}}}_{\perp }\cdot {{\rm{\nabla }}}_{\perp }){\boldsymbol{b}}\right],\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{b}}$ is the local direction of the magnetic field, s is the distance along the magnetic field. The first term in Equation (4) describes the effect of magnetic mirroring on the parallel motion. As was mentioned previously, and will become evident below, this term represents a redistribution between the parallel and perpendicular energies of the electron along its orbit. The second term contains both time- and space-dependent variations of the magnetic-field direction along the orbit. As demonstrated in detail by Northrop (1963), among other effects these terms give rise to the increase of parallel energy due to Fermi's mechanism.

The total rate of change of the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ velocity is

$$\begin{eqnarray}&&\displaystyle \frac{d{{\boldsymbol{U}}}_{\perp }}{{dt}}=\displaystyle \frac{\partial {{\boldsymbol{U}}}_{\perp }}{\partial t}+{v}_{\parallel }\displaystyle \frac{\partial {{\boldsymbol{U}}}_{\perp }}{\partial s}+({{\boldsymbol{U}}}_{\perp }\cdot {{\rm{\nabla }}}_{\perp }){{\boldsymbol{U}}}_{\perp }.\end{eqnarray} \tag{ 5 }$$

As will be shown below, direct contribution of ${{\boldsymbol{U}}}_{\perp }$ to the electron's kinetic energy is small, being quadratic in the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ velocity. Indirectly, however, the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ drift is extremely important, as it is responsible for advecting and compressing the magnetic-field lines (whose electric fields accelerate the electrons) within the evolving plasmoids.

From the results above, we now can compute the rate of change of the electron's energy,

$$\begin{eqnarray}\displaystyle \frac{{dW}}{{dt}} & = & \displaystyle \frac{{{dW}}_{\parallel }}{{dt}}+\displaystyle \frac{{{dW}}_{\perp }}{{dt}}+\displaystyle \frac{{{dW}}_{\mathrm{gyr}}}{{dt}}\\ & \approx & \displaystyle \frac{{{dW}}_{\parallel }}{{dt}}+\displaystyle \frac{{{dW}}_{\mathrm{gyr}}}{{dt}},\end{eqnarray} \tag{ 6 }$$

upon neglecting the contribution due to the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ drift. Using Equation (4), we find for the parallel energy ${W}_{\parallel }$ that

$$\begin{eqnarray}\displaystyle \frac{{{dW}}_{\parallel }}{{dt}} & = & {{mv}}_{\parallel }\displaystyle \frac{{{dv}}_{\parallel }}{{dt}}\\ & = & -\mu {v}_{\parallel }\displaystyle \frac{\partial B}{\partial s}+{{mv}}_{\parallel }{{\boldsymbol{U}}}_{\perp }\cdot \displaystyle \frac{d{\boldsymbol{b}}}{{dt}}.\end{eqnarray} \tag{ 7 }$$

For the gyration energy W_gyr we have, recalling the invariance of μ,

$$\begin{eqnarray}\displaystyle \frac{{{dW}}_{\mathrm{gyr}}}{{dt}} & = & \mu \displaystyle \frac{{dB}}{{dt}}\\ & = & \mu \left[\displaystyle \frac{\partial B}{\partial t}+{v}_{\parallel }\displaystyle \frac{\partial B}{\partial s}+({{\boldsymbol{U}}}_{\perp }\cdot {{\rm{\nabla }}}_{\perp })B\right].\end{eqnarray} \tag{ 8 }$$

Combining Equations (7) and (8), and noting the cancellation of the common magnetic mirroring terms (the first and second, respectively), for the total energy change we finally obtain

$$\begin{eqnarray}\displaystyle \frac{{dW}}{{dt}} & \approx & \mu \left[\displaystyle \frac{\partial B}{\partial t}+{{\boldsymbol{U}}}_{\perp }\cdot {{\rm{\nabla }}}_{\perp }B\right]\\ & & +{{mv}}_{\parallel }{{\boldsymbol{U}}}_{\perp }\cdot \left[\displaystyle \frac{\partial {\boldsymbol{b}}}{\partial t}+{v}_{\parallel }\displaystyle \frac{\partial {\boldsymbol{b}}}{\partial s}+({{\boldsymbol{U}}}_{\perp }\cdot {{\rm{\nabla }}}_{\perp }){\boldsymbol{b}}\right].\end{eqnarray} \tag{ 9 }$$

In this equation, the two terms on the first line come from the change in the gyration energy, and are commonly identified as betatron acceleration terms. Clearly, they arise from changes in the magnetic-field strength due either to a local compression or rarefaction of the flux density ($\partial B/\partial t$) or to an advection of the field line into a region of spatially varying flux density (${{\boldsymbol{U}}}_{\perp }\cdot {{\rm{\nabla }}}_{\perp }B$). The three terms on the second line come from the change in the parallel energy, and are commonly identified as Fermi acceleration terms. They arise from changes in the field direction due to a local reorientation of the field ($\partial {\boldsymbol{b}}/\partial t$), to parallel motion of the electron into a region of differing direction ($\partial {\boldsymbol{b}}/\partial s$), or to perpendicular advection of the field line into a region of differing direction (${{\boldsymbol{U}}}_{\perp }\cdot {{\rm{\nabla }}}_{\perp }{\boldsymbol{b}}$).

The great benefit of utilizing the guiding-center approximation in our simulation is the elimination of the enormous restriction on the time step. The scale of this restriction is demonstrated by the following estimation. The value of the magnetic field in the region of interest is $\sim {10}^{-5}$ T, which results in a gyration period of $m/({qB})\sim {10}^{-6}$ s, where q is the electron's charge. In order to simulate a process that lasts 10 minutes, one would need to perform at least half a trillion time steps per electron. Since the number of test electrons per field line needed to produce reliable statistics is of the order of a million, the problem would be virtually impossible to simulate on modern hardware. Thus, the adoption of guiding-center motion equations greatly facilitates simulation of the problem.

In order for their motion to be adiabatic, electrons need to move fast enough. The Alfvén speed in the domain is ${V}_{{\rm{A}}}\sim {10}^{5}$ m s⁻¹. For this reason, the speed range of electron injection was chosen to be between ${v}_{\min }^{\mathrm{inj}}=5\times {10}^{6}$ m s⁻¹ and ${v}_{\max }^{\mathrm{inj}}=1.2\times {10}^{7}$ m s⁻¹. The electrons in this energy range were randomly sampled from the Maxwellian distribution

$$\begin{eqnarray}&&f({\boldsymbol{v}})\propto \exp \left\{-\displaystyle \frac{m{({\boldsymbol{v}}-{\boldsymbol{U}})}^{2}}{2{k}_{{\rm{B}}}T}\right\},\end{eqnarray} \tag{ 10 }$$

where k_B is the Boltzmann constant and T and ${\boldsymbol{U}}$ are local plasma temperature and bulk velocity.

With the energy range set as above, one can estimate the contribution of the cross-field drift to the motion of electrons' guiding centers (Northrop 1963). Under ideal MHD, the ${\boldsymbol{E}}\times {\boldsymbol{B}}$ drift accounts for the advection of the magnetic-field lines and, therefore, for the cross-field motion of electrons attached to them. Its magnitude, ${U}_{\perp }\sim {10}^{5}\,{\rm{m}}\,{{\rm{s}}}^{-1}\ll {v}_{\min }^{\mathrm{inj}}$, has a negligible contribution to the total electron energy (see Equations (4) and (9)). However, the advection is responsible for both the betatron and Fermi acceleration of electrons as they sample regions with varying magnetic-field strength and direction. The gradient and curvature drifts both are on the order of $W/{eBL}\lesssim m{({v}_{\max }^{\mathrm{inj}})}^{2}/{eBL}\sim 10$ m s⁻¹, where $L\sim {10}^{7}\,{\rm{m}}$ (see Figure 2) is the characteristic spatial scale of the magnetic field's geometry. The drift associated with acceleration of the field lines' advection, ${\boldsymbol{a}}$, is $-m{\boldsymbol{a}}\times {\boldsymbol{B}}/{{eB}}^{2}$. With the characteristic time of change $\tau \sim {10}^{2}\,{\rm{s}}$ (see Figure 2), the estimate for this drift is $\sim {mU}/{eB}\tau \sim {10}^{-3}$ m s⁻¹. Therefore, the cross-field drift that may displace electrons from their original field lines is much less significant than ${\boldsymbol{E}}\times {\boldsymbol{B}}$ drift and can be neglected in the guiding-center equations. Overall, one may safely assume that each electron stays on the same field line throughout the simulation.

Zoom In Zoom Out Reset image size

It is also worth noting that the guiding-center approximation is naturally tied to the notion of the magnetic-field lines. Thus, the underlying physics of the process motivates improvements in the numerical approach. Since the motion of high energy electrons is essentially one-dimensional, as viewed on a larger scale, and follows the magnetic-field lines, it is most practical to adopt a data structure that takes advantage of it.

Rather than utilizing a Cartesian grid, a linear structure that follows a magnetic-field line has been developed. Under the assumption of axial symmetry used to simulate the background field and plasma (see Section 3.1), which reduces the dimensionality of the problem, we introduce a scalar potential for poloidal (out of the plane of Figure 2) components of the magnetic field as well as the magnetic flux function, Ψ. The latter quantity greatly simplifies the problem of field line extraction, while keeping high accuracy. The procedure is as follows.

The poloidal projection of a magnetic-field line is defined as a line of constant value of Ψ. If one chooses a starting point for a magnetic loop, it is possible to extract the whole loop point-by-point.

The initial point, X₀, is taken at an arbitrary location within the magnetic island. The value of the magnetic flux function, ${{\rm{\Psi }}}_{0}$, must be the same for all points associated with this magnetic loop. The candidate point is taken in the direction of the magnetic field, then it is corrected via gradient descent iterations until the value of the magnetic flux function is sufficiently close to ${{\rm{\Psi }}}_{0}$.

$$\begin{eqnarray}{\tilde{{\boldsymbol{X}}}}_{0} & = & {{\boldsymbol{X}}}_{n-1}+{\boldsymbol{B}}\ast \delta \\ {\tilde{{\boldsymbol{X}}}}_{k} & = & {\tilde{{\boldsymbol{X}}}}_{k-1}+\alpha \displaystyle \frac{({\tilde{{\rm{\Psi }}}}_{k}-{{\rm{\Psi }}}_{0})}{{B}^{2}}{\rm{\nabla }}{\tilde{{\rm{\Psi }}}}_{k},\quad k\geqslant 1\\ {{\boldsymbol{X}}}_{n} & = & {\tilde{{\boldsymbol{X}}}}_{\infty },\end{eqnarray} \tag{ 11 }$$

where ${{\boldsymbol{X}}}_{n-1}$ is the last extracted point on the field line, ${\boldsymbol{B}}$ is the magnetic field at this location, ${{\boldsymbol{X}}}_{n}$ is the point being extracted, k enumerates iterations, ${\tilde{{\boldsymbol{X}}}}_{k}$ and ${\tilde{{\rm{\Psi }}}}_{k}$ are values at iteration k, ${\tilde{{\boldsymbol{X}}}}_{\infty }$ is the limit of iterations, δ controls distance between consecutive extracted points, and α controls stability of gradient descent.

The effectiveness and accuracy of this procedure are demonstrated in Figure 1. As one can see, even complex features such as the stretched parts of field lines near an X-point additionally complicated with low values of the magnetic field are carefully resolved.

Each point thus extracted moves with the local plasma velocity. Under ideal MHD, frozen-in conditions apply, hence allowing us to trace the evolution of an individual field line via following motion of fluid elements extracted along it.

Although only the poloidal projections of field lines were extracted, the full 3D plasma velocity is used to update their state. The new locations of thus-obtained fluid elements are then projected to the original poloidal plane. The same procedure was applied to the electron motion.

The background plasma data from our ARMS simulation are imported into AMPS, with a data file cadence of 1 s. To maintain higher accuracy, linear time interpolation between files is used. Extracted field lines are regularly updated with each plasma fluid element being moved with the velocity of the background plasma and corrected via gradient descent to match the original value of Ψ0, which should be preserved under ideal MHD.

Another reason why our approach is better suited for particle tracking in the current case is the issue with electrons staying within the island. While the physical reason for electrons to escape the island is the cross-field drift, artificial electron leakage from the island is observed when a full 3D simulation on a regular grid is performed. Since the field evolves rapidly, an electron on the edge of the island commonly finds itself outside the island when the plasma state is updated. Such electrons escape the domain and information is lost. Our solution is to correct the position of an electron in order to put it back into the island (Equation (11)). With the adopted field-aligned data structure, this problem was effectively non-existent. Since all electrons are essentially attached to field lines, they cannot escape the island due to non-physical reasons.