NONTHERMALLY DOMINATED ELECTRON ACCELERATION DURING MAGNETIC RECONNECTION IN A LOW-β PLASMA

Magnetic reconnection is a fundamental plasma process during which the magnetic field restructures itself and converts its energy into plasma kinetic energies (e.g., Priest & Forbes 2000). It occurs ubiquitously in laboratory, space, and astrophysical magnetized plasmas. An important unsolved problem is the acceleration of nonthermal particles in the reconnection region. Magnetic reconnection has been suggested as a primary mechanism for accelerating nonthermal particles in solar flares (Masuda et al. 1994; Krucker et al. 2010; Lin 2011), Earth's magnetosphere (Øieroset et al. 2002; Fu et al. 2011; Huang et al. 2012), the sawtooth crash of tokamaks (Savrukhin 2001), and high-energy astrophysical systems (Colgate et al. 2001; Zhang & Yan 2011). In particular, observations of solar flares have revealed an efficient particle energization with 10%–50% of magnetic energy converted into energetic electrons and ions (Lin & Hudson 1976). The energetic particles usually develop a power-law energy distribution that contains energy on the same order of the dissipated magnetic energy (Krucker et al. 2010; Oka et al. 2015). Some observations find that the emission has no distinguishable thermal component, indicating that most of the electrons are accelerated to nonthermal energies (Krucker et al. 2010; Krucker & Battaglia 2014). This efficient production of energetic particles poses a challenge to current theories of particle acceleration.

Particle acceleration associated with reconnection has been studied in reconnection-driven turbulence (Miller et al. 1996), at shocks in the outflow region (Tsuneta & Naito 1998; Guo & Giacalone 2012), and in the reconnection layer (Drake et al. 2006; Fu et al. 2006; Oka et al. 2010; Kowal et al. 2012; Guo et al. 2014; Zank et al. 2014). Previous kinetic simulations have examined various acceleration mechanisms during reconnection, including the Fermi-type mechanism in magnetic islands (flux ropes in three-dimensional simulations; Drake et al. 2006; Guo et al. 2014) and direct acceleration in the diffusion region (Pritchett 2006; Huang et al. 2010). Most simulations focus on regimes with plasma β ≥ 0.1 with no obvious power-law distributions observed. It was argued that particle loss from the simulation domain is important for developing a power-law distribution (Drake et al. 2010). Early simulations of relativistic reconnection showed power-law distributions within the X-region (Zenitani & Hoshino 2001). Recent kinetic simulations with a highly magnetized ($\sigma =\frac{{B}^{2}}{4\pi {n}_{e}{m}_{e}{c}^{2}}\gg 1$) pair plasma found global power-law distributions without the particle loss, although the loss mechanism may be important in determining the spectral index (Guo et al. 2014, 2015). It is unknown whether or not this is valid for reconnection in a nonrelativistic proton–electron plasma since its property is different from the relativistic reconnection (Liu et al. 2015).

Motivated by the results of relativistic reconnection, here, we report fully kinetic simulations of magnetic reconnection in a nonrelativistic proton–electron plasma with a range of electron and ion β_e = β_i = 0.007–0.2. The low-β regime was previously relatively unexplored due to various numerical challenges. We find that reconnection in the low-β regime drives efficient energy conversion and accelerates electrons into a power-law distribution f(E) ∼ E⁻¹. At the end of the low-β cases, more than half of the electrons in number and 90% in energy are in the nonthermal electron population. This strong energy conversion and particle acceleration led to a post-reconnection region with the kinetic energy of energetic particles comparable to magnetic energy. Since most electrons are magnetized in the low-β plasma, we use a guiding-center drift description to demonstrate that the main acceleration process is a Fermi-type mechanism through the particle curvature drift motion along the electric field induced by fast plasma flows. The development of power-law distributions is consistent with the analytical model (Guo et al. 2014). The nonthermally dominated energization may help explain the efficient electron acceleration in the low-β plasma environments, such as solar flares and other astrophysical reconnection sites.

In Section 2, we describe the numerical simulations. In Section 3, we present simulation results and discuss the conditions for the development of power-law distributions. We discuss and conclude the results in Section 4.

The kinetic simulations are carried out using the VPIC code (Bowers et al. 2008), which solves Maxwell's equations and follows particles in a fully relativistic manner. The initial condition is a force-free current sheet with a magnetic field ${\boldsymbol{B}}={B}_{0}\mathrm{tanh}(z/\lambda )\hat{x}+{B}_{0}\mathrm{sech}(z/\lambda )\hat{y}$, where λ = d_i is the half thickness of the layer. Here, d_i is the ion inertial length. The plasma consists of protons and electrons with a mass ratio m_i/m_e = 25. The initial distributions for both electrons and protons are Maxwellian with uniform density n₀ and temperature ${{kT}}_{i}={{kT}}_{e}=0.01{m}_{e}{c}^{2}.$ A drift velocity for electrons U_e is added to represent the current density that satisfies Ampere's law. The initial electron and ion ${\beta }_{e}={\beta }_{i}=8\pi {n}_{0}{{kT}}_{e}/{B}_{0}^{2}$ are varied by changing ω_pe/Ω_ce, where ${\omega }_{{pe}}=\sqrt{4\pi {n}_{0}{e}^{2}/{m}_{e}}$ is the electron plasma frequency and ${{\rm{\Omega }}}_{{ce}}={{eB}}_{0}/({m}_{e}c)$ is the electron gyrofrequency. Quantities β_e = 0.007, 0.02, 0.06, and 0.2 correspond to ω_pe/Ω_ce = 0.6, 1, $\sqrt{3}$ and $\sqrt{10},$ respectively. The domain sizes are ${L}_{x}\times {L}_{z}=200{d}_{i}\times 100{d}_{i}.$ We use ${N}_{x}\times {N}_{z}=4096\times 2048$ cells with 200 particles per species per cell. The boundary conditions are periodic along the x-direction, perfectly conducting boundaries for fields and reflecting boundaries for particles along the z-direction. A long wavelength perturbation is added to induce reconnection (Birn et al. 2001).

Under the influence of the initial perturbation, the current sheet quickly thins down to a thickness of ∼d_e (electron inertial length c/ω_pe) that is unstable to the secondary tearing instability (Daughton et al. 2009; Liu et al. 2013b). Figures 1(a) and (b) show the evolution of the out-of-plane current density. The reconnection layer breaks and generates a chain of magnetic islands that interact and coalesce with each other. The largest island eventually grows comparable to the system size and the reconnection saturates at tΩ_ci ∼ 800. Figure 1(c) shows the time evolution of the magnetic energy in the x-direction (the reconnecting component) _bx and the kinetic energy of electrons K_e and ions K_i for the case with β_e = 0.02, respectively. Throughout the simulation, 40% of the initial _bx is converted into plasma kinetic energy. Of the converted energy, 38% goes into electrons and 62% goes into ions. We have carried out simulations with larger domains (not shown) to confirm that the energy conversion is still efficient and weakly depends on system size. Since the free magnetic energy overwhelms the initial kinetic energy, particles in the reconnection region are strongly energized. Eventually, K_e and K_i are 5.8 and 9.4 times their initial values, respectively. Figure 1(d) shows the ratio of the electron energy gain ΔK_e to the initial electron energy K_e(0) for different cases. While the β_e = 0.2 case shows only mild energization, cases with lower β_e give stronger energization as the free energy increases.

Zoom In Zoom Out Reset image size

The energy conversion drives strong nonthermal electron acceleration. Figure 2(a) shows the final electron energy spectra over the whole simulation domain for the four cases. More electrons are accelerated to high energies for lower-β cases, similar to earlier simulations (Bessho & Bhattacharjee 2010). More interestingly, in the cases with β_e = 0.02 and 0.007, the energy spectra develop a power-law-like tail f(E) ∼ E^−p with the spectral index p ∼ 1. This is similar to results from relativistic reconnection (Guo et al. 2014, 2015). We have carried out one simulation with m_i/m_e = 100 and β_e = 0.02 and find a similar electron spectrum. In contrast, the case with β_e = 0.2 does not show any obvious power-law tail, consistent with earlier simulations (Drake et al. 2010). The nonthermal population dominates the distribution in the low-β cases. For example, when we subtract the thermal population by fitting the low-energy distribution as Maxwellian, the nonthermal tail in the β_e = 0.02 case contains 55% of the electrons and 92% of the total electron energy. The power-law tail breaks at energy E_b ∼ 10E_th for β_e = 0.02 and extends to a higher energy for β_e = 0.007. Figure 2(b) shows the fraction of nonthermal electrons for different cases. For β_e = 0.007, the nonthermal fraction goes up to 66%, but it decreases to 17% for β_e = 0.2. Figures 2(c) and (d) show n_acc/n_e at tΩ_ci = 125 and 400 for the case with β_e = 0.02, where n_acc is the number density of accelerated electrons with energies larger than three times their initial thermal energy, and n_e is the total electron number density. The fraction of energetic electrons is over 40% and up to 80% inside the magnetic islands and reconnection exhausts, indicating a bulk energization for most of electrons in the reconnection layer. The energetic electrons will eventually be trapped inside the largest magnetic island. The nonthermally dominated distribution contains most of the converted magnetic energy, indicating that energy conversion and particle acceleration are intimately related.

Zoom In Zoom Out Reset image size

To study the energy conversion, Figure 3(a) shows the energy conversion rate $d{\varepsilon }_{c}/{dt}$ from the magnetic field to electrons through directions parallel and perpendicular to the local magnetic field. We define $d{\varepsilon }_{c}/{dt}={\displaystyle \int }_{{\mathcal{D}}}{\boldsymbol{j}}^{\prime} \cdot {\boldsymbol{E}}{dV},$ where ${\mathcal{D}}$ indicates the simulation domain and ${\boldsymbol{j}}^{\prime}$ is ${{\boldsymbol{j}}}_{\parallel }$ or ${{\boldsymbol{j}}}_{\perp }.$ We find that energy conversion from the perpendicular directions gives ∼90% of the electron energy gain. By tracking the trajectories (not shown) of a large number of accelerated electrons, we find various acceleration processes in the diffusion region, magnetic pile-up region, contracting islands, and island coalescence regions (Hoshino et al. 2001; Hoshino 2005; Drake et al. 2006; Fu et al. 2006; Huang et al. 2010; Oka et al. 2010; Dahlin et al. 2014; Guo et al. 2014). The dominant acceleration is by particles bouncing back and forth through a Fermi-like process accomplished by particle drift motions within magnetic structures (X. L. Li et al. 2015, in preparation). To reveal the role of particle drift motions, we use a guiding-center drift description to study the electron energization for the β_e = 0.02 case. The initial low β guarantees that this is a good approximation since the typical electron gyroradius ρ_e is smaller than the spatial scale of the field variation (∼d_i).

Zoom In Zoom Out Reset image size

By ensemble averaging the particle gyromotion and drift motions, the perpendicular current density for a single species can be expressed as (Parker 1957; Blandford et al. 2014)

$$\begin{eqnarray}{{\boldsymbol{j}}}_{\perp } & = & {P}_{\parallel }\displaystyle \frac{{\boldsymbol{B}}\times ({\boldsymbol{B}}\cdot {\rm{\nabla }}){\boldsymbol{B}}}{{B}^{4}}\\ & & +\;{P}_{\perp }\left(\displaystyle \frac{{\boldsymbol{B}}}{{B}^{3}}\right)\times {\rm{\nabla }}B-{\left[{\rm{\nabla }}\times \displaystyle \frac{{P}_{\perp }{\boldsymbol{B}}}{{B}^{2}}\right]}_{\perp }+\rho \displaystyle \frac{{\boldsymbol{E}}\times {\boldsymbol{B}}}{{B}^{2}}\\ & & +{\rho }_{m}\displaystyle \frac{{\boldsymbol{B}}}{{B}^{2}}\times \displaystyle \frac{d{{\boldsymbol{u}}}_{E}}{{dt}}\end{eqnarray} \tag{ 1 }$$

using a gyrotropic pressure tensor ${\mathcal{P}}={P}_{\perp }{\mathcal{I}}+({P}_{\parallel }-{P}_{\perp }){\boldsymbol{bb}},$ where ${P}_{\parallel }\equiv {m}_{e}\int {{fv}}_{\parallel }^{2}d{\boldsymbol{v}}$ and ${P}_{\perp }\equiv 0.5{m}_{e}\int {{fv}}_{\perp }^{2}d{\boldsymbol{v}},$ ρ is the particle charge density, and ρ_m is the particle mass density. The terms on the right are due to curvature drift, ∇B drift, magnetization, ${\boldsymbol{E}}\times {\boldsymbol{B}}$ drift, and polarization drift, respectively. The expression is simplified as ${{\boldsymbol{j}}}_{\perp }={{\boldsymbol{j}}}_{c}+{{\boldsymbol{j}}}_{g}+{{\boldsymbol{j}}}_{m}+{{\boldsymbol{j}}}_{{\boldsymbol{E}}\times {\boldsymbol{B}}}+{{\boldsymbol{j}}}_{{\boldsymbol{p}}},$ in which ${{\boldsymbol{j}}}_{{\boldsymbol{E}}\times {\boldsymbol{B}}}$ has no direct contribution to the energy conversion. This gives an accurate description for ${{\boldsymbol{j}}}_{\perp }$ if the pressure tensor is gyrotropic. To confirm this, we calculate the electron pressure agyrotropy $A\unicode{216}_{e}\equiv 2\frac{| {P}_{\perp e1}-{P}_{\perp e2}| }{{P}_{\perp e1}+{P}_{\perp e2}},$ where ${P}_{\perp e1}$ and ${P}_{\perp e2}$ are the two pressure eigenvalues associated with eigenvectors perpendicular to the mean magnetic field direction (Scudder & Daughton 2008). $A\unicode{216}_{e}$ measures the departure of the pressure tensor from cylindrical symmetry about the local magnetic field. It is zero when the local particle distribution is gyrotropic. Figure 3(b) shows that the regions with nonzero $A\unicode{216}_{e}$ are localized to X-points. The small $A\unicode{216}_{e}$ indicates that the electron distributions are nearly gyrotropic in most regions. Therefore, the drift description is a good approximation for electrons in our simulations even without an external guide field, which is required for this description in a high-β plasma (Dahlin et al. 2014).

Figures 4(a) and (b) show time-dependent d_c/dt and _c from different current terms, where ${\varepsilon }_{c}={\displaystyle \int }_{0}^{t}(d{\varepsilon }_{c}/{dt}){dt}.$ The contribution from polarization current and parallel current are small and not shown. The curvature drift term is a globally dominant term of ${{\boldsymbol{j}}}_{\perp }\cdot {\boldsymbol{E}},$ the ∇B term gives a net cooling, and the magnetization term is small compared to these two. Figure 4(c) shows the spatial distribution of ${{\boldsymbol{j}}}_{c}\cdot {\boldsymbol{E}}.$ When the flow velocity ${\boldsymbol{u}}$ is along the magnetic field curvature ${\boldsymbol{\kappa }}$ due to tension force, ${{\boldsymbol{j}}}_{c}\cdot {\boldsymbol{E}}\approx ({P}_{\parallel }{\boldsymbol{B}}\times {\boldsymbol{\kappa }}/{B}^{2})\cdot (-{\boldsymbol{u}}\times {\boldsymbol{B}})\gt 0.$ These regions are a few d_i along the z-direction but over 50d_i along the x-direction. The overall effect of ${{\boldsymbol{j}}}_{c}\cdot {\boldsymbol{E}}$ is a strong electron energization. Figure 4(d) shows that ${{\boldsymbol{j}}}_{g}\cdot {\boldsymbol{E}}$ is negative in most regions because the strong ∇B is along the direction out of the reconnection exhausts. Then, ${{\boldsymbol{j}}}_{g}\cdot {\boldsymbol{E}}\sim ({\boldsymbol{B}}\times {\boldsymbol{\nabla }}{\boldsymbol{B}})\cdot (-{\boldsymbol{u}}\times {\boldsymbol{B}})\lt 0.$ Note at some regions, ${{\boldsymbol{j}}}_{g}\cdot {\boldsymbol{E}}$ can give strong acceleration. Figure 4(e) shows the cumulation of ${{\boldsymbol{j}}}_{c}\cdot {\boldsymbol{E}}$ and ${{\boldsymbol{j}}}_{g}\cdot {\boldsymbol{E}}$ along the x-direction. In the reconnection exhaust region ($x=60-115{d}_{i}$), ${{\boldsymbol{j}}}_{c}\cdot {\boldsymbol{E}}$ is stronger than ${{\boldsymbol{j}}}_{g}\cdot {\boldsymbol{E}},$ so the electrons can be efficiently accelerated when going through these regions. In the pile-up region($x=120-140{d}_{i}$), ${\boldsymbol{\kappa }},$ ∇B and ${\boldsymbol{u}}$ are along the same direction, so both terms give electron energization. In the island coalescence region($x\sim 150{d}_{i}$), ${{\boldsymbol{j}}}_{c}\cdot {\boldsymbol{E}}$ gives electron heating, while ${{\boldsymbol{j}}}_{g}\cdot {\boldsymbol{E}}$ gives strong electron cooling. Although the net effect is electron cooling, island coalescence can be efficient in accelerating electrons to the highest energies (Oka et al. 2010).

Zoom In Zoom Out Reset image size

It has been shown that the curvature drift acceleration in the reconnection region corresponds to a Fermi-type mechanism (Dahlin et al. 2014; Guo et al. 2014, 2015). To develop a power-law energy distribution for the Fermi acceleration mechanism, the characteristic acceleration time τ_acc = 1/α needs to be smaller than the particle injection time τ_inj (Guo et al. 2014, 2015), where $\alpha =(1/\varepsilon )(\partial \varepsilon /\partial t),$ and $\partial \varepsilon /\partial t$ is the energy change rate of particles. To estimate the ordering of acceleration rate from the single-particle drift motion, consider the curvature drift velocity ${{\boldsymbol{v}}}_{c}={v}_{\parallel }^{2}{\boldsymbol{B}}\times {\boldsymbol{\kappa }}/({{\rm{\Omega }}}_{{ce}}B)$ in a curved field where ${R}_{c}={| {\boldsymbol{\kappa }}| }^{-1},$ so the time for a particle to cross this region is $\sim {R}_{c}/{v}_{\parallel }$ and the electric field is mostly induced by the Alfvénic plasma flow ${\boldsymbol{E}}\sim {-{\boldsymbol{v}}}_{{\rm{A}}}\times {\boldsymbol{B}}/c.$ The energy gain in one cycle is $\delta \varepsilon \sim {{mv}}_{{\rm{A}}}{v}_{\parallel }.$ The time for a particle to cross the island is ${L}_{\mathrm{island}}/{v}_{\parallel }.$ Then, the acceleration rate $\partial \varepsilon /\partial t\sim \varepsilon {v}_{{\rm{A}}}/{L}_{\mathrm{island}}$ for a nearly isotropic distribution. The characteristic acceleration time ${\tau }_{{acc}}\sim {L}_{\mathrm{island}}/{v}_{{\rm{A}}}.$ Taking ${L}_{\mathrm{island}}\sim 50{d}_{i}$ and ${v}_{{\rm{A}}}\sim 0.2c,$ the acceleration time ${\tau }_{{acc}}\sim 250{{\rm{\Omega }}}_{{ci}}^{-1}.$ The actual acceleration time may be longer because the outflow speed will decrease from v_A away from the X-points, and the ∇B term gives a non-negligible cooling effect. Our analysis has also found that pre-acceleration and trapping effects at the X-line region can lead to more efficient electron acceleration by the Fermi mechanism and are worthwhile to investigate further (Hoshino 2005; Egedal et al. 2015; Huang et al. 2015). Taking the main energy release phase as the injection time ${\tau }_{\mathrm{inj}}\sim 800{{\rm{\Omega }}}_{{ci}}^{-1},$ the estimated value of τ_inj/τ_acc ∼ 3.2, well above the threshold. For the case with β_e = 0.2, the ratio τ_inj/τ_acc ∼ 0.32 < 1, so there is no power-law energy distribution.

Nonthermal power-law distributions have rarely been found in previous kinetic simulations of nonrelativistic magnetic reconnection (Drake et al. 2010). We find that two essential conditions are required for producing power-law electron distribution. The first is that the domain should be large enough to sustain reconnection for a sufficient duration. A power-law tail develops as the acceleration accumulates long enough (τ_inj/τ_acc > 1). The second condition is that plasma β must be low to form a nonthermally dominated power-law distribution by providing enough free energy (∝1/β) for nonthermal electrons. Assuming 10% of magnetic energy is converted into nonthermal electrons with spectral index p = 1, one can estimate that β_e is about 0.02 for half of the electrons to be accelerated into a power law that extends to 10E_th. This agrees well with our simulation. We point out that a loss mechanism or radiation cooling can affect the final power-law index (Fermi 1949; Guo et al. 2014) of nonthermal electrons. Consequently, including loss mechanisms in a large three-dimensional open system is important, for example, to explain the observed power-law index in solar flares and other astrophysical processes. Another factor that may influence our results is the presence of an external guide field B_g. Our preliminary analysis has shown that the Fermi acceleration dominates when B_g ≲ B₀. The full discussion for the cases including the guide field will be reported in another publication (X. L. Li et al. 2015, in preparation). A potentially important issue is the three-dimensional instability, such as kink instability that may strongly influence the results. Unfortunately, the corresponding three-dimensional simulation is beyond the available computing resources. We note that results from three-dimensional simulations with pair plasmas have shown development of strong kink instability but appear to have no strong influence on particle acceleration (Guo et al. 2014; Sironi & Spitkovsky 2014). The growth rate of the kink instability can be much less than the tearing instability for a high mass ratio (Daughton 1999), and therefore the kink instability may be even less important for electron acceleration in a proton–electron plasma.

In our simulations, the low-β condition is achieved by increasing magnetic field strength (or equivalently decreasing density). We have carried out low-β simulations with the same magnetic field but lower temperature and found a similar power-law distribution (X. L. Li et al. 2015, in preparation).

The energy partition between electrons and protons shows that more magnetic energy is converted into protons. For simulations with a higher mass ratio m_i/m_e = 100, the energetic electrons still develop a power-law distribution, and the fraction of electron energy to the total plasma energy is about 33%, indicating that the energy conversion and electron acceleration are still efficient for higher mass ratios. Our results show that ions also develop a power-law energy spectrum for low-β cases and the curvature drift acceleration is the leading mechanism. However, the ion acceleration has a strong dependence on the mass ratio m_i/m_e for our relatively small simulation domain (∼100d_i). We therefore defer the study of ion acceleration to a future work (X. L. Li et al., 2015, in preparation).

The energetic electrons can generate observable X-ray emissions. As nonthermal electrons are mostly concentrated inside the magnetic islands, the generated hard X-ray flux can be strong enough to be observed during solar flares in the above-the-loop-top region (Masuda et al. 1994; Krucker et al. 2010) and the reconnection outflow region (Liu et al. 2013a). The nonthermal electrons may also account for the X-ray flares in the accretion disk corona (Galeev et al. 1979; Haardt et al. 1994; Li & Miller 1997).

In summary, we find that in a nonrelativistic low-β proton–electron plasma magnetic reconnection is highly efficient at converting the free energy stored in a magnetic shear into plasma kinetic energy and accelerating electrons into nonthermal energies. The nonthermal electrons contain more than half of the total electrons, and their distribution resembles power-law energy spectra with spectral index p ∼ 1 when particle loss is absent. This is in contrast to the high-β case, where no obvious power-law spectrum is observed. It is important to emphasize that the particle acceleration discussed here is distinct from the acceleration by shocks, where the nonthermal population contains only about 1% of particles (Neergaard Parker & Zank 2012).

We gratefully thank William Daughton for providing access to the VPIC code and for useful discussions. We also acknowledge the valuable discussions with Andrey Beresnyak and Yi-Hsin Liu. This work was supported by NASA Headquarters under the NASA Earth and Space Science Fellowship Program-Grant NNX13AM30H and by the DOE through the LDRD program at LANL and DOE/OFES support to LANL in collaboration with CMSO. Simulations were performed with LANL institutional computing.