© The Authors 2026

1. Introduction

In the solar wind, the velocity distributions (VDs) deviate from the idealized Maxwellian profile specific to thermal equilibrium, and exhibit several nonthermal features, such as temperature anisotropy relative to the magnetic field direction, and suprathermal populations as high-energy tails decreasing as Kappa (or κ) power laws of the velocity. Introduced as an empirical model to reproduce the electron VDs observed in the terrestrial magnetosphere (Olbert 1968; Vasyliunas 1968), and later the distributions of all plasma species in the solar wind (Christon et al. 1989, 1991; Collier et al. 1996; Maksimovic et al. 1997), Kappa models have become a revolutionary kinetic tool in the last decades, enabling nonequilibrium plasma states to be described and understood (Štverák et al. 2008; Viñas et al. 2015; Lazar et al. 2015; Dzifčáková et al. 2018; Yoon 2021; Shaaban et al. 2021; López et al. 2021; Moya et al. 2021).

Since space plasmas are weakly collisional, self-generated wave instabilities are expected to control the kinetic anisotropies of the VDs. Observations seem to confirm this hypothesis, such as in the case of temperature anisotropy A = T_⊥/T_∥ ≠ 0 (relative to the direction of the magnetic field), and especially for quasi-thermal core populations at low energies (Štverák et al. 2008; Bale et al. 2009; Shaaban et al. 2025). The observed distributions are sufficiently gyrotropic, and one can invoke bi-Maxwellian models to describe the core populations with temperature anisotropy, and bi-Kappa models for anisotropic distributions with suprathermal tails (Summers & Thorne 1991; Štverák et al. 2008; Lazar et al. 2017a). For bi-Kappa populations, a refined analysis is still necessary (Moya et al. 2021; Shaaban et al. 2025), to account for not only the relaxation of the anisotropy, but also the modifications suffered by the suprathermal tails via the parameter κ, as revealed by numerical simulations (Lazar et al. 2018; Shaaban et al. 2024). This paper proposes a new quasi-linear (QL) approach that generalizes the zero-order QL approach in Moya et al. (2021). The QL theory has the ability to explain the results observed in numerical simulations, in particular the mechanisms by which the growth and stabilization of instabilities are achieved, under the action of the wave fluctuations thus generated. We do this for the concrete case of whistler instability, also known as electromagnetic electron-cyclotron (EMEC) mode instability, triggered by anisotropic electrons with A > 1. However, the new approach proposed here can easily be extended to instabilities of other EM modes, such as the firehose instability (driven by A < 1), as well as to instabilities driven by anisotropic protons and ions.

Our paper is structured as follows. In Sect. 2 we construct the new first-order QL formalism of the whistler instability, which goes beyond the limits of a zero-order approach (Moya et al. 2021), realistically assuming that the quasithermal electron core is still dominant, but cannot evolve independently of the suprathermal component. The new QL model is based on the kurtosis components, introduced as fourth-order moments of the Kappa distribution, which allow us to treat the temporal variation of the parameter κ(t) in a self-consistent manner. The parametric analysis in Sect. 3 is performed on the results obtained for various initial conditions, by varying the key parameters, the temperature anisotropy (A), the plasma beta parameter (parallel component, β_∥), and the parameter κ. The last section, Sect. 4, draws the main conclusions and analyzes a series of perspectives opened by these results.

2. Quasi-linear model

For our model, we consider a quasi-neutral magnetized plasma composed of electrons and cold ions, in which electrons follow a bi-Kappa VD given by

$$\begin{array}{c}\hfill f_0(v_{\perp },v_{\Vert })=\frac{n_0}{{\pi }^{3/2}{\theta }_{\perp }^2{\theta }_{\Vert }}\frac{\mathrm{\Gamma }(\kappa )}{{\kappa }^{1/2}\mathrm{\Gamma }(\kappa -1/2)}(1+\frac{v_{\perp }^2}{\kappa {\theta }_{\perp }^2}+\frac{v_{\Vert }^2}{\kappa {\theta }_{\Vert }^2}{)}^{-\kappa -1},\end{array}$$(1)

where n₀ is the total number density and κ is the power-law index representing the presence of suprathermal particles. The κ index increases for decreasing high-energy tails, and κ → ∞ corresponds to a Maxwellian VD. In Eq. (1) the θ_⊥ and θ_∥ parameters define the parallel, T_∥^(κ), and perpendicular, T_⊥^(κ), kinetic temperatures by the second-order moment of the distribution function. Namely,

$$\begin{array}{cc}\hfill M_{2\Vert }=k_B{T}_{\Vert }^{(\kappa )}& {\displaystyle =\int m_e{v}_{\Vert }^2{f}_0(v_{\perp },v_{\Vert })\mathrm{d}\mathbf{v}=\frac{m_e{\theta }_{\Vert }^2}{2}\frac{\kappa }{\kappa -3/2}}\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill M_{2\perp }=k_B{T}_{\perp }^{(\kappa )}& {\displaystyle =\int \frac{1}{2}{m}_e{v}_{\perp }^2{f}_0(v_{\perp },v_{\Vert })\mathrm{d}\mathbf{v}=\frac{m_e{\theta }_{\perp }^2}{2}\frac{\kappa }{\kappa -3/2},}\hfill \end{array}$$(3)

where k_B is the Boltzmann constant and m_e the electron mass. We can also interpret the θ_j parameters as the thermal speeds of a bi-Maxwellian distribution with temperatures given by

$$\begin{array}{c}\hfill {\theta }_{\Vert }^2=\frac{2k_B{T}_{\Vert }}{m_e},{\theta }_{\perp }^2=\frac{2k_B{T}_{\perp }}{m_e}.\end{array}$$(4)

Thus, combining Eqs. (2), (3), and (4), the relation between the kinetic and Maxwellian temperatures is the following:

$$\begin{array}{c}\hfill T_{\Vert }^{(\kappa )}=F(\kappa )T_{\Vert },T_{\perp }^{(\kappa )}=F(\kappa )T_{\perp },\end{array}$$(5)

where

$$\begin{array}{c}\hfill F(\kappa )=\frac{\kappa }{\kappa -3/2}.\end{array}$$(6)

2.1. Whistler-cyclotron dispersion relation

Under this model, the Vlasov dispersion relation for right-hand polarized waves propagating in the direction of the background magnetic field, B₀, is given by (Viñas et al. 2015; Lazar et al. 2017b; Moya et al. 2021)

$$\begin{array}{c}\hfill 0=\frac{c^2{k}_{\Vert }^2}{{\omega }_{\mathit{pe}}^2}-(A-1)-\frac{A\omega -(A-1)|{\mathrm{\Omega }}_e|}{k_{\Vert }{\theta }_{\Vert }}{Z}_{\kappa }\left(\frac{\omega -|{\mathrm{\Omega }}_e|}{k_{\Vert }{\theta }_{\Vert }}\right),\end{array}$$(7)

where ω and k_∥ are the wave frequency and wavenumber, respectively, and ${\omega }_{\mathit{pe}}=\sqrt{4\pi n_0{e}^2/m_e}$ and Ω_e = eB₀/m_ec are the electron plasma frequency and gyrofrequency. Also, in Eq. (7) e is the elementary charge, c is the speed of light, and

$$\begin{array}{c}\hfill Z_{\kappa }(\xi )=\frac{1}{{\pi }^{1/2}{\kappa }^{1/2}}\frac{\mathrm{\Gamma }(\kappa )}{\mathrm{\Gamma }(\kappa -1/2)}{\int }_{-\infty }^{\infty }\mathrm{d}x\frac{{(1+x^2/\kappa )}^{-\kappa }}{x-\xi },\mathrm{Im}(\xi )>0\end{array}$$(8)

is the modified dispersion function (Hellberg & Mace 2002; Lazar et al. 2008). It is important to mention that in the Maxwellian limit, κ → ∞, Z_k becomes the well-known plasma dispersion function, Z, defined by Fried & Conte (1961).

For temperature anisotropy

$$\begin{array}{c}\hfill A=\frac{T_{\perp }^{(\kappa )}}{T_{\Vert }^{(\kappa )}}=\frac{{\theta }_{\perp }^2}{{\theta }_{\Vert }^2}=\frac{T_{\perp }}{T_{\Vert }}>1,\end{array}$$(9)

the dispersion relation Eq. (7) is unstable to whistler EMEC waves, with a maximum growth rate that increases with increasing temperature anisotropy, A > 1, and increasing parallel plasma beta, β_∥ = 8πn₀k_BT_∥/B₀² (see e.g. Gary 1993, and references therein). This behavior is similar to that of a collisionless plasma composed of anisotropic electrons, and, for Kappa-distributed electrons, several studies have already been conducted (Viñas et al. 2015; Lazar et al. 2017a, 2018). In the applications, we consider combinations of temperature anisotropy and plasma beta, and for each combination various values of κ. For comparison between Kappa distributions and the Maxwellian limit (κ → ∞), we consider the kinetic beta parameter

$$\begin{array}{c}\hfill {\beta }_{\Vert }^{(\kappa )}=\frac{8\pi n_0{k}_B{T}_{\Vert }^{(\kappa )}}{B_0^2}=F(\kappa )\frac{8\pi n_0{k}_B{T}_{\Vert }}{B_0^2}=F(\kappa ){\beta }_{\Vert },\end{array}$$(10)

which reduces to the lower Maxwellian limit ${\beta }_{\Vert }^{(\kappa )}\to {\beta }_{\Vert }<{\beta }_{\Vert }^{(\kappa )}$, when κ → ∞.

2.2. QL evolution of power-index κ

One of the most traditional ways to include nonlinear effects of kinetic instabilities is the use of QL theory (see e.g. Yoon 2017, and references therein), which provides expressions for the time evolution of the plasma VD and electromagnetic fluctuations due to nonlinear interactions

$$\begin{array}{cc}\hfill \frac{\partial f_0}{\partial t}& {\displaystyle =\frac{i{e}^2}{4m_e^2{c}^2}{\int }_{-\infty }^{\infty }\frac{d{k}_{\Vert }}{k_{\Vert }^2}[({\omega }^{\ast }-k_{\Vert }{v}_{\Vert })\frac{\partial }{\partial v_{\perp }}+k_{\Vert }{v}_{\perp }\frac{\partial }{\partial v_{\Vert }}]}\hfill \\ \hfill & \times \frac{v_{\perp }^2\delta B^2(k_{\Vert })}{\omega -k_{\Vert }{v}_{\Vert }+{\mathrm{\Omega }}_e}[(\omega -k_{\Vert }{v}_{\Vert })\frac{\partial f_0}{\partial v_{\perp }}+k_{\Vert }{v}_{\perp }\frac{\partial f_0}{\partial v_{\Vert }}],\hfill \end{array}$$(11)

where ω^* represents the complex conjugate of the complex solutions of the dispersion relation Eq. (7), ω = ω(k_∥)+iγ(k_∥), and δB²(k_∥) is the wave energy density whose temporal evolution is given by the wave kinetic equation

$$\begin{array}{c}\hfill \frac{\partial \delta B^2(k_{\Vert })}{\partial t}=2\gamma (k_{\Vert })\delta B^2(k_{\Vert }).\end{array}$$(12)

Then, taking the temporal derivative on both sides of Eqs. (2) and (3), and using Eq. (11), from QL theory we obtain

$$\begin{array}{cc}\hfill k_B\frac{\mathrm{d}{T}_{\Vert }^{(\kappa )}}{\mathrm{d}t}& =\int \frac{1}{2}{m}_e{v}_{\Vert }^2\frac{\partial f_0}{\partial t}\mathrm{d}\mathbf{v}\hfill \\ \hfill & =\frac{e^2}{m_e}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}{k}_{\Vert }}{c^2{k}_{\Vert }^2}\gamma (k_{\Vert })(\frac{c^2{k}_{\Vert }^2}{{\omega }_{\mathit{pe}}^2}+1)\delta B^2(k_{\Vert }),\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill k_B\frac{\mathrm{d}{T}_{\perp }^{(\kappa )}}{\mathrm{d}t}& =\int m_e{v}_{\perp }^2\frac{\partial f_0}{\partial t}\mathrm{d}\mathbf{v}\hfill \\ \hfill & =-\frac{e^2}{m_e}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}{k}_{\Vert }}{c^2{k}_{\Vert }^2}\gamma (k_{\Vert })(\frac{c^2{k}_{\Vert }^2}{{\omega }_{\mathit{pe}}^2}+\frac{1}{2})\delta B^2(k_{\Vert }).\hfill \end{array}$$(14)

Equations (14) and (13), together with the wave kinetic equation, Eq. (12), and the dispersion relation, Eq. (7), provide a closed system of differential equations that describe the time evolution and saturation of the EMEC instability in terms of the coupled variation of the plasma temperatures and the wave energy. From the expressions given in (5), their variation in time (i.e., $d{T}_{\Vert ,\perp }^{(\kappa )}/dt\ne 0$) implies variations in time not only for T_∥, ⊥ (as assumed in a simplified model) but also for the exponent κ. In previous studies, it has been customary to consider κ = κ₀ = constant, as is assumed in several QL approaches with κ distributions (see e.g Lazar et al. 2017a, 2018). Thus, in such cases, Eqs. (7), (12), and (14)–(13) provide a complete QL description of a temperature anisotropy instability and its relaxation.

On the other hand, if κ changes in time, i.e., κ = κ(t), as is suggested by the Vlasov and PIC simulations in Lazar et al. (2017b) and Lazar et al. (2018), respectively, we require an additional independent equation – different from but complementary to Eqs. (14) and (13) – to account for the time evolution of κ. Since f₀ is symmetric in the magnetic field direction, all odd moments are identically zero. Therefore, the additional equation can be obtained by considering the fourth-order moment (the smallest non-null moment besides the kinetic temperature and number density) of the Kappa VD, which is related to the kurtosis of the distribution. The kurtosis, K, is a well-known measure of the relative weight of the high-energy tails of the VD with respect to the core; namely,

$$\begin{array}{cc}\hfill K=\int v^4{f}_0(v_{\perp },v_{\Vert })d\mathbf{v}=& \left(\frac{\kappa }{\kappa -3/2}\right)\left(\frac{\kappa }{\kappa -5/2}\right)\hfill \\ \hfill & \times (2{\theta }_{\perp }^4+{\theta }_{\perp }^2{\theta }_{\Vert }^2+\frac{3}{4}{\theta }_{\Vert }^4).\hfill \end{array}$$(15)

Note that, for the convergence of the fourth-order moment, κ > 5/2 is required, a condition that we also take into account when obtaining and analyzing the results. Note that the adoption of the regularized Kappa distribution eliminates the need to restrict the Kappa index (Scherer et al. 2017; Thanh et al. 2022). However, for the sake of simplicity, the present study does not take such an additional feature into account.

In terms of temperatures, using Eq. (4) and (5), the kurtosis is given by

$$\begin{array}{cc}\hfill K(t)=\int v^4& f_0(v_{\perp },v_{\Vert })\mathrm{d}\mathbf{v}=\frac{4k_B^2}{m_e^2}\frac{\kappa (t)}{\kappa (t)-3/2}\frac{\kappa (t)}{\kappa (t)-5/2}\hfill \\ \hfill & \times [2{(T_{\perp }(t))}^2+T_{\perp }(t)T_{\Vert }(t)+\frac{3}{4}{(T_{\Vert }(t))}^2].\hfill \end{array}$$(16)

Therefore, in the case of Kappa VDs, the kurtosis is a measure or indicator of the shape of the distribution and is highly controlled by the κ parameter. Thus, considering the time derivative of kurtosis K provides a way to obtain a dynamic equation for the time evolution of κ. Indeed, taking the temporal derivative on both sides of Eqs. (2), (3), and (16), and using Eq. (11), from QL theory we obtain the following system of coupled differential equations:

$$\begin{array}{cc}\hfill \left(\frac{m_e}{k_B}\right)\frac{\mathrm{d}{M}_{2\Vert }}{\mathrm{d}t}& =F(\kappa )\frac{\mathrm{d}{T}_{\Vert }}{\mathrm{d}t}+T_{\Vert }\frac{\mathrm{d}F}{\mathrm{d}\kappa }\frac{\mathrm{d}\kappa }{\mathrm{d}t}\hfill \end{array}$$(17)

$$\begin{array}{cc}\hfill \left(\frac{m_e}{k_B}\right)\frac{\mathrm{d}{M}_{2\perp }}{\mathrm{d}t}& =F(\kappa )\frac{\mathrm{d}{T}_{\perp }}{\mathrm{d}t}+T_{\perp }\frac{\mathrm{d}F}{\mathrm{d}\kappa }\frac{\mathrm{d}\kappa }{\mathrm{d}t}\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill {\left(\frac{m_e}{2k_B}\right)}^2\frac{\mathrm{d}K}{\mathrm{d}t}& =G(\kappa )(4T_{\perp }+T_{\Vert })\frac{\mathrm{d}{T}_{\perp }}{\mathrm{d}t}\hfill \\ \hfill & +G(\kappa )(T_{\perp }+\frac{3}{2}{T}_{\Vert })\frac{\mathrm{d}{T}_{\Vert }}{\mathrm{d}t}\hfill \\ \hfill & +(2T_{\perp }^2+T_{\perp }{T}_{\Vert }+\frac{3}{4}{T}_{\Vert }^2)\frac{\mathrm{d}G}{\mathrm{d}\kappa }\frac{\mathrm{d}\kappa }{\mathrm{d}t},\hfill \end{array}$$(19)

where

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{M}_{2\Vert }}{\mathrm{d}t}& =\frac{{\omega }_{\mathit{pe}}^2}{m_e{n}_0}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}{k}_{\Vert }}{c^2{k}_{\Vert }^2}\gamma (k_{\Vert })(\frac{c^2{k}_{\Vert }^2}{{\omega }_{\mathit{pe}}^2}+1)\delta B^2(k_{\Vert })\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{M}_{2\perp }}{\mathrm{d}t}& =-\frac{{\omega }_{\mathit{pe}}^2}{m_e{n}_0}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}{k}_{\Vert }}{c^2{k}_{\Vert }^2}\gamma (k_{\Vert })(\frac{c^2{k}_{\Vert }^2}{{\omega }_{\mathit{pe}}^2}+\frac{1}{2})\delta B^2(k_{\Vert })\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}K}{\mathrm{d}t}& =\frac{4k_B{T}_{\Vert }{\omega }_{\mathit{pe}}^2}{m_e^2{n}_0}\frac{{\kappa }^{1/2}{(\kappa -1)}^{1/2}}{\kappa -3/2}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}{k}_{\Vert }}{c^2{k}_{\Vert }^2}\mathrm{Im}\{[\frac{{|\omega |}^2}{k_{\Vert }{\theta }_{\Vert }}\hfill \\ \hfill & +{\left(\frac{\kappa }{\kappa -1}\right)}^{1/2}(\frac{T_{\perp }}{T_{\Vert }}-1){\omega }^{\ast }{\xi }_{\kappa }]Z_{\kappa -1}({\xi }_{\kappa })\}\delta B^2(k_{\Vert })\hfill \\ \hfill & +\frac{4k_B{T}_{\Vert }{\omega }_{\mathit{pe}}^2}{m_e^2{n}_0}{\int }_{-\infty }^{\infty }\frac{\mathrm{d}{k}_{\Vert }}{c^2{k}_{\Vert }^2}\gamma (k_{\Vert })[\left(\frac{{\omega }^2(k_{\Vert })+{\gamma }^2(k_{\Vert })-{\mathrm{\Omega }}_e^2}{k_{\Vert }^2{\theta }_{\Vert }^2}\right)\hfill \\ \hfill & \times \frac{c^2{k}_{\Vert }^2}{{\omega }_{\mathit{pe}}^2}-\frac{1}{2}\left(\frac{\kappa }{\kappa -3/2}\right)(\frac{T_{\perp }}{T_{\Vert }}-1)(1+4\frac{T_{\perp }}{T_{\Vert }})\hfill \\ \hfill & +\frac{{\omega }^2(k_{\Vert })+{\gamma }^2(k_{\Vert })}{k_{\Vert }^2{\theta }_{\Vert }^2}]\delta B^2(k_{\Vert }).\hfill \end{array}$$(22)

Here,

$$\begin{array}{c}\hfill {\xi }_{\kappa }={\left(\frac{\kappa }{\kappa -1}\right)}^{1/2}\left(\frac{\omega -|{\mathrm{\Omega }}_e|}{k_{\Vert }{\theta }_{\Vert }}\right),\end{array}$$(23)

F(κ) is given by Eq. (6), and

$$\begin{array}{c}\hfill G(\kappa )=\left(\frac{\kappa }{\kappa -3/2}\right)\left(\frac{\kappa }{\kappa -5/2}\right).\end{array}$$(24)

In summary, Eq. (19) provides a new independent equation to close the system formed by Eqs. (7), (12), (17), and (18), allowing the self-consistent variation of the κ parameter due to the relaxation of electron VD under the action of EMEC instability. It is important to mention that energy conservation is guaranteed within our QL approach. Indeed, although Eqs. (17)–(22) go beyond the typical QL picture, the number of particles, momentum, and energy conservation are always satisfied in QL theory, regardless of the shape of the VD. A proof of these facts can be found in plasma physics textbooks (see e.g. Chapter 10 in Krall & Trivelpiece 1973).

3. Results

To solve the QL equations, we employed a discrete grid within the wavenumber space, normalized to the electron inertial length c/ω_pe, featuring 256 points ranging from 0.001 < |ck_∥|/ω_pe < 3. All cases were simulated up to |Ω_e|t_end = 655 (or 2¹⁵ time steps) using a time step of Δt = 0.02/|Ω_e|. The magnetic field spectrum was specified as δB_k²/B₀² = 10⁻⁵ at t = 0. Following this methodology, and knowing the magnetic field spectrum and the value of all parameters (A, β_∥, and κ) at any given time, t, we solved the dispersion relation to determined the complex frequency, ω + iγ, of whistler-cyclotron waves as a function of k_∥ at that specific moment. Subsequently, we evaluated the time derivative of each parameter. The entire system was then advanced to the next time, t + dt, using a second-order Runge-Kutta method. Furthermore, in all cases, we considered ω_pe/|Ω_e| = 20, m_p/m_e = 1836, and cold ions, but in each case, different initial conditions (t₀ = 0) were set in terms of the key electron parameters. This computational strategy was applied to solve the QL system using 560 combinations of κ, temperature anisotropy, and parallel beta as initial conditions, as is shown in Table 1.

Figures 1–3 present the QL solutions in detail, such as the time variations in the properties of the unstable EMEC waves, such as the maximum growth rate, γ_max/|Ω_e| (normalized by the electron gyrofrequency), and the wave magnetic energy density, δB²/B₀² (normalized by the regular magnetic energy density), as well as the variations in the key electron parameters on which these instabilities depend, such as the (parallel) beta parameter, β_∥(t), the temperature anisotropy, A(t) = T_⊥/T_∥(t), and of course the parameter κ(t). The new results obtained with a more realistic QL model, in which the κ parameter is allowed to vary self-consistently over time, are presented and compared with those from the simplified QL model, in which κ is held constant.

Fig. 1. QL runs with the same initial anisotropy, A(0) = 2, three values of initial β∥(0) = 0.1 (left), 0.5 (middle), and 8.0 (right), and for eight different (initial) values of κ(0) shown in various shades of brown. The first line (top) indicates the variation in κ, and the following ones indicate the variations in key parameters β∥(t) and A(t) = T⊥(t)/T∥(t), and variations in wave properties δB2(t)/B02 and γmax(t), for both cases in which κ varies in time and ones in which κ is constant.

Fig. 2. QL runs with the same initial β∥(0) = 0.1, three values of initial κ(0) = 2.6 (left), 4.0 (middle), and 8.0 (right), and for six different (initial) anisotropies, T⊥/T∥(0), shown in various shades of brown. The first line (top) indicates the variation in κ, and the following ones indicate the variations in key parameters β∥(t) and A(t) = T⊥(t)/T∥(t), and variations in wave properties δB2(t)/B02 and γmax(t), for both cases in which κ varies in time and ones in which κ is constant.

Fig. 3. QL runs with the same initial κ(0) = 8.0, three values of initial anisotropy A(0) = 2.0 (left), 3.9 (middle), and 6.0 (right), and for ten different (initial) values of β∥(0) shown in various shades of brown. The first line (top) indicates the variation in κ, and the following ones indicate the variations in key parameters β∥(t) and A(t) = T⊥/T∥(t), and variations in wave properties δB2(t)/B02 and γmax(t), for both cases in which κ varies in time and ones in which κ is constant.

Figure 1 displays QL runs with the same initial anisotropy A(0) = 2, and each column corresponds to three different values of β_∥(0) = 0.10 (left), 0.50 (middle), and 8.0 (right). The solutions represented in eight shades of brown correspond, from the darkest to the lightest, to different values of initial κ(0) = 2.6, 2.8, 3.0, 4.0, 5.0, 6.0, 7.0, and 8.0. The difference from simplified QL models is the variation in time of the parameter κ, which may eventually decrease if β_∥(0) is not very small (second and third columns), and also if the initial value of κ(0) is large enough. The larger the β_∥(0), the lower the initial value of κ(0) affected by the decrease in time. The QL decrease in κ(t) over time does not indicate a relaxation toward thermodynamic equilibrium, characterized by quasi-Maxwellian distributions with κ (much) larger than the initial one, but an increase in high-energy suprathermal tails. Whistler waves generated by the instability can contribute to electron energization, enhancing suprathermal tails (Ma & Summers 1999; Vocks et al. 2008), a sufficiently robust effect captured by simulations (Lazar et al. 2017b, 2018) and also indicated by the first QL models with zero-order approaches for the temporal variation in the kappa parameter (Moya et al. 2021).

Figure 2 shows QL runs with the same initial β_∥(0) = 0.1 and each column corresponds to three different initial values of κ(0) = 2.6 (left), 4.0 (middle), and 8.0 (right). The solutions in eight shades of brown correspond, from the darkest to the lightest, to different values of initial anisotropy A(0) = T_⊥/T_∥(0) = 2.0, 2.5, 3.1, 3.9, 4.8, and 6.0. In this case, the temporal evolutions obtained for the properties of the wave fluctuations and for the key parameters describing the anisotropic electrons show no significant differences between the models with κ= constant and those with variable κ. However, it should be noted that when κ can vary, it does so in both directions. When the initial value of κ(0) is not very high, as in the first two cases, κ(0) = 2.6 and 4, we obtain a decrease in κ; that is, an increased suprathermalization, as was also found in the previous cases. However, it should be noted that these variations in κ(t) are much smaller than in Fig. 1. In contrast, in the latter case, when the initial suprathermalization is reduced, with κ(0) = 8, this parameter tends to increase further, as evidence of the relaxation of the initial distribution through thermalization. This evolution is apparently more robust in the case of small initial anisotropies, T_⊥/T_∥(0) < 4, while for larger anisotropies the initial increase in κ saturates to a maximum (around κ = 10), after which it starts to decrease, most likely under the action of an increased level of fluctuations generated by instability in this case.

In Fig. 3 we plot QL runs with the same initial κ(0) = 8 and each column corresponds to three different initial anisotropies: A(0) = 2.0 (left), 3.9 (middle), and 6.0 (right). The solutions in eight shades of brown correspond, from the darkest to the lightest, to different values of initial β_∥(0) = 0.1, 0.13, 0.16, 0.20, 0.25, 0.50, 1.00, 2.00, 4.0, and 8.0. As in Fig. 2, the variations in κ(t) are similarly non-monotonic, but show different temporal profiles very sensitive to the initial values of the parameters. In this case, the initial κ(0) = 8 is already large, but κ(t) can still increase over time, even for larger anisotropies, if the initial β_∥(0) is low enough. If the initial anisotropy is not very low, for example A(t) = 3.9 and 6.0, the increase in κ(t) occurs in longer runs, after a short but significant drop. In general, these increases are not very steep, and the final values reached at saturation do not necessarily exceed the initial ones by much.

To these details, we add summarizing graphs in Figs. 4, 5, and 6 that provide insights from the long runs, including those specific to instability saturation at the end of the QL simulations. Thus, Fig. 4 displays diagrams with the dynamic paths of the complete QL evolutions as a function of the main kinetic parameters, T_⊥/T_∥ and β_∥. The quasi-stationary states reached at saturation correspond to an instability threshold of γ_max/|Ω_e| = 5 × 10⁻³. For the same dynamic paths, we present with different color codes the relative variation in the parameter κ(t) (top), the variation in |ξ_r = ω_r − |Ω_e|/(kθ_∥)|, and the temporal variation in the magnetic wave energy density δB²(t)/B₀². The top panels provide a comprehensive picture of the regimes in which the parameter κ(t) decreases or increases. Its increase, meaning the relaxation of the electron distribution toward less suprathermal states, is limited only to regimes with sufficiently small β_∥ < 0.3, and only when the initial value of κ is sufficiently large, such as in the third case with κ(0) = 8. In these regimes, the levels of wave fluctuations induced by instability remain low, as is shown in the corresponding bottom panel. Otherwise, for β_∥ > 0.3, the variation in κ is non-monotonic with an initial increase followed by decreases. ξ_r, as the real part of the argument of the plasma dispersion function, is an indicator of the resonant wave-electron interaction, which becomes maximum when |ξ_r|∼1 (see the narrow white bands), involving particles with a velocity close to the thermal velocity. |ξ_r|< 1 (bluish) implies the resonance of the core electrons that damp rather than amplify the waves, while |ξ_r|> 1 (reddish) signifies the resonance of electrons from the high-energy (suprathermal) tails that amplify the waves (Lazar et al. 2019, 2022). At small anisotropies near the threshold, the suprathermal populations are highly resonant. This is also the case for regimes with low β_∥ < 1, where an increase in the parameter κ may eventually occur. In addition, the bottom panels show that, as was expected, an increased wave energy density is obtained for large values of the anisotropy and beta parameters.

Fig. 4. QL dynamic paths in diagrams of T⊥/T∥ vs. β∥, showing for each run time variations (color coded) in κ(t) (top), the resonance factor (middle, see text), and the wave energy density (bottom). Each column groups runs with the same initial κ(0) = 2.6 (left), 4.0 (middle), and 8.0 (right).

Fig. 5. Diagram of T⊥/T∥ vs. β∥, showing final values of κ(tend) (color-coded) obtained for the final states (filled circles), which reach the instability thresholds (light gray, γmax/|Ωe| = 5 × 10−3) predicted by linear theory for different initial values of κ(0) (as is partially indicated).

Fig. 6. Relative variation in final κ(tend) with respect to the initial values, as a function of initial values for κ (panels a to c), T⊥/T∥(0) (panels d to f), and β∥(0) (all panels).

Figure 5 uses the same diagram of (T_⊥/T_∥ vs. β_∥) states to show only the final states, with filled circles of different colors depending on the final value of the parameter κ(t_end). These are aligned with the corresponding EMEC thresholds obtained for different initial values of κ(0), some of which are indicated to mark their increasing order. If large final values of κ(t_end) are obtained, it is clear that these are specific to regimes with low β < 1, but also to sufficiently high initial values of κ(0) > 5. Otherwise, we are dealing with an increased suprathermalization over time; that is, a decrease in κ to values of κ(t_end) < 5, also associated with the lowest anisotropies in the quasi-stationary final states. These results are in complete agreement with previous analyses, both in linear and QL theory and in numerical simulations (Lazar et al. 2019, 2022; Moya et al. 2021).

Figure 6 outlines the total relative variations obtained for the parameter κ after saturation of the EMEC instability. The six panels present these relative variations as a function of the initial conditions characterized by T_⊥/T_∥(0), β_∥(0), and κ(0). Positive relative variations are found only in panels (c), (e), and (f), which are specific to the low regimes β_∥(0) < 0.3. In all other cases, the final value of the kappa parameter is less than the initial value, κ(t_end) < κ(0), and this behavior strongly depends on the initial conditions (T_⊥/T_∥(0), β_∥(0)), which is subsequently related to the amount of free energy available to generate electromagnetic waves. Thus, it should be noted that, under the effect of the EMEC instability, which increases wave fluctuations, the predominance of regimes in which the parameter κ decreases also enhances the suprathermalization of electrons.

To further explore the connection between the decrease in kappa and the increase in the level of wave fluctuations in the system, Fig. 7 shows a scatter plot of κ(t_end) as a function of κ(0). In the figure, each dot corresponds to a given QL run, the color represents the level of magnetic energy at the end of each calculation t = t_end, and the dashed line corresponds to the identity. From the figure, we can see a clear correlation between the level of magnetic fluctuations and the evolution of the kappa index, quantifying the trends shown in Fig. 4. The larger the level of magnetic energy at the end of the QL run, the more pronounced the departure of the kappa parameter from its initial value. These results indicate that the relaxation of instabilities through the generation of magnetic fluctuations almost always affects the kappa value. For a large initial kappa value (κ(0) > 6), the generation of electromagnetic waves induces an increase in κ, suggesting that the relaxation of the EMEC instability induces a thermalization of the VD.

Fig. 7. Scatter plot showing the final value of the kappa parameter, κ(tend), as a function of the initial value. Each dot corresponds to a given QL run. The color bar represents the level of magnetic energy at the end of each calculation. The dashed line corresponds to the identity.

However, for κ < 6, the relaxation of the EMEC instability usually results in kappa decreasing, meaning that the VD moves away from the Maxwellian limit, which defies the intuition that equilibrium corresponds to Maxwellian distributions. These small kappa regimes are ubiquitous to the solar wind (Lazar et al. 2017a, 2020; Eyelade et al. 2025) and planetary environments, such as the Earth’s magnetosphere (Espinoza et al. 2018; Eyelade et al. 2021), where the beta plasma parameter for electrons is generally sufficiently large, β_e > 0.3 (Štverák et al. 2008; Wilson et al. 2019a,b; Abraham et al. 2022). In situ observations, particularly in the solar wind, reveal mainly quasi-stationary states in the vicinity of the lowered EMEC thresholds, whose suprathermalization increases with heliocentric distance as the parameter κ decreases toward 1 AU and beyond.

4. Conclusions

In this paper we have proposed a new dynamical model for the QL analysis of EMEC (or whistler) instabilities triggered by anisotropic electrons with bi-Kappa VDs, as observed in space plasmas. In contrast to the first QL models, which constrain the parameter κ (or the distribution’s functional form) to remain constant, the new approach allows κ to vary self-consistently over time. Moya et al. (2021) proposed a first zero-order model with variable κ in which they nevertheless constrained the low-energy core populations (described by the Maxwellian limit κ → ∞ of Kappa distributions) to remain isolated from the suprathermal ones (in the high-energy tails of Kappa distributions). In the new QL approach, we have also overcome this constraint by completing the QL equations with those for the fourth-order moment (kurtosis).

The results refine and extend the previous ones, this time making a clear distinction between alternative regimes that lead to either a decrease or an increase in the parameter κ during the stabilization of the instability. In agreement with the zero-order approach, a decrease in κ, i.e., an excess of suprathermalization, prevails. However, this time, the (initial) regimes in which the VDs can relax to a quasi-Maxwellian form due to the increase in κ are also identified. These (initial) conditions are generally restricted to low beta parameters, but not to very small initial values for κ ≳ 5. The main conclusion is that, in general, the bi-Kappa electron VDs undergo only a partial relaxation under the action of the EMEC instability. Only the temperature anisotropy decreases, while wave fluctuations result in excess and maintain but also enhance the suprathermal component, as is evidenced by the decrease in κ. In other words, for small κ regimes, during the relaxation of the EMEC instability, electromagnetic instabilities are produced as temperature anisotropy decreases due mainly to resonant interactions with the quasi-thermal core (see the central panel of Fig. 4). Subsequently, the self-generated waves can interact with suprathermal particles, generating energetic power-law tails represented by a small κ value. Additionally, the relations obtained between magnetic fluctuations and the kappa index suggest that a fixed value of κ is not an adequate strategy for studying kinetic instabilities in nonthermal plasmas. Thus, the consideration of a dynamic kappa should be customary for the study of nearly collisionless space plasmas.

That said, it is worth mentioning that the whistler mode can be driven by temperature anisotropies and/or heat flux. In the case of the heat-flux-driven instability, the free-energy source corresponds to the asymmetry of the electron velocity distribution function (VDF), quantified by the third moment of the VDF, as is shown in Tong et al. (2019). In the solar wind, this is usually provided by inherent asymmetry (Zenteno-Quinteros et al. 2021; Zenteno-Quinteros & Moya 2022; Zenteno-Quinteros et al. 2023) or the relative drift between the core and halo (Vasko et al. 2020) or strahl (López et al. 2020). In this case, to focus on the temperature-anisotropy-driven EMEC instability, our model assumes a single, symmetric Kappa distribution, such that all odd moments of the VDF are identically zero. It would be interesting to address the role of asymmetry and heat flux for the evolution of Kappa distributions during the relaxation of the whistler heat-flux instability. In such a case, a model based on an asymmetric VDF will allow for the third-order moment to be considered in order to complement the kurtosis equation, or replace it, under the assumption of a constant electron drift speed along the magnetic-field direction. However, analyzing such a model is beyond the scope of the current manuscript and will be left for future work.

In summary, our results suggest that the microscopic and macroscopic characteristics of the plasma, and the fine structure of electromagnetic field fluctuations, waves, and turbulence, strongly depend on the shape of the VD and its evolution, and vice versa, particularly after excitation and during relaxation processes associated with collisionless dissipation. Our results already show preliminary agreement with in situ observations in the solar wind and the magnetosphere, suggesting that the new QL model could provide a theoretical basis for explaining similar kinetic instabilities in natural plasmas with bi-Kappa-type distributions. This hypothesis may benefit in the future from confirmations through numerical simulations and a detailed analysis of anisotropic distributions during relaxation under the action of instabilities.

Acknowledgments

The authors acknowledge support from the University of Chile, Ruhr-University Bochum, and the Katholieke Universiteit Leuven. Portions of this research were supported by the International Space Science Institute (ISSI) in Bern, through the ISSI International Team project 24-612: Excitation and Dissipation of Kinetic-Scale Fluctuations in Space Plasmas. These results were also obtained in the framework of the projects C16/24/010 (C1 project Internal Funds KU Leuven), G002523N (FWO-Vlaanderen), 4000145223 SIDC Data Exploitation (SIDEX2), ESA Prodex, and Fondecyt ANID Chile No. 1240281 (PSM) and No. 1251712 (RAL). SP is funded by the European Union (ERC, Open SESAME, 101141362). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

All Tables

All Figures

Fig. 1. QL runs with the same initial anisotropy, A(0) = 2, three values of initial β∥(0) = 0.1 (left), 0.5 (middle), and 8.0 (right), and for eight different (initial) values of κ(0) shown in various shades of brown. The first line (top) indicates the variation in κ, and the following ones indicate the variations in key parameters β∥(t) and A(t) = T⊥(t)/T∥(t), and variations in wave properties δB2(t)/B02 and γmax(t), for both cases in which κ varies in time and ones in which κ is constant.

In the text

Fig. 2. QL runs with the same initial β∥(0) = 0.1, three values of initial κ(0) = 2.6 (left), 4.0 (middle), and 8.0 (right), and for six different (initial) anisotropies, T⊥/T∥(0), shown in various shades of brown. The first line (top) indicates the variation in κ, and the following ones indicate the variations in key parameters β∥(t) and A(t) = T⊥(t)/T∥(t), and variations in wave properties δB2(t)/B02 and γmax(t), for both cases in which κ varies in time and ones in which κ is constant.

In the text

Fig. 3. QL runs with the same initial κ(0) = 8.0, three values of initial anisotropy A(0) = 2.0 (left), 3.9 (middle), and 6.0 (right), and for ten different (initial) values of β∥(0) shown in various shades of brown. The first line (top) indicates the variation in κ, and the following ones indicate the variations in key parameters β∥(t) and A(t) = T⊥/T∥(t), and variations in wave properties δB2(t)/B02 and γmax(t), for both cases in which κ varies in time and ones in which κ is constant.

In the text

	Fig. 4. QL dynamic paths in diagrams of T⊥/T∥ vs. β∥, showing for each run time variations (color coded) in κ(t) (top), the resonance factor (middle, see text), and the wave energy density (bottom). Each column groups runs with the same initial κ(0) = 2.6 (left), 4.0 (middle), and 8.0 (right).
In the text

	Fig. 5. Diagram of T⊥/T∥ vs. β∥, showing final values of κ(tend) (color-coded) obtained for the final states (filled circles), which reach the instability thresholds (light gray, γmax/|Ωe| = 5 × 10−3) predicted by linear theory for different initial values of κ(0) (as is partially indicated).
In the text

	Fig. 6. Relative variation in final κ(tend) with respect to the initial values, as a function of initial values for κ (panels a to c), T⊥/T∥(0) (panels d to f), and β∥(0) (all panels).
In the text

	Fig. 7. Scatter plot showing the final value of the kappa parameter, κ(tend), as a function of the initial value. Each dot corresponds to a given QL run. The color bar represents the level of magnetic energy at the end of each calculation. The dashed line corresponds to the identity.
In the text