Scattering of energetic electrons by heat-flux-driven whistlers in flares

TL;DR

This study uses particle-in-cell simulations to investigate how heat-flux-driven whistler waves scatter energetic electrons in solar flares, revealing potential confinement mechanisms for high-energy particles.

Contribution

It demonstrates that heat-flux-driven whistlers can effectively scatter energetic electrons, contributing to their confinement in flare environments, a novel insight into flare particle transport.

Findings

Whistler waves are driven by heat flux and cold return currents.

Energetic electrons resonate with whistlers and are scattered within a hundred electron cyclotron times.

Electron scattering reduces heat flux and increases perpendicular energy, aiding confinement.

Abstract

The scattering of electrons by heat-flux-driven whistler waves is explored with a particle-in-cell (PIC) simulation relevant to the transport of energetic electrons in flares. The simulation is initiated with a large heat flux that is produced using a kappa distribution of electrons with positive velocity and a cold return current beam. This system represents energetic electrons escaping from a reconnection-driven energy release site. This heat flux system drives large amplitude oblique whistler waves propagating both along and against the heat flux, as well as electron acoustic waves. While the waves are dominantly driven by the low energy electrons, including the cold return current beam, the energetic electrons resonate with and are scattered by the whistlers on time scales of the order of a hundred electron cyclotron times. Peak whistler amplitudes of $\tilde{B} / B_{0} \sim 0.125$ and angles of $\sim 60 \degree$ with respect to the background magnetic field are observed. Electron perpendicular energy is increased while the field-aligned electron heat flux is suppressed. The resulting scattering mean-free-paths of energetic electrons are small compared with the typical scale size of energy release sites in flares, which might lead to the effective confinement of energetic electrons that is required for the production of very energetic particles.

Equations

$$\begin{array}{r}{f}_{h,\kappa }=\frac{n_0\Gamma (\kappa +1)}{{\pi }^{3/2}{\theta }_{h,\perp }^2{\theta }_{h,x}{\kappa }^{3/2}\Gamma (\kappa -1/2)}\times \\ {\left[1+\frac{v_x^2}{\kappa {\theta }_x^2}+\frac{v_{\perp }^2}{\kappa {\theta }_{\perp }^2}\right]}^{-(\kappa +1)}\Theta (v_x).\end{array}$$

$${\theta }_{hx,\perp }={\left[(\kappa -3/2)/\kappa \right]}^{1/2}{V}_{Thx,\perp }$$

$$\begin{array}{r}{f}_c=\frac{n_0}{{\pi }^{3/2}}\frac{e^{-[(v_x+v_d{)}^2+v_{\perp }^2]/v_{Tc}^2}}{v_{Tc}^3(1+\text{erf}(v_d/v_{Tc}))}\Theta (-v_x),\end{array}$$

$$\omega =\frac{|k_x|k{d}_e^2{\Omega }_e}{1+k^2{d}_e^2}$$

$$\omega -k_x{v}_{x,r}-\frac{n{\Omega }_{e0}}{\gamma }=0$$

$$\Delta v_{x,n}=2\sqrt{2}{\left|\frac{{\Omega }_e}{k_x}\frac{\tilde{B}}{B_0}\left[\frac{\omega }{|k_x|}\frac{k_y}{k}{k}^2{d}_e^2{J}_n+\frac{v_{\perp 0}}{2}\left(\left(\frac{k}{|k_x|}-1\right)J_{n+1}+\left(\frac{k}{|k_x|}+1\right)J_{n-1}\right)\right]\right|}^{1/2}$$

$$\Delta v_{x,0}=2\sqrt{2}\sqrt{\frac{{\Omega }_e}{k_x}\frac{\tilde{E}}{B_0}c}$$

$$v_{r,x}=\frac{v_{p,x}}{1+{\alpha }_n^2}\pm \sqrt{\frac{{\alpha }_n^2}{1+{\alpha }_n^2}\left(c^2-\frac{v_{p,x}^2}{1+{\alpha }_n^2}\right)}$$

$$\frac{{\left(v_x-v_{x,c}\right)}^2}{R/\left({\gamma }_0^2+v_{p,x}^2/c^2\right)}+\frac{v_y^2}{R/{\gamma }_0^2}=1$$

Full text

flares — particle-in-cell simulation — whistler waves — space plasma

††thanks: Work was carried out at the

University of Maryland, College Park

1 Introduction

Magnetic reconnection is the driver of explosive energy release in the sun’s corona (Forbes, 1988) and in plasma environments throughout the universe (Michel, 1994). In solar flares electrons can have energies exceeding a MeV (Lin et al., 2003; Gary et al., 2018) and in astrophysical environments can reach a PeV (Abdo et al., 2011). An important question that arises in building a model for electron acceleration in the sun and elsewhere is how electrons remain within the region where magnetic energy is being released for a sufficiently long time to reach the energies seen in observations.

The magnetic energy release rate in flares is controlled by magnetic reconnection in which the upper limit on the rate of reconnection is of order of $0.1\>V_{A}$ where $V_{A}$ is the local Alfvén speed. For typical parameters ($B\sim 50\>\text{G}$, $n\sim 10^{9}\>\text{cm}^{-3}$) $V_{A}$ is around $3000$ km/s. The scale size of magnetic energy release in large solar flares can reach $10^{4}$ km so energy release rates are tens of seconds. In the absence of a confinement mechanism the transit time of a relativistic electron out of this region is of the order of $0.03\>\text{s}$. Consistent with these estimates, the decay time of hard X-ray emission from flares exceeds the transit time of energetic electrons across the source by two orders of magnitude (Masuda et al., 1994; Krucker et al., 2007, 2010).

Thus, some mechanism for electron confinement is necessary to hold electrons in the energy release zone for a long enough time to reach relativistic energies. Because the Larmor radius of even relativistic electrons is small ($\simeq 7\>\text{cm}$ for the given numbers), one possibility is a form of magnetic confinement. However, magnetic field lines that link flare energy release sites end either on the chromosphere or the solar wind so confinement is not effective unless the electrons can mirror. Significant mirroring of energetic electrons would require that their velocity have a large component perpendicular to the ambient magnetic field. Recent models of electron acceleration during reconnection suggest that electron energy gain is mostly parallel to $\mathbf{B}$ and is dominated either by parallel electric fields or Fermi reflection (Drake et al., 2006; Dahlin et al., 2014, 2016). The consequence is that electron distribution functions are strongly anisotropic even in 3D reconnecting systems, which become turbulent with the development of multiple x-lines and chaotic magnetic fields (Dahlin et al., 2017).

There is evidence from RHESSI spacecraft observations that the electron energy flux in the source regions of flares can exceed that measured at the chromosphere by up to an order of magnitude (Simões et al., 2013). It has been suggested that double layers driven by the return current of cold electrons interacting with ambient ions would cause reflection of some hot electrons. However, since the resulting potential drop across the double layers scales like $T_{eh}\ll m_{e}c^{2}$ (Li et al., 2013, 2014) with $T_{eh}$ the temperature of hot escaping electrons, the amplitude of double layers is not sufficient to trap relativistic electrons. Further, the direct measurement of radio emission from gyro-synchrotron emission from energetic electrons in flares requires significant energy in the perpendicular motion of electrons (Gary et al., 2018). Thus, if the dominant acceleration mechanism of electrons in flares is parallel to the ambient magnetic field, a mechanism is required to scatter the parallel energy into perpendicular energy.

A potential scattering mechanism for energetic electrons is via oblique whistler waves (e.g. Artemyev et al. (2012, 2014)). The whistler resonance condition with electrons is given by $\omega-k_{\parallel}v_{\parallel}-n\Omega_{e}/\gamma=0$, where $\Omega_{e}=eB_{0}/(m_{e}c)$ is the electron cyclotron frequency, $n$ is an integer that can take on positive and negative values (Krall & Trivelpiece, 1986) and $\gamma$ is the relativistic Lorentz factor. For typical waves with $kd_{e}\sim 1$ the resonant velocities are given by $v_{\parallel}\sim nV_{Ae}$, where $d_{e}$ is the electron skin depth and $V_{Ae}$ is the electron Alfvén speed. Scattering by an oblique whistler at each resonance can be efficient. Furthermore when overlap of resonances occurs, scattering can be strongly increased for a large fraction of electrons. (Roberg-Clark et al., 2016; Karimabadi et al., 1992). Kinetic simulations with boundary conditions that impose a fixed temperature jump along an ambient magnetic field have established that in high $\beta$ systems the fluctuating magnetic field from whistlers is comparable to the initial ambient magnetic field, which is sufficient to strongly scatter electrons. This limits their effective streaming velocity to the whistler phase speed, which in a $\beta\sim 1$ system is of order $V_{Ae}$.

In these previous works treating high-$\beta$ systems (Roberg-Clark et al., 2018a, b; Komarov et al., 2018) the imposed electron distribution functions were Maxwellians with specified temperature jumps that drove a heat flux. Here we consider a system in which the heat flux arises from electrons with a $\kappa$ distribution propagating in one direction and a cold electron beam that produces a return current. In this system oblique whistlers develop that propagate both along and against the heat flux. Oblique whistlers propagating along the direction of the heat flux are driven at early time by the “fan” instability of the (anomalous) $n=-1$ cyclotron resonance (Kadomtsev & Pogutse, 1968; Haber et al., 1978; Fülöp et al., 2006; Krafft & Volokitin, 2010; Vasko et al., 2019; Verscharen et al., 2019). At later time the growth of parallel-propagating electron acoustic waves facilitates the extraction of energy from the bulk of the hot electrons in tandem with the oblique whistlers. Whistlers propagating against the heat flux are driven by the Landau resonance with the return current electrons. Both classes of whistlers resonate with and scatter the most energetic electrons in the tail of the $\kappa$ distribution, reducing their heat flux and substantially increasing their velocity perpendicular to the magnetic field.

2 Simulation Method

We carry out a two-dimensional (2D) simulation using the PIC code $\tt{p3d}$ (Zeiler et al., 2002) to model energetic electrons and a cold return current electron beam. $\tt{p3d}$ calculates particle trajectories using the relativistic Newton-Lorentz equations and the electromagnetic fields are advanced using Maxwell’s equations. An initially uniform magnetic field $\mathbf{B_{0}}=B_{0}\mathbf{\hat{x}}$ threads the plasma. $v_{x}$ is therefore the parallel velocity and $v_{y}$ and $v_{z}$ are the perpendicular velocities, while $y$ is the perpendicular spatial coordinate. Boundary conditions are periodic in $x$ and $y$. The initial particle distribution function has two components. The first is a parallel ($v_{x}>0$) bi-$\kappa$ distribution with temperatures $T_{h,x}\gg T_{h,\perp}$,

[TABLE]

$n_{0}$ is the initial density of each of the electron components, $\Gamma$ is the gamma function,

[TABLE]

are the effective thermal speeds, $V_{Thx,\perp}=\sqrt{2T_{hx,\perp}/m_{e}}$ are the regular thermal speeds, $\Theta(v_{x})$ is the Heaviside step function and $\kappa$ is a parameter that tunes the steepness of the nonthermal tail of the distribution.

The second electron component is the cold return current beam (moving against $B_{0}$) which takes the form of a drifting isotropic Maxwellian,

[TABLE]

where $V_{Tc}=\sqrt{2T_{c}/m_{e}}$ is the cold thermal speed and $v_{d}$ is a drift speed that ensures zero net current ($\langle v_{\parallel}\rangle=0$) in the initial state while the error function $\text{erf}(v_{d}/v_{Tc})$ makes the density of hot and cold particles equal. This choice is motivated by observations of flares (Krucker et al., 2010; Oka et al., 2013) suggesting that the population of energetic electrons is large.

For the simulation presented we chose $\kappa=4$, $T_{x}=20\>T_{\perp}=20\>T_{c}$, $\beta_{e0h}=8\pi nT/B^{2}_{0}\sim 8\pi n_{0}T_{hx}/B_{0}^{2}=2$. While this is a relatively high $\beta$ for the corona , a $\beta\sim 1$ system is consistent with specific flare observations (Krucker et al., 2010), is inferred more generally from rough equipartition of energy release between heated and energetic ions and electrons in flares (Emslie et al., 2005, 2012), and is seen in reconnection simulations (Dahlin et al., 2017). The large value of $T_{hx}/T_{h\perp}$ is motivated by 3D PIC simulations of reconnection in Dahlin et al. (2017) showing $P_{\parallel}/P_{\perp}\sim 100$. This system is marginally stable with respect to the fluid firehose criterion, i.e. $\beta_{\parallel}-\beta_{\perp}\lesssim 2$.

The simulation domain lengths are $L_{x}=L_{0}=163.84\ d_{e}$ and $L_{y}=L_{0}/2$, where $d_{e}=c/\omega_{pe}$ is the electron skin depth and $\omega_{pe}=(4\pi n_{0}e^{2}/m_{e})^{1/2}$ is the electron plasma frequency. Other parameters in the simulation include $\omega_{pe}/\Omega_{e0}=5\sqrt{2}$, and $T_{hx}/(m_{e}c^{2})=0.02$, which sets $v_{Thx}/c=1/5$. The characteristic velocity of whistlers depends on the wavelength but has an upper limit that scales with the electron Alfvén speed $V_{A,e}=d_{e}\Omega_{e0}$, where $\Omega_{e0}=eB_{0}/m_{e}c$. Ions, with mass ratio $m_{i}/m_{e}=1600$, are initialized with a Maxwellian distribution of temperature $T_{i0}=T_{eh}/2$ and do not play a significant role in the simulation. The simulation uses $560$ particles per species per cell, has a grid of $4096$ by $2048$ cells, and is run to the time $t\Omega_{e0}=1332$.

3 Simulation Results

The initial distribution function [the total of Eqs. (1) and (3)] drives magnetic fluctuations unstable in the system. Shown in Fig. 1(a) is the time evolution of the box-averaged magnetic fluctuation energies $\langle B^{2}_{z}\rangle$ and $\langle\delta B^{2}_{x}\rangle$, where $\delta B_{x}=B_{x}-B_{0}$. The energy in $B_{z}$ peaks at around $t\Omega_{e0}=190$ after a fast growth phase while a second peak occurs at $t\Omega_{e0}\sim 375$. The magnetic fluctuations damp significantly by the end of the simulation, approaching early-time noise levels.

To establish that the growth of the fluctuations is associated with the scattering of electrons, we track the time dependence of the electron distribution function $f_{e}(v_{x},v_{y})$, with $v_{z}$ as well as the two spatial coordinates averaged out. We use $v_{x}$ and $v_{y}$ as proxies for the parallel and perpendicular velocities since $\mathbf{B}\sim B_{0}\hat{x}$ and the distribution remains gyrotropic $(v_{y}\sim v_{z}\sim v_{\perp})$ for the duration of the simulation (not shown). We divide velocity space into ranges of width $\delta v=1.75V_{Ae0}$, treating $v=\sqrt{v^{2}_{x}+v^{2}_{y}}$ as the total velocity. In Fig. 1(b) the time dependence of $\langle v^{2}_{y}\rangle$ is shown for each velocity range with labels indicating the velocity intervals. For each of the energy ranges $\langle v^{2}_{y}\rangle$ increases with time. Figure 2 reveals that the distribution of energy $E=(\gamma-1)m_{e}c^{2}$ is almost unchanged during the course of the simulation since the initial and late-time distributions (at $t\Omega_{e0}=0$ and $710$) basically overlap for most energies (and only vary slightly at high energy where there are few particles). Thus, the increase in $\langle v^{2}_{y}\rangle$ in time corresponds to a decrease in $\langle v^{2}_{x}\rangle$ so the dominant scattering is in pitch angle. The sharpest increase in $\langle v^{2}_{y}\rangle$ occurs at around $t\Omega_{e0}=300$. The characteristic time over which $\langle v^{2}_{y}\rangle$ increases is roughly $100\>\Omega^{-1}_{e0}$, which we take to be the scattering time of the electrons. Saturation of $\langle v^{2}_{y}\rangle$ (the “rollover”) takes place for a long period of time starting around $t\Omega_{e0}=400$ and continues until the end of the simulation.

Scattering of the energy from $v_{x}$ to $v_{y}$ also reduces the energy flux, which is shown in Fig. 1(c) where we plot the time dependence of $\langle v_{x}v^{2}\rangle$ for the same energy bins as in Fig. 1b. Note that the data is presented on a log scale and we do not divide by the number of particles in each velocity range for this quantity. Curves representing most of the velocity ranges show a monotonic decrease of $\langle v_{x}v^{2}\rangle$ with time, in some cases resulting in a drop by a factor of $2$. An exception is the curve corresponding to $0<v/V_{ae0}<1.75$ which increases until it peaks at $t\Omega_{e0}\sim 190$ and then drops off.

To establish the nature of the fluctuations that develop in the simulation we show in Fig. 3 the structure and motion of the out-of-plane magnetic fluctuations $B_{z}$ and the parallel electric field fluctuations $E_{x}$. Figure 3a is a 2D plot of $B_{z}$ at early time, $t\Omega_{e0}=175$, during the first growth phase. The fluctuations travel at angles of roughly $60$ degrees relative to $B_{0}$ as predicted in Verscharen et al. (2019) for whistler scattering of energetic particles in a low-$\beta$ system. The Verscharen et al. result is constrained by the location of the $n=-1$ resonance and a minimization of Landau damping, which is ensured if the whistler phase speed is large enough compared to the core electron thermal speed. In our simulation, however, the particle distribution is flat in $v_{\parallel}$ near the whistler phase speed and so the whistler does not suffer from large Landau damping in the $\beta\sim 1$ case. Rather, the angle of $60\degree$ maximizes the wave-particle interaction for the $n=-1$ resonance (equation 6). The crossed interference pattern that emerge in the simulation seem to be a consequence of the symmetry between the plus and minus $y$ directions. For every oblique wave generated with a velocity in the plus $y$ direction, there will be a companion wave with a velocity in the negative $y$ direction, which leads to a standing-wave-like pattern. Spatial inhomogeneities in the wave pattern are probably due to variations in noise levels in the system at early time.

A spacetime diagram of $B_{z}$ at a cut at $y=64.32\>d_{e}$ is shown in Fig. 3c. Waves moving in the $+\hat{x}$ direction become visible around $t\Omega_{e0}=200$ (recall the first peak in the fluctuating $B_{z}$ in Fig. 1b at $\Omega_{e0}t\sim 175$) but are beginning to slow down and reverse their direction by this time. The parallel phase speed of these waves at early time ($t\Omega_{e0}\sim 100$) is $v_{p,x}/V_{Ae0}=(\omega/k_{x})/V_{Ae0}\sim+0.5$ while at later time the leftward-propagating waves move at roughly $v_{p,x}/V_{Ae0}\sim-0.5$ with noticeably longer parallel wavelengths.

In Fig. 3d are plots of $B_{y}$ and $B_{z}$ along a cut at $y/d_{e}=20.48$ revealing what appears to be a $90$ degree phase shift between the two components. We also show a hodogram of these quantities from $x/d_{e}=40.96$ to $81.92$ in Fig. 3e demonstrating right-handed elliptical polarization, confirming that these are indeed whistlers.

The cold plasma dispersion relation for whistlers (neglecting the displacement current) is

[TABLE]

with $k_{x}>0$ ($<0$) corresponding to a “rightward” (“leftward”) whistler propagating along (against) $B_{0}$. From Eq. (4) the phase speed of $0.5\>V_{Ae0}$ seen at early time in the simulation corresponds to a wave with $kd_{e}\simeq 1$, a value that is consistent with the wavelength of fluctuations in Fig. 3a. The fluctuations in the electric field at $t\Omega_{e0}=175$ (Fig. 3b) show crossed patterns similar to those in 3a. These are the parallel electric fields of the rightward moving oblique whistlers, which are significant for large oblique angles and for $kd_{e}\gtrsim 1$. In Fig. 3b there are shorter-wavelength, parallel-propagating modes with $k\lambda_{De}\sim 10$ and $\omega\simeq 0.4\>\omega_{pe}$ which are electron acoustic waves (EAWs) [$\lambda_{De}=V_{Te}/(\sqrt{2}\omega_{pe})$ is the electron Debye length]. The waves grow through the Landau resonance between the thermal speeds of the hot and cold electron components and reduce the relative drift between those components (Gary & Tokar, 1985; Agapitov et al., 2018; Vasko et al., 2018). These modes are of large amplitude, $\tilde{E}/((V_{Ae0}/c)B_{0})\sim 0.3$ (note Fig. 3b has its color scale capped at $0.15$ to bring out the whistler signal).

Figures 4a-c illustrate the scattering of the electron distribution function as time proceeds in the simulation. The color plot shows the contours of constant $f(v_{x},v_{y})$. At $t=0$ (Fig. 4a) the large discontinuity in $f$ at $t=0$ separates the return current beam with $v_{x}<0$ and the $\kappa$ distribution with $v_{x}>0$. At $t\Omega_{e0}=177$ (Fig. 4b) the distribution develops horn-like structures near $v_{x}/V_{Ae0}\simeq-1,1.8,3,$ and $3.5$ that demonstrate that particles from the initial distribution near the $v_{y}$=0 axis have been scattered to higher $v_{y}$ and lower $v_{x}$. The largest number of scattered particles is in the structure at $v_{x}/V_{Ae0}=1.8$. The discontinuity in Fig. 4a has been filled in and the contours of the distribution are fairly flat in the vicinity of the whistler phase speed $v_{p,x}/V_{Ae0}\sim 0.5$. We attribute the flattening to large-amplitude electrostatic fluctuations that quickly grow up and damp in the simulation (not shown). When $t\Omega_{e0}=355$ (4c), the distribution is significantly more isotropic in the $v_{x}>0$ half-plane. While some particles have been scattered to $v_{x}<0$, most of the scattering seems to be limited to $v_{x}>0$, suggesting that $\langle v_{y}^{2}\rangle$ saturates in Fig. 1a because the $v_{x}>0$ half-plane has become nearly isotropic. Figure 4c is therefore representative of the late-time structure of the distribution function.

4 Resonances

To explain scattering in the simulation we invoke the basic theory of resonant interaction between oblique whistlers and electrons (see e.g. Roberg-Clark et al. (2016) and references therein). We write the resonance condition as

[TABLE]

where $n=0,\pm 1,\pm 2,...$, $v_{x,r}$ is the parallel resonant velocity, and $\gamma=(1-v_{r}^{2}/c^{2})^{-1/2}$ with $v_{r}$ the total electron velocity. We first discuss the non-relativistic case ($\gamma=1$) for which the resonance is represented by a vertical line $v_{x}=v_{x,r}$ in the $v_{x}-v_{y}$ space. In the long-wavelength limit, $k^{2}d^{2}_{e}\ll 1$, the electric field of a single whistler is eliminated in a frame moving along $B_{0}$ at the speed $(v_{p,x}=\omega/k_{x})$. Energy conservation, $(v_{x}-v_{p,x})^{2}+v_{y}^{2}=\text{const}$, then requires that particle orbits lie on circles centered around $v_{x}=v_{p,x}$ as they oscillate in the fields of the whistler. When $kd_{e}\gtrsim 1$ the whistler retains a finite electric field in its frame and particle energy is not exactly conserved (Karimabadi et al., 1990, 1992). The nonlinear trapping width associated with a resonance $n$ can be calculated from the electron equation of motion in the whistler frame. Using the linearized cold plasma dispersion relation to obtain the whistler field components (neglecting the displacement current) we find the parallel trapping width of the $n$th resonance to be {widetext}

[TABLE]

where $J_{n}(k_{y}v_{y0}/\Omega_{e0})$ is the Bessel function of order n, $v_{y0}$ is the initial perpendicular velocity of the particle as it becomes resonant, and the whistler eigenvector amplitude $\tilde{B}$ is that of $\tilde{B_{y}}$. Figure 5 shows a schematic of the trapping widths for the $n=0,\pm 1$ resonances with $k_{x}d_{e}=0.6$, $k_{y}d_{e}=1$, which is the location of the peak in the spectrum obtained from the FFT of $B_{z}$ at $t\Omega_{e0}=177$ (not shown). The peak whistler amplitude of $\tilde{B}/B_{0}=0.125$ is chosen to evaluate the trapping widths. The trapping widths in $v_{x}$ are comparable to $V_{Ae0}$ for $v_{y0}/V_{Ae0}\gtrsim 1.5$.

Figure 5 establishes that the $n=-1$ resonance strongly scatters particles with $v_{x}>0$, producing the large horn-like feature near $v_{x}/V_{Ae0}\lesssim 2$ in Fig. 4b and transferring energy from the particles to the wave. The horn-like feature could also explain why the heat flux of the lowest-energy particles in Fig. 1 peaks for $t\Omega_{e0}\sim 100$ since initially higher-energy particles could be scattered into the lowest velocity bin. The $n=1$ resonance scatters some of the particles in the return current beam with $v_{x}/V_{Ae0}\sim-1$ to higher $v_{y}$ and smaller $v_{x}$ (note the small horn-like feature at this location in Fig. 4b). This interaction has a damping effect on the wave but the drive from $n=-1$ still dominates since the number of upscattered cold return current electrons is small. Initially the distribution function is flat near the whistler phase speed (Fig. 4) so any effect on the energy of the whistler by the $n=0$ resonance is small. As a result the fan instability associated with the $n=-1$ resonance is what drives the wave at early time.

The rightward-moving electron acoustic waves (EAWs) in Fig. 3b are parallel-propagating, so only the Landau resonance is significant. It produces a trapping width

[TABLE]

where $\tilde{E}$ is the EAW amplitude. Setting $k_{x}d_{e}=3.92$ and $\tilde{E}=0.3$ (the peak amplitude observed) we find that the EAW Landau resonance width (indicated by green vertical lines) overlaps with that of the whistler (orange curves) for a region of the phase space near $v_{y}/V_{Ae0}=2.4$ and nearly overlaps that of the $n=-1$ resonance (blue curves) near $v_{y}/V_{Ae0}\lesssim 3$. Such overlap should lead to irreversible diffusion in phase space, dragging particles towards $v_{x}=0$ as seen over time in Fig. 4(a-c). If the overlap is not quite reached between the two whistler resonances and that of the EAW, chaotic orbits will still set in as overlap is approached. A small spread in $k$ can also shift the centroids of the traps and bridge the gap between the resonances. Importantly, a particle’s energy can be reduced by a whistler and an EAW acting in tandem. The $n=-1$ hornlike feature in Fig. 4b has been shifted to the EAW phase speed $v_{p,x}/V_{Ae0}=1.17$ in Fig. 4c, pulling some particles into the whistler $n=0$ resonance at late time. This is a simple, graphical way to describe the nonlinear interaction of whistlers and EAWs (Agapitov et al., 2018; Vasko et al., 2018; Drake et al., 2015). The interaction between the whistler $n=-1$ resonance and the EAW resonance is analogous to the upscattering of high-energy particles and subsequent driving of plasma waves via a bump-on-tail instability as discussed in Haber et al. (1978).

As particles lose energy and diffuse towards $v_{x}=0$ they also lose momentum. To maintain zero net current the system generates an inductive field that slows the average drift speed of the return current beam. The leftward moving whistler grows by drawing energy from the return current beam via the Landau resonance at the phase speed $v_{p,x}/V_{Ae0}\sim-0.5$. The flattening of the distribution across $v_{x}$ in Fig. 4c is evidence that this is taking place. This leftward whistler does not grow at early time because of damping of this wave by particles at the $n=1$ resonance (Verscharen et al., 2019). At late time the electron distribution at this resonance has been flattened by the action of the rightward propagating whistler so growth of the negative propagating wave through the Landau resonance takes place. Note that the large number of particles in the return current means the dominant energy transfer to the leftward wave is through the Landau resonance, so it is a beam-driven mode.

5 Relativistic effects

Since there are very few particles at high energy for the initial distribution function chosen for the simulation, the high-energy component basically act as test particles. Their feedback on the wave dynamics is small. For these particles relativistic effects are important. In the relativistic case ($\gamma\gg 1$) the resonant surfaces and constant energy surfaces become elliptical in the $v_{x}-v_{\perp}$ plane (Omidi & Gurnett, 1982; Karimabadi et al., 1990). Trapping occurs around the intersection point of these two ellipses (Karimabadi et al., 1990), assuming that to lowest order in wave amplitude the particle energy is conserved in the wave frame.

Intersection of the resonance ellipse with the $v_{x}$ axis is given by

[TABLE]

where $\alpha_{n}=n\Omega_{e0}/(k_{x}c)$ and $v_{p,x}$ is calculated from Eq. (4). The constant-energy ellipses (Karimabadi et al., 1990) are given by

[TABLE]

with $R=c^{2}(\gamma_{0}^{2}-1)+v_{p,x}^{2}/(\gamma_{0}^{2}+v_{p,x}^{2}/c^{2})$, $v_{x,c}=v_{p,x}/(\gamma_{0}^{2}+v_{p,x}^{2}/c^{2})$ and $\gamma_{0}$ the initial Lorentz factor of a particle in the lab frame. For the whistlers in the simulation with $v_{p,x}\ll c$ the $v_{x}$ and $v_{y}$ axes of the ellipses are nearly equal. In the highly relativistic limit ($\gamma_{0}\rightarrow\infty$) the surface is a circle centered around $v=0$ with radius $c$ as expected.

Figure 4b displays the location of the resonance intersections for $n=1$ through $n=-5$ using (8) and nearby energy surfaces using (9) assuming $k_{x}d_{e}=0.6$ for the rightward moving wave. The energy surfaces, the horn-like structures, and the resonance intersections line up surprisingly well, justifying the use of the relativistic theory and further confirming that the waves are whistlers. The $n=-1$ through $n=5$ resonances are also shown for the leftward wave ($k_{x}d_{e}=-0.2$, $k_{y}d_{e}=0.6$) in Fig. 4c, along with two energy surfaces. Since the leftward wave is long-wavelength, only the Landau resonance acts on low-energy particles. Scattering by the normal resonances $n=1,2$ likely aids in the diffusion of high-energy particles although the bulk of the scattering has already taken place by the time the leftward wave has a large amplitude. The most significant relativistic effect is the location of the resonances in Fig. 4. Expressions for the relativistic trapping widths, which we do not invoke here, can be found in Karimabadi et al. (1990).

6 Discussion

We have found that energetic electrons escaping from a flare-like system with $\beta\sim 1$, which is expected from magnetic reconnection (Dahlin et al., 2017), efficiently drive whistler waves that scatter the escaping electrons. The whistlers pitch-angle scatter high-energy electrons on a rapid time scale of hundreds of cyclotron periods by means of cyclotron resonances, suppressing energy flux and increasing the perpendicular velocities of electrons. This is a local mechanism which can operate under the generic conditions of a reconnection-driven flare. Moderate reduction (up to a factor of two) of the field-aligned electron energy flux occurs. Since this scattering tends to increase perpendicular velocity, electrons will more effectively mirror when they encounter small-scale magnetic fluctuations in the corona or if they are accelerated towards the sun where the ambient magnetic field is stronger. Thus, this scattering mechanism will facilitate the confinement of energetic electrons in energy release sites in flares which is required for electrons to reach the relativistic velocities seen in observations (Lin et al., 2003; Krucker et al., 2010; Gary et al., 2018).

Although our model was designed to study flares, we point out some possible implications for transport into the outer corona and ultimately the solar wind. The fluctuations in our simulation damp out after scattering is complete. Nevertheless, the heat flux associated with the initial electron distribution is permanently reduced. We suggest on the basis of the present simulations that oblique waves could have grown to large amplitude either in the corona or in the solar wind in the outer reaches of the corona, scattered the electrons to reduce the heat flux and then died away, leaving a remnant reduced heat flux. We thus propose that the electron energy flux produced in the corona as a result of reconnection or other mechanisms could be suppressed by oblique whistlers, leading to marginally stable electron distributions that then propagate outward, leaving no trace of the self-generated turbulence that limited the heat flux. It has been shown that distribution functions can stream ballistically over several mean free paths before collisions dominate and the distribution becomes Maxwellian-like (Malkov, 2017). Since the mean free path of the solar wind is roughly $1\text{AU}$, the streaming of heat-flux-carrying distributions out of the corona is thus feasible.

There is mounting evidence that the electron heat flux in the solar wind is limited by whistlers in regions with large plasma $\beta$ (Tong et al., 2018). However, the large-amplitude, oblique whistlers that would limit the heat flux have not been measured in the solar wind at $1$ AU (Wilson et al., 2013; Lacombe et al., 2014; Stansby et al., 2016; Tong et al., 2019a, b; Kuzichev et al., 2019), although some observations have implied effective pitch-angle scattering by whistlers with small angles relative to $B_{0}$ (Kajdič et al., 2016). Oblique whistler scattering of the solar wind strahl into the halo is currently under active investigation (e.g. Verscharen et al. (2019), Boldyrev & Horaites (2019), Vasko et al. (2019)).

Our simulation addresses the $\beta\sim 1$ limit of flares. However, for a lower $\beta$ simulation, whistler growth would be suppressed since less free energy would be available (Roberg-Clark et al., 2018b). Scattering would still be expected to occur and the mode would retain a finite angle to $B_{0}$. Since the electron Alfvén speed would increase relative to the thermal speed ($\beta=(V_{T}/V_{Ae})^{1/2}$) the whistler phase speed might lie farther out in the tail of the distribution, leading to scattering of the highest energy particles.

We caution that our simulation suffers from the usual constraints of particle-in-cell simulations such as small simulated domains and short time scales compared to relevant scales in the corona. The model heat flux distribution function we use also contains a very sharp gradient near $v_{x}=0$. However, sharp gradients are likely to form during flares, during which particles are rapidly accelerated and our intent is to explore the possible mechanisms at play that would limit electron escape. Our model is relevant to $\beta\sim 1$ weakly collisional plasmas and could also be applied to transport in astrophysical coronae and low-luminosity accretion flows.

The authors acknowledge support from NSF Grant No. PHY1805829, from NASA grant NNX17AG27G and from the FIELDS team of the Parker Solar Probe (NASA Contract No. NNN06AA01C). GTRC was funded by an Ann G. Wylie Dissertation Fellowship from the University of Maryland, College Park. OA and JFD were supported by NASA grant 80NNSC19K0848. OA was partially supported by NSF grant number 1914670. JFD acknowledges partial support from NSF Grant No. PHY1748958 at the Kavli Institute for Theoretical Physics at UCSB. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No.DE-AC02-05CH11231. Simulation data is available upon request.