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.

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/T1 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>T can be limited by the firehose instability, while an excess temperature in the perpendicular direction T>T 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., β<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 wind7 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 Gary11 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 κ), 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/T=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) 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.

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=ẑB0, where B0 is the magnitude of the magnetic field and ẑ 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

f(v,v)=1π3/2θ2θΓ(κ+1)κ3/2Γ(κ1/2)(1+v2κθ2+v2κθ2)κ1,
(1)

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

T=meθ22kBκκ3/2,
(2)

and

T=meθ22kBκκ3/2.
(3)

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

k2c2ωp2=TT1+(ωωce)T/T+ωcek2kBT/me13/(2κ)Zκ(ωωcek2kBT/me13/(2κ)),
(4)

in terms of the dispersion function for a bi-Kappa distributed plasma28 

Zκ(ξ)=1π1/2κ1/2Γ(κ)Γ(κ1/2)dx(1+x2/κ)κxξ,(ξ)>0.
(5)

In the limit κ=, the bi-Kappa distribution converges to the bi-Maxwellian distribution function, and the dispersion relation can be written29 

k2c2ωpe2=TT1+(ωωce)T/T+ωcek2kBT/meZ(ωωcek2kBT/me),
(6)

where

Z(ξ)=iπexp(ξ2)[1+erf(iξ)],
(7)

is the plasma dispersion function30 and erf(z) is the error function of complex arguments.

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

f(x,v,t)=f̂(x,η,t)eiη·vd3η,
(8)
f̂(x,η,t)=1(2π)3f(x,v,t)eiη·vd3v,
(9)

which yields the Fourier transformed bi-Kappa distribution function

f̂(η,η)=2(2π)3(n2)κ1/2Kκ1/2(n)Γ(κ1/2),
(10)

where n=κθ2η2+κθ2η2, η is the Fourier transformed velocity variable, Γ(ξ) 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 fd3v=1 corresponds to f̂η=0=1/(2π)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 ωce/ωpe1.

We have carried out three sets of simulations, (i) where the initial electron number density and magnetic field are chosen such that ωpe/ωce=12.65, using a relatively large temperature anisotropy T/T=8, (ii) for a larger density relative to the magnetic field such that ωpe/ωce=40 and T/T=8, and finally (iii) ωpe/ωce=40 and a lower temperature anisotropy T/T=4. In all cases, we use an electron temperature such that kBT/mec2=1/2500, corresponding to T=2.37MK. 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

β=2n0kBTμ0B02=2kBTmec2ωpe2ωce2,
(11)

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 distributions5 have shown that the instability has the largest growth-rate in the parallel direction for β>0.025, while for β<0.025 the instability has larger growth-rates at large angles giving rise to waves with large electrostatic components. For the selected parameters, with β>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 κ= (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.

FIG. 1.

The spatiotemporal evolution of the wave magnetic field, showing the nonlinear saturation of the EMEC instability by large amplitude waves for ωpe/ωce=12.65,T/T=8, and (a) κ=, (b) κ = 8, (c) κ = 4, and (d) κ = 2. (e) The time evolutions of the root mean square magnetic field, and the same data are plotted in (f) with logarithmic scale on the vertical axis.

FIG. 1.

The spatiotemporal evolution of the wave magnetic field, showing the nonlinear saturation of the EMEC instability by large amplitude waves for ωpe/ωce=12.65,T/T=8, and (a) κ=, (b) κ = 8, (c) κ = 4, and (d) κ = 2. (e) The time evolutions of the root mean square magnetic field, and the same data are plotted in (f) with logarithmic scale on the vertical axis.

Close modal
FIG. 2.

The spatiotemporal evolution of the wave magnetic field, showing the nonlinear saturation of the EMEC instability by large amplitude waves for ωpe/ωce=40,T/T=8, and (a) κ=, (b) κ = 8, (c) κ = 4, and (d) κ = 2. (e) The time evolutions of the root mean square magnetic field, and the same data are plotted in (f) with logarithmic scale on the vertical axis.

FIG. 2.

The spatiotemporal evolution of the wave magnetic field, showing the nonlinear saturation of the EMEC instability by large amplitude waves for ωpe/ωce=40,T/T=8, and (a) κ=, (b) κ = 8, (c) κ = 4, and (d) κ = 2. (e) The time evolutions of the root mean square magnetic field, and the same data are plotted in (f) with logarithmic scale on the vertical axis.

Close modal
FIG. 3.

The spatiotemporal evolution of the wave magnetic field, showing the nonlinear saturation of the EMEC instability by large amplitude waves for ωpe/ωce=40,T/T=4, and (a) κ=, (b) κ = 8, (c) κ = 4, and (d) κ = 2. (e) The time evolutions of the root mean square magnetic field, and the same data are plotted in (f) with logarithmic scale on the vertical axis.

FIG. 3.

The spatiotemporal evolution of the wave magnetic field, showing the nonlinear saturation of the EMEC instability by large amplitude waves for ωpe/ωce=40,T/T=4, and (a) κ=, (b) κ = 8, (c) κ = 4, and (d) κ = 2. (e) The time evolutions of the root mean square magnetic field, and the same data are plotted in (f) with logarithmic scale on the vertical axis.

Close modal

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/ωpe0.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/ωpe(T/T1)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, ωR0.700.75ω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.

FIG. 4.

Growth-rates (top panels) and corresponding real frequencies (bottom panels) of the EMEC instability obtained in the Vlasov simulation (“+” markers) and the theoretical growth values (solid lines) for ωpe/ωce=12.65,T/T=8,κ= [(a) and (b)], κ = 8 [(c) and (d)], κ = 4 [(e) and (f)], and κ = 2 [(g) and (h)]. The maximum growth-rate is ωI/ωce=0.32 at kc/ωpe=1.2 for κ=, ωI/ωce=0.31 at kc/ωpe=1.2 for κ = 8, ωI/ωce=0.29 at kc/ωpe=1.2 for κ = 4, and ωI/ωce=0.22 at kc/ωpe=1.2 for κ = 2.

FIG. 4.

Growth-rates (top panels) and corresponding real frequencies (bottom panels) of the EMEC instability obtained in the Vlasov simulation (“+” markers) and the theoretical growth values (solid lines) for ωpe/ωce=12.65,T/T=8,κ= [(a) and (b)], κ = 8 [(c) and (d)], κ = 4 [(e) and (f)], and κ = 2 [(g) and (h)]. The maximum growth-rate is ωI/ωce=0.32 at kc/ωpe=1.2 for κ=, ωI/ωce=0.31 at kc/ωpe=1.2 for κ = 8, ωI/ωce=0.29 at kc/ωpe=1.2 for κ = 4, and ωI/ωce=0.22 at kc/ωpe=1.2 for κ = 2.

Close modal
FIG. 5.

Growth-rates (top panels) and corresponding real frequencies (bottom panels) of the EMEC instability obtained in the Vlasov simulation (“+” markers) and the theoretical growth values (solid lines) for ωpe/ωce=40,T/T=8,κ= [(a) and (b)], κ = 8 [(c) and (d)], κ = 4 [(e) and (f)], and κ = 2 [(g) and (h)]. The maximum growth-rate is ωI/ωce=1.09 at kc/ωpe=1.1 for κ=, ωI/ωce=1.06 at kc/ωpe=1.1 for κ = 8, ωI/ωce=1.01 at kc/ωpe=1.1 for κ = 4, and ωI/ωce=0.79 at kc/ωpe=1.0 for κ = 2.

FIG. 5.

Growth-rates (top panels) and corresponding real frequencies (bottom panels) of the EMEC instability obtained in the Vlasov simulation (“+” markers) and the theoretical growth values (solid lines) for ωpe/ωce=40,T/T=8,κ= [(a) and (b)], κ = 8 [(c) and (d)], κ = 4 [(e) and (f)], and κ = 2 [(g) and (h)]. The maximum growth-rate is ωI/ωce=1.09 at kc/ωpe=1.1 for κ=, ωI/ωce=1.06 at kc/ωpe=1.1 for κ = 8, ωI/ωce=1.01 at kc/ωpe=1.1 for κ = 4, and ωI/ωce=0.79 at kc/ωpe=1.0 for κ = 2.

Close modal
FIG. 6.

Growth-rates (top panels) and corresponding real frequencies (bottom panels) of the EMEC instability obtained in the Vlasov simulation (“+” markers) and the theoretical growth values (solid lines) for ωpe/ωce=40,T/T=4,κ= [(a) and (b)], κ = 8 [(c) and (d)], κ = 4 [(e) and (f)], and κ = 2 [(g) and (h)]. The maximum growth-rate is ωI/ωce=0.5 at kc/ωpe=0.85 for κ=, ωI/ωce=0.48 at kc/ωpe=0.84 for κ = 8, ωI/ωce=0.45 at kc/ωpe=0.80 for κ = 4, and ωI/ωce=0.33 at kc/ωpe=0.80 for κ = 2.

FIG. 6.

Growth-rates (top panels) and corresponding real frequencies (bottom panels) of the EMEC instability obtained in the Vlasov simulation (“+” markers) and the theoretical growth values (solid lines) for ωpe/ωce=40,T/T=4,κ= [(a) and (b)], κ = 8 [(c) and (d)], κ = 4 [(e) and (f)], and κ = 2 [(g) and (h)]. The maximum growth-rate is ωI/ωce=0.5 at kc/ωpe=0.85 for κ=, ωI/ωce=0.48 at kc/ωpe=0.84 for κ = 8, ωI/ωce=0.45 at kc/ωpe=0.80 for κ = 4, and ωI/ωce=0.33 at kc/ωpe=0.80 for κ = 2.

Close modal
FIG. 7.

The electron distribution function (parallel component) for ωpe/ωce=12.65, T/T=8 at x = 0 and t = 0 (solid line), and (a) κ= at t=65ωce1 (dashed line), t=118ωce1 (dashed-dotted line), and t=172ωce1 (dotted line); (b) for κ = 8 at t=65ωce1 (dashed line), t=121ωce1 (dashed-dotted line), and t=178ωce1 (dotted line); (c) for κ = 4 at t=68ωce1 (dashed line), t=124ωce1 (dashed-dotted line), and t=180ωce1 (dotted line); and (d) for κ = 2 at t=122ωce1 (dashed line), t=163ωce1 (dashed-dotted line), and t=203ωce1 (dotted line).

FIG. 7.

The electron distribution function (parallel component) for ωpe/ωce=12.65, T/T=8 at x = 0 and t = 0 (solid line), and (a) κ= at t=65ωce1 (dashed line), t=118ωce1 (dashed-dotted line), and t=172ωce1 (dotted line); (b) for κ = 8 at t=65ωce1 (dashed line), t=121ωce1 (dashed-dotted line), and t=178ωce1 (dotted line); (c) for κ = 4 at t=68ωce1 (dashed line), t=124ωce1 (dashed-dotted line), and t=180ωce1 (dotted line); and (d) for κ = 2 at t=122ωce1 (dashed line), t=163ωce1 (dashed-dotted line), and t=203ωce1 (dotted line).

Close modal

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=fd2v. At t = 0 (the solid lines in Figs. 7–9), we have from Eq. (1)

f(v)=1π1/2θΓ(κ+1)κ3/2Γ(κ1/2)(1+v2κθ2)κ.
(12)

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 quantities23,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 simulations26 using a dense core (95% density) Maxwellian distribution and a low-density (5% density) bi-Kappa distributed halo with T/T=25, but in the latter case a decrease of κ was instead observed.

FIG. 8.

The electron distribution function (parallel component) for ωpe/ωce=40, T/T=8 at x = 0 and t = 0 (solid line), and (a) κ= at t=28ωce1 (dashed line), t=58ωce1 (dashed-dotted line), and t=96ωce1 (dotted line); (b) for κ = 8 at t=29ωce1 (dashed line), t=59ωce1 (dashed-dotted line), and t=96ωce1 (dotted line); (c) for κ = 4 at t=30ωce1 (dashed line), t=62ωce1 (dashed-dotted line), and t=103ωce1 (dotted line); and (d) for κ = 2 at t=46ωce1 (dashed line), t=70ωce1 (dashed-dotted line), and t=111ωce1 (dotted line).

FIG. 8.

The electron distribution function (parallel component) for ωpe/ωce=40, T/T=8 at x = 0 and t = 0 (solid line), and (a) κ= at t=28ωce1 (dashed line), t=58ωce1 (dashed-dotted line), and t=96ωce1 (dotted line); (b) for κ = 8 at t=29ωce1 (dashed line), t=59ωce1 (dashed-dotted line), and t=96ωce1 (dotted line); (c) for κ = 4 at t=30ωce1 (dashed line), t=62ωce1 (dashed-dotted line), and t=103ωce1 (dotted line); and (d) for κ = 2 at t=46ωce1 (dashed line), t=70ωce1 (dashed-dotted line), and t=111ωce1 (dotted line).

Close modal
FIG. 9.

The electron distribution function (parallel component) for ωpe/ωce=40, T/T=4 at x = 0 and t = 0 (solid line), and (a) κ= at t=41ωce1 (dashed line), t=77ωce1 (dashed-dotted line), and t=119ωce1 (dotted line); (b) for κ = 8 at t=41ωce1 (dashed line), t=80ωce1 (dashed-dotted line), and t=124ωce1 (dotted line); (c) for κ = 4 at t=42ωce1 (dashed line), t=83ωce1 (dashed-dotted line), and t=128ωce1 (dotted line); and (d) for κ = 2 at t=62ωce1 (dashed line), t=92ωce1 (dashed-dotted line), and t=139ωce1 (dotted line).

FIG. 9.

The electron distribution function (parallel component) for ωpe/ωce=40, T/T=4 at x = 0 and t = 0 (solid line), and (a) κ= at t=41ωce1 (dashed line), t=77ωce1 (dashed-dotted line), and t=119ωce1 (dotted line); (b) for κ = 8 at t=41ωce1 (dashed line), t=80ωce1 (dashed-dotted line), and t=124ωce1 (dotted line); (c) for κ = 4 at t=42ωce1 (dashed line), t=83ωce1 (dashed-dotted line), and t=128ωce1 (dotted line); and (d) for κ = 2 at t=62ωce1 (dashed line), t=92ωce1 (dashed-dotted line), and t=139ωce1 (dotted line).

Close modal

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/T, the instability has largest growth-rate for large spectral indices (κ). 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<β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.

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.

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π/0.15)πc/ωpe41.9c/ω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×76×76 intervals, using a Fourier method31,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 ηz,max=6vT1 for κ=, ηz,max=7vT1 for κ = 8, ηz,max=8vT1 for κ = 4, and ηz,max=10vT1 for κ = 2, and half these values, ηx,max=ηy,max=ηz,max/2, in the perpendicular dimensions, where we denoted vT=kBT/me. The maximum represented velocity is given by vmax=π/Δη where Δη is the grid size in the Fourier transformed velocity space, and the represented grid spacing in velocity space is given by Δv=π/ηmax. Therefore, the maximum represented parallel velocities are vz,max20vT for κ=,vz,max17vT for κ = 8, vz,max15vT for κ = 4, and vz,max12vT for κ = 2, and twice as large represented velocities in the perpendicular vx and vy 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.

1.
E.
Marsch
, “
Kinetic physics of the solar corona and solar wind
,”
Living Rev. Sol. Phys.
3
,
1
(
2006
).
2.
S. P.
Gary
and
J.
Wang
, “
Whistler instability: Electron anisotropy upper bound
,”
J. Geophys. Res.
101
,
10749
10754
, doi: (
1996
).
3.
S.
Stverak
,
P.
Travnicek
,
M.
Maksimovic
,
E.
Marsch
,
A. N.
Fazakerley
, and
E. E.
Scime
, “
Electron temperature anisotropy constraints in the solar wind
,”
J. Geophys. Res.
113
,
A03103
, doi: (
2008
).
4.
S. P.
Gary
and
I. H.
Cairns
, “
Electron temperature anisotropy instabilities: Whistler, electrostatic and z mode
,”
J. Geophys. Res.
104
(
A9
),
19835
19842
, doi: (
1999
).
5.
S. P.
Gary
,
K.
Liu
, and
D.
Winske
, “
Whistler anisotropy instability at low electron β: Particle-in-cell simulations
,”
Phys. Plasmas
18
,
082902
(
2011
).
6.
R. A.
Helliwell
,
Whistlers and Related Ionospheric Phenomena
(
Stanford University Press
,
Stanford, CA
,
1965
).
7.
C.
Lacombe
,
O.
Alexandrova
,
L.
Matteini
,
O.
Santolik
,
N.
Cornilleau-Wehrlin
,
A.
Mangeney
,
Y.
de Conchy
, and
M.
Maksimovic
, “
Whistler mode waves and the electron heat flux in the solar wind: CLUSTER observations
,”
Astrophys. J.
796
,
5
(
2014
).
8.
C.
Cattell
,
J. R.
Wygant
,
K.
Goetz
,
K.
Kersten
,
P. J.
Kellogg
,
T.
von Rosenvinge
,
S. D.
Bale
,
I.
Roth
,
M.
Temerin
,
M. K.
Hudson
,
R. A.
Mewaldt
,
M.
Wiedenbeck
,
M.
Maksimovic
,
R.
Ergun
,
M.
Acuna
, and
C. T.
Russell
, “
Discovery of very large amplitude whistler-mode waves in Earths radiation belts
,”
Geophys. Res. Lett.
35
,
L01105
, doi: (
2008
).
9.
A.
Breneman
,
C.
Cattell
,
J.
Wygant
,
K.
Kersten
,
L. B.
Wilson
 III
,
S.
Schreiner
,
P. J.
Kellogg
, and
K.
Goetz
, “
Large-amplitude transmitter-associated and lightning-associated whistler waves in the Earths inner plasmasphere at L < 2
,”
J. Geophys. Res.
116
,
A06310
, doi: (
2011
).
10.
J. M.
Urrutia
,
R. L.
Stenzel
, and
K. D.
Strohmaier
, “
Nonlinear electron magnetohydrodynamics physics. IV. Whistler instabilities
,”
Phys. Plasmas
15
,
062109
(
2008
).
11.
S. P.
Gary
,
Theory of Space Plasma Microinstabilities
(
University Press
,
Cambridge
,
1993
).
12.
R. L.
Mace
, “
Whistler instability enhanced by suprathermal electrons within the Earth's foreshock
,”
J. Geophys. Res.
103
(A
7
),
14643
14654
, doi: (
1998
).
13.
R. L.
Mace
and
R. D.
Sydora
, “
Parallel whistler instability in a plasma with an anisotropic bi-Kappa distribution
,”
J. Geophys. Res.
115
,
A07206
, doi: (
2010
).
14.
M.
Lazar
,
S.
Poedts
, and
R.
Schlickeiser
, “
Instability of the parallel electromagnetic modes in Kappa distributed plasmas—I. Electron whistler-cyclotron modes
,”
Mon. Not. R. Astron. Soc.
410
,
663
670
(
2011
).
15.
M.
Lazar
,
S.
Poedts
, and
M. J.
Michno
, “
Electromagnetic electron whistler-cyclotron instability in bi-Kappa distributed plasmas
,”
Astron. Astrophys.
554
,
A64
(
2013
).
16.
M.
Lazar
,
S.
Poedts
, and
R.
Schlickeiser
, “
The interplay of Kappa and core populations in the solar wind: Electromagnetic electron cyclotron instability
,”
J. Geophys. Res.: Space Phys.
119
,
9395
9406
(
2014
).
17.
M.
Lazar
,
S.
Poedts
,
R.
Schlickeiser
, and
C.
Dumitrache
, “
Towards realistic parametrization of the kinetic anisotropy and the resulting instabilities in space plasmas. Electromagnetic electron-cyclotron instability in the solar wind
,”
MNRAS
446
,
3022
3033
(
2015
).
18.
S. P.
Gary
,
R. S.
Hughes
,
J.
Wang
, and
O.
Chang
, “
Whistler anisotropy instability: Spectral transfer in a three-dimensional particle-in-cell simulation
,”
J. Geophys. Res.: Space Phys.
119
,
1429
1434
(
2014
).
19.
A. A.
Galeev
and
R. Z.
Sagdeev
,
Nonlinear Plasma Theory
(
Benjamin
,
New York
,
1966
).
20.
V. I.
Karpman
, “
Nonlinear effects in the ELF waves propagating along the magnetic field in the magnetosphere
,”
Space Sci. Rev.
16
,
361
388
(
1974
).
21.
P. H.
Yoon
,
J. J.
Seough
,
K. H.
Kim
, and
D. H.
Lee
, “
Empirical versus exact numerical quasilinear analysis of electromagnetic instabilities driven by temperature anisotropy
,”
J. Plasma Phys.
78
,
47
54
(
2012
).
22.
R. C.
Davidson
,
D. A.
Hammer
,
I.
Haber
, and
C. E.
Wagner
, “
Nonlinear development of electromagnetic instabilities in anisotropic plasmas
,”
Phys. Fluids
15
,
317
333
(
1972
).
23.
R. C.
Davidson
and
D. A.
Hammer
, “
Nonequilibrium energy constants associated with large-amplitude electron whistlers
,”
Phys. Fluids
15
,
1282
1284
(
1972
).
24.
S. L.
Ossakow
,
I.
Haber
, and
E.
Ott
, “
Simulation of whistler instabilities in anisotropic plasmas
,”
Phys. Fluids
15
,
1538
1540
(
1972
).
25.
S. L.
Ossakow
,
E.
Ott
, and
I.
Haber
, “
Nonlinear evolution of whistler instabilities
,”
Phys. Fluids
15
,
2314
2326
(
1972
).
26.
Q.
Lu
,
L.
Zhou
, and
S.
Wang
, “
Particle-in-cell simulations of whistler waves excited by an electron κ distribution in space plasma
,”
J. Geophys. Res.
115
,
A02213
, doi: (
2010
).
27.
C. F.
Kennel
and
H. E.
Petschek
, “
Limit on stably trapped particle fluxes
,”
J. Geophys. Res.
71
,
1
28
, doi: (
1966
).
28.
D.
Summers
and
R. M.
Thorne
, “
The modified plasma dispersion function
,”
Phys. Fluids
B
3
,
1835
1847
(
1991
).
29.
H.
Stix
,
Waves in Plasmas
(
Springer-Verlag
,
New York
,
1992
).
30.
B. D.
Fried
and
S. D.
Conte
,
The Plasma Dispersion Function/The Hilbert transform of the Gaussian
(
Academic Press
,
New York
,
1961
).
31.
B.
Eliasson
, “
Outflow boundary conditions for the Fourier transformed three-dimensional Vlasov-Maxwell system
,”
J. Comput. Phys.
225
,
1508
1532
(
2007
).
32.
B.
Eliasson
, “
Numerical simulations of the Fourier transformed Vlasov-Maxwell system in higher dimensions—Theory and applications
,”
Transp. Theory Stat. Phys.
39
(
5&7
),
387
465
(
2010
).
33.
The modified Bessel function of the second kind for non-negative argument and non-negative order is calculated numerically using the Fortran 77 routine RKBESL, see http://www.netlib.org/specfun/rkbesl.
34.
R. A.
Treumann
and
W.
Baumjohann
,
Advanced Space Plasma Physics
(
Imperial College Press
,
London
,
1997
).
35.
M. J.
Aschwanden
,
Physics of the Solar Corona
(
Springer
,
New York
,
2004
).