3.1 Conversion of kinetic energy

Before discussing the generation of Langmuir and electron cyclotron waves, let us consider the temporal variation of the kinetic energy of the beam electrons (see figure 6$a1$–$a3$) and of the background electrons (see figure 6$b1$–$b3$). Here the total kinetic energy is calculated as $\mathcal {E}_t=\frac {1}{2}m_e\overline {v^2}\cdot N_e\cdot M$, the bulk flow kinetic energy as $\mathcal {E}_b=\frac {1}{2}m_e\bar {v}^2\cdot N_e\cdot M$ and the electron thermal energy (heat) as $\mathcal {E}_{\textrm {th}}=\mathcal {E}_t-\mathcal {E}_b$. The averages are over all the particles, e.g. $\bar {v}=(\sum _i v_i)/N_e$, where $v_i=\sqrt {v_{ix}^2 + v_{iy}^2 + v_{iz}^2}$ is the speed of the $i$th particle. Here $N_e$ is the number of macro-electrons of the beam or background plasma and $M$ is the ratio of physical to numerical particles as defined before. The total electric field energy is calculated as $\mathcal {E}_{E}=1/(8{\rm \pi} )\cdot \sum _{i=1}^{N_x} \sum _{j=1}^{N_y} E_{ij}^2\Delta V$, where $\Delta V=\Delta x\cdot \Delta x$ is the cell area (or ‘volume’) for a two-dimensional simulation. The total magnetic energy is calculated in a similar fashion by replacing electric field $E$ with magnetic field $B$.

Figure 6.

Temporal energy variations for all simulations. ($a1$–$a3$) Variation of total kinetic energy, bulk flow energy and thermal energy of beam electrons, respectively. ($b1$–$b3$) The same for background electrons. ($c1$–$c3$) Variations of the total kinetic energy, electric energy and magnetic energy of the beam–plasma system. All the energy variations are normalized by the initial bulk flow kinetic energy of the electron beam $\mathcal {E}_0$ of Run1.

In order to compare the electron energy gain or loss, we evaluate the normalized variation of each kind of energy as $\Delta \mathcal {E}(t)/\mathcal {E}_0=(\mathcal {E}(t)-\mathcal {E}(t=0))/\mathcal {E}_0$. All energy variations are normalized by the corresponding initial bulk flow kinetic energy of the electron beam $\mathcal {E}_0=\frac {1}{2}m_eu_{d\parallel }^2\cdot N_{e,\textrm {bm}}\cdot M$. Since the number of macro-electrons of the beam in all PIC code simulations is unchanged, and the electron beams are initialized with the same parallel drift speed $u_{d\parallel }$, as a result, the initial energy $\mathcal {E}_0$ is the same for all simulations. Figure 6($c1$–$c3$) also shows evolutions of the total kinetic energy and the electric and magnetic energy fluctuations of beam–plasma systems. The results displayed in figure 6 show that:

1. (i) most of the bulk flow kinetic energy of the beam electrons is transferred into thermal energy of the background electrons while only a small part of it is converted into electric field energy;

2. (ii) the background electrons are heated while their bulk kinetic energy remains nearly zero; and

3. (iii) the electric energy increases at the expense of the bulk flow kinetic energy of the beam electrons.

Comparing the energy variations obtained for different parameter regimes displayed in figure 6, we found that the density gradient of the background plasma significantly influences the energy conversion process for both Maxwellian and ring beams: the larger the density gradient $\eta$, the more efficiently is the kinetic energy converted into other forms of energy. Comparing energy variations for different cases with the same density gradient $\eta$, e.g. Run1 and Run2, or Run3 and Run4, or Run5 and Run6, we found that when $\eta \le 50$ the energy conversion of the Maxwellian-beam case is more efficient than that of the ring-beam case. When $\eta > 50$, however, the energy conversion in ring-beam cases is more efficient than that in Maxwellian-beam cases, e.g. Run7 and Run8.

3.2 Langmuir waves

Our simulations revealed the generation of electrostatic Langmuir waves as well as of electromagnetic electron cyclotron waves. Let us first discuss the formation of Langmuir waves and their consequences for possible radio wave emissions.

Langmuir waves are generated by the bump-on-tail instabilities for both Maxwellian- beam (Runs 1, 3, 5, 7) and ring-beam (Runs 2, 4, 6, 8) cases. The results of our simulations show that the bump-on-tail EVDFs of Maxwellian- and ring-beam plasma systems with the same density gradient have nearly the same influence on the generation of Langmuir waves. In this part, only results associated with Maxwellian beams are discussed in detail, because the same results also hold for the ring-beam cases.

How do the electron populations relax in the parallel velocity space? Figure 7 depicts the temporal evolution of the bump-on-tail EVDF $f(v_{\parallel })$ of Maxwellian beam plasma cases (i.e. Runs 1, 3, 5, 7) with increasing density gradients (i.e. $\eta =1, 5, 50, 200$) at $t=0, 100, 200$ and $275 \omega _{\textrm {pe}}^{-1}$, respectively. Since the beam density is much smaller than the background density, we normalize the one-dimensional EVDFs of background and beam independently in order to visualize them in the same plot. Note that $t=275 \omega _{\textrm {pe}}^{-1}$ is beyond the maximum time $t=200 \omega _{\textrm {pe}}^{-1}$ of the duration in which we analysed plasma waves.

Figure 7.

Temporal evolution of the parallel EVDF for Runs 1, 3, 5, 7 at time $t=0$, $100\omega _{\textrm {pe}}^{-1}$, $200\omega _{\textrm {pe}}^{-1}$ and $275 \omega _{\textrm {pe}}^{-1}$. The background (centred at $v_{\parallel }=0$) and beam (centred at $u_{d\parallel }$) EVDFs are probability density functions but independently normalized for visualization purposes, so that their respective integrations in the whole $v_{\parallel }$ space yield 1.

The beam of Run7 (red curves in figure 7), for example, is relaxing to form a characteristic plateau between $t=100 \omega _{\textrm {pe}}^{-1}$ and $200 \omega _{\textrm {pe}}^{-1}$, while the plateau is already formed around $t=200\text {--}275 \omega _{\textrm {pe}}^{-1}$. As the density gradient of background plasma increases (from Run1 to Run 7), this beam relaxation process occurs faster. This behaviour is correlated to the more efficient transfer of the bulk flow kinetic energy of the electron beam into energy of plasma waves (see figure 6). Note that the beam EVDF has not yet fully relaxed in each simulation while harmonics of Langmuir waves are generated even at earlier stages of the beam relaxation (see figure 8). This means that the free energy released from the electron beam causes beam–plasma instabilities, thus generating Langmuir waves and their harmonics from the very beginning of the beam relaxation. For example, figure 8 shows that harmonics of Langmuir waves are generated between $t=100 \omega _{\textrm {pe}}^{-1}$ and $164 \omega _{\textrm {pe}}^{-1}$.

Figure 8.

Illustration of fundamental (F) and harmonics (H) of Langmuir wave (LW) modes. The PSDs in the $k_{\parallel }$–$\omega$ plane for Runs 1, 3, 5, 7 (Maxwellian beam) during $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$ are displayed. The PSD is evaluated as $\log _{10}|E_{\iota }/B_0|^2$. Dispersion relation curves of beam mode (red dashed line), fundamental (blue dashed curve) and harmonics (magenta dashed curves) of Langmuir waves are overlaid.

Note that a full beam relaxation, in particular for a small relative density gradient, would require a much larger simulation box. This would be computationally very demanding. In more precise terms, let us consider the relative gradient $G=(n_{\max }-n_{\min })/(L_x/2)=n_{\max }(1-\eta ^{-1})/(L_x/2)$, where $L_x$ is the simulation box size along the $x$ direction, so that the denominator is the distance between the location of the maximum density at the left-hand boundary and the minimum density at the centre of the simulation box. Note that the quantity $\eta =n_{\max }/n_{\min }$ in the numerator is the density gradient defined before. In general, in order to keep $G$ constant, a larger simulation domain $L_x$ would require a larger $\eta$. However, in order to avoid numerical artefacts, like numerical heating and strong electric field fluctuations due to charge separation, the minimum number of particles per cell $N_{\textrm {bg}}\cdot \eta ^{-1}$ should not be smaller than approximately 5–10 particles per cell, assuming a constant macro-factor or weight (ratio of physical to numerical particles) for each macro-particle. This situation could be partially avoided by using a variable macro-factor so that regions with low physical density can be represented by more numerical particles, but we do not utilize such a feature in the present simulations. On the other hand, in our simulations the maximum number of particles per cell is constrained to be smaller than 1000 due to computational reasons. Thus the upper bound of our density gradient is $\eta \le N_{\textrm {bg}}/5\approx 200$. Therefore, if we want to have the same gradient $G$ but we are constrained by an upper bound on $\eta$, we can only choose a maximum simulation domain size $L_x$ for the given $\eta$. A larger simulation domain size $L_x$ for a given $\eta$ would imply a smaller gradient $G$.

In the following we will only concentrate on the electric field $\boldsymbol {E}$. We refer to $E_{\iota }$, $E_{\tau }$ and $E$ as the longitudinal, transverse and total electric field in Fourier space, i.e.

(3.1)\begin{equation} \left. \begin{gathered} E_{\iota}=E_x,\\ E_{\tau}=\sqrt{E_y^2+E_z^2},\\ E=\sqrt{E_x^2+E_y^2+E_z^2}. \end{gathered} \right\} \end{equation}

The PSD of each electric field component is simply evaluated by squaring each of those components and normalizing appropriately.

For a time window spanning $64 \omega _{\textrm {pe}}^{-1}$, the variation of the local plasma frequency along the beam trajectory can be neglected. This implies that the the effective plasma frequency can be treated as a constant when performing the fast Fourier transform on the electromagnetic fields. Figure 8 shows the resulting PSD derived from longitudinal electric field $E_{\iota }$ for Runs 1, 3, 5 and 7 (Maxwellian beams) in the $(k_{\parallel },\omega )$ domain during $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$ as well as the corresponding analytical dispersion curves of fundamental Langmuir waves and up to their sixth harmonic. A PSD enhancement along these dispersion relation curves indicates that those harmonics of Langmuir waves are indeed generated.

The dispersion relation of the fundamental Langmuir waves $L$ can be expressed by the following expression:

(3.2)\begin{equation} \omega_{L}=\sqrt{\omega_{\textrm{loc}}^2+3v_{\textrm{the}}^2k^2}, \end{equation}

where $\omega _{\textrm {loc}}$ is the effective local plasma frequency, which is dependent on the location of the electron beam as it propagates through the background plasma, and $v_{\textrm {the}}$ is the electron thermal speed of background plasma. Note that in a homogeneous background plasma the local plasma frequency remains constant, i.e. $\omega _{\textrm {loc}}=\omega _{\textrm {pe}}$, and then (3.2) becomes the standard Bohm–Gross dispersion relation of fundamental Langmuir waves (Yoon et al. Reference Yoon, Gaelzer, Umeda, Omura and Matsumoto2003; Melrose Reference Melrose2017).

For the higher harmonic of Langmuir waves, we generalized (3.2) to the dispersion relation in the following form:

(3.3)\begin{equation} \omega_{L_n}(k)=\sqrt{n^2\omega_{\textrm{loc}}^2+3v_{\textrm{the}}^2k^2}. \end{equation}

Here $L_n (n\ge 2)$ denotes the $n$th harmonic of Langmuir waves. This is just an empirical relation that fits well with harmonics of the Langmuir waves in our simulations.

In order to compare our formula (3.3) with previous works, we can take its long-wavelength limit, i.e. $\omega _{\textrm {loc}}^2\gg 3v_{\textrm {the}}^2k^2$, so that

(3.4)\begin{equation} \omega_{L_n}(k)\approx n\omega_{\textrm{loc}}\left(1+\frac{3}{2n^2}\lambda_{D}^2k^2\right). \end{equation}

This is similar to but not the same as previously reported results (e.g. Yoon Reference Yoon2000; Gaelzer et al. Reference Gaelzer, Ziebell and Yoon2002; Yoon et al. Reference Yoon, Gaelzer, Umeda, Omura and Matsumoto2003). For the second harmonic ($n=2$), (3.3) as well as (3.4) are equivalent to (34) of Yoon (Reference Yoon2000), while for higher harmonics of Langmuir waves, i.e. $n>2$, (3.3) and (3.4) are more accurate due to a factor $1/n$ in their second term within the parentheses on the right-hand side of (3.4) (comparing also to (18) in Yoon et al. (Reference Yoon, Gaelzer, Umeda, Omura and Matsumoto2003)). This discrepancy with previous theoretical works is an indication of a different physical process causing those harmonic waves in our simulations. A theoretical investigation of the dispersion relation of harmonics of Langmuir waves based on kinetic theory is needed but beyond the scope of this study. Note that (3.4) is still valid for $\omega _{\textrm {loc}}^2\approx 3v_{\textrm {the}}^2k^2$, as confirmed by our simulations.

As the electron beam propagates through the background plasma, the local plasma frequency changes due to the density gradient. Table 3 shows an algorithm to estimate the effective local plasma frequency $\omega _{\textrm {loc}}$ based on the power spectrum derived from the longitudinal electric field $E_{\iota }$. Once the effective local plasma frequency is determined, the dispersion relation can be solved for the harmonics of Langmuir waves by using (3.3).

Table 3.

Algorithm to determine the effective local plasma frequency $\omega _{\textrm {loc}}$.

Figure 8 shows that the steeper the background density gradient (larger $\eta$), the greater the number of Langmuir harmonics that appear. For Run1 ($\eta =1$; see figure 8$a1$) only three harmonics are generated, while for Run5 ($\eta =50$; see figure 8$b1$) at least six harmonics are generated. As figure 8 shows, the wavenumber $k_{\parallel }$ of the power spectrum peak of each harmonic mode also increases with the harmonic order $n$.

In order to quantitatively analyse the frequencies of harmonic Langmuir waves, we calculate the power spectrum in the frequency domain by integrating the PSD of the electrostatic (longitudinal) electric field over $k_{\parallel }$, i.e.

(3.5)\begin{equation} \mathcal{P}(\omega)=\int |E_{\iota}\left(\omega,k_{{\parallel}},k_{{\perp}}=0\right)|^2 \,\textrm{d}k_{{\parallel}}. \end{equation}

Figure 9 shows the power spectra $\mathcal {P}(\omega )$ of all simulations during the time period $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$. We conclude that:

1. (i) For beam–plasma systems with the same background density gradient $\eta$, we observe practically the same power spectra, in particular at the characteristic frequencies of Langmuir waves and their harmonics. This is because we initialized both Maxwellian and ring beams with the same parallel drift speed $u_{d\parallel }$ and thus the same bump-on-tail EVDF.

2. (ii) For beam–plasma systems with different background density gradients, as density gradients $\eta$ increase, the characteristic frequencies of fundamental and harmonics of Langmuir waves shift to lower frequencies. This is because for larger density gradients $\eta$, the effective local plasma frequency at a beam location becomes smaller.

3. (iii) As $\eta$ increases, more harmonics of Langmuir waves are generated. For example, for Maxwellian-beam plasma system Run1 ($\eta =1:F, H_2, H_3$; see figure 9$a$), Run3 ($\eta =5:F, H_2, \ldots , H_6$; see figure 9$b$), Run5 and Run7 ($\eta =50$ and $200: F, H_2, \ldots , H_7$; see figure 9$c$,$d$).

4. (iv) The local maximum power of each harmonic of Langmuir waves significantly decreases as the harmonic number $n$ increases in each beam–plasma system.

5. (v) Nearly no difference is found between the power spectra of the beam–plasma systems with $\eta =50$ and $\eta =200$ (i.e. Runs 5, 6, 7, 8; see figure 9$c$,$d$). This means that when $\eta \ge 50$, background density gradients affect the Langmuir waves to nearly the same extent. Figure 9 shows that characteristic frequencies of the third harmonic of Run5 and Run6 (see figure 9$c$) and of Run7 and Run8 (see figure 9$d$) slightly deviate from each other. This might be due to the frequency resolution $\Delta \omega$ of the fast Fourier transform. Such deviation is expected to disappear when performing fast Fourier transform on a longer time series allowing a higher frequency resolution.

Figure 9.

Power spectrum $\mathcal {P}(\omega )$ derived from longitudinal electric field $E_{\iota }$ for all simulations in $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$. The vertical dashed lines denote characteristic frequencies of the fundamental and harmonics of Langmuir waves. The unit of the spectrum is $\textrm {statV}^2~\textrm {cm}^{-2}~\textrm {s}$. Here ‘M’ indicates a Maxwellian beam, while ‘R’ denotes a ring beam.

The characteristic frequencies and wavenumbers of fundamental and harmonics of Langmuir waves in our simulations satisfy the wave–wave conditions, also known as beat or Manley–Rowe conditions (see e.g. (1) in Melrose Reference Melrose2017). Indeed, by virtue of the dispersion relation (3.4), one can prove the beat condition of the frequencies $\omega _{L_{n-1}}+\omega _{L} \approx \omega _{L_n}$, where $n$ is the harmonic number. The wavenumber of each harmonic of Langmuir waves also satisfies the corresponding beat condition $k_{n}=nk_{1}$, where $k_n$ corresponds to the $k_{\parallel }$ location of the PSD peak along the dispersion relation curve of the $n$th harmonic of Langmuir waves. The method to determine such $k_n$ will be explained below. This provides evidence that the most likely generation mechanism of harmonics of Langmuir waves in our simulations is the coalescence of beam-generated (fundamental) Langmuir waves $L$ and adjacent Langmuir harmonic $L_{n-1}$, which leads to the $n$th harmonic of Langmuir waves, i.e.

(3.6)\begin{equation} L_{n-1}+L\to L_n. \end{equation}

Here $n\ge 2$. In particular, if $n=2$, (3.6) indicates that the second harmonic of Langmuir waves is caused by a coalescence of beam-generated Langmuir waves $L$. This mechanism was already predicted by previous theoretical studies (see e.g. Yoon et al. Reference Yoon, Gaelzer, Umeda, Omura and Matsumoto2003; Yi et al. Reference Yi, Yoon and Ryu2007; Rhee et al. Reference Rhee, Ryu, Woo, Kaang, Yi and Yoon2009). It is important to emphasize that the specific harmonic wave dispersion relation of those previous works is different from ours (see (3.3)), probably because the underlying physical processes are not exactly the same. Nevertheless, the harmonics of Langmuir waves observed in our simulations also satisfy the coalescence process (3.6).

To analyse the characteristic wavenumbers $k_n$ ($n=1,2,3,\ldots$) of harmonics of Langmuir waves, we calculate the power distribution $\mathcal {P}_i(k_{\parallel })$ of each harmonic based on (2.10) after appropriately ruling out the influence of the beam mode:

(3.7)\begin{equation} \mathcal{P}_i(k_{{\parallel}})=\int \frac{1}{\sqrt{2{\rm \pi}}\sigma}\exp\left[-\frac{(\omega-\omega_{L_i})^2}{2\sigma^2}\right]\cdot |E_l(k_{{\parallel}},k_{{\perp}}=0,\omega)|^2\,\textrm{d}\omega. \end{equation}

Here dispersion relation $\omega _{L_i}$ is determined by (3.3). Then $k_n$ of each Langmuir wave harmonic is determined by the wavenumber $k_{\parallel }$ of the local maximum of $\mathcal {P}(k_{\parallel })$.

Note that this characteristic wavenumber $k_n$ is time-dependent. Figure 10 shows the temporal evolution of the PSD of longitudinal electric field $E_{\iota }$ of Run5 ($\eta =50$) at time periods $0\text {--}64, 50\text {--}114, 100\text {--}164, 150\text {--}214 \omega _{\textrm {pe}}^{-1}$. As time evolves, the $k_{\parallel }$ location of the PSD peak along its dispersion relation curve of each Langmuir harmonic moves to a higher-wavenumber region.

Figure 10.

The PSD in the $k_{\parallel }$–$\omega$ plane for Run5 ($\eta =50$, Maxwellian beam) in four consecutive time periods $t=0\text {--}64, 50\text {--}114, 100\text {--}164, 150\text {--}214 \omega _{\textrm {pe}}^{-1}$. The PSD is evaluated by $\log _{10}|E_{\iota }/B_0|^2$. Dispersion relation curves of beam mode (red dashed lines), fundamental (blue dashed curves) and harmonics (magenta dashed curves) of Langmuir waves are overlaid.

Linear regression analysis of the characteristic wavenumbers $k_n$ of Langmuir modes is carried out at different time stages. Figure 11 shows the results of Runs 3, 4, 5, 6 during three consecutive time periods $t = 50\text {--}114, 100\text {--}164, 150\text {--}214 \omega _{\textrm {pe}}^{-1}$. We find that not only does the characteristic wavenumber of the (fundamental) Langmuir mode increase as the harmonic number increases, but also that of each Langmuir harmonic slightly increases over time. The latter is already confirmed in figure 10. It always agrees very well with the beat conditions of the wavenumber $k_n=nk_1$ and thus supports (3.6).

Figure 11.

Linear regression analysis of characteristic wavenumbers of harmonics of Langmuir waves for Runs 3, 4, 5 and 6 in three consecutive time periods $50\text {--}114 \omega _{\textrm {pe}}^{-1}$ (red dots), $100\text {--}164 \omega _{\textrm {pe}}^{-1}$ (blue squares) and $150\text {--}214 \omega _{\textrm {pe}}^{-1}$ (magenta triangles). The horizontal axis denotes the harmonic number $n$.

We also find evidence of waves near $k_{\parallel }\to 0$ in our simulations, i.e. for long wavelengths. Waves with $k_{\parallel }=0$ cannot of course take place in a numerical simulation with a finite-size simulation box, but it is nonetheless interesting to analyse the behaviour of the spectral power for those small $k$ wavenumbers. Figure 12 shows the power distribution of the fundamental Langmuir wave mode calculated by (3.7) for RunA, Run1 and Run3 at different time windows. Here RunA was carried out as a convergence test, where no electron beam is present (for details, see the Appendix). The power $\mathcal {P}(k_{\parallel })$ of the fundamental Langmuir mode for long wavelengths ($k_{\parallel }\to 0$) of Run1 (with beam, $\eta =1$) is several orders of magnitude larger than that of RunA (without beam), which is comparable to the noise level (see figure 12$a$). The power of the fundamental Langmuir mode of Run1 is consistently larger than the noise level of RunA for all time intervals. It shows a well-defined local maximum in the long-wavelength limit (see figure 12$b$). However, Run3 (with beam, $\eta =5$) has a clearly defined local maximum near $k_{\parallel }\to 0$ only during the initial time period (i.e. $0\text {--}64 \omega _{\textrm {pe}}^{-1}$). Indeed, figure 12$(c)$ shows that density inhomogeneities of background plasma inhibit these long-wavelength waves. The wavevector $k_F$ in figure 12 is the characteristic wavenumber of the beam-generated Langmuir mode, being determined by the location where the dispersion relation curves of the beam mode and of the fundamental Langmuir mode meet. In our simulations, $k_F\approx 2.18 \omega _{\textrm {pe}}/c$. The power distribution $\mathcal {P}(k_{\parallel })$ at $k_F$ is similar for both Run1 and Run3. The peaks of these power enhancements move slightly to larger $k_{\parallel }$ over time.

Figure 12.

Power distribution $\mathcal {P}(k_{\parallel })$ (in unit of $\textrm {statV}^2\cdot\textrm {cm}^{-1}$) of fundamental Langmuir mode for Runs A, 1 and 3 in different time windows $t=0\text {--}64, 50\text {--}114, 100\text {--}164, 150\text {--}214 \omega _{\textrm {pe}}^{-1}$. Here $k_F=2.18 \omega _{\textrm {pe}}/c$ (denoted by the vertical dashed line) is the characteristic wavenumber of the beam-generated Langmuir mode.

The enhancement of spectral power in that long-wavelength region may be due to waves with a wavenumber equal to or smaller than (i.e. unresolved) the spectral resolution, which is $k_{\parallel }=2{\rm \pi} /L_x=0.03 d_e^{-1}$. Another source for waves in this region can be due to the non-zero initial net current carried by the electron beam, causing a standing wave oscillating at the plasma frequency (i.e. plasma oscillations). This is in turn a direct consequence of the displacement current in Ampere's law and the periodic boundary conditions. Although it could be possible to suppress this initial net current by subtracting it at each time step of the simulation, the physical consequences of such a procedure are even more unclear in our plasma system with a localized beam.

3.3 Electron cyclotron waves

The EVDF of a ring beam can offer a perpendicular source of free energy (positive gradient in the perpendicular EVDF). This can cause ECMIs and generate electron cyclotron waves via the ECME mechanism. Simulation Runs 2, 4, 6, 8 address the nonlinear evolution of localized ring beams propagating in background plasmas with varying density gradients $\eta$. Figure 13 depicts the perpendicular EVDF for Runs 2, 4, 6, 8 at three different moments of time $t=0, 150$ and $275 \omega _{\textrm {pe}}^{-1}$ (see figure 13$a1$–$a3$). Results show that as time evolves the electron beams relax, i.e. the positive gradients in these perpendicular EVDFs are reduced, with the consequent transfer of kinetic energy to plasma heating and wave energy. This relaxation occurs faster as the background density gradient becomes steeper (e.g. compare green curves for $\eta =1$ with the red curves for $\eta =200$), indicating a more efficient energy transfer. Figure 13($b1$–$b3$) shows the time evolution of the two-dimensional EVDF in the $v_{\perp 1}$–$v_{\perp 2}$ plane for Run6 ($\eta = 50$). The two-dimensional EVDFs clearly show the characteristic broadening of the ring beam over time, although still far from its complete relaxation. A considerable part of the free energy of the ring beams is transferred to growing plasma waves.

Figure 13.

($a1$–$a3$) Temporal evolution of perpendicular EVDF for Runs 2, 4, 6, 8 (ring beam) at time: $t=0, 150, 275 \omega _{\textrm {pe}}^{-1}$, respectively. ($b1$–$b3$) Normalized two-dimensional beam EVDF in the $v_{\perp 1}$–$v_{\perp 2}$ plane for Run6 at the same three moments of time.

Figure 14 depicts the PSD of the transverse electric field in the ($\omega , k_{\perp }$) domain during $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$ for simulation Runs 1 and 2 ($\eta =1$) and 5 and 6 ($\eta =200$). Note that Runs 1 and 5 are carried out with a Maxwellian beam and Runs 2 and 6 with a ring beam, so that we can compare the effect of perpendicular sources of free energy (positive gradients in the EVDFs). Figure 14 shows that fundamental and harmonics of electron cyclotron waves due to ECRs are observed in all the simulations (e.g. Runs 1, 2, 5, 6; see red dashed lines in figure 14). On the other hand, fundamental and harmonics of electron cyclotron waves due to ECMIs are only observed in simulations with a ring beam (e.g. Runs 2, 6; see blue dashed lines in figure 14$b1$,$b2$). The dispersion relations of various electron cyclotron waves are discussed below. Results show that those waves, including high-order harmonics, are already generated at the early stage of the instability growth (e.g. $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$) and well before the relaxation of the velocity gradient in perpendicular EVDF starts (see figure 13). By comparing results between Runs 1 and 5 and between Runs 2 and 6, we notice that the background density gradients practically have no influence on the generation of electron cyclotron waves.

Figure 14.

Illustration of fundamental (F) and harmonics (H) of electron cyclotron wave (ECW) modes. The PSD in the $k_{\perp }$–$\omega$ plane for Runs 1, 2, 5, 6 in the time window $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$ are displayed. The PSD is evaluated by $\log _{10}|E_{\tau }/B_0|^2$. Dispersion relation curves of fundamental and harmonics of electron cyclotron waves due to ECRs (denoted by red dashed lines) and ECMIs (denoted by blue dashed lines), X (green dashed curve) and Z (magenta dashed curve) modes are overlaid.

As mentioned above, two mechanisms contribute to the generation of the fundamental and harmonics of electron cyclotron waves. The first one is due to the ECRs at $\omega =n\varOmega _{\textrm {ce}}\ (n=1,2,3,\ldots )$ (see red dashed lines in figure 14) and the other is attributed to the ECMIs caused by the ring beam (see blue dashed lines in figure 14$b1$,$b2$). We discuss now their corresponding dispersion relations.

The dispersion relation of the well-known electron cyclotron Bernstein modes is expressed in the following form (see e.g. Melrose Reference Melrose1986):

(3.8)\begin{equation} 1-\frac{2\omega_{\textrm{pe}}^2}{\lambda_e}\textrm{e}^{-\lambda_e}\sum_{n=1}^{\infty}\frac{n^2\textrm{I}_n(\lambda_e)}{\omega^2-n^2\varOmega_{\textrm{ce}}^2}=0, \end{equation}

where $\textrm {I}_n(\lambda )$ is the modified Bessel function of the first kind with argument $\lambda _e=k^2v_{\textrm {the}}^2/\varOmega _{\textrm {ce}}^2$ and $v_{\textrm {the}}$ is the electron thermal speed of the background plasma. In the limit $\lambda _e\ll 1$, $\textrm {e}^{-\lambda _e}\rightarrow 1$ and the asymptotic expression of the modified Bessel function $\textrm {I}_n(\lambda )$ becomes $\textrm {I}_n(\lambda _e)\approx ({1}/{n!})({\lambda _e}/{2})^n$. Using this approximation of $\textrm {I}_n(\lambda _e)$ in (3.8) we can get approximate solutions for the frequency $\omega \approx \sqrt {\omega _{\textrm {pe}}^2+\varOmega _{\textrm {ce}}^2}$, $2\varOmega _{\textrm {ce}}$, $3\varOmega _{\textrm {ce}},\ldots$ While in the opposite limit of $\lambda _e\gg 1$, the modified Bessel function $\textrm {I}_n(\lambda )$ can be approximated as $\textrm {e}^{-\lambda _e}\textrm {I}_n(\lambda _e)\approx ({1}/{\sqrt {2{\rm \pi} \lambda _e}})\textrm {e}^{-n^2/2\lambda _e}$. By solving (3.8) for the frequency with this approximated expression, we get solutions for the harmonics of the electron cyclotron frequency as $\omega =n\varOmega _{\textrm {ce}}$, $n=1,2,3,\ldots$

The dispersion relation of electron cyclotron wave modes by ECMIs is obtained by a small modification related to relativistic effects in (3.8), yielding

(3.9)\begin{equation} 1-\frac{2\omega_{\textrm{pe}}^2}{\lambda_e^2\gamma_{{\perp}}^2}\textrm{e}^{-\lambda_e^2\gamma_{{\perp}}^2}\sum_{n=1}^{\infty}\frac{n^2\textrm{I}_n(\lambda_e^2\gamma_{{\perp}}^2)}{\omega^2-n^2\varOmega_{\textrm{ce}}^2/\gamma_{{\perp}}^2}=0. \end{equation}

The reason is that those electron cyclotron wave modes are associated with velocity gradient in the perpendicular EVDF of the ring beam. Relativistic effects related to the perpendicular drift speed of the beam are essential for this maser mechanism (Melrose Reference Melrose2017). Specifically, the difference between the standard (3.8) and the electron cyclotron maser counterpart (3.9) has to do with the the perpendicular Lorentz factor $\gamma _{\perp }$ derived from the perpendicular drift speed of the ring beam, i.e.

(3.10)\begin{equation} \gamma_{{\perp}}=\sqrt{1+u_{d\perp}^2/c^2}. \end{equation}

Accordingly, the electron cyclotron frequency term $\varOmega _{\textrm {ce}}$ in (3.8) has to be replaced by $\varOmega _{\textrm {ce}}/\gamma _{\perp }$, consequently, to get (3.9). The consequence of this difference is a downward shift of the resonance frequencies, i.e. the ECMI-related dispersion curves (see blue dashed lines in figure 14$b1$,$b2$) are located always below the ECR-related dispersion curves (see red dashed lines in figure 14$b1$,$b2$).

Note that a numerical aliasing effect appears in the form of line-like PSD enhancements in figure 14($a1$,$a2$) between frequencies $\omega =3\varOmega _{\textrm {ce}}$ and $4 \varOmega _{\textrm {ce}}$, between $4\varOmega _{\textrm {ce}}$ and $5 \varOmega _{\textrm {ce}}$ as well as between $5\varOmega _{\textrm {ce}}$ and $6 \varOmega _{\textrm {ce}}$. This is an artefact due to our maximum sampling frequency of $2{\rm \pi} \varOmega _{\textrm {ce}}$. Hence, for any wave mode with wave frequency above this frequency, i.e. $\omega >2{\rm \pi} \varOmega _{\textrm {ce}}\approx 6 \varOmega _{\textrm {ce}}$, a folding to lower-frequency region (i.e. in which $\omega < 2{\rm \pi} \varOmega _{\textrm {ce}}$) takes place. We also notice that in the large-wavevector regime (e.g. $k_{\perp }>4.5 \omega _{\textrm {ce}}/c$) the analytical dispersion relation curve of X mode deviates corresponding enhancement in PSD slightly for each simulation (see figure 14). The deviations between analytical dispersion relation curves and simulation results, in particular in the larger wavevector and frequency regime, are inherent in numerical simulations due to the usage of finite grid cell size and finite time step.

To compare the characteristic frequencies of electron cyclotron wave modes due to the different mechanisms mentioned above, we calculated the power spectrum by integrating the PSD of the transverse electric field over perpendicular wavevector $k_{\perp }$ in the $(\omega ,k_{\perp })$ domain, i.e.

(3.11)\begin{equation} \mathcal{P}_{\tau}(\omega)=\int |E_{\tau}\left(\omega,k_{{\parallel}}=0,k_{{\perp}}\right)|^2\,\textrm{d}k_{{\perp}}. \end{equation}

Figure 15 shows the resulting power spectra for all eight simulations in the time window $100\text {--}164 \omega _{\textrm {pe}}^{-1}$. The numbered peaks indicate the harmonics of electron cyclotron wave modes. The blue/red squares correspond to harmonics of electron cyclotron waves due to ECRs while red circles indicate harmonics due to the ECMIs. For electron cyclotron waves caused by ECRs, the power spectra generated by Maxwellian beam and by ring beam are practically the same at the characteristic frequencies of these wave modes. While for electron cyclotron waves produced by ECMIs, the power spectra generated by the ring beam significantly exceed those generated by the Maxwellian beam at their respective frequencies. Figure 15 further confirms that the density inhomogeneities of the background plasma have no influence on the generation of electron cyclotron waves.

Figure 15.

Power spectrum derived from transverse electric field $E_{\tau }$ for all simulations during $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$. The coloured squares denote harmonics of electron cyclotron mode due to ECRs and coloured circles indicate those due to ECMIs. Other labels and notations are used in the same way as in figure 9.

3.4 Doppler frequency shift

Figure 16 shows the PSD of transverse electric field $E_{\tau }$ for the simulations with homogeneous background plasma (i.e. Runs 1 and 2, $\eta =1$; see figure 16$a1{,}b1$) as well as with inhomogeneous background plasma (i.e. Runs 5 and 6, $\eta =50$; see figure 16$a2{,}b2$). As shown in figure 16, electron thermal fluctuations develop near the electron cyclotron frequency inside triangular regions delimited by black dashed lines, which are given by

(3.12)\begin{equation} \omega=\varOmega_{\textrm{ce}}\pm 3v_{\textrm{the}}k_{{\parallel}}. \end{equation}

For small $|k|\leq 5 \omega _{\textrm {pe}}/c$, the Z mode (indicated by a red dashed curve in figure 16) does not enter the triangular regions. In the short-wavelength regime $|k|>5 \omega _{\textrm {pe}}/c$, the Z mode asymptotically approaches the electron cyclotron frequency $\varOmega _{\textrm {ce}}=4 \omega _{\textrm {pe}}$ and enters into the triangular regions, becoming heavily damped.

Figure 16.

The PSD in the $k_{\parallel }$–$\omega$ plane for Runs 1, 2, 5, 6 in the time interval $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$. The PSD is evaluated as $\log _{10}|E_{\tau }/B_0|^2$. The dispersion relation curves of the R (magenta dashed curve), L (green dashed curve) and Z (red dashed curve) modes as well as those predicted by (3.13) (blue solid curve) and (3.14) (blue dashed line) are overlaid.

For the simulations with a ring beam, Runs 2 and 6, a significant enhancement of power near the Doppler-shifted ECR region is observed (see figure 16$b1{,}b2$). This region with enhanced power crosses the Z and L wave modes near frequencies $\omega \approx 2\text {--}3 \omega _{\textrm {pe}}$ at $k_{\parallel }\approx -4 \omega _{\textrm {pe}}/c$, and crosses the L and R wave modes above the ECR frequencies for $\omega \approx 4\text {--}6 \omega _{\textrm {pe}}$ at $k_{\parallel }\approx 4 \omega _{\textrm {pe}}/c$. These PSD enhancements can be fitted by an empirical expression (see blue solid curves denoted by ‘Doppler frequency shift’ in figure 16$b1{,}b2$):

(3.13)\begin{equation} \omega=\frac{1}{\gamma_{{\perp}}}\sqrt{\varOmega_{\textrm{ce}}\left(\varOmega_{\textrm{ce}}+u_{d\perp}k_{{\parallel}}\right)}. \end{equation}

In the long-wavelength limit, i.e. for $k_{\parallel }\to 0$, the fitting (3.13) reveals another linear expression (see blue dashed oblique lines denoted by ‘Doppler frequency shift (linear)’ in figure 16$b1{,}b2$):

(3.14)\begin{equation} \omega\approx \frac{1}{\gamma_{{\perp}}}\left(\varOmega_{\textrm{ce}}+\frac{1}{2}u_{d\perp}k_{{\parallel}}\right). \end{equation}

The fittings (3.13) and (3.14) indicate that the power enhancements are attributed to a Doppler frequency shift of waves at the electron cyclotron frequency associated with the beam drift speed.

Figure 17 shows the power spectra of the R, L and Z wave modes. These power spectra are calculated from the transverse electric field according to (3.11) in the time window $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$. Results show that, for the ring-beam cases (i.e. Runs 2, 4, 6, 8; see dashed curves in figure 17), the power of the R mode wave is enhanced mainly in the frequency range $\omega =4\text {--}6 \varOmega _{\textrm {ce}}$, while the L mode is enhanced for both $2\text {--}3 \varOmega _{\textrm {ce}}$ and $4\text {--}6 \varOmega _{\textrm {ce}}$. The Z wave mode is enhanced mainly for $2\text {--}3 \varOmega _{\textrm {ce}}$. The density inhomogeneities also influence the power spectra of wave modes. With an increase of the density gradient $\eta$, the power of the R mode decreases from approximately $10^{-1}$ to $10^{-6}$ in the frequency range $\omega =4\text {--}6 \varOmega _{\textrm {ce}}$. The same happens to the power of the Z mode: its power decreases from approximately $10^{-1}$ to $10^{-2}$ in the frequency range $\omega =2\text {--}3 \varOmega _{\textrm {ce}}$ as the density gradient increases. The L mode shows the opposite behaviour: its power increases from approximately $10^{-6}$ to $10^{-3}$ for $\omega =2\text {--}3 \varOmega _{\textrm {ce}}$ and from approximately $10^{-6}$ to $10^{-1}$ in the range $\omega =4\text {--}6 \varOmega _{\textrm {ce}}$ as the density gradient increases.

Figure 17.

Power spectrum (in unit of $\textrm {statV}^2\cdot\textrm {cm}^{-2} \cdot \textrm {s}$) of R $(a)$, L $(b)$ and Z $(c)$ modes for all simulations in $t=100\text {--}164 \omega _{\textrm {pe}}^{-1}$.