# High-harmonic Plasma Emission Induced by Electron Beams in Weakly Magnetized
Plasmas

Chuanyang Li (李传洋) Institute of Frontier and Interdisciplinary Science,
Shandong University, Qingdao, Shandong, 266237, People’s Republic of China
College of Physics and Electronic information, Dezhou University, Dezhou,
Shandong, 253023, People’s Republic of China Yao Chen (陈耀) Institute of
Frontier and Interdisciplinary Science, Shandong University, Qingdao,
Shandong, 266237, People’s Republic of China Institute of Space Sciences,
Shandong University, Shandong, 264209, People’s Republic of China Zilong Zhang
(张子龙) Institute of Space Sciences, Shandong University, Shandong, 264209,
People’s Republic of China Hao Ning (宁昊) Institute of Frontier and
Interdisciplinary Science, Shandong University, Qingdao, Shandong, 266237,
People’s Republic of China TangMu Li (李汤姆) Institute of Frontier and
Interdisciplinary Science, Shandong University, Qingdao, Shandong, 266237,
People’s Republic of China

(Received ***; Revised ***; Accepted ***)

###### Abstract

Electromagnetic radiation at higher harmonics of the plasma frequency
($\omega\sim n\omega_{pe},n>2$) has been occasionally observed in type II and
type III solar radio bursts, yet the underlying mechanism remains
undetermined. Here we present two-dimensional fully kinetic electromagnetic
particle-in-cell simulations with high spectral resolution to investigate the
beam-driven plasma emission process in weakly magnetized plasmas of typical
coronal conditions. We focused on the generation mechanisms of high-harmonic
emission. We found that a larger beam velocity ($u_{d}$) favors the generation
of the higher-harmonic emission. The emissions grow later for higher harmonics
and decrease in intensity by $\sim$2 orders of magnitude for each jump of the
harmonic number. The second and third harmonic ($\rm H_{2}$ and $\rm H_{3}$)
emissions get closer in intensity with larger $u_{d}$. We also show that (1)
the $\rm H_{3}$ emission is mainly generated via the coalescence of the $\rm
H_{2}$ emission with the Langmuir waves, i.e., $\rm H_{2}+L\rightarrow H_{3}$,
wherein the coalescence with the forward-propagating beam-Langmuir wave leads
to the forward-propagating $\rm H_{3}$, and coalescence with the backward-
propagating Langmuir wave leads to the backward-propagating $\rm H_{3}$; and
(2) the $\rm H_{4}$ emission mainly arises from the coalescence of the $\rm
H_{3}$ emission with the forward- (backward-) propagating Langmuir wave, in
terms of $\rm H_{3}+L\rightarrow H_{4}$.

Plasma astrophysics (1261); Solar electromagnetic emission(1490); Solar radio
emission(1522); Radio bursts (1339); Space plasmas(1544); Solar coronal radio
emission (1993);

##  1 Introduction

Plasma emission (PE) is defined as electromagnetic (EM) radiation at
frequencies close to the plasma frequency $\omega_{pe}$ (F) and its second
harmonic ($\rm H_{2}$). Several types of solar radio burst, the radio
radiation from the outer heliospheric boundary, and some radio emissions from
the planetary magnetosphere have been understood as PE (e.g., Etcheto &
Faucheux, 1984; Cairns, 1995, 1998; Gurnett et al., 1998; Cairns & Zank, 2002;
Kuncic & Cairns, 2005; Chen et al., 2014; Vasanth et al., 2016, 2019; Lv et
al., 2017; Li et al., 2017; Píša et al., 2017; Tasnim et al., 2022). Emission,
probably at the third harmonic ($\rm H_{3}$), has been occasionally reported
for type II (Bakunin et al., 1990; Kliem et al., 1992; Zlotnik et al., 1998;
Brazhenko et al., 2012) and type III solar radio bursts (Kundu, 1965; Takakura
& Yousef, 1974; Zlotnik, 1978; Cairns, 1986; Reiner et al., 1992; Reiner &
MacDowall, 2019). For other examples, Cairns (1986) reported harmonics at
$n\omega_{pe}$ ($n$ = 3–5) in the foreshock region of the bow shock with the
ISEE1 data, while Reiner & MacDowall (2019) detected possible $\rm H_{3}$
emissions in interplanetary type III bursts based on the Wind data.

Ginzburg & Zhelezniakov (1958) proposed the original framework of the theory
of PE, which has been developed into the standard PE model (e.g., Melrose,
1970a, b; Zheleznyakov & Zaitsev, 1970a, b; Melrose, 1980, 1987; Cairns, 1987;
Robinson et al., 1994). The model involves a multistep nonlinear process of
wave-particle and wave-wave interactions: (1) efficient excitation of Langmuir
(L) turbulence by electron beams through the kinetic bump-on-tail instability,
(2) scattering of L waves by ion-acoustic (IA) waves or ion density
inhomogeneities to generate the fundamental emission and/or backward-
propagating Langmuir ($\rm L^{\prime}$) waves ($\rm L\pm IA\rightarrow F$ and
$\rm L\pm IA\rightarrow L^{\prime}$), and (3) resonant coupling of forward-
and backward-propagating Langmuir turbulence to generate the harmonic emission
($\rm L+L^{\prime}\rightarrow H_{2}$). The standard PE theory has been widely
used to explain radio bursts in space and astrophysical plasmas, such as solar
radio bursts in terms of types I–V (see, e.g., Chen et al., 2014; Vasanth et
al., 2016; Lv et al., 2017; Li et al., 2017; Vasanth et al., 2019). Numerical
verifications based on fully kinetic EM particle-in-cell (PIC) simulations
(e.g., Thurgood & Tsiklauri, 2015; Che et al., 2017; Henri et al., 2019; Ni et
al., 2020; Chen et al., 2022; Zhang et al., 2022) or weak turbulence theory
have been performed (e.g., Yoon, 2000, 2005, 2006; Yoon et al., 2012; Ziebell
et al., 2015, 2016; Lee et al., 2019, 2022).

However, the generation mechanism of high-harmonic emissions ($\rm
H_{3},H_{4},\cdots$) remains inconclusive (Zheleznyakov & Zlotnik, 1974;
Cairns, 1987; Yin et al., 1998; Zlotnik et al., 1998; Ziebell et al., 2015).
Three different mechanisms have been proposed: (1) coalescence of three L
waves in terms of $\rm L+L^{\prime}+L^{\prime\prime}\rightarrow H_{3}$ (Kliem
et al., 1992); (2) coalescence of an L wave with $\rm H_{2}$ in terms of $\rm
L+H_{2}\rightarrow H_{3}$ (Zlotnik, 1978), where this mechanism has been
generalized by Cairns (1988) as $\rm L+H_{n}\rightarrow H_{n+1}$; and (3) $\rm
L+L_{n}\rightarrow H_{n+1}$ (Yi et al., 2007), i.e., coalescence of the
primary L wave with the $n$th-harmonic L mode ($\rm L_{n}$).

Earlier numerical studies mainly focused on the generation of the F and $\rm
H_{2}$ emissions (e.g., Pritchett & Dawson, 1983; Yin et al., 1998; Kasaba et
al., 2001), with few reports on high-harmonic emission. Using the two-
dimensional relativistic EM PIC code to simulate the beam-plasma interaction
process (Matsumoto & Omura , 1993), Rhee et al. (2009) found that high-
harmonic emission (up to the fifth) can be generated with the beam velocity
exceeding $0.5c$. By examining the radiation patterns and the temporal
correlation between the electrostatic (ES) and EM waves, they suggested that
the $\rm H_{3}$ ($\rm H_{4}$) emission resulted from the merging of the $\rm
H_{2}$ ($\rm H_{3}$) emission with the L waves, in agreement with the second
process mentioned above. Note that Rhee et al. (2009) utilized a relatively
small domain ($512\ \lambda_{De}\ \times\ 512\ \lambda_{De}$, where
$\lambda_{De}$ is the electron Debye length) and a limited duration (328
$\omega_{pe}^{-1}$) with a relatively low number of macroparticles (the number
of macroparticles per cell per species, NPPCPS, which, for background
electrons and protons and beam electrons, was 80, 80, and 8, respectively).
These limitations led to excessive numerical noise along dispersion curves and
a limited resolution. Thus, further validation is required.

Recently, Krafft & Savoini (2022) conducted two-dimensional PIC simulations in
weakly magnetized plasmas using larger-domain ($1448\ \lambda_{De}\ \times\
1448\ \lambda_{De}$), longer-duration (9000 $\omega_{pe}^{-1}$), and larger
NPPCPS values (1800, 1800, 1800). Their study focused on the coalescence
processes to generate $\rm H_{3}$ and $\rm H_{4}$ in the beam-plasma system
with or without density fluctuations. They suggested the following processes
for $\rm H_{3}$: $\rm H_{2}+L^{\prime}(L^{\prime\prime})\rightarrow H_{3}$,
where $\rm L^{\prime}(L^{\prime\prime})$ is supposed to be the product of the
first (second) cascade of the ES decay of the beam-Langmuir (BL) waves ($\rm
L\rightarrow L^{\prime}+S^{\prime}$ and $\rm L^{\prime}\rightarrow
L^{\prime\prime}+S^{\prime\prime}$); a very similar process has been suggested
for $\rm H_{4}$. They did not reject the possibility of other processes such
as $\rm L_{2}^{\prime}+L\rightarrow H_{3}$ and $\rm
L_{3}^{\prime}+L\rightarrow H_{4}$.

With careful evaluations, we suggest that the result of Krafft & Savoini
(2022) remains questionable and should be further tested due to the following
concerns.

(1) Regarding the Fourier analysis used to obtain the high-harmonic spectra in
the wavevector space, the authors restricted the frequency range to a small
window centered around $\omega_{n}\approx n\omega_{pe}$ (refer to their
Figures 1(c) and (d) and 5(c) and (d)). This means only signals within the
prescribed narrow range have been processed. Such analysis could not tell the
existence of high-harmonic radiation since, even in the thermal case, such
analysis can give enhanced signals around $\omega_{n}$, or any frequency
($>\omega_{pe}$). This is due to the local enhancement of numerical noise
along dispersion curves with PIC simulations. We have confirmed this with our
result, following the same procedure as employed by Krafft & Savoini (2022).

(2) Regarding the uncertainties of the energy evaluation of $\rm H_{4}$, the
$\rm H_{4}$ intensity is only marginally larger than the initial value in
their homogeneous case, with no discernable increase of $\rm H_{4}$ intensity
in the other case with density fluctuations (see Figures 2(a) and 4(a) in
Krafft & Savoini, 2022). This raises doubts as to whether the obtained $\rm
H_{4}$ emission is really excited by the beam-plasma interaction or is simply
a part of the numerical noise.

Another point is that the $\rm H_{3}$ and $\rm H_{4}$ emissions can be excited
at a beam velocity of $0.25c$ according to Krafft & Savoini (2022). This
contradicts the abovementioned result of Rhee et al. (2009). To clarify these
issues, we repeated the simulations of the beam interaction with a uniform
plasma, with basically the same configurations as used by Krafft & Savoini
(2022), such as the simulation domain, grid spacing, number of particles,
etc.. This is done with the Vector-PIC (VPIC) code (described in the next
section) and the same computational time (9000 $\omega_{pe}^{-1}$). Our
simulations only yielded weak $\rm H_{2}$ emission, without significant high-
harmonic emissions whose generation mechanism thus remains to be revealed.
This presents the main purpose of the present study.

##  2 Numerical Setup and Parameters

The simulations are performed using the VPIC code developed and released by
Los Alamos National Labs, which is two-dimensional in space and three-
dimensional for particle velocity and EM fields. VPIC employs a second-order,
explicit, leapfrog algorithm to update charged particle positions and
velocities in order to solve the relativistic kinetic equation for each
species, along with a full Maxwell description for electric and magnetic
fields evolved via a second-order finite-difference time-domain solver (Bowers
et al., 2008a, b, 2009). Periodic boundary conditions are used. The background
magnetic field is set to be $\vec{B}_{0}(=B_{0}\hat{e}_{z})$, and the
wavevector $\vec{k}$ is in the $xOz$ plane. The plasmas consist of three
components, including background electrons and protons with a Maxwellian
distribution and an electron beam with the following velocity distribution
function:

$f_{e}=A_{e}\exp\left(-\frac{u_{\perp}^{2}}{2u_{0}^{2}}-\frac{(u_{\parallel}-u_{d})^{2}}{2u_{0}^{2}}\right)$
(1)

where $u_{\parallel}$ and $u_{\perp}$ are the parallel and perpendicular
components of the momentum per mass, $u_{d}$ is the average drift momentum per
mass of the beam electrons, $u_{0}$ is the thermal velocity of energetic
electrons, and $A_{e}$ is the normalization factor.

Based on the observations of solar radio bursts (e.g., Wild & McCready, 1950;
Wild et al., 1959; Alvarez & Haddock, 1973; Reid et al., 2014) and in situ
data of energetic electrons (e.g., Lin et al., 1973, 1981, 1986), the drift
speed of the electron beam is set to be $u_{d}$ = 0.3–0.7$c$ and the thermal
velocity of background and beam electrons with a fixed value of
$u_{0}=0.028c$. The ratio of plasma oscillation frequency to electron
gyrofrequency ($\omega_{pe}/\Omega_{\rm ce}$) is set to be 100. The density
ratio of beam-background electrons ($n_{b}/n_{0}$) is set to be $0.5\%$. All
particles initially distribute homogeneously in space.

To facilitate the comparison with Krafft & Savoini (2022), we use the same
numerical configurations, such as the simulation domain ($L_{x}=L_{z}=1024\
\Delta$, where $\Delta=\sqrt{2}\ \lambda_{De}$ is the grid spacing and
$\lambda_{De}$ is the Debye length of the background electrons), the unit of
time ($\omega_{pe}^{-1}$), the simulation time 3000 $\omega_{pe}^{-1}$ and the
time step $\Delta t=0.01$ $\omega_{pe}^{-1}$. The corresponding resolvable
range of $|k|$ is [0.15, 79] $\omega_{pe}/c$, and the range of $\omega$ is
[0.002, 6.4] $\omega_{pe}$ (for the time interval of 3000 $\omega_{pe}^{-1}$).
The NPPCPS is 1800 for background electrons and 900 for both background
protons and beam electrons. The ratio of the ion to the electron plasma
temperatures is set to be $T_{i}/T_{e}=0.1$ to avoid strong Landau damping. At
the start of the simulation, we maintain charge neutrality by setting proper
weights to each species of macroparticle and neutralize the electron beam
current by setting a proper bulk speed to the background electrons (see, e.g.,
Henri et al., 2019; Chen et al., 2022; Zhang et al., 2022).

##  3 Numerical Results

We first present the reference case (Case R) with beam velocity $u_{d}=0.5c$
and the realistic proton-electron mass ratio ($m_{i}/m_{e}=1836$) that will be
compared with the corresponding thermal case (Case T) to tell the excitations
of waves. Then, we vary $u_{d}$ to investigate its effect. One experiment with
a larger mass ratio ($m_{i}/m_{e}=18360$) is conducted to expand the
investigation.

###  3.1 Wave Analysis for Case R

We first present the energy curves of the six field components and the
negative change of the total kinetic energy of all electrons ($-\Delta E_{k}$)
for Cases R and T (Figure 1(a) and (b)). The energies of all field components
in Case R are significantly stronger than those in Case T by at least 2–4
orders of magnitude. This means these fields arise due to the beam-plasma
interaction. Based on the energy profiles, the initial stage (0–100
$\omega_{pe}^{-1}$) is characterized by the rapid rise of $E_{x}$ and $E_{z}$
in energy, corresponding to the growth of the primary BL mode. At the end of
this stage, $-\Delta E_{k}$ reaches $\sim$4$\%$ of $E_{k0}$. The maximum of
$-\Delta E_{k}$ is $\sim$6$\%$ of $E_{k0}$, obtained at $\sim$300
$\omega_{pe}^{-1}$. From 100 to 900 $\omega_{pe}^{-1}$, the BL mode saturates
at an energy level of 0.03 $E_{k0}$. This is followed by the damping stage,
within which the wave energy is returned to electrons. Between 0 and 2000
$\omega_{pe}^{-1}$, the intensities of the other three components ($E_{y}$,
$B_{x}$, and $B_{z}$) slowly rise and saturate. The $E_{y}$ component is
mainly associated with the Z mode and the harmonic radiation, and the $B_{x}$
and $B_{z}$ components mainly carry the W mode, as well as the harmonic
radiation.

Figure 2 presents wave-energy maps in the wavevector ($\vec{k}$) space for the
ES modes around $\omega_{pe}$, 2 $\omega_{pe}$, and 3 $\omega_{pe}$ and the
harmonic radiations. The maps present the intensity maxima of waves at the
corresponding $\vec{k}$, regardless of their frequencies. The mode nature can
be identified with the analytic dispersion curves plotted in Figure 3 (and the
accompanying movies). The modes excited here include the $\rm L_{3}$ mode, the
$\rm L_{2}$ and $\rm L_{2b}$ modes, the BL and $\rm L_{b}$ modes, and the $\rm
H_{2}$ and $\rm H_{3}$ emissions. Here, the subscript “f” (“b”) denotes the
forward (backward)-propagating portion of the corresponding mode, while the
number in the subscript indicates the harmonics. The dominant $E_{x}$ and
$E_{z}$ components correspond to the primary forward-propagating ES-BL mode
(Figure 2(c)). The BL mode extends $\sim$$\pm 77^{\circ}$ away from the
parallel direction in $E_{z}$, with a significant perpendicular electric field
component of $E_{x}$. The movie accompanying Figure 2 shows its evolution in
the $\vec{k}$ space. During the rapid-excitation stage, the $E_{z}$ component
presents a rather narrow range in the $\vec{k}$ space, with $1.5\leq
k_{\parallel}(\omega_{pe}/c)\leq 3.5$ and $k_{\perp}(\omega_{pe}/c)\leq 8$;
later, it expands toward larger $k_{\parallel}$. It is within the range of
$2\leq k_{\parallel}(\omega_{pe}/c)\leq 15$ at $\sim$900 $\omega_{pe}^{-1}$,
followed by the decline of intensity and the shrinkage of the $\vec{k}$ range.

Figures 2(c) and 3(c) reveal the presence of two additional L components on
the left (backward) side of the BL mode, referred to as the forward- and
backward-propagating generalized Langmuir (GL) modes, represented by $\rm
GL_{f}$ and $\rm L_{b}$, respectively (see also Chen et al., 2022; Zhang et
al., 2022). As depicted in Figure 2 and its accompanying movie, both the $\rm
GL_{f}$ and $\rm L_{b}$ modes experience delayed growth with weaker
intensities compared to the BL mode. The $k_{\parallel}$ range of the $\rm
GL_{f}$ mode remains relatively narrow ($0\leq
k_{\parallel}(\omega_{pe}/c)\leq 2$), while the $k_{\parallel}$ of the $\rm
L_{b}$ mode expands gradually and reaches $k_{\parallel}(\omega_{pe}/c)\sim$
-15 at $\sim$900 $\omega_{pe}^{-1}$, with an angular pattern resembling that
of the primary BL mode. Subsequently, both modes damp with decreasing ranges
of $\vec{k}$.

The intensities of the BL, $\rm L_{2}$, and $\rm L_{3}$ modes decrease
gradually with increasing ranges of $k_{\parallel}$. The $k_{\parallel}$ range
of $\rm L_{2}$ expands from [3, 6] $\omega_{pe}/c$ initially to [5, 15]
$\omega_{pe}/c$ around 900 $\omega_{pe}^{-1}$. For $\rm L_{3}$, its
$k_{\parallel}$ range extends from [6, 8] to [9, 13] $\omega_{pe}/c$ during
[100, 900] $\omega_{pe}^{-1}$. The $\rm L_{2b}$ mode appears around 650
$\omega_{pe}^{-1}$, within $-10\leq k_{\parallel}(\omega_{pe}/c)\leq-4$. After
$\sim$900 $\omega_{pe}^{-1}$, these ES waves begin to decay with the shrinkage
of the $k_{\parallel}$ range.

The circular patterns shown in Figure 2(d) represent the $\rm H_{2}$ and $\rm
H_{3}$ radiation. We observe no signatures of higher-harmonic emissions.
According to the movie accompanying Figure 2, the $\rm H_{3}$ emission grows
after $\rm H_{2}$. The $\rm H_{2}$ radiation is stronger than $\rm H_{3}$, and
the forward-propagating portion of $\rm H_{3}$ is stronger than the backward-
propagating portion.

In Figure 3, we present the $\omega-k$ dispersion analysis along the
propagation angle $\theta_{kB}=30^{\circ}$. The accompanying movie shows the
variation of the dispersion relation at a step of $10^{\circ}$. Note that the
frequency ranges of $\rm H_{2}$ and $\rm H_{3}$ vary similarly with
$\theta_{kB}$. The overplotted dashed lines represent the dispersion relations
of the L wave with thermal effects
($\omega^{2}=\omega_{pe}^{2}(1+3k^{2}\lambda_{De}^{2})$) and wave modes given
by the classical cold plasma magnetoionic theory. Table 1 summarizes the
$\omega$ and $k$ ranges for each mode just before the rapid decay of the ES
modes ($\sim$900 $\omega_{pe}^{-1}$).

We examine the temporal development of various wave modes with energy profiles
for specific field component(s) (see Figures 1(c) and (e)), which are
calculated by integrating the field energy within a specified spectral range
(Figure 3(c)) along the corresponding dispersion curve according to Parseval’s
theorem. The forward-propagating BL, $\rm L_{2}$, and $\rm L_{3}$ modes
display the strongest and fastest growth; their intensities saturate around
100 $\omega_{pe}^{-1}$, then remain at an almost constant level before
declining gradually after $\sim$900 $\omega_{pe}^{-1}$. The $\rm L_{b}$ mode
is slightly stronger than $\rm GL_{f}$, both with similar energy profiles.
Both modes grow gradually over time, reaching the maximum level near 900
$\omega_{pe}^{-1}$. The $\rm L_{2b}$ mode shows a delayed onset of growth at
$\sim$600 $\omega_{pe}^{-1}$, followed by further enhancement and subsequent
decay to the thermal noise level.

The comparison with the corresponding thermal case (Figures 1(e) and (f))
confirms the significant enhancement of the $\rm H_{2}$ and $\rm H_{3}$
emissions over thermal noise, while emissions at the fourth and higher
harmonics are absent.

The energy profiles of $\rm H_{2}$ and $\rm H_{3}$ are similar. Both modes get
enhanced gradually over time, reaching the saturation level near 900
$\omega_{pe}^{-1}$ and remaining at a nearly constant level thereafter.
However, the onset of $\rm H_{3}$ is later by $\sim$ 400 $\omega_{pe}^{-1}$
compared to $\rm H_{2}$. The $\rm H_{3}$ intensity gets weaker than $\rm
H_{2}$, and their intensity ratio ($\rm W_{H_{2}}/W_{H_{3}}$) at saturation is
$\sim$120.

The forward- (f) and backward- (b) propagating portions of $\rm H_{2}$ show
little difference, and the intensity ratio ($\rm W_{H_{2f}}/W_{H_{2b}}$) at
saturation is $\sim$1.4. However, the forward-propagating portion of $\rm
H_{3}$ is much stronger than its backward-propagating portion, with $\rm
W_{H_{3f}}/W_{H_{3b}}$ $\sim$ 4 at saturation, agreeing with the wave-energy
maps shown in Figure 2(d).

###  3.2 Effect of the Beam Velocity and the Mass Ratio ($m_{i}/m_{e}$)

Two sets of numerical experiments are carried out. We first investigate the
effect of the beam velocity ($u_{d}$) with the physical mass ratio
($m_{i}/m_{e}$ = 1836). This gives four additional cases (A, B, C, and D),
corresponding to $u_{d}=0.3c,\ 0.4c,\ 0.6c$, and $0.7c$, respectively. We then
investigate the effect of $m_{i}/m_{e}$, since this parameter can affect the
intensity of the ES modes and the radiation process (Chen et al., 2022; Zhang
et al., 2022). This gives Case E with $u_{d}=0.5c$ and $m_{i}/m_{e}$ = 18360.

Figures 4 and 5 present the wave-energy maps in the $\vec{k}$ space for Cases
A ($u_{d}=0.3c$) and D ($u_{d}=0.7c$), respectively. The accompanying movies
show the evolution. In comparison with Case R ($u_{d}=0.5c$), the forward-
propagating ES modes ($\rm BL$, $\rm L_{2}$ and $\rm L_{3}$) have larger
$k_{\parallel}$ values in Case A, but smaller values in Case D. For instance,
the $k_{\parallel}$ range of $\rm BL$ in Case A extends from [3, 6.5]
$\omega_{pe}/c$ initially to [4, 11] $\omega_{pe}/c$, while it expands from
[1, 2] $\omega_{pe}/c$ to [1.5, 15] $\omega_{pe}/c$ in Case D.

For EM harmonics, the $\rm H_{2}$ and $\rm H_{3}$ emissions exhibit weaker
intensities in Case A than in Case R, with $\rm H_{3}$ only presenting a faint
backward-propagating part ($\rm H_{3b}$; see Figure 4(d) and the accompanying
movie). Figure 5(d) shows the significant enhancements of $\rm H_{2}$ and $\rm
H_{3}$, with the additional appearance of $\rm H_{4}$ emission in the forward-
propagation direction ($\rm H_{4f}$). According to the movie accompanying
Figure 5, the $\rm H_{3}$ emission grows after $\rm H_{2}$, and $\rm H_{4}$
grows after $\rm H_{3}$.

Figure 6 illustrates the temporal energy profiles of the ES modes for
different cases. As $u_{d}$ increases from $0.3c$ to $0.7c$, the intensity
maximum of each ES mode increases correspondingly. Among them, the BL, $\rm
L_{2}$, and $\rm L_{3}$ modes exhibit a similar trend, with the onset of the
decay being earlier for larger $u_{d}$. The $\rm L_{b}$ and $\rm GL_{f}$ modes
arise simultaneously and both grow and decay faster with increasing $u_{d}$.
Both modes reach their intensity maximum earlier with increasing $u_{d}$. The
times of their intensity maxima are indicated by circles in Figure 6 at 3000,
1282, 933, 735, and 677 $\omega_{pe}^{-1}$ for $\rm L_{b}$ mode and 1622,
1294, 993, 848, and 766 $\omega_{pe}^{-1}$ for $\rm GL_{f}$ mode. These times
correspond to the start of the decay of the respective BL mode.

Figure 7 illustrates the energy profiles of $\rm H_{2}$ and $\rm H_{3}$. Both
show considerable enhancements in Cases B–E and R, while significant $\rm
H_{4}$ emission appears only in Cases C and D, corresponding to $u_{d}=0.6c$
and $0.7c$. For different cases, the energy profiles of $\rm H_{2}$, $\rm
H_{3}$, and $\rm H_{4}$ are quite similar, with comparable growth rates and
saturation times. The $\rm H_{2}/H_{3}$ emissions saturate earlier with
increasing $u_{d}$. The $\rm H_{2}$ exhibits higher intensity and grows
earlier than $\rm H_{3}$. The $\rm H_{2}/H_{3}$ energy ratio ($\rm
W_{H_{2}}/W_{H_{3}}$) decreases with increasing $u_{d}$. For instance, as
$u_{d}$ increases from $0.3c$ to $0.7c$, $\rm W_{H_{2}}/W_{H_{3}}$ decreases
from 206.8 to 44.1; see Table 2 for details. Likewise, $\rm H_{3}$ exhibits
higher intensity and grows earlier than $\rm H_{4}$. In Cases C and D, the
$\rm H_{3}/H_{4}$ energy ratios ($\rm W_{H_{3}}/W_{H_{4}}$) are 93.6 and 88.5,
respectively.

The $\rm H_{2f}$ portion is stronger than $\rm H_{2b}$ for smaller $u_{d}$ (=
$0.3c$ and $0.4c$), while $\rm H_{2f}$ is closer to $\rm H_{2b}$ in intensity
with larger $u_{d}$ (from $0.5c$ to $0.7c$). Regarding $\rm H_{3}$, the $\rm
H_{3b}$ portion is stronger than $\rm H_{3f}$ with $u_{d}=0.3c$ (Case A). As
$u_{d}$ increases to $0.4c$ (Case B), $\rm H_{3b}$ is initially weaker, yet it
gets stronger than $\rm H_{3f}$ later. For $u_{d}\geq 0.5c$, $\rm H_{3f}$ is
always stronger than $\rm H_{3b}$. In Cases C and D, $\rm H_{4f}$ is always
stronger than $\rm H_{4b}$.

###  3.3 On the Generation Mechanism of $\rm H_{3}$

As introduced, three nonlinear wave-wave coupling processes have been proposed
to account for the $\rm H_{3}$ emission. The modes involved must satisfy the
matching conditions: $\vec{k}_{1}+\vec{k}_{2}=\vec{k}_{3}$ and
$\omega_{1}+\omega_{2}=\omega_{3}$. For Case R with $u_{d}=0.5c$, according to
the $\omega$ and $k$ ranges of various modes listed in Table 1, three
processes can satisfy the matching conditions: (a) $\rm BL+L_{2b}\rightarrow
H_{3}$, (b) $\rm L_{b}+L_{2}\rightarrow H_{3}$, and (c) $\rm
BL/L_{b}+H_{2}\rightarrow H_{3}$.

These possibilities can be further evaluated by analyzing the energy profiles
of relevant modes. From Figures 1(c) and (e), we have (1) the energy curves of
$\rm L_{b}$, $\rm H_{2}$, and $\rm H_{3}$ are correlated with similar growth
rates and consistent saturation times ($\sim$900 $\omega_{pe}^{-1}$); (2) $\rm
H_{2}$ exhibits higher intensity and grows earlier than $\rm H_{3}$; (3) $\rm
L_{2b}$ is much weaker than other modes and damps to the noise level shortly
after the brief enhancement; and (4) $\rm H_{3}$ grows earlier than $\rm
L_{2b}$.

The first point suggests that $\rm L_{b}$ may play a role in the generation of
both $\rm H_{2}$ and $\rm H_{3}$. Combining the first two points, $\rm H_{2}$
may get involved in generating $\rm H_{3}$, which supports process (c).

The third and forth points rule out the participation of $\rm L_{2b}$ in
generating $\rm H_{3}$, i.e., process (a).

Process (b) ($\rm L_{b}+L_{2}\rightarrow H_{3}$) can be largely rejected by
considering the late-stage (say, after $t$ = 2000 $\omega_{pe}^{-1}$)
characteristics of $\rm H_{3}$ and $\rm L_{2}$ in Cases A, B, and E. We
observe from Figures 7(a), (b), and (f) that $\rm H_{3}$ presents persistent
enhancement during this late stage. Yet, according to Figures 6(a), (b), and
(f) and the movies accompanying Figure 4, in the late stage, (1) in Cases A
and B, $\rm L_{2}$ damps significantly and approaches its initial level of
numerical noise; and (2) in Case A the $k_{\parallel}$ of $\rm L_{2}$ shifts
from the earlier range of $\sim$11–17 $\omega_{pe}/c$ to $\sim$13–17
$\omega_{pe}/c$, while that of $\rm L_{b}$ remains in the range of -6–0
$\omega_{pe}/c$, making it impossible to meet the corresponding matching
condition of process (b). A similar shift occurs to Case E in which the
$k_{\parallel}$ of $\rm L_{2}$ shifts from the earlier range of $\sim$6–11 to
$\sim$9–12 $\omega_{pe}/c$ while that of $\rm L_{b}$ remains in the range of
-7–0 $\omega_{pe}/c$, making it difficult to meet the corresponding matching
condition for the generation of $\rm H_{3b}$ with process (b). We conclude
that these observations do not favor the significance of process (b) in
generating $\rm H_{3}$.

Thus, the most likely process to generate $\rm H_{3}$ is (c), $\rm
BL/L_{b}+H_{2}\rightarrow H_{3}$. Based on the matching conditions of three-
wave interaction, we can further separate the process into two subprocesses:
$\rm BL+H_{2b}\rightarrow H_{3f}$ and $\rm L_{b}+H_{2f}\rightarrow H_{3b}$.
Note that $\rm H_{3f}$ is stronger and rises earlier than $\rm H_{3b}$,
indicating that the process $\rm BL+H_{2b}\rightarrow H_{3f}$ is more
effective.

We observe a different directional pattern of $\rm H_{3}$ in different cases
(A, B, and R). This can be used to constrain the generation mechanism of $\rm
H_{3f}$ and $\rm H_{3b}$. According to Figure 4 and the accompanying movie,
the $k$ value of the BL mode ($k_{\rm BL}$) is too large in Case A
($u_{d}=0.3c$) to coalesce with $\rm H_{2}$ so as to generate $\rm H_{3f}$;
this explains the absence of $\rm H_{3f}$ in this case. On the other hand,
$\rm L_{b}$ has the appropriate wavenumber to coalesce with $\rm H_{2}$ and
generate $\rm H_{3b}$. For Case B with $u_{d}=0.4c$ (see the second part of
the accompanying movie), the value of $k_{\rm BL}$ allows the occurrence of
the process $\rm BL+H_{2b}\rightarrow H_{3f}$. Note that $k_{\rm BL}$
increases with increasing time; this makes it more difficult to meet the
corresponding matching conditions, while the conditions for $\rm
L_{b}+H_{2f}\rightarrow H_{3b}$ can still be satisfied. This can explain the
relative intensity variation of $\rm H_{3b}$ and $\rm H_{3f}$. For Case R with
$u_{d}=0.5c$ (see Figure 2 and the accompanying movie), both the BL and $\rm
L_{b}$ modes can coalesce with $\rm H_{2}$ during the simulation. We get a
stronger $\rm H_{3f}$ than $\rm H_{3b}$ since BL is always stronger than $\rm
L_{b}$.

Note that in Case A we observed the simultaneous presence of $\rm GL_{f}$ and
$\rm H_{2}$, yet the $\rm H_{3f}$ radiation is always at the noise level. This
indicates that the process of $\rm GL_{f}+H_{2}\rightarrow H_{3}$ is not
important here.

###  3.4 On the Generation Mechanism of $\rm H_{4}$

In the following, we discuss the generation mechanism of $\rm H_{4}$. Based on
the matching conditions, the following processes are possible: (a) $\rm
BL+L_{3b}\rightarrow H_{4}$, (b) $\rm L_{b}+L_{3}\rightarrow H_{4}$, (c) $\rm
BL/L_{b}+H_{3}\rightarrow H_{4}$, and (d) $\rm L_{2}/L_{2b}+H_{2}\rightarrow
H_{4}$.

First, we can exclude process (a) involving the $\rm L_{3b}$ mode. According
to Figures 6(d) and (e) and 7(d) and (e), during the initial growth of $\rm
H_{4}$ $\rm L_{3b}$ is at the noise level before 600 $\omega_{pe}^{-1}$ and
then presents a local peak from 600 to 800 $\omega_{pe}^{-1}$. Therefore, it
is unlikely to be important for the radiation of $\rm H_{4}$.

We can also exclude process (b) based on the following two aspects. First,
according to Figures 6(d) and (e) and 7(d) and (e), $\rm L_{b}$ first
increases and then decreases in intensity along with the $\rm H_{4}$ growth,
while $\rm L_{3}$ exhibits a decreasing trend in intensity. This is different
from the $\rm H_{4}$ intensity profile. Second, $\rm L_{3}$ does not change
much in intensity when $u_{d}$ increases from $0.6c$ to $0.7c$, while $\rm
H_{4}$ reveals significant enhancement.

We can further exclude process (d). From Figure 7(e), $\rm H_{3}$ starts to
grow earlier than $\rm H_{4}$. As stated earlier, $\rm H_{2}$ contributes to
the $\rm H_{3}$ radiation. If it also contributes to the $\rm H_{4}$
radiation, e.g., via process (d), then $\rm H_{3}$ and $\rm H_{4}$ shall rise
at a close time. This is inconsistent with our simulation. In addition, when
$u_{d}$ increases from $0.6c$ to $0.7c$, $\rm L_{2b}$ does not change much in
intensity, while $\rm H_{4}$ presents substantial enhancement (see Figures
7(d) and (e)). This indicates that the process involving $\rm L_{2b}$ in (d)
can be ruled out.

Thus, we are left with process (c), with $\rm BL/L_{b}+H_{3}\rightarrow
H_{4}$. This is supported by the similar intensity profiles of $\rm L_{b}$,
$\rm H_{3}$, and $\rm H_{4}$. If further taking the matching conditions into
account, we find that the following subprocesses may act in the system: $\rm
BL+H_{3b}\rightarrow H_{4f}$ and $\rm L_{b}+H_{3f}\rightarrow H_{4b}$.

##  4 Conclusions and Discussion

In this study, we conduct fully EM PIC simulations to investigate the
interaction between beam electrons and a weakly magnetized plasma with the
solar coronal-interplanetary conditions. The purpose is to reveal the
generation mechanism of high-harmonic plasma radiations. We obtained
significant second, third, and fourth harmonic ($\rm H_{2}$, $\rm H_{3}$, and
$\rm H_{4}$) radiations. The neighboring harmonics are different in intensity
by about 2 orders of magnitude, with higher harmonics being weaker. Cases with
larger beam velocities tend to favor the generation of higher harmonics and
stronger radiation for the specific harmonic, which is consistent with the
statement of Rhee et al. (2009), “by increasing the beam velocity, the present
case allows for the excitation of higher-harmonic modes.” We suggest that the
$\rm H_{3}$ radiation is mainly generated via the process of $\rm
H_{2}+BL/L_{b}\rightarrow H_{3}$ and the $\rm H_{4}$ radiation is mainly
generated via the process of $\rm H_{3}+BL/L_{b}\rightarrow H_{4}$, with the
forward portion given by the coalescence of $\rm H_{2}$ ($\rm H_{3}$) with BL
mode and the backward portion given by its coalescence with $\rm L_{b}$.

The obtained results are helpful to understand the data of solar radio bursts
with high-harmonic emission. For instance, the presence of $\rm H_{3}$
emission in observations suggests that the energetic electron beam has a
higher-than-usual velocity (say, $>\sim 0.5c$), and the intensity ratio of
$\rm H_{3}$ and $\rm H_{2}$ depends on the beam velocity, with a larger
intensity ratio corresponding to a faster beam.

According to our earlier simulations (Chen et al., 2022; Zhang et al., 2022)
that investigated the PE process for weakly magnetized plasmas with
$\omega_{pe}/\Omega_{\rm ce}=10$ and nonmagnetized plasmas (i.e.,
$\omega_{pe}/\Omega_{ce}\sim\infty$), the excitation of primary modes and
escaping emissions is quite similar for these two cases. Thus, we expect that
the major conclusions deduced here still hold for other weakly magnetized
plasmas as long as the values of $\omega_{pe}/\Omega_{\rm ce}$ are much larger
than unity.

## Acknowledgements

This study is supported by NNSFC grants (Nos. 12103029, 1973031, 12203031), a
NSFSP (Natural Science Foundation of Shandong Province) grant (No.
ZR2021QA033) and the China Postdoctoral Science Foundation (No. 2021M691904).
The authors acknowledge the National Supercomputer Centers in Tianjin and the
Beijing Super Cloud Computing Center (BSCC; http://www.blsc.cn/) for providing
high-performance computing (HPC) resources, and the open-source Vector-PIC
(VPIC) code provided by Los Alamos National Labs (LANL). The authors are
grateful to the anonymous referee for valuable comments.

## References

  * Alvarez & Haddock (1973) Alvarez, H. & Haddock, F. T. 1973, Sol. Phys., 30, 175. doi:10.1007/BF00156186
  * Bakunin et al. (1990) Bakunin, L. M., Ledenev, V. G., Kosugi, T., et al. 1990, Sol. Phys., 129, 379. doi:10.1007/BF00159048
  * Bowers et al. (2008a) Bowers, K. J., Albright, B. J., Yin, L., et al. 2008a, Physics of Plasmas, 15, 055703. doi:10.1063/1.2840133
  * Bowers et al. (2008b) Bowers, K. J., Albright, B. J., Bergen, B., et al. 2008b, in SC ’08: Proc. 2008ACM/IEEE Conf. on Supercomputing No. 63 (Piscataway, NJ: IEEE press), 1 doi: 10.1109/SC.2008.5222734.
  * Bowers et al. (2009) Bowers, K. J., Albright, B. J., Yin, L., et al. 2009, Journal of Physics Conference Series, 180, 012055. doi:10.1088/1742-6596/180/1/012055
  * Brazhenko et al. (2012) Brazhenko, A. I., Melnik, V. N., Konovalenko, A. A., et al. 2012, Odessa Astronomical Publications, 25, 181
  * Cairns (1986) Cairns, I. H. 1986, J. Geophys. Res., 91, 2975. doi:10.1029/JA091iA03p02975
  * Cairns (1987) Cairns, I. H. 1987, Journal of Plasma Physics, 38, 169. doi:10.1017/S0022377800012496
  * Cairns (1995) Cairns, I. H. 1995, Geophys. Res. Lett., 22, 3433. doi:10.1029/95GL03331
  * Cairns (1988) Cairns, I. H. 1988, J. Geophys. Res., 93, 858. doi:10.1029/JA093iA02p00858
  * Cairns (1998) Cairns, I. H. 1998, ApJ, 506, 456. doi:10.1086/306215
  * Cairns & Zank (2002) Cairns, I. H. & Zank, G. P. 2002, Geophys. Res. Lett., 29, 1143. doi:10.1029/2001GL014112
  * Che et al. (2017) Che, H., Goldstein, M. L., Diamond, P. H., et al. 2017, Proceedings of the National Academy of Science, 114, 1502. doi:10.1073/pnas.1614055114
  * Chen et al. (2022) Chen, Y., Zhang, Z., Ni, S., et al. 2022, ApJ, 924, L34. doi:10.3847/2041-8213/ac47fa
  * Chen et al. (2014) Chen, Y., Du, G., Feng, L., et al. 2014, ApJ, 787, 59. doi:10.1088/0004-637X/787/1/59
  * Etcheto & Faucheux (1984) Etcheto, J. & Faucheux, M. 1984, J. Geophys. Res., 89, 6631. doi:10.1029/JA089iA08p06631
  * Ginzburg & Zhelezniakov (1958) Ginzburg, V. L. & Zhelezniakov, V. V. 1958, Soviet Ast., 2, 653
  * Gurnett et al. (1998) Gurnett, D. A., Allendorf, S. C., & Kurth, W. S. 1998, Geophys. Res. Lett., 25, 4433. doi:10.1029/1998GL900201
  * Henri et al. (2019) Henri, P., Sgattoni, A., Briand, C., et al. 2019, Journal of Geophysical Research (Space Physics), 124, 1475. doi:10.1029/2018JA025707
  * Kasaba et al. (2001) Kasaba, Y., Matsumoto, H., & Omura, Y. 2001, J. Geophys. Res., 106, 18693. doi:10.1029/2000JA000329
  * Kliem et al. (1992) Kliem, B., Krueger, A., & Treumann, R. A. 1992, Sol. Phys., 140, 149. doi:10.1007/BF00148435
  * Krafft & Savoini (2022) Krafft, C. & Savoini, P. 2022, ApJ, 934, L28. doi:10.3847/2041-8213/ac7f28
  * Kuncic & Cairns (2005) Kuncic, Z. & Cairns, I. H. 2005, Journal of Geophysical Research (Space Physics), 110, A07107. doi:10.1029/2004JA010953
  * Kundu (1965) Kundu, M. R. 1965, New York: Interscience Publication, 1965
  * Lee et al. (2019) Lee, S.-Y., Ziebell, L. F., Yoon, P. H., et al. 2019, ApJ, 871, 74. doi:10.3847/1538-4357/aaf476
  * Lee et al. (2022) Lee, S.-Y., Yoon, P. H., Lee, E., et al. 2022, ApJ, 924, 36. doi:10.3847/1538-4357/ac32bb
  * Li et al. (2017) Li, C. Y., Chen, Y., Wang, B., et al. 2017, Sol. Phys., 292, 82. doi:10.1007/s11207-017-1108-1
  * Lin et al. (1981) Lin, R. P., Potter, D. W., Gurnett, D. A., et al. 1981, ApJ, 251, 364. doi:10.1086/159471
  * Lin et al. (1986) Lin, R. P., Levedahl, W. K., Lotko, W., et al. 1986, ApJ, 308, 954. doi:10.1086/164563
  * Lin et al. (1973) Lin, R. P., Evans, L. G., & Fainberg, J. 1973, Astrophys. Lett., 14, 191
  * Lv et al. (2017) Lv, M. S., Chen, Y., Li, C. Y., et al. 2017, Sol. Phys., 292, 194. doi:10.1007/s11207-017-1218-9
  * Matsumoto & Omura (1993) Matsumoto, H., & Omura, Y. 1993, in Computer Space Plasma Physics: Simulation Techniques and Software, ed. H. Matsumoto & Y. Omura (Tokyo: Terra Sci. Pub.)
  * Melrose (1980) Melrose, D. B. 1980, Space Sci. Rev., 26, 3. doi:10.1007/BF00212597
  * Melrose (1987) Melrose, D. B. 1987, Sol. Phys., 111, 89. doi:10.1007/BF00145443
  * Melrose (1970a) Melrose, D. B. 1970a, Australian Journal of Physics, 23, 871. doi:10.1071/PH700871
  * Melrose (1970b) Melrose, D. B. 1970b, Australian Journal of Physics, 23, 885. doi:10.1071/PH700885
  * Ni et al. (2020) Ni, S., Chen, Y., Li, C., et al. 2020, ApJ, 891, L25. doi:10.3847/2041-8213/ab7750
  * Píša et al. (2017) Píša, D., Kurth, W. S., Hospodarsky, G. B., et al. 2017, European Planetary Science Congress
  * Reid et al. (2014) Reid, H. A. S., Vilmer, N., & Kontar, E. P. 2014, A&A, 567, A85. doi:10.1051/0004-6361/201321973
  * Reiner & MacDowall (2019) Reiner, M. J. & MacDowall, R. J. 2019, Sol. Phys., 294, 91. doi:10.1007/s11207-019-1476-9
  * Reiner et al. (1992) Reiner, M. J., Stone, R. G., & Fainberg, J. 1992, ApJ, 394, 340. doi:10.1086/171586
  * Rhee et al. (2009) Rhee, T., Ryu, C.-M., Woo, M., et al. 2009, ApJ, 694, 618. doi:10.1088/0004-637X/694/1/618
  * Robinson et al. (1994) Robinson, P. A., Cairns, I. H., & Willes, A. J. 1994, ApJ, 422, 870. doi:10.1086/173779
  * Pritchett & Dawson (1983) Pritchett, P. L. & Dawson, J. M. 1983, Physics of Fluids, 26, 1114. doi:10.1063/1.864222
  * Takakura & Yousef (1974) Takakura, T. & Yousef, S. 1974, Sol. Phys., 36, 451. doi:10.1007/BF00151214
  * Tasnim et al. (2022) Tasnim, S., Zank, G. P., Cairns, I. H., et al. 2022, ApJ, 928, 125. doi:10.3847/1538-4357/ac5031
  * Thurgood & Tsiklauri (2015) Thurgood, J. O. & Tsiklauri, D. 2015, A&A, 584, A83. doi:10.1051/0004-6361/201527079
  * Vasanth et al. (2019) Vasanth, V., Chen, Y., Lv, M., et al. 2019, ApJ, 870, 30. doi:10.3847/1538-4357/aaeffd
  * Vasanth et al. (2016) Vasanth, V., Chen, Y., Feng, S., et al. 2016, ApJ, 830, L2. doi:10.3847/2041-8205/830/1/L2
  * Wild & McCready (1950) Wild, J. P. & McCready, L. L. 1950, Australian Journal of Scientific Research A Physical Sciences, 3, 387. doi:10.1071/CH9500387
  * Wild et al. (1959) Wild, J. P., Sheridan, K. V., & Neylan, A. A. 1959, Australian Journal of Physics, 12, 369. doi:10.1071/PH590369
  * Yi et al. (2007) Yi, S., Yoon, P. H., & Ryu, C.-M. 2007, Physics of Plasmas, 14, 013301. doi:10.1063/1.2424556
  * Yin et al. (1998) Yin, L., Ashour-Abdalla, M., El-Alaoui, M., et al. 1998, Geophys. Res. Lett., 25, 2609. doi:10.1029/98GL01989
  * Yoon et al. (2012) Yoon, P. H., Ziebell, L. F., Gaelzer, R., et al. 2012, Physics of Plasmas, 19, 102303. doi:10.1063/1.4757224
  * Yoon (2006) Yoon, P. H. 2006, Physics of Plasmas, 13, 022302. doi:10.1063/1.2167587
  * Yoon (2000) Yoon, P. H. 2000, Physics of Plasmas, 7, 4858. doi:10.1063/1.1318358
  * Yoon (2005) Yoon, P. H. 2005, Physics of Plasmas, 12, 042306. doi:10.1063/1.1864073
  * Zhang et al. (2022) Zhang, Z., Chen, Y., Ni, S., et al. 2022, ApJ, 939, 63. doi:10.3847/1538-4357/ac94c6
  * Zheleznyakov & Zlotnik (1974) Zheleznyakov, V. V. & Zlotnik, E. Y. 1974, Sol. Phys., 36, 443. doi:10.1007/BF00151213
  * Zheleznyakov & Zaitsev (1970a) Zheleznyakov, V. V. & Zaitsev, V. V. 1970a, Soviet Ast., 14, 250
  * Zheleznyakov & Zaitsev (1970b) Zheleznyakov, V. V. & Zaitsev, V. V. 1970b, Soviet Ast., 14, 47
  * Ziebell et al. (2016) Ziebell, L. F., Petruzzellis, L. T., Yoon, P. H., et al. 2016, ApJ, 818, 61. doi:10.3847/0004-637X/818/1/61
  * Ziebell et al. (2015) Ziebell, L. F., Yoon, P. H., Petruzzellis, L. T., et al. 2015, ApJ, 806, 237. doi:10.1088/0004-637X/806/2/237
  * Zlotnik (1978) Zlotnik, E. I. 1978, Soviet Ast., 22, 228
  * Zlotnik et al. (1998) Zlotnik, E. Y., Klassen, A., Klein, K.-L., et al. 1998, A&A, 331, 1087

Figure 1:  (a) and (b) Temporal profiles of energies of various field
components ($E_{x},\ E_{y},\ E_{z},\ B_{x},\ B_{y}$, and $B_{z}$) and the
negative change of the total electron energy ($-\Delta E_{k}$). Temporal
energy profiles of different modes, where (c) and (d) are for the ES modes and
(e) and (f) are for the EM radiation. The subscripts “f” and “b” indicate the
forward- and backward-propagating parts of the respective modes, while the
subscript numbers represent the corresponding harmonics. The left column
presents the results of Case R, and the right column displays the
corresponding thermal Case T. The values in the plots are normalized by the
respective initial kinetic energy of total electrons ($E_{k0}$).

Figure 2:  Maximum intensity of the six field components in the $\omega$
domain ((a) 2.5–3.5, (b) 1.5–2.5, (c) 0.5–1.5, and (d) 1.5–6 $\omega_{pe}$) as
a function of $k_{\parallel}$ and $k_{\perp}$ over the interval of
$900<\omega_{pe}t<950$ for Case R. The upper three rows illustrate the
corresponding ES modes, while the fourth row represents the EM radiation
displayed by the magnetic field components. An animation of this figure is
available. The animation begins with the map for the interval of [0, 50]
$\omega_{pe}^{-1}$ and ends with that for the interval of [2950, 3000]
$\omega_{pe}^{-1}$. The real-time duration of the video is 12 s.

(An animation of this figure is available.)

Figure 3:  Wave dispersion diagrams of the six field components along
$\theta_{kB}=0^{\circ}$ ($\theta_{kB}$ is the angle between the wavevector,
$\vec{k}$, and the background magnetic field, $\vec{B}_{0}$). The curve near
$\omega_{pe}$ in panel (c) represents the dispersion relation of L wave with
thermal effects ($\omega^{2}=\omega_{pe}^{2}(1+3k^{2}\lambda_{De}^{2})$),
while the other white dashed lines represent the dispersion curves given by
the magnetoionic theory. “Z” denotes the Z mode, and “W” denotes the whistler
mode. The yellow boxes mark the regions selected for the evaluation of mode
energy. The white boxes in panels (b) and (e) zoom in on the $\rm H_{2}$ and
$\rm H_{3}$ spectra. An animation of this figure is available. The animation
begins at $\theta_{kB}=0^{\circ}$ and advances $10^{\circ}$ at a time up to
$\theta_{kB}=90^{\circ}$. The real-time duration of the video is 3 s.

(An animation of this figure is available.)

Table 1: Ranges of $\omega$ and $k$ for each mode at the moment before the rapid decay of each ES harmonic mode ($\sim$900 $\omega_{pe}^{-1}$) for Case R. | BL | $\rm L_{b}$ | $\rm GL_{f}$ | $\rm L_{2}$ | $\rm L_{2b}$ | $\rm L_{3}$ | $\rm H_{2}$ | $\rm H_{3}$
---|---|---|---|---|---|---|---|---
$\omega\ (\omega_{pe})$ | 0.8–1.18 | 0.88–1.04 | 0.9–1 | 1.8–2.06 | 1.88–2 | 2.9–3.06 | 1.88–2.06 | 2.89–3.01
$k\ (\omega_{pe}/c)$ | 2–15 | -15–0 | 0–2 | 5–15 |  -10$\sim$-4 | 9–13 | -1.8–1.8 | -2.8–2.8

Figure 4:  Maximum intensity of the six field components in the $\omega$
domain ((a) 2.5–3.5, (b) 1.5–2.5, (c) 0.5–1.5, and (d) 1.5–6 $\omega_{pe}$) as
a function of $k_{\parallel}$ and $k_{\perp}$ over the interval of
$2950<\omega_{pe}t<3000$ for Case A. The upper three rows illustrate the
corresponding ES modes, while the fourth row represents the EM radiation
displayed by the magnetic field components. An animation of this figure is
available, and it presents the mode evolution for Cases A, B and E, with each
part beginning with the map for the interval of [0, 50] $\omega_{pe}^{-1}$ and
ending with that for the interval of [2950, 3000] $\omega_{pe}^{-1}$. The
real-time duration of the video is 36 s.

(An animation of this figure is available.)

Figure 5:  Maximum intensity of the six field components in the $\omega$
domain ((a) 2.5–3.5, (b) 1.5–2.5, (c) 0.5–1.5, and (d) 1.5–6 $\omega_{pe}$) as
a function of $k_{\parallel}$ and $k_{\perp}$ over the interval of
$800<\omega_{pe}t<850$ for Case D. The upper three rows illustrate the
corresponding ES modes, while the fourth row represents the EM radiation
displayed by the magnetic field components. The video begins with the map for
the interval of [0, 50] $\omega_{pe}^{-1}$, ending with that for the interval
of [2950, 3000] $\omega_{pe}^{-1}$. The real-time duration of the video is 12
s.

(An animation of this figure is available.)

Figure 6:  The temporal energy profiles of the ES harmonic modes (BL, $\rm
L_{b},\ GL_{f},\ L_{2},\ L_{2b},\ L_{3},\ L_{3b}$), normalized by the
respective initial kinetic energy of total electrons ($E_{k0}$). (a)–(e) are
for cases with $m_{i}/m_{e}$ = 1836 and $u_{d}$ = 0.3–0.7$c$, and (f) is for
Case E with $u_{d}=0.5c$ and $m_{i}/m_{e}=18360$. The magenta and cyan circles
in panels (a)–(e) mark the intensity maxima of the $\rm L_{b}$ and $\rm
GL_{f}$ modes, respectively.

Figure 7:  The temporal energy profiles of the second, third, and fourth EM harmonic emission. The solid lines represent the sum of the energy in each direction of the corresponding harmonic emission, and the dotted and dashed lines represent the forward (f; $\ k_{\parallel}>0$) and backward (b; $\ k_{\parallel}<0$) propagation part of the energy of the corresponding radiation, respectively. The values in the figures are normalized by the respective initial kinetic energy of total electrons ($E_{k0}$).  Table 2: The intensity ratios of $\rm H_{2}/H_{3}$ and $\rm H_{3}/H_{4}$ at saturation for Cases R and A–E. The energy ratios of the forward- (f) and backward- (b) propagating portions for $\rm H_{2}$, $\rm H_{3}$, and $\rm H_{4}$ are also listed. | Case A | Case B | Case R | Case C | Case D | Case E
---|---|---|---|---|---|---
$\rm W_{H2}/W_{H3}$ | 206.8 | 169.3 | 120.0 | 58.3 | 44.1 | 124.8
$\rm W_{H3}/W_{H4}$ | – | – | – | 93.6 | 88.5 | –
$\rm W_{H2f}/W_{H2b}$ | 2.8 | 3.0 | 1.4 | 1.6 | 0.8 | 4.3
$\rm W_{H3f}/W_{H3b}$ | 0.3 | 0.7 | 4.0 | 7.5 | 8.5 | 1.1
$\rm W_{H4f}/W_{H4b}$ | – | – | – | 1.4 | 3.7 | –