Preferential Heating and Acceleration of Heavy Ions in Impulsive Solar Flares

TL;DR

This study uses 3D PIC simulations to demonstrate that in impulsive solar flares, heavy ions are preferentially heated and accelerated perpendicular to the magnetic field, aligning with observed non-thermal features.

Contribution

It provides the first detailed simulation-based evidence that heavy ions are preferentially energized via cyclotron resonance in turbulent solar flare environments.

Findings

Heavy ions gain energy faster than lighter ions.

Ion heating is predominantly perpendicular to the magnetic field.

Simulation results match analytical estimates of cyclotron resonance heating.

Abstract

We simulate decaying turbulence in a homogeneous pair plasma using three dimensional electromagnetic particle-in-cell (PIC) method. A uniform background magnetic field permeates the plasma such that the magnetic pressure is three times larger than the thermal pressure and the turbulence is generated by counter-propagating shear Alfv\'en waves. The energy predominately cascades transverse to the background magnetic field, rendering the turbulence anisotropic at smaller scales. We simultaneously move several ion species of varying charge to mass ratios in our simulation and show that the particles of smaller charge to mass ratios are heated and accelerated to non-thermal energies at a faster rate, in accordance with the enhancement of heavy ions and non-thermal tail in their energy spectrum observed in the impulsive solar flares. We further show that the heavy ions are energized mostly in the direction perpendicular to the background magnetic field with a rate consistent with our analytical estimate of the rate of heating due to cyclotron resonance with the Alfv\'en waves of which a large fraction is due to obliquely propagating waves.

Equations

$$\omega \pm k_{\parallel }{v}_{\parallel }=n{\omega }_c$$

$$\Delta v^2=2Q{\vec{E}}_{\perp }\cdot {\vec{v}}_{\perp }\Delta t/M=2Q(\delta {\vec{u}}_{\perp }\times {\vec{B}}_0)\cdot {\vec{v}}_{\perp }/M$$

$$\Delta v^2=2Q/M(\sum _{\vec{k}}\delta {\vec{u}}_{\perp }(k_x,k_y,k_z)\times {\vec{B}}_0)\cdot {\vec{v}}_{\perp }$$

$$\begin{array}{rl}D\simeq & \frac{2\pi Q(1+v_{\parallel }/V_A)}{M}\times \sum _{k_y,k_z}\frac{J_1(M{k}_r{v}_{\perp }/Q{V}_A)}{k_r^2}\\ & \times (B_y^2(k_{\parallel }^{res},k_y,k_z)+B_z^2(k_{\parallel }^{res},k_y,k_z))/B_0^2\end{array}$$

Full text

11affiliationtext: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA22affiliationtext: Department of Physics, Ben-Gurion University, Be’er-Sheba 84105, Israel33affiliationtext: Einstein and Spitzer Fellow

1. Introduction

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 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 the flare events (Reames, 2015). Impulsive SEPs are often characterized by greatly enhanced ${}^{3}He/^{4}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 which form a non-thermal tail in energy distribution (Mason et al., 2002; Mason, 2007). Additionally, ${}^{3}He$ rich SEPs show enhancement of other heavy ions as well, 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 which disfavors processes such as the diffusive acceleration in shocks where the rate of acceleration depends on rigidity of charged particles (see, however Lee & Wu (2000)). It has been suggested that the cyclotron damping of Alfvenic 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 gyro-frequency and would therefore resonate with waves of 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-Alfvenic particles can resonate with even smaller frequency waves since in the frame of particles the frequency of Alfven waves is Doppler shifted due to their parallel velocity along the magnetic field lines. Therefore, for super-Alfvenic particles energy gain in one gyration time is an increasing function of particle velocity which 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 anisotropic nature of the turbulence cascade which creates turbulent eddies elongated along the large scale magnetic field (see e.g., Strauss, 1976; Montgomery & Turner, 1981; Shebalin et al., 1983; Chandran, 2000b, a; Eichler, 2014). Nevertheless, it has been argued that there may be enough power in the Alfven waves propagating along or obliquely to the mean magnetic field to produce observed enhancement of heavier nuclei.

In this paper we numerically study relative heating and acceleration rate of plasma species of varying Q/M in a Alfvenic 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 Alfvenic turbulence is a plausible scenario of preferential heating and acceleration of heavy ions as observed in He3 rich SEPs.

2. Numerical Simulation

We use fully electromagnetic particle-in-cell (PIC) method (Buneman, 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 Cartesian grid with a background magnetic field $B_{0}$ along the x-axis. Periodic boundary conditions are imposed along all three cartesian axes. The simulation is started by exciting six shear Alfven waves with wave vector $\vec{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 $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.18B_{0}$ for each wave. The value of initial electromagnetic field (other than the background magnetic field $B_{0}$) on the grid is superposition of local transverse perturbations corresponding to each initial Alfven modes. The velocity distribution of the charged particles is a drifting Maxwellian where the drift of a particles is the superposition of local $E\times B$ and Polarization drift due to each Alfven modes, where E and B are local value 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 (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’s 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 number density of heavy ions is small or damping of turbulence due to heavy ions is small enough that 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 ratio 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 dimension 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^{r}_{e}}$, where $q_{e}$ is the charge of an electron, $n_{0}$ 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 Aflven speed $V_{A}=B_{0}/\sqrt{4\pi\rho}=0.1c$, where $\rho$ is total initial mass density. The thermal speed of particle is such that the plasma $\beta$, which is 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 3D simulations requires larger simulation box size and computational time to capture cyclotron resonance. However, we substantiate our results with analytical estimates which can be extended to more realistic physical parameters. In other coronal and heliospheric environments the beta parameter may vary over even larger range (at least 0.01 to 10). A numerical study of the dependence of heating and acceleration mechanism on the beta parameter would be attempted in future simulations.

3. Turbulence Cascade

The non-linear interactions among Alfven waves propagating in opposite directions distorts each wave packet, hence generating waves of larger wave numbers (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 wave numbers, reaching a near maximum at $k\sim\omega_{pe}/c$ in one Alfven crossing time $\tau_{A}=L_{x}/V_{A}=128/\Omega_{i}$, where $\Omega_{i}$ is the ion (positron) cyclotron frequency ($\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). In units of As the turbulent energy cascades reaches smaller spatial scales electromagnetic field fluctuations are eventually converted into the plasma heat. The amplitude of the large scale Alfven 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 one-dimensional spectrum of transverse magnetic field 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.

4. Heating and Acceleration

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 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 dependence of heating rate on charge to mass ratios of the ions. Figure 2 shows the evolution of mean thermal speed of particles of differing charge to mass ratio as a function of time in units of ion cyclotron time $1/\Omega_{i}$. The thermal velocity of each particle is defined as $\vec{v}_{th}=\vec{v}-\vec{E}\times{B}_{0}\hat{x}/B_{0}^{2}$, where electric field $\vec{E}$ is the local value of the electric field interpolated from the grid and $\vec{B}_{0}=B_{0}\hat{x}$. As evident from the figure, the heating rate strongly depends on the the Q/M of 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 rate for smaller Q/M particles persists. Additionally, equal thermal velocity criteria clearly identifies the differences in the heating rate solely due to differences in the cyclotron frequency of the ion species.

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 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 Li et al., 1999; Hollweg, 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 specie.

We suggest that cyclotron resonance with Alfven waves is the dominant mechanism for the heating of ions of smaller Q/M and hence of smaller gyro-frequency. 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 field of Alfven waves are in the transverse direction and power in the resonant wave is a steeply increasing function of the M/Q. In the following we estimate the rate of heating of heavy ions due to cyclotron resonance with Alfven 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 large scale background magnetic field $B_{0}$ can resonantly be accelerated or decelerated by the the similarly rotating transverse electric field $E_{\perp}$ of Alfven waves of frequency $\omega$ propagating along the magnetic field. The resonance condition is given by

[TABLE]

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 make random walk in the momentum space and the net energy of the ensemble of particles following Maxwellian distribution would diffusively increase with time.

For the purpose of estimating the rate of heating we consider interaction of individual particles with turbulent electric and magnetic field which is considered to be an ensemble of shear Alfven waves, although the properties of the ensemble itself is 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 amount or barely altering the turbulent dynamics, which is the case for the heavy ions in the simulations which 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 $\Delta t$ is given by

[TABLE]

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

[TABLE]

where sum is over all discrete three dimensional wave vectors 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 a net change in the particle’s thermal speed. Therefore, we restrict the sum to wave vectors satisfying the resonant condition, i.e. $k^{res}_{\parallel}(V_{A}\pm v_{\parallel})=n\omega_{c}$. We now define a dimensionless diffusion rate as, $D=<(\Delta v^{2}_{th}/V^{2}_{A})^{2}>/\Omega_{i}\Delta t$, where the brackets refer to an average over all particles constituting an ion species. The average is estimated as in the Section 3 of Eichler (2014), which gives

[TABLE]

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 turn over 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 radius is 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 Alfven 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 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^{res}_{\parallel}$, $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 $\Delta t$ we get $<\Delta v_{th}^{2}/V^{2}_{A}>\propto\Delta t$. The instantaneous heating rate thus obtained at $t=45/\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 on species of decreasing charge to mass ratio as compared to the simulation ions. Up until about one Alfven crossing time the spectrum of magnetic fluctuation is steep and ions of larger M/Q are preferentially heated at a faster rate. The Alfven waves that ions can resonate with have a parallel wavelength of $\lambda_{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.5q_{i}/m_{i}$ ) resonating with the decaying flatter part of the spectrum are heated at slower rate and cease to enjoy a preferential heating.

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 Alfven waves of smaller frequency or smaller wave number. Therefore, rate of gain in energy is an increasing function of particle’s speed which can lead to 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.

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

5. Discussion and Conclusions

We have simulated turbulence cascade in a magnetized plasma using fully electromagnetic PIC method which 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 is studied by simultaneously evolving all particles in the phase-space. We find that the plasma species of 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 Alfven waves.

The physical size of the simulation box limits the extent of inertial scale in the simulation which determines the range of value of 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 SEPs 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 rather extended inertial range of the turbulence the cyclotron resonance would imply 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 suggest 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 ${}^{3}He$ rich SEPs appear to be consistent with the generic nature of the heating and acceleration in anisotropic Aflvenic turbulence. However, the rate of heating via cyclotron resonance strongly depends on the spectrum of Alfvenic fluctuations and the extent of inertial range of the turbulence which determines total power in the Alfven 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 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}K$ at the source of impulsive SEPs, preferential acceleration is observed for nuclei with $M/Q\gtrsim 10$. In case of acceleration due to the turbulence it would imply that the scale at which turbulence is isotopic is at least hundred times larger than the mean gyroradius of protons. The future observations, numerical simulations, and a 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 Solar corona and fast solar wind(Kohl et al., 2006; Landi & Cranmer, 2009; von Steiger et al., 1995; Kasper et al., 2013; Tracy et al., 2015). While some of 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 Alfvenic 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 helpful discussions. We acknowledge 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).