The Wave Energy Density and Growth Rate for the Resonant Instability in Relativistic Plasmas

We start our derivation from analysing the dielectric tensor $\overleftrightarrow{\mathbf{K}}$ in relativistic plasmas. The derivation of $\overleftrightarrow{\mathbf{K}}$ is explained by many plasma physics textbooks such as Stix (1992) and Gurnett & Bhattacharjee (2017), and the dielectric tensor in relativistic plasmas is then

where

and $\overleftrightarrow{\mathbf{1}}$ is the unit tensor. The subscript j indicates the particle species. We denote the relativistic momentum coordinates perpendicular and parallel with respect to $\mathbf{B}_{0}$ as $p_{\perp}=m_{j}v_{\perp}$ and $p_{\parallel}=m_{j}v_{\parallel}$ where $m_{j}=\gamma m_{0j}$, $m_{0j}$ is the rest mass of particle species j, $v_{\perp}$ and $v_{\parallel}$ are the perpendicular and parallel velocity coordinates, and $\gamma=\sqrt{1+(p_{\parallel}^{2}+p_{\perp}^{2})/m_{0j}^{2}c^{2}}$. The normalised and gyrotropic particle MDF of species j is denoted as $f_{j}\equiv f_{j}(p_{\perp},p_{\parallel})$ so that $\int d^{3}\mathbf{p}f_{j}=1$. Note that, because $f_{j}$ is gyrotropic, we are free to select the value of $\phi_{k}$ in Eqs. (8) to (15).

The integer $n$ determines the order of the resonance, where $n=0$ corresponds to the Landau resonance and $n\neq 0$ corresponds to cyclotron resonances. The plasma frequency of species j is defined as $\omega_{pj}^{2}=\omega_{0pj}^{2}/\gamma=4\pi n_{j}q_{j}^{2}/m_{j}$ where $q_{j}$ and $n_{j}$ are the charge and density of particle species j. We denote the $n^{th}$-order Bessel function as $J_{n}=J_{n}(\rho_{j})$ where $\rho_{j}\equiv k_{\perp}v_{\perp}/\Omega_{j}$. The cyclotron frequency of species j is defined as $\Omega_{j}=\Omega_{0j}/\gamma=q_{j}B_{0}/m_{j}c$ and $c$ is the light speed. In this paper, we express that $\sum_{j}\sum_{n}\equiv\sum_{j,n}$.

The homogeneous wave equation is then expressed as

The condition that the non-trivial solution of Eq. (19) exists is that the determinant of the $3\times 3$ matrix in Eq. (19) must be zero. This determinant corresponds to the plasma dispersion relation. Note that Eq. (19) reduces to the homogeneous wave equation leading to the cold-plasma dispersion relation if we apply the cold-particle MDF, expressed in the plasma-rest frame as

We denote Dirac’s $\delta$-function as $\delta$. The superscript $C$ indicates the cold particle population.

To calculate the Hermitian part $\overleftrightarrow{\mathbf{K}}^{h}$, and anti-Hermitian part $\overleftrightarrow{\mathbf{K}}^{a}$, which appear in Eqs. (1) and (2), we find the real and imaginary part of the dielectric tensor Eq. (6) based on the weak-growth rate approximation (i.e., $|\gamma_{k}|\ll|\omega_{k}|$). ¹¹1Vocks, C. & Mann, G. (2006) discuss details about the real and imaginary parts of the dielectric tensor. To do so, we use the Sokhotski-Plemelj relation for the $p_{\parallel}$-integral which faces the singularity problem in Eq. (6). However, the $p_{\parallel}$-singularity in Eq. (6), representing the resonant parallel momentum, is not yet in an explicit form because of the dependence of $\gamma$ on $p_{\parallel}$. Therefore, we rearrange the integrand of Eq. (6) as

We define the index of refraction in the direction parallel with respect to $\mathbf{B}_{0}$ as

and, in the present paper, we only consider wave modes with

Then, the $p_{\parallel}$-singularities of the right-hand side of Eq. (21) are

where

A careful analysis for Eq. (24) with Eq. (25) verifies that

and

Based on the assumption that $|\gamma_{k}|\ll|\omega_{k}|$, the Sokhotski-Plemelj relation for the $p_{\parallel}$-integral, therefore, approximates Eq. (6) as

where

and

At the second term of the right-hand side of Eq. (28), the integral with a dash indicates the principal value integral, defined as

Note that, after applying the Sokhotski-Plemelj relation, $\omega$ becomes $\omega_{k}$ in all following equations. By applying Eq. (28) to Eqs. (3) and (4), we then find the Hermitian part $\overleftrightarrow{\mathbf{K}}^{h}$ and anti-Hermitian part $\overleftrightarrow{\mathbf{K}}^{a}$, given as

and

where

Note that $\hat{G}\equiv\hat{H}$ at $p_{\parallel}=p_{\parallel res}^{+}$.

We calculate the energy density of the resonant wave Eq. (1) by using Eq. (32). In calculating the second term of Eq. (1), we obtain

where $\overleftrightarrow{\mathbf{U}}=\overleftrightarrow{\mathbf{T}}\!f_{j}/(\omega_{k}\gamma\Lambda_{\alpha})$ and $\overleftrightarrow{\mathbf{U}}_{res}\equiv\overleftrightarrow{\mathbf{U}}\big{|}_{p_{\parallel}=p_{\parallel res}^{+}}$. Eq. (35) is obtained through the definition of the principal value integral given as Eq. (31), Leibniz integral rule, and the following relation

Then, by using Eq. (35), Eq. (1) is expressed as

The first term of Eq. (37) is the pure energy density of the wave fields, expressed as

The second term of Eq. (37) is the energy density of the wave fields associated with the kinetic energy of particles through Landau and cyclotron resonances, expressed as

where

and $(\hat{\Pi}f_{j})_{res}\equiv\hat{\Pi}f_{j}\big{|}_{p_{\parallel}=p_{\parallel res}^{+}}$. Re indicates the real part of a complex function. We denote the right- and left-circularly polarised components of the electric field as $\mathit{{E}}_{k}^{R}\equiv(\mathit{{E}}_{k}^{x}-i\mathit{{E}}_{k}^{y})/\sqrt{2}$ and $\mathit{{E}}_{k}^{L}\equiv(\mathit{{E}}_{k}^{x}+i\mathit{{E}}_{k}^{y})/\sqrt{2}$. The longitudinal component of the electric field is $\mathit{{E}}_{k}^{z}$. Note that, if the plasma is cold (i.e., if $f_{j}=\overline{n}_{j}^{C}f_{j}^{C}$ where $\overline{n}_{j}^{C}$ is the relative density of the cold particle population of species $j$ in total particles), Eq. (39) reduces to

Finally, by substituting Eq. (33) into Eq. (2), the growth rate for a resonant wave in relativistic plasmas is expressed as

where

In the resonance with the wave mode of $\eta_{\parallel}^{2}>1$, a non-relativistic limit is achievable by setting $c\rightarrow\infty$ in Eqs. (39) and (48). By doing so, Eq. (48) reduces to the growth rate given by Kennel & Wong (1967). Note that the only approximation, which the derivation of Eqs. (39) and (48) is based on, is that $|\gamma_{k}|\ll|\omega_{k}|$.

Earth’s inner magnetosphere includes a relatively cold and dense ($\sim 1$ eV) electron population originating from the Earth’s ionosphere, alongside more sparse, but equally important, energetic electron populations ($>1$ keV) which are sourced from the magnetotail plasma sheet and further processed by wave-particle interactions in the outer radiation belt (Borovsky & Valdivia, 2018). These energetic electron populations have been studied statistically from geosynchronous orbit (Denton et al., 2010), demonstrating that there are typically two energetic electron populations: a hotter population with a number density of $\sim 5\times 10^{-4}$ cm^-3 and a temperature of $\sim 150$ keV, and a less-hot population with a number density of $\sim 10^{-2}$ cm^-3 and a temperature of $\sim 30$ keV. If these energetic electron populations exhibit the loss-cone MDF, they are unstable to whistler-mode wave in the inner magnetosphere (Ratcliffe & Watt, 2017; Li et al., 2019).

In our numerical calculations for Eqs. (39) and (48), we first assume that dominant cold electron and proton populations, and a minor relativistic electron population exist in the radiation belt. In the plasma-rest frame, for the cold populations, we use Eq. (20), while, for the relativistic electron population, we use the loss-cone MDF typically expressed in the plasma-rest frame as (Xiao et al., 1998)

where

$v_{th,e}\equiv\sqrt{2k_{B}T_{e}^{R}/m_{0e}}$ and $\Gamma(x)$ is the $\Gamma$-function. The superscript $R$ indicates the relativistic particle population. The free parameter $q$ and $s$ determine the anisotropy and steepness of the loss-cone MDF, $T_{e}^{R}$ is the temperature of the relativistic electron population, and $k_{B}$ is the Boltzmann constant. For our calculation, we then express the proton and electron MDFs of a quasi-neutral plasma as

and

where $\overline{n}_{e}^{C}$ and $\overline{n}_{e}^{R}$ indicate the relative density of the cold and relativistic electron population in total electrons, and $\overline{n}_{e}^{C}+\overline{n}_{e}^{R}=1$. The subscript $p$ and $e$ indicate the proton and electron.

In our calculation, we set the plasma parameters, which are representative in the radiation belt of the Earth, notwithstanding the wide range of natural variation (Borovsky & Denton, 2011; Watt et al., 2019). We calculate two different cases for the relativistic electron population, set based on the in-situ data (Denton et al., 2010; Gao et al., 2014) and defined as

Case 1, 1

Soft population: $T_{e}^{R}=25$ keV and $\overline{n}_{e}^{R}=0.05$,

Case 2, 2

Hard population: $T_{e}^{R}=140$ keV and $\overline{n}_{e}^{R}=0.0025$.

Fig. 1 illustrates the normalised $f_{e}^{R}$, which are the initially unstable loss-cone MDFs for the relativistic electron populations and plotted through Eq. (50) with our plasma parameters of the (a) soft and (b) hard populations. In Eq. (50), we arbitrarily set to $q=4$ and $s=2$ for our calculation due to the lack of in-situ data, which need to be determined by fitting with in-situ data at future studies. We set $\omega_{0pe}=3|\Omega_{0e}|$ (Horne et al., 2013).

To calculate Eqs. (39) and (48), we preferentially find the frequency and polarisation properties of the whistler-mode wave as a function of $k_{\parallel}$ by numerically solving a dispersion relation from Eq. (19). We confirm that the effect of our relativistic electron populations on the solution of a full relativistic-plasma dispersion relation is very subtle due to its small density (i.e., $5\%$ of soft population or $0.25\%$ of hard population in total electrons) compared to the cold plasma. Thus, to reduce the computational amount and time in finding its solution, we only apply $\overline{n}_{e}^{C}f_{e}^{C}$ of Eq. (53) to Eq. (19), which means that we numerically solve the cold-plasma dispersion relation in our calculation. ²²2In case that the effect of the relativistic electron population is not negligible, a full relativistic-plasma dispersion relation must be solved. To do so, the ALPS (Arbitrary Linear Plasma Solver) is an appropriate tool, which is an advanced and parallelised code numerically solving a full dispersion relation in relativistic regime (Verscharen et al., 2018).

Without loss of generality, we set $\omega_{k}>0$ in this calculation. We set $\phi_{k}=0^{\circ}$ because Eqs. (52) and (53) are gyrotropic, and denote the angle of the wave propagation with respect to $\mathbf{B}_{0}$ as $\theta$ (i.e., $k_{x}=k_{\perp}=k_{\parallel}\tan\theta$ and $k_{y}=0$). We choose our coordinate system so that the whistler-mode waves with an arbitrary angle $\theta$ have $k_{\parallel}>0$. We then find $\omega_{k}\equiv\omega_{k}(k_{\parallel})$, which makes the determinant of the $3\times 3$ matrix in Eq. (19) zero, in cold-plasma regime and with a given $\theta$. In doing so, we use the fsolve function in our Python code, which finds the root of the non-linear equation. Once we find $\omega_{k}\equiv\omega_{k}(k_{\parallel})$, we use it to find the polarisation properties of the whistler-mode wave from Eq. (19) as the ratios $|\mathit{{E}}_{k}^{y}|/|\mathit{{E}}_{k}^{x}|$ and $|\mathit{{E}}_{k}^{z}|/|\mathit{{E}}_{k}^{x}|$, and as a function of $k_{\parallel}$. In this way, because we are at liberty to set $E_{k}^{x}$, we set $E_{k}^{x}$ as a constant.

By using the profiles of $\omega_{k}$, $|\mathit{{E}}_{k}^{y}|/|\mathit{{E}}_{k}^{x}|$ and $|\mathit{{E}}_{k}^{z}|/|\mathit{{E}}_{k}^{x}|$, we calculate Eqs. (39) and (48) as a function of $k_{\parallel}$ based on plasma parameters of soft and hard populations, respectively. We numerically evaluate the principal value integral in Eq. (39), defined as Eq. (31), through the following relation,

where $h(X)$ is an arbitrary function of the variable $X$, and $X_{s}$ indicates a singularity. We assume that the integration of the first term on the right-hand side of Eq. (54) converges. We normalise $p_{\perp}$ and $p_{\parallel}$ by $m_{0e}c$, and frequency by $|\Omega_{0e}|$. For numerical reasons in our calculations, we set the finite ranges of $p_{\perp}/m_{0e}c$ and $p_{\parallel}/m_{0e}c$ as $0.01<p_{\perp}/m_{0e}c\leq 3.8$ and $-3.8\leq p_{\parallel}/m_{0e}c\leq 3.8$ for integrals in Eqs. (39) and (48), confirming that these ranges are wide enough to make the integration converge. We numerically calculate the $p_{\parallel}$- and $p_{\perp}$-integration based on the Simpson’s $3/8$ rule. We set the range of $ck_{\parallel}/|\Omega_{0e}|$ as $0.3<ck_{\parallel}/|\Omega_{0e}|\leq 10$, and the step size of $ck_{\parallel}/|\Omega_{0e}|$ is $0.065$ in our calculation.

Fig. 2 shows the profiles of $\omega_{k}$, $|\mathit{{E}}_{k}^{y}|/|\mathit{{E}}_{k}^{x}|$, and $|\mathit{{E}}_{k}^{z}|/|\mathit{{E}}_{k}^{x}|$ as a function of $k_{\parallel}$ with different propagation angles, ranging from $\theta=0.01^{\circ}$ to $\theta=89^{\circ}$. In each figure, the solid curves correspond to the numerical results from the cold-plasma dispersion relation. With the assumptions that the proton contribution is neglected, $\omega_{k}^{2}\ll\omega_{0pe}^{2}$ and $\Omega_{0e}^{2}\ll\omega_{0pe}^{2}$, the cold-plasma dispersion relation also provides the analytic expressions for $\omega_{k}$, $|\mathit{{E}}_{k}^{y}|/|\mathit{{E}}_{k}^{x}|$, and $|\mathit{{E}}_{k}^{z}|/|\mathit{{E}}_{k}^{x}|$ of the whistler-mode wave as (Stix, 1992)

and

The dashed curves in Fig. 2 with the $\theta$-colour scheme same as the $\theta$-colour scheme of solid curves correspond to the analytical results from Eqs. (55) to (57). The numerical results for $\omega_{k}$, $|\mathit{{E}}_{k}^{y}|/|\mathit{{E}}_{k}^{x}|$ and $|\mathit{{E}}_{k}^{z}|/|\mathit{{E}}_{k}^{x}|$ closely follow the analytical results. This verifies that our numerical solutions of the cold-plasma dispersion relation are correct.

By using our numerical results for $\omega_{k}$, $|\mathit{{E}}_{k}^{y}|/|\mathit{{E}}_{k}^{x}|$, and $|\mathit{{E}}_{k}^{z}|/|\mathit{{E}}_{k}^{x}|$, we calculate Eq. (48) as a function of $k_{\parallel}$ with different $\theta$ based on Eq. (53). Note that the cold-particle MDF provides no contribution to Eq. (48), which means that Eq. (48) is only estimated at $f_{e}=\overline{n}_{e}^{R}f_{e}^{R}$, illustrated in Fig. 1. Fig. 3 shows the profiles of $\gamma_{k}$, $\gamma_{k,e}^{n=+1}$, $\gamma_{k,e}^{n=-1}$ and $\gamma_{k,e}^{n=0}$ estimated at $f_{e}=\overline{n}_{e}^{R}f_{e}^{R}$ for the soft population of relativistic electrons, ranging from $\theta=0.01^{\circ}$ to $\theta=89^{\circ}$. In each plot, the solid black curve, and dashed red, blue, and green curves correspond to $\gamma_{k}$, $\gamma_{k,e}^{n=+1}$, $\gamma_{k,e}^{n=-1}$ and $\gamma_{k,e}^{n=0}$, respectively. In Fig. 3, $\gamma_{k,e}^{n=-1}$ only contributes to the wave instability (i.e., the positive contribution to $\gamma_{k}$). On the other hand, $\gamma_{k,e}^{n=+1}$ and $\gamma_{k,e}^{n=0}$ contribute to the wave damping (i.e., the negative contribution to $\gamma_{k}$), and counteract the driving of the wave instability. Therefore, the $n=-1$ resonance must dominate over the $n=+1$ and $n=0$ resonances, so that the loss-cone driven instability occurs. When $\theta=0^{\circ}$, the loss-cone driven instability is the most efficient. But, when increasing $\theta$, it becomes inefficient due to an increase of the wave damping through the $n=0$ resonance, and, from around $\theta=30^{\circ}$, the loss-cone driven instability disappears.