Preferential Heating and Acceleration of Heavy Ions in Impulsive Solar Flares

The ceaseless stream of ions and electrons originating from the Sun often shows short-term enhancement in the flux of energetic charged particles known as the solar energetic particles (SEPs). While long duration gradual SEP events lasting about a few days are associated with shock waves due to coronal mass ejections (CMEs), impulsive SEP events lasting several hours appear to have originated from compact regions near the Sun, which are associated with flare events (Reames 2015). Impulsive SEPs are often characterized by a greatly enhanced ${}^{3}\mathrm{He}{/}^{4}\mathrm{He}$ ratio over the solar abundance ratio. Measurements of the flux of particles in differential energy bins suggest that a fraction of particles are accelerated to very high energies that form a non-thermal tail in energy distribution (Mason et al. 2002; Mason 2007). Additionally, ³He rich SEPs show enhancement of other heavy ions and the enhancement factor is generally larger for the heavier ions (Mason et al. 2004).

The physical mechanisms that give rise to preferential acceleration of heavier ions in impulsive SEPs are not well understood, although a wide range of observations have set some constraints on the environment where the acceleration takes place. The strong dependence of enhancement on charge to mass ratio, Q/M, of ions indicates that a Q/M dependent acceleration process is at play that disfavors processes such as the diffusive acceleration in shocks where the rate of acceleration depends on the rigidity of charged particles (see, however Lee & Wu 2000). It has been suggested that the cyclotron damping of Alfvénic turbulence by the tail of particle distribution may lead to preferential acceleration of particles with smaller Q/M and larger velocity (Eichler 1979, 2014). Particles of smaller Q/M have smaller gyrofrequency and would therefore resonate with waves of a smaller wavenumber. If the power spectrum is a sufficiently steep function of parallel wavenumber, it would imply a larger heating and acceleration rate for smaller Q/M particles. Moreover, faster moving super-Alfvénic particles can resonate with even smaller frequency waves, since the frequency of Alfvén waves is Doppler shifted in the frame of particles due to their parallel velocity along the magnetic field lines. Therefore, for super-Alfvénic particles, energy gain in one gyration time is an increasing function of particle velocity that can potentially give rise to a power-law tail in velocity distribution (Eichler 1979, 2014). On the other hand, acceleration due to cyclotron resonance is believed to be inefficient due to the anisotropic nature of the turbulence cascade that creates turbulent eddies elongated along the large-scale magnetic field (see, e.g., Strauss 1976; Montgomery & Turner 1981; Shebalin et al. 1983; Chandran 2000a, 2000b; Eichler 2014). Nevertheless, it has been argued that there may be enough power in the Alfvén waves propagating along or obliquely to the mean magnetic field to produce observed enhancement of heavier ions.

In this paper we numerically study relative heating and acceleration rates of plasma species of varying Q/M in Alfvénic turbulence. The ion species heavier than simulation ions (positrons in the case of pair plasma) are treated as test particles in our simulations (see, e.g., Dmitruk & Matthaeus 2006; Lehe et al. 2009; Markovskii & Vasquez 2010; Lynn et al. 2012; Servidio et al. 2015; González et al. 2016) and our electromagnetic PIC simulations resolve kinetic scales of all species involved. We show that the rate of energy gain is strongly dependent on Q/M. We then suggest that resonant damping of Alfvénic turbulence is a plausible scenario of preferential heating and acceleration of heavy ions as observed in ³He rich SEPs.

We use a fully electromagnetic particle-in-cell (PIC) method (Buneman et al. 1993; Kumar et al. 2015) to simulate the development of turbulence in a magnetized collisonless plasma in three spatial dimensions. The simulation box contains uniformly distributed thermal plasma embedded in a uniform grid with a background magnetic field B₀ along the x-axis. Periodic boundary conditions are imposed along all three Cartesian axes. The simulation is started by exciting six shear Alfvén waves with a wave-vector ${\boldsymbol{k}}=({k}_{x},{k}_{y},{k}_{z})$. Specifically, the initial wave vectors are chosen such that the triplet $({k}_{x}/(2\pi /{L}_{x}),{k}_{y}/(2\pi /{L}_{y}),{k}_{z}/(2\pi /{L}_{z}))$, where the L_x, L_y, and L_z are size of the box along x, y, and z axes, respectively, takes the values (1, 1, 0); (1, 2, 0); (−2, 1, 0); (−1, 0, −1); (−1, 0, −2); (2, 0, −1). The energy is equally partitioned among these waves and the amplitude of magnetic fluctuation is taken to be $0.18{B}_{0}$ for each wave. The value of initial electromagnetic field (other than the background magnetic field B₀) on the grid is the superposition of local transverse perturbations corresponding to each initial Alfvén modes. The velocity distribution of the charged particles is a drifting Maxwellian where the drift of particles is the superposition of local E × B and polarization drift due to each Alfvén modes, where E and B are local values of electric and magnetic field at the particle's initial location.

In addition to the simulation particles, we inject heavier ion species in our simulation (a few percent of the simulation particles by number) at the beginning of the simulation to study simultaneous evolution of velocity distribution of ion species characterized by differing Q/M and initial temperature. Heavier ions are treated as test particles. That is to say, heavy ions do not contribute any current in the Maxwell equations, but they are subject to the Lorentz force due to electric and magnetic fields, like the simulation particles. The assumption is justified as long as the number density of heavy ions is small or damping of turbulence due to heavy ions is small enough that the energy budget is not monopolized by them. Like simulation particles, the heavy ions are initially homogeneously distributed in the simulation box and their velocity distribution is a drifting Maxwellian. We initialize several populations of heavy ions with different charge-to-mass ratios and initial temperatures.

In order to accommodate large resonant scales of heavy ions in the simulation box while simultaneously resolving the kinetic scales of plasma, we choose to simulate electron-positron plasma (although we still refer to positrons as simulation ions). The physical dimensions of the box L_x, L_y, and L_z in units of electron skin depth are 128, 64, and 64, respectively. Electron skin depth $c/{\omega }_{{pe}}$ is resolved by 8 grid cells along each axes, where c is the speed of light and ${\omega }_{{pe}}=\sqrt{4\pi {n}_{0}{q}_{e}^{2}/{m}_{e}^{r}}$, where q_e is the charge of an electron, n₀ is the initial number density of electrons, and m^r_e is the reduced mass of an electron (for electron-positron plasma considered here reduced mass of each particle is half of its physical mass). On average, there are 12 macro-particles of electrons and ions per unit cell. The Alfvén speed, ${V}_{A}={B}_{0}/\sqrt{4\pi \rho }=0.1c$, where ρ is total initial mass density. The thermal speed of particle is such that the plasma β, which is the ratio of plasma thermal pressure to magnetic pressure, is 1/3. The plasma beta in the environments where impulsive SEPs are produced is likely much larger than 1. A large plasma beta in three-dimensional (3D) simulations requires larger simulation box size and computational time to capture cyclotron resonance. However, we substantiate our results with analytical estimates that can be extended to more realistic physical parameters. In other coronal and heliospheric environments the beta parameter may vary over an even larger range (at least 0.01–10). A numerical study of the dependence of heating and acceleration mechanisms on the beta parameter will be attempted in future simulations.

The nonlinear interactions among Alfvén waves propagating in opposite directions distort each wave packet, hence generating waves of larger wavenumbers (see, e.g., Kraichnan 1965; Howes & Nielson 2013, and the references therein). The cascade of energy from larger to smaller length scales rapidly increases power in the waves of larger wavenumbers, reaching a near maximum at $k\sim {\omega }_{{pe}}/c$ in one Alfvén crossing time ${\tau }_{A}={L}_{x}/{V}_{A}=128/{{\rm{\Omega }}}_{i}$, where ${{\rm{\Omega }}}_{i}$ is the ion (positron) cyclotron frequency (${{\rm{\Omega }}}_{i}={q}_{i}{B}_{0}/{m}_{i}c$, where q_i and m_i are the charge and mass of the simulation ion, respectively). As the turbulent energy reaches smaller spatial scales electromagnetic field fluctuations are eventually converted into the plasma heat. The amplitude of the large scale Alfvén waves decreases with time, since the turbulence is not driven (the total energy in the simulation remains almost constant and does not change by more than 0.1 percent until $1.5{\tau }_{A}$). Therefore, the turbulence simulated here does not achieve a steady state turbulence cascade.

In Figure 1 we show the temporal evolution of a one-dimensional spectra of transverse magnetic fields along parallel and perpendicular directions with respect to the background magnetic field. As evident from the spectra, the turbulent cascade of magnetic energy is anisotropic (Goldreich & Sridhar 1995; Schekochihin et al. 2009; Wan et al. 2015). The energy cascades preferentially in the direction perpendicular to the background magnetic field and creates elongated eddies.

Zoom In Zoom Out Reset image size

The turbulent cascade generates electric and magnetic field perturbations from injection to kinetic scales of the plasma. Perturbations at certain dissipative scales are absorbed by the ions and electrons and are converted into heat. The dissipation process may depend on several physical properties of the ion species, such as their charge, mass, and temperature. We have therefore placed several ion species in the simulation box to study the dependence of heating rate on charge-to-mass ratios of the ions. Figure 2 shows the evolution of the mean thermal speeds of particles of differing charge-to-mass ratios as a function of time in units of ion cyclotron time $1/{{\rm{\Omega }}}_{i}$. The thermal velocity of each particle is defined as ${{\boldsymbol{v}}}_{\mathrm{th}}={\boldsymbol{v}}-{\boldsymbol{E}}\times {B}_{0}\hat{x}/{B}_{0}^{2}$, where the electric field ${\boldsymbol{E}}$ is the local value of the electric field interpolated from the grid and ${{\boldsymbol{B}}}_{0}={B}_{0}\hat{x}$. As evident from the figure, the heating rate strongly depends on the Q/M of the ion species. The thermal energy per nucleon for the ions of smaller Q/M increases at a much faster rate than that of the ions of larger Q/M. All the ion species shown in Figure 2 were initially at the same thermal velocity. Although the rate of heating marginally depends on the initial temperature of the ions, the trend of stronger heating rates for smaller Q/M particles persists. Additionally, equal thermal velocity criteria clearly identify the differences in the heating rate solely due to differences in the cyclotron frequency of the ion species.

Zoom In Zoom Out Reset image size

The heating of heavier ion species is anisotropic in momentum space (see also Chandran et al. 2010; Xia et al. 2013; Vasquez 2015). In Figure 3, we show the instantaneous net work done on all particles by parallel and perpendicular components of the electric field. It shows that most of the heating of heavy ions is due to the perpendicular component of the electric field. Notably, the net work done by the parallel component of the electric field is marginally negative, which could possibly be due to time-dependent compressible modes with electric field along the magnetic field (see also Hollweg 1999; Li et al. 1999). Anisotropic heating of the heavy ions leads to an anisotropic temperature of the ions where the perpendicular temperature of the ions is larger than the parallel temperature. The development of temperature anisotropy is illustrated in Figure 4 for a heavy ion species.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We suggest that cyclotron resonance with Alfvén waves is the dominant mechanism for the heating of ions of smaller Q/M and hence, of smaller gyrofrequency. The perpendicular heating of the ions and a strong dependence of heating rate on Q/M are both consistent with the resonant heating, since the electric fields of Alfvén waves are in the transverse direction, and power in the resonant wave is a steeply increasing function of the M/Q. In the following section, we estimate the rate of heating of heavy ions due to cyclotron resonance with Alfvén waves.

4.1. Estimate of Cyclotron Resonant Heating Rate

The magnetic perturbations lead to random scattering of the pitch angle while the electric fluctuations stochastically alter the energy of the particles. Individual ions of charge-to-mass ratios Q/M gyrating in a large-scale background magnetic field B₀ can resonantly be accelerated or decelerated by the similarly rotating transverse electric field ${E}_{\perp }$ of Alfvén waves of frequency ω propagating along the magnetic field. The resonance condition is given by:

$$\begin{eqnarray}&&\omega \pm {k}_{\parallel }{v}_{\parallel }=n{\omega }_{c},\end{eqnarray} \tag{ 1 }$$

where n is a positive integer and ${\omega }_{c}$ is the cyclotron frequency (${\omega }_{c}={{QB}}_{0}/{Mc}$). If the resonance takes place at a uniformly random angle between the electric field and the transverse component of the velocity of the particle, the particle would perform a random walk in the momentum space and the net energy of the ensemble of particles following a Maxwellian distribution will diffusively increase with time.

For the purpose of estimating the rate of heating we consider the interaction of individual particles with turbulent electric and magnetic fields, which is considered to be an ensemble of shear Alfvén waves, although the properties of the ensemble itself are governed by the collective motion of charged particles in the plasma. The wave particle interaction picture is more justified for the species occurring in trace amounts or barely altering the turbulent dynamics, which is the case for the heavy ions in the simulations that are treated as the test particles. The change in the squared velocity of a particle (a measure of energy per unit mass) during a time interval ${\rm{\Delta }}t$ is given by:

$$\begin{eqnarray}&&{\rm{\Delta }}{v}^{2}=2Q{{\boldsymbol{E}}}_{\perp }\cdot {{\boldsymbol{v}}}_{\perp }{\rm{\Delta }}t/M=2Q(\delta {{\boldsymbol{u}}}_{\perp }\times {{\boldsymbol{B}}}_{0})\cdot {{\boldsymbol{v}}}_{\perp }/M,\end{eqnarray} \tag{ 2 }$$

where $\delta {{\boldsymbol{u}}}_{\perp }$ is the local transverse drift velocity of all simulation particles. Considering a Fourier series representation of $\delta {u}_{\perp }$ we write:

$$\begin{eqnarray}&&{\rm{\Delta }}{v}^{2}=2Q/M\,\left(\displaystyle \sum _{{\boldsymbol{k}}}\delta {{\boldsymbol{u}}}_{\perp }({k}_{x},{k}_{y},{k}_{z})\times {{\boldsymbol{B}}}_{0}\right)\cdot {{\boldsymbol{v}}}_{\perp },\end{eqnarray} \tag{ 3 }$$

where the sum is over all discrete three-dimensional wavevectors in the simulation box. In the case of a sharp resonance, only the modes satisfying the resonant condition described by Equation (1) contribute to the net change in the particle's thermal speed. Therefore, we restrict the sum to wavevectors satisfying the resonant condition, i.e., ${k}_{\parallel }^{\mathrm{res}}({V}_{A}\pm {v}_{\parallel })=n{\omega }_{c}$. We now define a dimensionless diffusion rate as, $D=\langle {({\rm{\Delta }}{v}_{\mathrm{th}}^{2}/{V}_{A}^{2})}^{2}\rangle /{{\rm{\Omega }}}_{i}{\rm{\Delta }}t$, where the brackets refer to an average over all particles constituting an ion species. The average is estimated following Section 3 of Eichler (2014), which gives:

$$\begin{eqnarray}D & \simeq & \displaystyle \frac{2\pi Q(1+{v}_{\parallel }/{V}_{A})}{M}\times \displaystyle \sum _{{k}_{y},{k}_{z}}\displaystyle \frac{{J}_{1}({{Mk}}_{r}{v}_{\perp }/{{QV}}_{A})}{{k}_{r}^{2}}\\ & & \times ({B}_{y}^{2}({k}_{\parallel }^{\mathrm{res}},{k}_{y},{k}_{z})+{B}_{z}^{2}({k}_{\parallel }^{\mathrm{res}},{k}_{y},{k}_{z}))/{B}_{0}^{2},\end{eqnarray} \tag{ 4 }$$

where ${k}_{r}=c\sqrt{{k}_{y}^{2}+{k}_{z}^{2}}/{\omega }_{{pe}}$ and J stands for Bessel functions of the first kind. Here we have considered only n = 1 resonance and eddy turnover time at the length scale $1/{k}_{\parallel }$ is taken to be ${k}_{\parallel }{V}_{A}$. The first term in the summation indicates that only a fraction of particles whose Larmor radii are smaller than the perpendicular wavelength of the oblique modes are efficiently energized due to that particular oblique mode.

The above expression for the diffusion rate can be used to obtain an estimate of the heating rate due to cyclotron resonance for any species in the simulation. However, until one Alfvén crossing time, the net magnetic power in the resonance modes for the ions species shown in Figure 2 is increasing with time as the turbulent energy cascades down from the larger length scales. Additionally, increasing the parallel velocity of the ions due to heating also enhances the rate of heating at later times. Therefore, the diffusion rate is increasing with time and the rate of heating in the simulation must be obtained from a time-dependent diffusion in momentum space. We use the mean value of ${k}_{\parallel }^{\mathrm{res}}$, ${v}_{\parallel }$, and ${v}_{\perp }$ to estimate the value D for each species as a function of time. Taking $D={D}_{0}t$ within a time interval ${\rm{\Delta }}t$ we get $\langle {\rm{\Delta }}{v}_{\mathrm{th}}^{2}/{V}_{A}^{2}\rangle \propto {\rm{\Delta }}t$. The instantaneous heating rate thus obtained at $t=45/{{\rm{\Omega }}}_{i}$ is shown in Figure 2 as straight lines. The analytically estimated rate of heating due to cyclotron resonance appears to be consistent with the observed rate of heating of the heavy ions.

In Figure 5 we show the rate of heating at different stages of the simulation for heavier ion species of decreasing charge-to-mass ratio as compared to the simulation ions. Up until about one Alfvén crossing time, the spectrum of magnetic fluctuation is steep and ions of larger M/Q are preferentially heated at a faster rate. The Alfvén waves that ions can resonate with have a parallel wavelength of ${\lambda }_{\mathrm{res}}\sim 2\pi (c/{\omega }_{{pe}})\,\times (1+{v}_{\parallel }/{V}_{A})\times M/Q$. At later stages the magnetic power in large-scale waves declines and the spectrum tends to be flatter (see Figure 1). Consequently, the ions ($Q/M\lesssim 0.5{q}_{i}/{m}_{i}$) resonating with the decaying flatter part of the spectrum are heated at slower rates and cease to enjoy a preferential heating.

Zoom In Zoom Out Reset image size

4.2. Non-thermal Acceleration

The resonance condition in Equation (1) implies that particles moving faster along the magnetic field lines are able to resonate with Alfvén waves of smaller frequency or smaller wavenumber. Therefore, rate of gain in energy is an increasing function of a particle's speed, which can lead to the acceleration of particles to very high energies. The time evolution of the distribution of speed is illustrated in Figure 6. It is clear from the left panel that the velocity distribution, which is initially a drifting Maxwellian, develops a non-thermal tail at later times.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 7, we show the history of the mean parallel and perpendicular velocities of a  few selected particles that become significantly more energetic than the particles constituting the bulk of the thermal population. It is clear that these selected particles had significantly larger speeds to begin with and were at the tail of initial Maxwellian distribution.

We have simulated a turbulence cascade in a magnetized plasma using a fully electromagnetic PIC method that is suitable for modeling kinetic physics of all plasma species. The relative rate of absorption of turbulent energy trickling down from larger scales by particles of different charge-to-mass ratios was studied by simultaneously evolving all particles in the phase-space. We find that the plasma species with smaller charge-to-mass ratios are heated and accelerated at a significantly higher rate. We have suggested that the heating and acceleration of these heavier ions are predominantly due to cyclotron resonance with Alfvén waves.

The physical size of the simulation box limits the extent of inertial scale in the simulation that determines the range of value of the Q/M of heavy ions that are preferentially heated and accelerated. The range of mass-to-charge ratios and the temporal and spatial scales involved in impulsive SEP is much larger than the scales simulated here. However, the rate of heating of heavy ions observed in the simulation is consistent with our analytical estimates of the heating expected due to the cyclotron resonance. In impulsive SEPs where the resonance can take place over a rather extended inertial range of the turbulence, the cyclotron resonance implies large enhancement of even ultra-heavy ions.

The sources of impulsive SEPs are not entirely known but their association with Type III ratio bursts (Wang et al. 2012) has led to identification of their origin as extended regions in the solar corona associated with coronal jets (Nitta et al. 2015). The variation in the energy spectrum of particles and enhancement ratios, which does not show any apparent correlation with the properties of the sources, suggests that the observed heating may take place over a wide range of plasma parameters, which favors a rather robust physical process viable over varying source conditions. The observed features of ³He rich SEPs appear to be consistent with the generic nature of the heating and acceleration in anisotropic Alfvénic turbulence. However, the rate of heating via cyclotron resonance strongly depends on the spectrum of Alfvénic fluctuations and the extent of inertial range of the turbulence that determines total power in the Alfvén waves at the resonate scale of the particles. A steeper spectrum enhances differentiation in the resonant heating of particles of varying M/Q. On the other hand, a steeper spectrum implies smaller power in the resonant modes at the kinetic scales and therefore a larger heating and acceleration time for all plasma species. For an estimated best-fit temperature of $\,{10}^{7}\,{\rm{K}}$ at the source of impulsive SEPs, preferential acceleration is observed for nuclei with $M/Q\gtrsim 10$. In the case of acceleration due to turbulence, it implies that the scale at which turbulence is isotopic is at least a hundred times larger than the mean gyroradius of protons. The future observations, numerical simulations, and better statistics of events might shed more light on the environments and mechanisms responsible for acceleration and heating of ions in impulsive SEPs.

The preferential heating of heavy ions is evident in other heliospheric environments as well, such as in the solar corona and the fast solar wind (von Steiger et al. 1995; Kohl et al. 2006; Landi & Cranmer 2009; Kasper et al. 2013; Tracy et al. 2015). While some of the enhancement of heavy nuclei could be due to other physical mechanisms, such as reacceleration in shocks, magnetic reconnections (Drake et al. 2009), and resonance with ion cyclotron waves (Fisk 1978; Knizhnik et al. 2011), cyclotron damping of Alfvénic turbulence may still be playing an important or possibly a dominant role (Cranmer 2014) in heating and acceleration process.

We thank M. Kunz, A. Philippov, and A. Schekochihin for their helpful discussions. We acknowledge the support from the Israel-U.S. Binational Science Foundation, the Israeli Science Foundation (ISF), the ISF-University Grant Commission (India), NSF grant AST-1517638, and the Joan and Robert Arnow Chair of Theoretical Astrophysics (D.E.). M.G. is supported by NASA through Einstein Postdoctoral Fellowship Award Number PF-160137 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-03060. Simulations in this paper used computational resources supported by PICSciE-OIT High Performance Computing Center and the NASA/Ames HEC Program (SMD-16-6751; Pleiades).