This paper presents a numerical study of the linear and nonlinear evolution of the electromagnetic electron-cyclotron (EMEC) instability in a bi-Kappa distributed plasma. Distributions with high energy tails described by the Kappa power-laws are often observed in collision-less plasmas (e.g., solar wind and accelerators), where wave-particle interactions control the plasma thermodynamics and keep the particle distributions out of Maxwellian equilibrium. Under certain conditions, the anisotropic bi-Kappa distribution gives rise to plasma instabilities creating low-frequency EMEC waves in the whistler branch. The instability saturates nonlinearly by reducing the temperature anisotropy until marginal stability is reached. Numerical simulations of the Vlasov-Maxwell system of equations show excellent agreement with the growth-rate and real frequency of the unstable modes predicted by linear theory. The wave-amplitude of the EMEC waves at nonlinear saturation is consistent with magnetic trapping of the electrons.

## I. INTRODUCTION

Direct *in-situ* measurements in space plasmas reveal particle distributions far from thermal (Maxwellian) equilibrium by the presence of suprathermal populations that enhance the high-energy tails of the distribution functions, as well as kinetic anisotropies manifested as a temperature anisotropy $T\u22a5/T\u2225\u22601$ of the electron and proton velocity distributions (where the perpendicular and parallel directions are with respect to the direction of the ambient magnetic field), and sometimes as beam-like features.^{1} Particle-particle collisions are indeed rare in space plasmas and cannot establish thermal equilibrium, but wave-particle interactions can resonantly accelerate part of the particle population and explain the deviations from Maxwellian equilibrium. However, these anisotropies are not as high as expected from exospheric models, and hence there seems to be a limiting factor on the temperature anisotropy.^{2,3} In the absence of collisions, linear theory predicts that temperature anisotropy driven instabilities could lead to a partial thermalization of the particles until marginal stability is reached. For electrons, an excess of temperature in the parallel direction $T\u2225>T\u22a5$ can be limited by the firehose instability, while an excess temperature in the perpendicular direction $T\u22a5>T\u2225$ is limited by the whistler instability, also known as the electromagnetic electron cyclotron (EMEC) instability.^{2} For a sufficiently small electron plasma beta, i.e., $\beta <0.1$, the EMEC instability is competed by electrostatic instabilities propagating obliquely to the magnetic field,^{4,5} or by the slow extraordinary mode (Z mode) that becomes unstable when the electron gyro-frequency is larger than the electron plasma frequency.^{4}

Whistler waves propagating along the magnetic field lines are right-hand circularly polarized electromagnetic waves with frequencies between the ion and electron cyclotron frequencies and were first observed in ionospheric physics.^{6} Low-frequency whistler waves triggered by lightnings in the Earth's equatorial region have frequencies in the kilohertz range and propagate between the hemispheres with a descending whistler-like tone. On the other hand, whistler waves triggered by electron anisotropies and suprathermal populations of electrons are expected to play an important role at dissipative scales in the solar wind^{7} and planetary magnetospheres.^{8,9} Whistler instabilities have also been observed in laboratory experiments.^{10}

In the first investigations of the EMEC instability, the anisotropic electrons were described by idealized bi-Maxwellian models (see the textbook of Gary^{11} and references therein), while more recently, the EMEC instability has been studied using generalized power-laws, such as bi-Kappa or product-bi-Kappa models, which incorporate suprathermal electrons in the distribution functions.^{12–17} A variety of plasma parameters have also been invoked, ranging from the weakly magnetized plasma of the solar wind, via the moderately magnetized plasmas of the coronal outflows, flares, and the plasma sheet, to the highly magnetized plasmas in the tail lobe and the radiation belts in the magnetosphere. Compared to the case with bi-Maxwellian electrons (recovered in the limit of a very large power-index $\kappa \u2192\u221e$), the influence of suprathermal populations on the EMEC instability was found to be highly dependent on the shape of the velocity distribution as well as the temperature anisotropy. For bi-Kappa distributed electrons with perpendicular temperature much higher than parallel temperature, the growth-rate tends to increase with increasing the kappa-index,^{13,14} while the opposite holds at a low temperature anisotropy.^{12} On the other hand, the lowest thresholds relevant for the marginal instability condition decrease with increasing the density of suprathermal populations.^{15}

Large amplitude whistlers can undergo dual cascades to both larger and smaller wavenumbers involving parametric processes.^{18} A spectrum of waves can change the background electron distribution function via quasilinear and other nonlinear processes.^{19,20} Introduced first for an arbitrary distribution of plasma particles,^{22,23} the quasilinear approach has been specialized and tested for bi-Maxwellian distributed plasmas,^{21,24,25} showing that the increase of the EMEC magnetic fields tends to decrease the temperature anisotropy to a point of marginal stability. Similar results were found in particle-in-cell simulations using a Maxwellian dense core and highly temperature anisotropic ($T\u22a5/T\u2225=25$) bi-Kappa distributed tails.^{26} Using a set of conserved quantities for whistler waves propagating parallel to the magnetic field,^{23} it is possible to show that an increase of the wave magnetic field energy leads to an increase of the total parallel kinetic energy and a decrease of the perpendicular kinetic energy of the electrons, while conserving the total energy of the system.

The aim of this paper is to numerically investigate linear growth and nonlinear saturation of the EMEC instability in a bi-Kappa distributed plasma, and to compare the results with the existing theoretical predictions and simulations. In Sec. II, we provide a brief mathematical formulation for the bi-Kappa model invoked and the resulting dispersion relations of the EMEC modes. We restrict to EMEC modes propagating parallel ($k=k\u2225$) to the magnetic field, because in parallel direction these modes exhibit maximum growth rates and decouple from the electrostatic modes.^{27} The simulation setup and the principal results are presented to a large detail in Sec. III. The simulations are based on a previously developed Vlasov-Maxwell solver that uses a Fourier transform technique in velocity space to solve the Vlasov equation.^{31,32} The numerical code, which showed excellent agreement with theory for the whistler instability with a bi-Maxwellian electron distribution,^{31,32} is here extended to study the more general case of a bi-Kappa distribution function for the electrons. Similarities and differences between the simulation results and previous results using the bi-Maxwellian electron distributions are discussed in detail. The main conclusions of our study are iterated in Sec. IV.

## II. MATHEMATICAL FORMULATION

In preparation for the numerical simulations, we assume a geometry of the problem in which the background magnetic field is directed along the *z*-axis, $B0=z\u0302B0$, where *B*_{0} is the magnitude of the magnetic field and $z\u0302$ is a unit vector along the *z*-axis. The fluctuating electric and magnetic fields, $E$ and $B$, respectively, may have components in all three dimensions.

To model an anisotropic electron distribution with high energy tails, we use the bi-Kappa distribution function

where *κ* is the power index, the perpendicular and parallel (to the magnetic field) velocity components are denoted $v\u22a5=vx2+vy2$ and $v\u2225=vz$, respectively, and the parallel and perpendicular thermal speeds $\theta \u2225$ and $\theta \u22a5$ are defined in terms of the respective kinetic temperatures $T\u2225$ and $T\u22a5$ via the relations (see, e.g., Ref. 15)

and

The linear dispersion relation for the instability of EMEC-whistler waves propagating along the background magnetic field for a bi-Kappa distribution (omitting subscript $\u2225$ for the wavenumber *k*) reads (see, e.g., Ref. 15)

in terms of the dispersion function for a bi-Kappa distributed plasma^{28}

In the limit $\kappa =\u221e$, the bi-Kappa distribution converges to the bi-Maxwellian distribution function, and the dispersion relation can be written^{29}

where

is the plasma dispersion function^{30} and $erf(z)$ is the error function of complex arguments.

## III. SIMULATION SETUP AND NUMERICAL RESULTS

Fully electromagnetic and kinetic simulations are carried out by solving numerically the Maxwell equations for the electromagnetic field and the Vlasov equation for the electron dynamics, using a Fourier method in velocity space,^{31,32} while the ions are assumed to be immobile. The initial condition for the electron distribution function is obtained by using the Fourier transform pair in velocity space

which yields the Fourier transformed bi-Kappa distribution function

where $n=\kappa \theta \u22252\eta \u22252+\kappa \theta \u22a52\eta \u22a52$, $\eta $ is the Fourier transformed velocity variable, $\Gamma (\xi )$ is the Gamma function of argument *ξ*, and $Kv(n)$ is the modified Bessel function of second kind of non-integer order *v* and argument *n*.^{33} The normalization $\u222bf\u2009d3v=1$ corresponds to $f\u0302\eta =0=1/(2\pi )3$ in the Fourier transformed velocity space.^{31,32} More details of the simulation set-up are given in the Appendix. It has been confirmed that the simulations conserve the total energy to a high degree in the nonlinear evolution of the system: While the magnetic wave energy increases to about 3%–5% of the total (kinetic plus electromagnetic) energy and the electron kinetic energy decreases with approximately the same amount, the total energy typically decreases by about 0.1%. The electric energy is only about 0.01% as expected for low-frequency whistlers with $\omega ce/\omega pe\u226a1$.

We have carried out three sets of simulations, (i) where the initial electron number density and magnetic field are chosen such that $\omega pe/\omega ce=12.65$, using a relatively large temperature anisotropy $T\u22a5/T\u2225=8$, (ii) for a larger density relative to the magnetic field such that $\omega pe/\omega ce=40$ and $T\u22a5/T\u2225=8$, and finally (iii) $\omega pe/\omega ce=40$ and a lower temperature anisotropy $T\u22a5/T\u2225=4$. In all cases, we use an electron temperature such that $kBT\u2225/mec2=1/2500$, corresponding to $T\u2225=2.37\u2009MK$. These plasma parameters are typical for the conditions in the solar corona, solar wind, and solar flares; see, e.g., Refs. 14, 34, and 35. The parallel electron plasma beta

takes the value 0.128 in case (i) and 1.28 in cases (ii) and (iii). Theoretical and numerical investigations using bi-Maxwellian electron velocity distributions^{5} have shown that the instability has the largest growth-rate in the parallel direction for $\beta \u2225>0.025$, while for $\beta \u2225<0.025$ the instability has larger growth-rates at large angles giving rise to waves with large electrostatic components. For the selected parameters, with $\beta \u2225>0.1$, it is therefore expected that the instability has a maximum growth rate in the parallel direction, consistent with the simulation geometry.

Figures 1–3 show the spatiotemporal evolution of the wave magnetic field for $\kappa =\u221e$ (i.e., a bi-Maxwellian distribution), *κ* = 8, *κ* = 4, and *κ* = 2 for the three cases (i)–(iii). The unstable EMEC wave modes initially grow exponentially with time and saturate nonlinearly after the initial small-amplitude perturbations have increased about 8 orders of magnitude in amplitude. The resulting large amplitude right-hand polarized whistler-like EMEC waves are visible in panels (a)–(d). As seen in panels (e) and (f), the most unstable mode grows faster and the saturation amplitude is larger for larger values of *κ*. The amplitude of the wave field typically reaches a maximum as the linear growth saturates, and for the cases of higher plasma beta shown in Figs. 2 and 3 [panels (a)–(d)], the wave field at later times slowly decreases in amplitude as the wave spectrum cascades to longer wavelengths. A similar effect was observed in Ref. 25 in the form of a saturation and switching off of the growth of shorter wavelength waves and a shift of the unstable wave spectrum to longer wavelengths. The case of smaller temperature anisotropy shown in Fig. 3 gives lower growth-rates and saturation amplitudes compared to the higher anisotropy case in Fig. 2.

Figures 4–6 show the growth-rates and real frequencies obtained in the Vlasov simulations, and a comparison with the theoretical values obtained from the wave dispersion relations (4) and (6). The numerical values of the real frequency were obtained by studying the periodicity of the exponentially growing wave mode in time for each value of the wavenumber, while the growth-rates were obtained by linear fits (linear regression) in time to the logarithm of the amplitude of each mode during its linear growth phase. As seen in Figs. 4–6, the growth-rates and real frequencies obtained in the Vlasov simulations agree very well with the theoretical dispersion curves, which give confidence in the correctness of both the theoretical and numerical results. The modes with the highest growth-rates in Figs. 4–6 are the dominating mode first visible and saturating nonlinearly in Figs. 1–3. Also in agreement with Figs. 1–3, the growth-rates of the fastest growing modes are larger for larger values of *κ*, while the wavenumber and hence the wavelength of the fastest growing modes are only weakly dependent on *κ*. The ambient magnetic field has in general a stabilizing effect on the EMEC instability, and the spectrum of unstable wave modes shrink with increasing magnetic field. This effect can be observed in the low-beta case in Fig. 4, where the smallest wavenumbers (e.g., $kc/\omega pe\u22720.5$) are stabilized compared to the high-beta cases in Figs. 5 and 6. On the other hand, for a smaller spectral index *κ* = 2, the smaller wavenumber modes also become unstable, as seen in Fig. 4(g). As expected, the case of smaller temperature anisotropy in Fig. 7 also gives smaller growth-rates and a narrower range of unstable modes of the EMEC instability compared to Fig. 6. This is consistent with the condition $kc/\omega pe\u2264(T\u22a5/T\u2225\u22121)1/2$ for a bi-Maxwellian plasma,^{22,25} which also leads to a stabilization of the higher wavenumber modes as the temperature anisotropy decreases. For all cases, the most unstable wave modes have real frequencies somewhat below the electron cyclotron frequency, $\omega R\u22480.70$–$0.75\u2009\omega ce$. By comparing the results in Figs. 1–3 with Figs. 4–6, it is interesting to notice that the maximum amplitude of the wave magnetic field as the instability saturates nonlinearly increases with the maximum linear growth-rate. Similar results were obtained in Refs. 22 and 24 and were explained by magnetic trapping of electrons.

Figures 7–9 show the time evolution of the parallel electron distribution function, defined as the electron distribution function integrated over the perpendicular velocity space, $f\u2225=\u222bfd2v\u22a5$. At *t* = 0 (the solid lines in Figs. 7–9), we have from Eq. (1)

In general, the electron distribution function shows a broadening in the parallel direction until marginal stability is reached. The broadening in the parallel direction is indicative of wave-particle interactions via quasi-linear processes. For parallel propagation of the whistler waves, two conserved quantities^{23,25} ensure that an increase of magnetic wave energy leads to the simultaneous increase of parallel kinetic energy and decrease of perpendicular kinetic energy. Hence, the growth of whistler waves due to the EMEC instability leads to a broadening of the electron velocity distribution in the parallel direction (as seen in Figs. 7–9) and narrowing in the perpendicular direction until marginal stability is achieved. In a similar way as in Ref. 25 for a bi-Maxwellian plasma, when the anisotropy decreases, the electro distribution does not change the shape, remaining in this case approximately Kappa distributed. Thermal spread increases in parallel direction but, apparently, the high-energy tails are diminished (increasing of the power index *κ*). In general, the present results are also in line with more recent simulations^{26} using a dense core (95% density) Maxwellian distribution and a low-density (5% density) bi-Kappa distributed halo with $T\u22a5/T\u2225=25$, but in the latter case a decrease of *κ* was instead observed.

## IV. CONCLUSIONS

We have studied the linear and nonlinear evolution of the EMEC instability for a collisionless anisotropic plasma by means of numerical simulations of the Vlasov-Maxwell system, using a bi-Kappa distribution for the electrons to model high-energy tails in their distribution. The numerically obtained growth-rates and real frequencies show excellent agreement with theoretical models. For a relatively large temperature anisotropy $T\u22a5/T\u2225$, the instability has largest growth-rate for large spectral indices ($\kappa \u2192\u221e$). The EMEC instability saturates nonlinearly by exciting large-amplitude right-hand circularly polarized electromagnetic waves with frequencies near the electron cyclotron frequency. The EMEC instability saturates by magnetic trapping of electrons,^{22,24,25} leading to an increase of the saturated amplitude for higher linear growth-rates of the instability.

Values chosen for the plasma parameters in our numerical setup are relevant for the hot plasma in the solar corona and flares, where the electron temperature typically is a few million Kelvin and the electron plasma frequency is about ten times higher than the electron cyclotron frequency. But extensions of the present results to the outer corona conditions in the solar wind and terrestrial magnetosphere are straightforward because the plasma beta in these environments remains in the same range (e.g., $0.1<\beta e<10$ (Ref. 3)) as relevant in our simulations. Thus, the results obtained in the present paper provide major premises to approach more complex models of the solar wind particles, which combine a bi-Maxwellian core population at low energies and a bi-Kappa distribution modeling the suprathermal component.^{16,17} The anisotropy of each component contributes to different peaks in the unstable spectrum, and the nonlinear evolution of these cases of high complexity may now be addressed in future works. Moreover, the numerical technique presented here may easily be adapted to enable further studies of other instabilities driven by kinetic anisotropies in Kappa distributed plasmas, for example, oblique generation of whistlers that would include quasi-electrostatic modes with a component of the electric field along the wave vector that could lead to efficient acceleration of the tail electrons to high energies.

## ACKNOWLEDGMENTS

M.L. acknowledges support from the Katholieke Universiteit Leuven, Grant No. SF/12/003, and from the Ruhr-Universität Bochum. These results were obtained in the framework of the Project Nos. GOA/2015-014 (KU Leuven) and C 90347 (ESA Prodex 9). The research leading to these results has also received funding from the European Commission's Seventh Framework Programme (FP7/2007-2013) under the Grant Agreements SOLSPANET (Project No. 269299, www.solspanet.eu) and eHeroes (Project No. 284461, www.eheroes.eu). B.E. acknowledges support from the Engineering and Physical Sciences Research Council (EPSRC), UK, Grant No. EP/M009386/1. The simulation data associated with this paper are available at: http://dx.doi.org//10.15129/cd90b629-a2b3-4c75-8c22-58abe72ad03c.

### APPENDIX: SIMULATION PARAMETERS AND SETUP

The simulation geometry is a one-dimensional simulation box in the *z*-direction, parallel to the ambient magnetic field. The simulation box is of size $Lx=(2\pi /0.15)\pi c/\omega pe\u224841.9c/\omega pe$ and is resolved with 80 intervals with periodic boundary conditions. A pseudospectral method to calculate spatial derivatives accurately, and numerical dissipation is used at high wavenumbers to reduce aliasing effects. The three-dimensional velocity space is resolved by $76\xd776\xd776$ intervals, using a Fourier method^{31,32} and high-order difference schemes to calculate derivatives in the Fourier transformed velocity space. To accurately resolve the evolution in the Fourier transformed velocity space, we use a maximum $\eta z,max=6\u2009vT\u2225\u22121$ for $\kappa =\u221e$, $\eta z,max=7\u2009vT\u2225\u22121$ for *κ* = 8, $\eta z,max=8\u2009vT\u2225\u22121$ for *κ* = 4, and $\eta z,max=10\u2009vT\u2225\u22121$ for *κ* = 2, and half these values, $\eta x,max=\eta y,max=\eta z,max/2$, in the perpendicular dimensions, where we denoted $vT\u2225=kBT\u2225/me$. The maximum represented velocity is given by $vmax=\pi /\Delta \eta $ where $\Delta \eta $ is the grid size in the Fourier transformed velocity space, and the represented grid spacing in velocity space is given by $\Delta v=\pi /\eta max$. Therefore, the maximum represented parallel velocities are $vz,max\u224820\u2009vT\u2225$ for $\kappa =\u221e,\u2009vz,max\u224817\u2009vT\u2225$ for *κ* = 8, $vz,max\u224815\u2009vT\u2225$ for *κ* = 4, and $vz,max\u224812\u2009vT\u2225$ for *κ* = 2, and twice as large represented velocities in the perpendicular *v _{x}* and

*v*velocity dimensions where the anisotropic velocity distributions are wider. To seed the instability, a quasi-random wave magnetic field is introduced in the initial conditions with an amplitude about 8 orders of magnitude smaller than the ambient magnetic field. Identical initial conditions for the wave magnetic field are used in all simulations. A 4th-order Runge-Kutta scheme is used to advance the solution in time, with adaptive time-stepping to maintain numerical stability in the nonlinear evolution of the system. The number of time-steps used is 16 000.

_{y}