Wave emission of non-thermal electron beams generated by magnetic reconnection

Magnetic reconnection in solar flares can efficiently generate non-thermal electron beams. The energetic electrons can, in turn, cause radio waves through microscopic plasma instabilities as they propagate through the ambient plasma along the magnetic field lines. We aim at investigating the wave emission caused by fast moving electron beams (FEBs) with characteristic non-thermal electron velocity distribution functions (EVDFs) generated by kinetic magnetic reconnection: two-streaming EVDFs along the separatrices and in the diffusion region, and perpendicular crescent-shaped EVDFs closer to the diffusion region. For this purpose, we utilized 2.5D fully kinetic Particle-In-Cell (PIC) code simulations in this study. We found that:(1) the two-streaming EVDFs plus the background ions are unstable to electron/ion (streaming) instabilities which cause ion acoustic waves and Langmuir waves due to the net current. This can lead to multiple harmonic plasma emission in the diffusion region and the separatrices of reconnection. (2) The perpendicular crescent-shaped EVDFs can cause multiple harmonic electromagnetic electron cyclotron waves through the electron cyclotron maser instabilities in the diffusion region of reconnection. Our results are applicable to diagnose the plasma parameters which control magnetic reconnection in solar flares by means of radio waves observations.

monic electromagnetic electron cyclotron waves through the electron cyclotron maser instabilities in the diffusion region of reconnection. Our results are applicable to diagnose the plasma parameters which are associated to magnetic reconnection in solar flares by means of radio waves observations.
An important consequence of the acceleration of particles in magnetic reconnection is the formation of non-thermal electron velocity distribution functions (EVDFs) (Yao et al. 2022). Non-thermal EVDFs are prone to decay unstably (Lee et al. 2007), which might cause non-linear phenomena like double layers (Lee et al. 2008), anomalous transport and anisotropic heating of the coronal flare-loop plasma (Lee & Büchner 2010, 2011a. The non-thermal EVDFs can cause microscopic plasma instabilities and generate electromagnetic emission across the entire spectrum including radio waves. Previous radio and hard X-ray observations of non-thermal radiation have confirmed the existence of such non-thermal sources of wave emissions from reconnection regions (Glesener et al. 2012;Narukage et al. 2014;Fletcher & Hudson 2008). For example, solar radio bursts (SRBs) are proposed to be caused by electron beams with non-thermal velocity distribution functions energized by magnetic reconnection in solar flare regions (Benz 2002;Mann et al. 2009).
In general, several kinds of waves have been associated to magnetic reconnection (see reviews in Fujimoto et al. 2011;Khotyaintsev et al. 2019). Some of those waves can be electromagnetic and in the radio-band, which are thought to be caused by electron beams. Observations of those waves from, e.g., the flaring solar corona, can serve as a remote diagnostics tool to probe the plasma conditions of astrophysical reconnection processes. Based on observations of coronal Type III SRBs by the upgraded Karl G. Jansky Very Large Array (VLA), Chen et al. (2018) reported waves at the fundamental plasma frequency and its harmonics, emitted from two separate source regions near magnetic null points in the solar corona.
Other waves due to reconnection have been observed in situ in the Earth's magnetosphere. Not all of them are electromagnetic and therefore can escape the plasma, but they are part of certain wave-wave interactions which can finally lead to electromagnetic radiation. Large-amplitude and high-frequency Langmuir waves were confirmed in magnetopause reconnection (Burch et al. 2016aGraham et al. 2018). High-frequency electrostatic upper-hybrid (UH) waves, generated due to energetic electrons of perpendicular crescent-shaped EVDFs, were observed by the high spatial and temporal resolution multi-spacecraft Magnetospheric Multiscale (MMS) mission not only in the diffusion region of asymmetric magnetopause reconnection (Graham et al. 2017) but also in central reconnection regions of the Earth's magnetotail . Menietti (2002) observed electrostatic electron cyclotron waves (ECWs) due to energetic electron beams in Earth's polar regions detected by the Plasma Wave Instrument on the Polar spacecraft.
Harmonics of electrostatic ECWs have also been observed in an over-dense plasma in the diffusion region of asymmetric reconnection in the Earth's magnetopause (Rönnmark & Christiansen 1981;Tang et al. 2013;Zhou et al. 2016;Li et al. 2020). By using MMS data, Li et al. (2020) found that the generation of multiple harmonic ECWs stands in close relation to the observation of crescentshaped EVDFs perpendicular to the local magnetic field in the Earth's magnetopause.
There are several possible emission mechanisms that can be responsible for the generation of radio waves by non-thermal electrons accelerated by magnetic reconnection. The most widely accepted mechanism is the so-called plasma emission mechanism that relies on the wave-wave interaction of Langmuir waves induced by electron beam instabilities with ion-acoustic waves (e.g., see Ginzburg & Zhelezniakov 1958;Melrose 1970a,b;Reid & Ratcliffe 2014;Melrose 2017, and references therein).
According to the plasma emission mechanism, an electron beam first generates electrostatic Langmuir waves (L) via streaming-like instabilities, which are due to a positive velocity gradient in the EVDF in the direction parallel to the magnetic field: v · ∂f /∂v > 0. EVDFs prone to those instabilities are known to be generated by magnetic reconnection (Yao et al. 2022). The energy density of the Langmuir waves usually exceeds the thermal fluctuation level of the ambient plasma by several orders of magnitude. The beam-generated Langmuir waves L can decay into ion-acoustic waves S and fundamental emission F (transverse electromagnetic waves at the plasma frequency), i.e., L → S + F (Melrose 1987(Melrose , 2017. The fundamental emission can also be generated by a coalescence process of backward-scattered Langmuir waves L and ion-acoustic waves S, i.e., L + S → F (Melrose 1987(Melrose , 2017. The (second) harmonic emission H (= H 2 ), i.e., electromagnetic waves at harmonic frequencies ω = 2ω pe , can be generated by a coalescence of Langmuir waves L and back-scattered Langmuir waves L through L + L → H (Willes et al. 1996;Yoon 2006). These usually circularly polarized transverse waves can escape from the ambient plasma, if their frequencies exceed the local plasma frequency which depends on the plasma density in the source region (Melrose 2017;Melrose et al. 1978). The phase velocities of these waves must also approximately exceed the speed of light, namely, they have to be superluminal. Other mechanisms have also been recently developed to solve some problems of the standard plasma emission mechanism. For example, Che et al. (2017) proposed a wave-wave interaction process where the modulation of Langmuir waves by whistler waves during the non-linear stage of an electron two-stream instability allows a longer emission period (continuous coherent emission) than in previous models.
Not only electron beams can generate Langmuir waves but also a variety of other processes, like a current density possibly driven by the relative drift between electrons and ions, as often observed in current sheets. A relative electron/ion streaming can offer a source of free energy for a family of instabilities collectively called electron/ion (streaming) instabilities. There are at least two wellstudied instabilities. The high-frequency one occurs when the relative streaming speed between the electron and ion populations significantly exceeds the electron thermal speed, and it is called Buneman instability (Buneman 1958(Buneman , 1959, leading to the generation of Langmuir waves (Gary 1993;Treumann & Baumjohann 2001;Jain et al. 2011). The low-frequency one is called ion-acoustic instability, which requires a relative drift speed much lower than the Buneman instability, but well above the ion thermal speed and mainly occurs for T e T i , leading to the generation of ion acoustic waves (Gary 1993;Treumann & Baumjohann 2001;Hellinger et al. 2004). This way, a plasma featuring an electron/ion streaming could also provide the ingredients for the typical wave-wave interaction of the plasma emission mechanism to generate (fundamental) electromagnetic waves at the electron plasma frequency.
The third and higher harmonic plasma emissions, at the frequency of multiples of the local plasma frequency, were also observed in SRBs (Takakura & Yousef 1974;Reiner & MacDowall 2019;Cairns 1986). Various schemes have been proposed to explain the formation of higher harmonic plasma emission. Zlotnik (1978) put forward that the third harmonic emission H 3 is generated by the coalescence process of the Langmuir waves and the harmonic emission H, i.e., L + H → H 3 . Cairns (1988) generalized this idea to explain higher harmonics of transverse wave emissions, claiming that nth harmonic waves H n could be generated by a coalescence of beam-generated Langmuir waves L and adjacent transverse waves H n−1 , i.e., L + H n−1 → H n . An alternative way of generation of harmonic transverse waves would be an interaction of Langmuir waves L and harmonic electrostatic waves L n−1 , i.e., L + L n−1 → H n (Gaelzer et al. 2003;Yi et al. 2007). Ziebell et al. (2015) numerically solved the weak-turbulence equations and theoretically demonstrate the process of multiple harmonic plasma emission, including the non-linear conversion from Langmuir turbulence to electromagnetic radiation. Rhee et al. (2009);Ziebell et al. (2015); Yao et al. (2021) (Thurgood & Tsiklauri 2015;Henri et al. 2019). The PIC simulations of ring-beam EVDFs by Zhou et al. (2020) also provided evidence for the generation of second and third harmonic emission via wave-wave interactions. Annenkov et al. (2019) used PIC code simulations to describe localized beams and found that emissions at harmonics of the plasma frequency can be generated due to an antenna mechanism.
Electron cyclotron waves (ECWs) at the relativistic electron cyclotron frequency and their harmonics have been often observed in a variety of environments. They can be generated by the so-called electron cyclotron maser instability (ECMI). The source of free energy of this instability is due to a positive velocity gradient in EVDFs in the direction perpendicular to the local magnetic field: ∂f /∂v ⊥ > 0. ECMI could be caused by EVDFs including, e.g., cup-like distribution functions (Büchner & Kuska 1996), horseshoe-like EVDFs (Bingham & Cairns 2000;Melrose & Wheatland 2016), ring- (Pritchett 1984;Lee et al. 2011;Yao et al. 2021) and crescent-shaped EVDFs in the velocity space perpendicular to magnetic field (Burch et al. 2016b;Egedal et al. 2016;Chen et al. 2016b). The ECMI produces electromagnetic waves at the local electron cyclotron frequency and its harmonics as well as UH waves (Benáček & Karlický 2019;Ni et al. 2020). The electron cyclotron maser emission (ECME) mechanism predicts electromagnetic emission amplified by a quasi-linear wave-particle interaction in the magnetic field (Twiss 1958;Melrose et al. 1984;Melrose 2017).
Observations show that ECWs are usually X-mode polarized (Ellis 1962;Ellis & Mcculloch 1963;Melrose 2017). The ECME theory of X polarized ECWs at fundamental electron cyclotron frequency leads to the frequency condition Ω ce > ω pe (Ellis 1962;Ellis & Mcculloch 1963;Hewitt & Melrose 1985;Melrose et al. 1984;Melrose 1986Melrose , 2017, which is also necessary for electromagnetic ECWs to escape from source region. The above frequency condition implies a sufficiently low plasma density (ω pe ∝ √ n e , with n e the electron plasma density) in regions of strong magnetic fields (Ω ce ∝ B with B the background magnetic field strength). In the solar corona this condition can be fulfilled in some localized regions near active regions and solar flares (Régnier 2015;Morosan et al. 2016).
However, the plasma-to-cyclotron frequency ratio is usually the opposite one in typical space and astrophysical plasmas, in particular, in the solar corona, where the electron cyclotron frequency is often smaller than the local plasma frequency Ω ce < ω pe . On the other hand, observations of X polarized ECWs in dynamic spectra have shown that ECWs can possibly not start with the fundamental but with a higher harmonic of the electron cyclotron frequency. This implies its source region may possibly be located in a plasma where Ω ce /ω pe < 1 (Treumann et al. 2011;Treumann & Baumjohann 2017). Treumann & Baumjohann (2017) proposed that higher harmonic of ECWs can be generated in a source region with frequency condition Ω ce /ω pe < 1, they may escape out and be remotely observed if they are excited with high intensities. Treumann & Baumjohann (2017) pointed out that the n = 5 harmonic ECWs can serve as cyclotron harmonic seed of the X mode polarized electromagnetic fluctuations after being amplified. By means of PIC code simulations for the frequency condition Ω ce /ω pe = 0.1, Ni et al. (2020) found X-O polarized fundamental and harmonic plasma emissions induced by electrostatic UH waves and other electromagnetic modes such as Z and whistler modes.
A large number of numerical studies have been carried out to study the EVDFs formed by kinetic magnetic reconnection utilizing PIC code simulations (see, e.g., Hoshino et al. 2001;Pritchett & Coroniti 2004;Ng et al. 2012;Shuster et al. 2014;Muñoz & Büchner 2016;Yao et al. 2022, and references therein). The most typical found distributions were field-aligned electron beams (Drake et al. 2003;Pritchett & Coroniti 2004;Che et al. 2010Che et al. , 2011, which can cause electromagnetic wave emission. The generation mechanisms of those (parallel) beam distributions are well understood.
They are generated due to the reconnection electric field near the X-point and eventually parallel electric fields near the separatrices (Treumann & Baumjohann 2013).
Several formation mechanisms of those perpendicular crescent-shaped EVDFs have been proposed.
The redistribution of the kinetic energy of the electron motion from the direction parallel to the magnetic field into the motion in the direction perpendicular to the local magnetic field can cause the formation of the perpendicular ring-and crescent-shaped EVDFs (Voitcu & Echim 2012. This process always takes place in kinetic magnetic reconnection since this process naturally occurs in steep magnetic gradients. Voitcu & Echim (2018) demonstrated that as electron beams propagate through a tangential magnetic discontinuity, magnetic gradient drifts can redistribute the energy of the field-aligned electron motion into that of the motion perpendicular to the local magnetic field.
The formation of perpendicular crescent-shaped EVDFs in the diffusion region of reconnection can also be attributed to the meandering motion of non-gyrotropic and non-adiabatic electrons across reconnecting current sheets (Hesse et al. 2014;Bessho et al. 2016;Shay et al. 2016;Lapenta et al. 2017;Bessho et al. 2019). The exact mechanisms are still under debate, but they mostly have to do with the effects of either the E × B or gradient-B drifts (Shay et al. 2016;Bessho et al. 2016;Lapenta et al. 2017).
Many studies have been carried out to explain radio waves/emission generated by unstable EVDFs, which are not necessarily related to magnetic reconnection (Lee et al. 2011;Ganse et al. 2012;Reid & Ratcliffe 2014;Thurgood & Tsiklauri 2015;Zhou et al. 2020;Yao et al. 2021). There have also been studies analyzing radio wave emissions due to magnetic reconnection without focusing on specific features of the EVDFs caused by reconnection or their associated kinetic plasma instabilities (Sakai et al. 2005;Karlický & Bárta 2007). A direct relation of the generation of electromagnetic waves to non-thermal EVDFs generated by kinetic magnetic reconnection is, to the best of our knowledge, still missing. In this paper, we intend to bridge this gap.
Our approach is to directly use the non-thermal EVDFs generated by numerical simulation of kinetic magnetic reconnection, which feature sources of free energy for micro-instabilities, as initial conditions of beam-plasma simulations with the objective to investigate their consequent radio wave emission. This two-step approach allows a higher frequency and wavenumber resolution than that of the immediate diagnostics of radio waves in reconnection simulations. This approach provides thus a better understanding not only of the instabilities caused by non-thermal EVDFs but also of the properties of the resulting radio waves.
For this purpose, we used the results of PIC code simulations of 3D kinetic magnetic reconnection performed in a previous study (Yao et al. 2022). They identified non-thermal EVDFs as potential sources of free energy for the generation of radio waves in those reconnection simulations based on an unsupervised machine learning technique (Dupuis et al. 2020). According to Yao et al. (2022), two characteristic EVDFs that are prone to instabilities relevant to radio emission were identified: two-streaming EVDFs near the separatrices and diffusion region, while perpendicular crescent-shaped EVDFs in the diffusion and the outflow region of reconnection. We used those EVDFs in a parametrized way as initial conditions for separate beam-plasma simulations.
Previous works have also extracted and parametrized EVDFs from reconnection simulation to assess their stability properties via independent 2D simulations. This is because a high frequency and wavenumber resolution are needed to resolve the unstable waves, which is usually not feasible in 3D reconnection simulations. The usual method is to simulate that parametrized EVDF in a homogeneous background plasma in a 2D or 1D geometry. Indeed, Goldman et al. (2008) carried out 1D and 2D kinetic Vlasov simulations to understand the electron holes caused by the Buneman instability, which is triggered by a EVDF taken from 2D PIC simulations of magnetic reconnection.
They used a spatially varying parallel velocity profile to mimic the current sheet width but otherwise constant density. Divin et al. (2012) simulated separatrix instabilities driven by unstable EVDFs from 2D PIC simulations of magnetic reconnection with realistic mass ratio. A step forward of this work is the inclusion of the force balance at the separatrices into the initialization of the EVDF simulations. They neglected, however, gradients parallel to the magnetic field. Note that none of those works aimed at the analysis of the resulting radio emission or velocity distribution functions (VDFs) prone to the electron cyclotron maser instability, so our work attempts to bridge that gap. Hesse et al. (2018) numerically solved the electrostatic dispersion relations of typical EVDFs found in the separatrices of reconnection, with the purpose to assess their stability and role in generating the resulting electrostatic soliton-like waves.
The main objective of this paper is thus to relate the non-thermal EVDFs caused by magnetic reconnection with the resulting electromagnetic radiation due to the kinetic instabilities triggered by those EVDFs and also the ion distribution functions. This is of fundamental importance to understand, for example, the plasma properties of solar coronal reconnection based on the observed radio waves caused by the accelerated electron beams. We address only the initial stage of this process at its source region; a more complete model would require to take into account the transport of electrons as they propagate through the solar corona and solar wind.
The content of this paper is organized as follows: in Section 2 we first describe the numerical simulations of 3D kinetic magnetic reconnection and the resulting characteristic non-thermal EVDFs.
We then describe the numerical setup of simulations with a few of those non-thermal EVDFs as initial conditions. In Section 3 we presented our analysis of the resulting radio waves by these unstable EVDFs. Our conclusions are summarized in Section 4.

3D magnetic reconnection simulation
We first briefly summarized the 3D magnetic reconnection simulation (Yao et al. 2022), from which we obtained the EVDFs to be analyzed in this study. Those simulations were performed with the fully-kinetic 3D Particle-in-Cell (PIC) code ACRONYM , which is appropriate to model the collisionless plasmas of the solar corona as it solves the fully-kinetic Vlasov-Maxwell system of equations.
The initial equilibrium was a double Harris current sheet equilibrium without an external guide field (also known as antiparallel reconnection). Both electrons and ions are initialized following this equilibrium. From now on we will denote those species as Harris population. In addition, we also impose a homogeneous background consisting also of ions and electrons with the same temperature as the Harris population and with a density 0.1625 times that of the peak Harris current sheet density.
A small perturbation was applied to the equilibrium magnetic field in order to accelerate the reconnection onset. The number of grid points of the simulation box was 256 × 512 × 1024, and the physical simulation domain was 4d i × 8d i × 16d i , where d i is the ion skin depth. We apply periodic boundary conditions in all directions. The ion-to-electron mass ratio was m i /m e = 100. Other parameters of this simulation can be seen in the central column of Table 1.
As described in Yao et al. (2022), the evolution of this magnetic reconnection simulation is characterized by reconnection rates that keep increasing almost until the end of the simulated time.
The same can be observed by means of other diagnostics, like the root mean square (RMS) values of the magnetic field, magnetic energy, etc. Our system does not reach a steady state. The reconnection rates reach a value of 0.1 in normalized units (i.e., v A B ∞ /c in CGS units) at t = 5.25 Ω −1 ci . This is the standard value for fast reconnection (Cassak et al. 2017) as well as the value reported by some turbulent reconnection simulations (Wendel et al. 2013;Liu et al. 2018), although it can be higher in some cases, specially when Buneman-like turbulence is triggered in magnetic reconnection (Che et al. 2011;Che 2017;Muñoz & Büchner 2018). Simulations of kinetic turbulence where many sites of magnetic reconnection occur show a distribution of reconnection rates below and above 0.1 (Haggerty et al. 2017). In our case, the time after the reconnection rates reach values higher than 0.1 is dominated by unphysical effects, which are the interaction with the second current sheet in magnetic reconnection (because of periodic boundary conditions) and the depletion of the available magnetic flux (which is due to the reduced simulation box size). Previous 3D magnetic reconnection simulations also usually focused on diagnostics at times when the reconnection rates are near 0.1 (Che et al. 2011;Muñoz & Büchner 2018). In this study, we decided to analyzed EVDFs at t = 5.25 Ω −1 ci when the reconnection rates reach 0.1, same as the analysis of Yao et al. (2022). Table 1. A comparison of parameters of the 3D simulations of magnetic reconnection and 2.5D simulations of a beam-plasma system. For reconnection simulations, some parameters have different values depending on their distance from the current sheet center (as per the Harris current sheet equilibrium). In this case, the first indicated value is the asymptotic value, while the values to the right of the arrows are calculated at the current sheet center. Here λ D is the Debye length, d i the ion skin depth, ω pe the electron plasma frequency calculated at the current sheet center. Ω ce and Ω ci are the electron and ion cyclotron frequencies calculated asymptotically away from the current sheet center, respectively. v A is the Alfvén speed, β e the electron plasma beta, n e the electron number density, ∆x the grid cell size, ∆t the timestep. In this study, we choose two characteristic EVDFs formed at two different locations on a given reconnection plane. The locations are denoted by the points A and B (see red dots) in Figure 1(a0). The characteristic non-thermal EVDFs for FEBs at these locations are: (A) two-streaming EVDF generated in the separatrix region (see Figure 1(a1-a2)) and (B) perpendicular crescent-shaped EVDF formed near the diffusion region (see Figure 1(b1-b2)). Those EVDFs were calculated with the electrons that belong to the so-called Harris population (which initially establishes the initial equilibrium) which, for the purposes of the present work, can be considered as the electron beams. They can provide the source of free energy for electron/electron streaming instabilities depending on their positive velocity gradients, and also for electron/ion streaming instabilities which depend on the relative drift speed between the bulk flow speed of the entire EVDF and the ion distribution. A positive velocity gradient in the parallel EVDFs is a necessary condition for the streaming instabilities, eventually causing Langmuir waves due to inverse Landau damping, as per the so-called Penrose criterion (Penrose 1960). The rigorous necessary and sufficient conditions requires to not only determine the regions with positive velocity gradients, but also to take into account any negative contribution to the instability from the rest of the distribution function.
There is also a background population not shown in Figure 1, which tends to reduce positive velocity gradient(s) due to the beam population and thus suppress streaming-like instabilities. However, the electron beam can propagate through a background with different densities and so this effect will vary point to point within the reconnection region, so we focus mainly on the beam population.
In this 3D antiparallel magnetic reconnection simulation, the magnetic field lines roughly are onto the reconnection plane. Electron beams then move along the magnetic field lines (see white curves in Figure 1(a0)) from the X-point toward the separatrices.
Note that the third dimension of this magnetic reconnection study (in contrast to other 2D reconnection simulations) allows the release of energy of unstable distributions in this direction. More precisely, field-aligned beams along the third direction will generate unstable waves with a wavenumber in this direction. In a pure 2D configuration that process is prohibited, and therefore, a distribution that could be unstable in 3D will be stable in a 2D geometry. As a result, there should be very different waves due to those unstable EVDFs in a 3D simulation in comparison to its 2D counterpart. The unstable EVDFs in 3D magnetic reconnection are observed for relatively long time-scales (i.e. ion time-scales), and thus they are always continuously generated even though their free energy should be depleted due to streaming-like instabilities.
For further details about this simulation and a more comprehensive analysis of EVDFs, see (Yao et al. 2022).

2.5D beam-plasma simulations
We used two characteristic EVDFs described above as initial conditions for independent simulations performed with the 2.5D version of ACRONYM, i.e., a 2D mesh grid in space and the full 3D velocity coordinates. We considered an electron-ion plasma with an ion-to-electron mass ratio m p /m e = 100, same as the reconnection simulations. The plasma consists of three species of particles: background electrons and ions, and beam electrons streaming at a given drift speed. The initial electron plasma frequency is set to be ω pe = 5.0 × 10 9 rad/s, which corresponds to an electron number density of n 0 = 7.9 × 10 9 cm −3 , typical for the solar corona (Aschwanden 2005), and also the same values as the reconnection simulations. The ratio of the electron cyclotron frequency to the electron plasma frequency is Ω ce /ω pe = 0.45. Such frequency ratio is generally found in the separatrices in our 3D kinetic collisionless magnetic reconnection simulation. It can also be found at some locations in the solar corona (Morosan et al. 2016). This implies an electron plasma beta of β e = 0.047.
The size of the simulation box is (L x , L y ) = (N x , N y ) × ∆x = (2048, 2048) × ∆x along the x and y directions, respectively. We set the grid cell size to be ∆x = λ D (with the Debye length of λ D = 0.42 cm). In order to satisfy the Courant-Friedrichs-Lewy (CFL) condition, we imposed the condition c∆t/∆x = 1/2 < 1/ √ 2. Periodic boundary conditions are applied in both directions of the simulation box. The time step is ∆t = 0.035 ω −1 pe . The background magnetic field is assumed to be constant throughout the box, and the direction of the magnetic field defines the x direction of the simulation box, namely, B 0 = B 0 e x . From here on, we refer to the x direction as the parallel direction. We set the thermal speed of the background electron plasma as v the = 0.07c and of   Yao et al. (2022). Nevertheless, all the results obtained in this work are relatively insensitive to small variations in the thermal speed (or temperature) within the range v the /c = 0.05 ∼ 0.1. Table 1 shows a comparison between the parameters of the 3D simulation of magnetic reconnection and those of the 2.5D simulation of beam-plasma system. In contrast to magnetic reconnection simulations, these beam-plasma simulations allow to reach a higher frequency and wavenumber resolution and provide a better understanding of the stability properties and the resulting radio waves. For 2.5D simulation of beam-plasma systems, the size of the simulation box allows a wavenumber resolution of ∆k x = ∆k y = 0.044 ω pe /c. Based on a sampling period of 10∆t, the simulation allows to obtain a frequency up to ω max ≈ 3πω pe in the frequency domain. The presented dispersion relation analysis is based on a time window which allows a frequency resolution of ∆ω = 0.035 ω pe (by 512 samplings) or even higher frequency resolution such as ∆ω = 0.0175 ω pe (by 1024 samplings) or ∆ω = 0.0044 ω pe (by 4096 samplings). For the 3D simulation of magnetic reconnection, the size of the simulation box allows a wavenumber resolution of ∆k x = 0.08ω pe /c, ∆k y = 0.04ω pe /c, ∆k z = 0.02ω pe /c. However, the 3D simulation allows maximum frequency only up to ω max = 0.025ω pe , which is not enough to analyze the high frequency waves such as Langmuir waves and electron cyclotron waves caused by micro-instabilities. Due to the large output period comparable to the reconnection time scale at ion cyclotron time (i.e., 1429∆t = 0.025Ω −1 ci ), only a few outputs are available. This is because of computational limitations: the parallel output (in HDF5 format) of those simulations that utilize thousands of cores tends to slow down the entire simulation in most supercomputers file systems, so that a too often output significantly hinders the speed of the code. Therefore, a 2.5D simulation with high frequency and wavenumber resolutions is necessary to investigate the waves caused by the relevant micro-instabilities that are studied here.
For the simulations in this study, we choose N bg = 100 macro-particles per cell for the background plasma and N bm = 50 macro-particles per cell for the electron beam. This represents a beam-tobackground density ratio of n bm /n bg = 1/2 since the number of macro-particles per cell is proportional to the physical number density, provided a constant ratio of macro-particles to physical particles.
Such a relatively high beam-to-background density ratio can be generated in a local region in the above described 3D magnetic reconnection simulation (see locations denoted by solid red circles in Figure 1(a0)). Several examples of distributions with this beam-to-background density ratio can be found in Yao et al. (2022). We noted that our convergence tests of similar beam-plasma simulations verified that simulations with a larger number of macro-particles per cell for both background and beam plasma basically produce similar results (Yao et al. 2021). In order to reduce the level of numerical noise, a second order shape function was used.
From here on we used the term "momentum" for the momentum per unit mass, which is equivalent to the relativistic velocity in its four-vector form, namely, p = v , p ⊥ = v ⊥ . Note that the fourvelocities v and v ⊥ differ from the ordinary three-velocity by a factor of γ, where the Lorentz factor For the non-thermal EVDFs generated in our kinetic magnetic reconnection simulation, the averaged kinetic bulk flow speed (or drift speed) of beam electrons along the direction of the local magnetic field (i.e., the parallel drift speed) is |u d | < 0.2c. The maximum value 0.2c corresponds to the typical drifts speeds obtained by fitting Maxwellian distributions to our two selected EVDFs found in the separatrices and the diffusion region (see Figure 1(a1,b1)). We also found the excitation of harmonic(s) of plasma emission is sensitively dependent on the parallel drift speed of electron beam and the beam-to-background number density ratio (this will be discussed below). Therefore our beam-plasma simulations are initialized with a beam drift speed of u d = 0.2c. In observations, the drift speed of electron beams is found to be either non-relativistic, 0.2 − 0.5c (Wild et al. 1959;Alvarez & Haddock 1973) or mildly relativistic with > 0.6c (Poquerusse 1994;Klassen et al. 2003).
The drift speed of the electron beams can be deduced from observations of Type III radio bursts as summarized in Reid & Ratcliffe (2014); Reid & Kontar (2018). Meanwhile, the kinetic bulk flow speed perpendicular to the direction of the local magnetic field (i.e., the perpendicular drift speed) of the electron beam with perpendicular crescent-shaped EVDF is set to the same value obtained by the fitting of the EVDFs from Figure 1(b1,b2), i.e., u d⊥ = 0.2c. Both parallel and perpendicular drift speeds u d = u d⊥ = 0.2 c correspond to a Lorentz factor γ = 1.02 and kinetic energy E k = 10.13 keV .
We conducted three simulation runs that are summarized in the Table 2. Among them, Run1 represents a control case with only background plasma. This thermal and homogeneous plasma can excite most of the normal plasma modes (surfaces or curves in the frequency-wavenumber space) due to the fluctuations caused by the thermal noise. In a PIC simulation the physical thermal noise is provided by the macro-particles, where each of them represents a bulk of physical particles.
Therefore the thermal noise level in a PIC simulation is higher than in a real plasma (see, e.g., section 5.3 of Melzani et al. (2013)). In either case the normal plasma modes are excited above the thermal level, but only in the Fourier space (frequency-wavenumber) region where they are either undamped or weakly damped. The excitation of all plasma normal modes is routinely used as a test problem to verify the correctness of PIC-code algorithms. For example, Kilian et al. (2017) showed the presence of most plasma modes in simulations of a homogeneous thermal plasma (i.e., without sources of free energy) with a variety of kinetic codes and plasma models, including a fully-kinetic PIC code.
This way, waves are excited due to the electron beams in Run2 and Run3, which mostly (but not always) occur onto the surfaces or curves (in the frequency-wavenumber space) of the normal plasma modes, can be compared to the waves excited by thermal fluctuations onto those same normal modes in Run1. Note that any wave excited by a plasma instability, like in our Runs2 and Run3, features a spectral power that is orders of magnitude larger than the normal plasma modes excited by the thermal noise in Run1. So the waves by those different processes can be easily distinguished.
We initialized the beam-plasma system by prescribing a particle distribution function The distribution is spatially homogeneously distributed and thus independent on the (x, y) coordinates. The background particle distribution function is expressed as follows, here n 0,bg = N bg (M/∆V ) = n 0 is the background number density, N bg is number of macro-particles per cell of the background plasma, M is the ratio of physical to numerical particles and ∆V is the cell volume, v the is the thermal speed of the background electrons.
The distribution function of the beam electrons is expressed in the following form: here n 0,bm = N bm (M/∆V ) is the beam number density, N bm is number of beam macro-electrons per cell.
The two kinds of EVDFs of FEB for Run2 and Run3 are: The two-streaming EVDF is generated by the product of Eq. (5) and Eq. (6) and the perpendicular crescentshaped EVDF by the product of Eq. (7) and Eq. (8). See details about the specific parameters and the fitting of the reconnection simulations in the text.
where v th and v th⊥ are the thermal speeds along and perpendicular to the direction of the local magnetic field. We consider here only the case of an initially isotropic beam plasma, i.e, v th = v th⊥ = v th,bm .
2. perpendicular crescent-shaped EVDF is the polar angle with the two orthogonal components of the per- Here φ th determines the angular width spread of velocity distribution function about the angle φ 0 in the perpendicular velocity space v ⊥1 − v ⊥2 . The quantity N is the normalization factor and α = √ 2v th⊥ , β = √ 2φ th . In our simulations, φ 0 = 0 and φ th = 0.6π, which empirically correspond to the typical perpendicular crescent-shaped EVDF found in diffusion region (see Figure 1(b2)). Note that our simulations are implemented in Cartesian coordinates: the parallel direction is field-aligned as mentioned For the components of perpendicular velocity, we simply assign v ⊥1 = v y Both beam electrons, with the two-streaming and perpendicular crescent-shaped EVDFs, are initialized by using Eq. (4) (Figure 2). Note that both of these non-thermal EVDFs for FEBs have the two-streaming EVDF along the parallel direction in velocity space (see Figure 2(a1,b1) in v − v ⊥ plane) because fast moving electron beams inherently have a positive parallel drift speed (bulk flow kinetic speed), namely, u d > 0.
The main three beam-plasma simulations to be investigated in this paper and the associated kinetic processes that are expected in them can be summarized as follows: • Run1: There is no beam at all. This is used as a control case. Normal plasma wave modes are generated as a result of the background thermal plasma noise (See response to the Remark 3 above • Run2: We use an electron beam featuring a streaming-like EVDF. The kinetic processes are dominated by instabilities due to the field-aligned Maxwellian electron beam, which in general leads to a relative electron/ion streaming. This relative streaming causes a net current and thus oscillations at the electron plasma frequency. Because of thermal effects, waves at those frequencies experience dispersion and so they are commonly called Langmuir waves. Later during the evolution of the system the same electron/ion streaming leads to an electron/ion streaming instability, causing ion-acoustic waves at much lower frequencies. The existing Langmuir waves can then interact with the generated ion-acoustic waves through a wave-wave interaction process known as the plasma emission mechanism. The resulting electromagnetic waves at the electron plasma frequency (fundamental) and their harmonics can escape the plasma and be remotely observed.
• Run3: We use an electron beam featuring a crescent-shaped EVDF. This beam is also characterized by the same drift speed in the parallel direction to the local magnetic field as the EVDF for Run2. Therefore all kinetic processes taking place in Run2 also occur in this simulation.
In addition, this crescent-shaped EVDF features a positive velocity gradient in the direction perpendicular to the local magnetic field. The crescent-shaped EVDF can cause electron cyclotron maser instabilities and generate multiple harmonic electron cyclotron waves due to a wave-particle interaction at the relativistic resonance condition ω = nΩ ce /γ d + u d k . Here the Lorentz factor γ d is due to the motion of the electron beam, u d the parallel drift speed of the electron beam.

RESULTS
We presented now the results of the simulations described above, which were initialized with either a two-streaming EVDF (see Figure 2(a1,a2)) or a perpendicular crescent-shaped EVDF (see Figure 2(b1,b2)). Those initial EVDFs lead to the following instabilities, respectively, • Electron/ion streaming instabilities. They cause (electrostatic) ion-acoustic waves. In addition, the relative streaming also generates electrostatic Langmuir waves. Transverse emission in the form of harmonics of Langmuir waves is generated via a non-linear wave-wave interaction process.
• Electron cyclotron maser instabilities (ECMIs). They cause multiple harmonic electromagnetic ECWs due to a wave-particle process.

Evolution of the instabilities
Here we briefly describe the evolution of both instabilities to understand the evolution of the beamplasma system and the resulting radiation properties from magnetic reconnection EVDFs. More As time evolves, the kinetic energy of the electron beam is transferred into the thermal energy of the beam-plasma system and electric field energy. Meanwhile, the variation of magnetic energy is negligible (Yao et al. 2021). This way the signatures of the unstable waves due to this instability can be seen mostly in the electric field energy. Figure 3 (a) shows the temporal evolution of ln |max(E x )| 2 (the electric field component parallel to the background magnetic field) for three runs. For the control case Run1 (with electron beam), the fluctuations of the electric field component E x always slightly oscillate in the thermal noise level. The other two runs with a FEB, namely, Run2 and Run3, have initially larger electric field oscillations than those of Run1, and they are also similar to each other.
That indicates that this enhanced initial level of electric field oscillations is only due to the parallel drift speed of the beam, which contributes to the relative streaming between the total EVDF and the ion VDF (initially at rest), as we explain later.
The Maxwellian electron beam of Run2 causes longitudinal electric field oscillations ln |max(E x )| 2 which exponentially grow between t ≈ 400ω −1 pe and t ≈ 600ω −1 pe , with a estimated linear growth rate of 0.012ω pe . This linear growth phase is followed by a saturation, where the electric field energy remains roughly constant. This implies that the source of free energy is depleted.
In order to investigate the specific source of free energy for this observed electric field growth in Run2, and so the instabilities themselves, it is necessary to analyze the distribution functions. Run2. The ion VDF (see the solid grey curve in Figure 4 (a)) is also present at v x = 0. We verified that the ion VDF practically does not vary during the entire evolution of the system.
The initial total EVDF (sum of background and beam electrons) at t = 0.0 does not feature a positive velocity gradient and therefore is stable. But it has a net drift or bulk flow speed v x > 0.
That, combined with the fact that the ion VDF is located at v x = 0, results in a net relative drift speed between the ion VDF and EVDF, equal to V ei = 0.95v the . This relative drift represents a net current and is a possible source of free energy for instabilities and waves. Due to Ampere's law, this causes the finite initial electric field amplitude for Run2 in Figure 3, which represents strong electric field oscillations accompanied by associated oscillatory motions of the EVDF. The initial EVDF only slightly changes during the time interval in which the electric field energy remains constant in Figure 3, i.e., up to t ≈ 400ω −1 pe . The electric field oscillations occur at the plasma frequency modified by thermal effects, which can be identified as Langmuir waves.
After t ≈ 400ω −1 pe the EVDF shifts to the left, i.e., the bulk flow speed of the whole EVDF decreases and adopts negative values in order to decrease the net current (relative electron-ion streaming). As a consequence, the relative drift speed between electrons and ions becomes negative as well, reaching values up to V ei = −0.71v the at t ≈ 525.0ω −1 pe . This makes the electron-ion distribution unstable (see discussion below), since the bulk of the ion VDF falls into the monotonously decreasing part (or to the right side) of the EVDF. This unstable distribution is correlated with the exponential growth phase denoted by the red dashed line of Figure 3. As a consequence the source of free energy, the relative electron-ion drift speed, is exhausted quickly, so that the final EVDF shown in Figure 4(a) for t ≈ 612.5ω −1 pe moves to a location closer to the initial EVDF. This final EVDF also features a larger temperature (broader distribution).
In order to numerically quantify the growth rates of the instability driven by the relative electron-ion drift, we linearize the 1D electrostatic Vlasov equation for three plasma species, ion VDF (subscript i), background EVDF (subscript ec) and beam EVDF (subscript eb). We assume for each one a 1D drifting Maxwellian in the following way for ions and the total electron population, respectively. The parameters are obtained from the fitting of the actual simulation distribution functions at different times. The resulting dispersion relation where s = i, ec or eb. Z(ζ s ) is the plasma dispersion function, with arguments: Here ω ps is the plasma frequency, v th,s the thermal speed and u d,s the drift speed of the s-th species, respectively. In this subsection we normalize our results to the electron plasma frequency ω pe , electron thermal speed v the = k B T e /m e and electron number density n e of the initial (t = 0) electron background plasma (otherwise indicated with the subscript ec).
Note that we solve Eq. (12) for a complex ω + iγ, where ω is the real frequency and γ the growth rate (or damping rate) and for a given k. We search for the complex roots (ω + iγ,k) that satisfy the dispersion relation Eq. (12) by means of a multidimensional root finding algorithm. This method can find all infinite roots of Eq. (12), which are mostly heavily damped. We are only interested in the solutions with positive γ, namely, growing waves.
For this purpose, we first fitted the simulation VDFs shown in Figure 4 by the ion VDF Eq. (10) and the EVDF Eq. (11). Note that the ion VDF does not practically vary during the evolution of the system. The solutions of Eq. (12) show that most of the EVDFs shown in Figure 4 lead to stable solutions, i.e., negative or zero γ. Only two (fitted) EVDFs, namely, at t = 437.5ω −1 pe and t = 525.0ω −1 pe (see dashed curves shown in Figure 4 (a)), lead to a positive growth rate γ. The EVDF at t = 525.0ω −1 pe features a larger source of free energy, namely, a larger electron-to-ion drift speed and so a much larger growth rate than for t = 437.5ω −1 pe . The fitting parameters for this EVDF at t = 525.0ω −1 pe are drift speeds u d,ec = −1.87v the , u d,eb = 0.4v the , electron beam-to-background density ratio n eb /n ec = 1.04, thermal speeds v th,ec = 1.007v the and v th,eb = 1.573v the . The solution of Eq. (12) with those parameters (see skyblue solid curve in Figure 4(b2)) shows a positive growth rate γ > 0 in the wavenumber range about k ∈ [0, 6d −1 e ] with a maximum γ ≈ 0.01ω pe , which is in rough agreement with the value 0.012ω pe obtained from the linear fitting of the electric field amplitude in Figure 3 (a). The difference can be attributed to the dynamical evolution of the EVDF during the linear growth phase; namely, it is hard to choose a representative set of parameters to solve Eq. (12) for the linear theory analysis within this time period.
Figure 4 (b1) shows the real frequency vs the wavenumber of the modes shown in Figure 4 (b2). The unstable branches at t = 437.5ω −1 pe (solid light-green curve) and t = 525.0ω −1 pe (solid skyblue curve) belong to the same wave mode with a mostly linear dependence between frequency and wavenumber.
Its phase speed (or slope) is negative and equal to ω/k ≈ −0.25v the = −0.0175c (at t = 525.0ω −1 pe ). This implies the existence of a resonant interaction between the left tail of the ion VDF and the EVDF at such a phase speed. The source of free energy comes from the positive velocity gradient of the ion VDF combined with the negative velocity gradient of the EVDF. The ion acoustic speed, defined as (Baumjohann & Treumann 1997;Swanson 2003): is c s ≈ ±0.014c, which in roughly agreement with the theoretical phase speed of the unstable wave mode. The differences can be attributed to the non-Maxwellian nature of the EVDF, which in practice leads to a larger effective electron temperature and thus increases the effective sound speed.
Therefore, we can identify the unstable wave mode as ion-acoustic waves (see Section 3.2): strong ion density fluctuations that start to develop during the unstable period after t ≈ 400ω −1 pe , with frequency and wavenumber matching the linear dispersion relation of ion acoustic waves. Note that ion acoustic waves are the essential ingredient, besides Langmuir waves, to generate the fundamental plasma emission via wave-wave interactions (see Section 1).
Those agreements between linear theory calculations and simulation results are evidence that supports the development of an electron/ion instability with the ion-acoustic branch as the unstable mode in our system. This is despite the fact that our system does not satisfy the standard threshold for such an instability, which usually requires T e T i and it is based on a single drifting Maxwellian for the whole population (Gary 1993;Treumann & Baumjohann 2001). For a plasma with equal electron and ion temperatures, like the initial conditions of our simulations, ion-acoustic waves are actually heavily Landau damped, so they are not observed in the early time period (see Figure 5(a)).
It is only later in the evolution of the system (after t ≈ 400ω −1 pe ) that the strong non-Maxwellian EVDF, which is hotter than the initial one, with a heavy skewed tail and interacts with the ion EVDF, can allow the excitation of ion-acoustic waves with an amplitude that exponentially grows.
The linear growth phase of the longitudinal E x component of Run2 is also associated with an exponential growth of the transverse component E y . This is the component responsible for the electromagnetic radiation, which is observed as harmonics of the electron plasma frequency in the electromagnetic wave branch (see Figures 6 and 11). That means the excitation of ion-acoustic waves during the linear growth phase of E x is correlated with the excitation of transverse electromagnetic waves (harmonics) as seen in E y . Note, however, that the estimated linear growth rate is 0.07ω pe and therefore smaller than for E x , indicating that those waves are only indirectly caused by the relative electron-ion drift speed.
Run3 with a crescent-shaped EVDF has a somewhat different evolution. The evolution of the longitudinal electric field component E x also has a linear growth phase, similar to Run2. It starts, however, a bit earlier t ≈ 300ω −1 pe and ends near t ≈ 400ω −1 pe with the same saturation level as that of Run2 and with a similar growth rate as well: 0.013ω pe . That implies that the linear growth phase is also driven by the relative electron-ion streaming, since the parallel EVDFs for both Run2 and Run3 have the same parallel drift speed. The earlier start of the linear phase for the evolution of the E x component of Run3 can be attributed to the perpendicular EVDF, which is different from Run2.
Indeed, the transverse component E y for Run3 grows from the very beginning, saturating at the same time (t ≈ 400ω −1 pe ) as the longitudinal component. This implies that the positive perpendicular velocity gradients generate unstable waves from the very beginning due to the electron cyclotron maser instability, as shown below.

Langmuir waves
The introduction of beam plasma to the background plasma has a significant influence on the local plasma frequency for Run2 and Run3, because the local electron number density is heavily affected by the density fluctuations as the electron beam propagates through the ambient plasma.
By integrating the power spectral density (PSD) of E x over the wavenumber k , yields the power spectra with respect to frequency, namely The local plasma frequency ω loc can then be numerically determined by the corresponding frequencies of those maxima of the P(ω) (see more details in Yao et al. 2021). In this study, the effective local plasma frequency is about ω loc ≈ 0.95 ω pe both for Run2 and Run3, while ω loc = ω pe for Run1.
Harmonics of Langmuir waves at the frequency of multiples of the local plasma frequency, i.e., ω = nω loc (n = 1, 2, 3, ...), due to two-streaming EVDFs are observed both in Run2 and Run3. In this part, we analyze the results of Run2 as an example to explain the generation of Langmuir waves.
Similar phenomena are also observed in Run3.
Due to the small mass ratio µ = m i /m e = 100 used in this study, ion-acoustic waves can be generated within the time scale of our simulations. They start to be visible from t · ω pe ≈ 100 on (see Figure 5). Figure 5 (a) shows the ion number density fluctuation ∆n i = n i − n i within the time interval t · ω pe = 0 − 800 in the x − t plane. Figure 5 (b) shows the power spectral density (PSD) derived from the ion number density n i (k , ω) in Fourier space. The PSD of n i (k , ω) in the low-frequency region is broadly thermally spread near the dispersion relation curve of ion-acoustic waves (Baumjohann & Treumann 1997;Swanson 2003): here c s is the ion-sound speed defined in Eq. (15). Note that this expression contains short-wavelength corrections (contained in the denominator 1 + k 2 λ 2 D ). In the long-wavelength limit (kλ D 1, or kd e 14.28) this expression reduces to ω = c s k. All calculations and results in this work satisfy this approximation very well, since the spectral power is confined to k · d e < 4 as seen in Figure 5(b).
This way, spectral power near this branch implies the existence of low-frequency ion-acoustic waves S. , ω = ω 2 pe + 3v 2 the k 2 (Bohm & Gross 1949). The dispersion relations of R, L, and whistler modes are solved by using the dispersion relation of waves in the cold plasma approximation (Stix 1992;Swanson 2003).
For Run2, we found harmonics of electrostatic Langmuir waves in the parallel direction, i.e., k , (Figure 6(b1)) and harmonics of the transverse waves in the perpendicular direction i.e., k ⊥ , (Figure 6(b2)). In Figure 6(b2), the dispersion relation curves of R, L and Whistler modes are also overlaid. The dispersion relation curves of the electromagnetic R and L modes approach and merge when ω ≥ 2ω loc . The multiple harmonic wave modes with a parallel phase speed equal to or faster than the speed of light in the long-wavelength regime, i.e., the superluminal regime ω ≥ ck , are the transverse waves.
By assuming that the power spectrum of a wave mode follows a Gaussian power distribution along its dispersion relation curve in the k −ω domain, the integrated power distribution for each harmonic of Langmuir waves can be quantitatively estimated as follows: here σ = 0.02ω pe is the frequency width spread along the theoretical dispersion curve.
Figure 6 (c1) shows this integrated power spectra P i (k) versus wavenumber k of harmonics of the electrostatic Langmuir waves for Run2. The power spectrum of the fundamental Langmuir waves decrease from a magnitude of 10 −1 at k = 0 to 10 −7 at k · d e = 12. But most of the power is concentrated between kd e = 0 and kd e = 0.96, both of which represent a local maxima. Meanwhile, the power spectra of the second to fifth harmonic of electrostatic Langmuir waves increase to their maximum and then decrease to a magnitude of 10 −7 at k · d e = 12. As the harmonic number increases, the magnitude of power spectrum of each harmonic electrostatic Langmuir waves decreases significantly, roughly by two orders of magnitude. Similar phenomena are found in transverse emission at frequency of multiple ω loc . Figure 6 (c2) shows the power spectra P(k ) of harmonics of transverse waves at the frequency of multiple local plasma frequencies. They increase to their maximum and then decrease to a magnitude of 10 −11 at k · d e = 12. For fundamental waves, a "hump" in the power distribution curve corresponds to a wave interaction with the L mode at k · d e ≈ 0.48 (see the red curve in Figure 6 (c2)). In general, the humps represent an enhancement of the spectral power at the intersection of wave modes indicating a wave-wave interaction process (see more details at the end of Section 3.4). While for the higher harmonic transverse waves, humps form at k · d e ≈ 1.67, 2.67, 3.68, 4.69 respectively. They are related to a wave interaction with the R mode.
This implies the wave interactions with L and R modes have a significant influence on the power spectra of harmonics of transverse waves.   Table 3. The dashed oblique lines in each panel of Figure 7 indicate the symmetry axis of the angle range where the maximum radiation intensity is theoretically expected to take place. respectively. While for the fifth transverse emission, its radiation intensity is weak in the transverse direction comparing to that of lower harmonic transverse emission, and it mainly appears in the angle We found that the excitation of multiple harmonic plasma emission sensitively depends on the background-to-beam plasma particle number density ratio r = n bg /n bm (or its inverse n bm /n bg ) and (ω = ω loc ) to fifth harmonic (ω = 5ω loc ) of transverse waves, respectively. The radiation intensity is estimated by log 10 |E i (ω, k , k ⊥ )/B 0 | 2 (i = x, y) at given ω = nω loc for Run2 in the time window t · ω pe = 455 − 831, with n = 1, 2, ..., 5 the harmonic number. Here θ = arctan k ⊥ /k is the wave propagation angle. These dashed oblique lines indicate the symmetry axis of the angle region where the maximum radiation intensity of electrostatic Langmuir waves and transverse waves take place (see Appendix A). to-beam plasma particle number density ratio can be up to r = n bg /n bm ≥ 2. The observations show that the drift speed of the non-thermal electron beam is possibly between [0.1c, 0.5c] (Wild et al. 1959;Alvarez & Haddock 1973). Figure 8 shows an empirical parametric regime of parameter pairs (r, u d ) based on 20 additional simulations similar to Run2. Due to limitations in computational resources, those runs are implemented at several discrete parameter pairs (r, u d ) (denoted by orange and green dots in Figure 8), while other parameters are the same as those of Run2. Based on the formation of harmonics of transverse Langmuir waves, we separate the r − u d plane into three regimes: Regime I -where harmonics of transverse emission can be generated, Regime III -where transverse emission is not generated or at most the fundamental emission can be generated, as well as a transition Regime II.
As mentioned above, two-streaming EVDFs with |u d | ≤ 0.2c can be found in the diffusion region and separatrices. In these regions, the background-to-beam particle density ratio is usually n bg /n bm ≥ 2. Figure 8 implies that only when the linear relation u d /c ≥ 0.1n bg /n bm − 0.05 (denoted by the blue dashed line) is roughly satisfied, two-streaming EVDFs generated in magnetic reconnection can cause electron/ion streaming instabilities causing ion acoustic and Langmuir waves and, in this way, generate at least the harmonic plasma emission.
An analytical estimation of the threshold (r, u d ) for the generation of plasma emission, i.e., harmonics of electrostatic/transverse Langmuir waves, can in principle be constructed based on the understanding of energy conversion in the mode conversion process of plasma emission, which is beyond the scope of this study.
The dispersion relations of Bernstein modes are solved from the following equation (Bernstein 1958): where I n (λ) is the modified Bessel function of the first kind with argument λ e = k 2 ⊥ v 2 the /Ω 2 ce , v the is the electron thermal speed of the background plasma. In the short wavelength regime, the frequency of the nth harmonic of Bernstein mode asymptotically approaches to ω → nΩ ce .
The dispersion relation of ECWs due to electron gyroresonance conditions is ω = nΩ ce , this kind of gyroresonance-related electromagnetic waves are named as electron gyroresonance emission if their Electron beam with non-thermal EVDFs offering perpendicular sources of free energy can drive ECMIs, and generate electromagnetic ECWs at relativistic electron gyrofrequency and its harmonics (Melrose 1986), namely, where the Lorentz factor γ d = 1 + u 2 d + u 2 d⊥ /c 2 with u d and u d⊥ the drifting speed of the electron beam. For Run3, γ d ≈ 1.05. Figure 9 (a1-a2,b1-b2) show the power spectral density, i.e., log 10 |E i /B 0 | 2 (i = y, z), in the k ⊥ − ω plane for Run1 (without beam) and Run3 in the time interval t · ω pe = 455 − 634, respectively.
EBWs are observed in PSDs for Run1 and Run3. Multiple harmonic gyroresonance-related ECWs are observed in PSD of E z for Run1, since the local plasma frequency is ω pe ≈ 2.2Ω ce (denoted by white dashed lines in panels of Figure 9), only the third and higher ECWs may escape from the source region where they generated (see Figure 9(a2)). For Run3, a perpendicular crescent-shaped EVDF with a positive gradient in its 1D perpendicular EVDF, i.e., ∂f (v ⊥ )/∂v ⊥ > 0, can offer sources of free energy to cause ECMIs and generate harmonics of electromagnetic ECWs (see Figure 9(b1-b2)).
Similarly, since the effective local plasma frequency is about ω loc ≈ 2.1Ω ce for Run3, only the third and higher harmonic of ECMI-generated ECWs may escape from the source region.
In order to quantitatively estimate the power spectral density of harmonics of ECWs, Eq. (18) is applied to calculate the power spectra for each harmonic of gyroresonance-related for Run1 and ECMI-generated ECWs for Run3 along their dispersion relation curves separately. Figure 9 (c1) shows that the integrated power spectra of each harmonic of the ECMI-generated ECWs for Run3, while ECWs is absent in PSD of E y for Run1. For Run3, the magnitude of power spectra of ECWs is about 10 −7 − 10 −2 , the power spectra of the second to fourth harmonic ECWs are higher than that of the fundamental ECWs due to a wave coupling between EBWs and ECWs. In the long-wavelength regime, interaction between Z mode with the harmonic ECWs and X mode with the third harmonic ECWs enhance corresponding the power spectra of ECWs. The power spectra of the fifth ECWs is comparable to that of the fundamental ECWs. If its power spectra are strong enough, the third and higher harmonic ECWs may escape the source region where they are generated. and ECWs at frequency ω = nΩ ce /γ d (black dashed curves) are overlaid. Note the equivalent Lorentz factor for Run1 (without beam) is γ d = 1. The white dashed line corresponds to local plasma frequency ω = ω loc .
Integrated power spectra P(k ⊥ ) corresponding to the harmonics of ECWs ω = nΩ ce /γ d are calculated in (c1) E y and (c2) E z for Run1 (dashed curves) and Run3 (solid curves) separately. The power spectrum P is normalized by initial magnetic field power density.
Figure 9 (c2) shows power spectra of ECWs derived from PSD of E z . Power spectra of harmonics of the gyroresonance-related ECWs for Run1 are generally by order of 10 −4 less than them of ECMIgenerated ECWs for Run3. In particular, in the long-wavelength regime, power spectra of the second and higher harmonic are higher than that of the fundamental ECWs due to wave-wave interaction between Z and X modes with these ECWs for Run3. While in the short-wavelength regime, power spectra of the harmonic and higher harmonics of ECWs are slightly lower but generally comparable to them of the fundamental ECWs for Run3.  Figure 10. Power spectral density in the k − k ⊥ plane of harmonics of ECWs at ω = 3, 4, 5, 6Ω ce /γ d , respectively. The PSD is estimated by log 10 |E y (ω, k , k ⊥ )/B 0 | 2 at given ω for Run3 in the time interval t · ω pe = 455 − 831. When ω ≥ 4Ω ce /γ d , the L-O mode approaches and merges with the R-X mode (denoted by the red circle at each panel). Figure 10 show PSDs of E y at ω = 3, 4, 5, 6Ω ce /γ d , respectively, in k − k ⊥ plane within time duration t · ω pe = 455 − 831 for Run3. The radiation intensity of the third and fourth harmonic ECWs are nearly isotropic, while that of the fourth harmonic ECWs is slightly more intense along the perpendicular direction at k = 0 (see Figure 10(b)). For the fifth harmonic ECWs, its radiation intensity is slightly stronger along the parallel direction about k ⊥ = 0 at k ·d e = ±2 (see Figure 10(c)), but radiation on the top-half and bottom-half plane are still significant. The radiation intensity of the sixth harmonic ECW is weaker (see Figure 10(d)) than that of the lower harmonics of ECWs.
3.4. Transverse radio waves Figure 6 (b1-b2) show the dispersion relation of transverse waves due to the electron/ion streaming of Run2. Similarly, Figure 9 show the dispersion relation of transverse waves due to the ECMIs of Run3. However, not all those waves can escape the plasma and be eventually remotely detected as radio waves. Only a small region of those diagrams represents radio emission. Figure 11 shows dispersion diagrams of transverse (electromagnetic) waves but constrained to two conditions: (1) waves' frequencies higher than the plasma frequency, i.e. ω > ω pe , (2) superluminal waves, namely, whose phase speeds are larger than the speed of light, i.e., ω/k > c. Those conditions are commonly known as "escaping waves" in the radio emission literature (Melrose 1986). Figure 11 shows the resulting dispersion relation of escaping electromagnetic waves of our simulations, which mostly follow the dispersion relation of the R-X mode with the presence of some harmonics of the plasma frequency for Run 2.
Another common diagnostic for electromagnetic emission in the radio wave range from the Sun is spectrograms (frequency-time diagrams). In particular for solar radio bursts. So in the following we present this diagnostics from our simulations. The caveat of our simulated spectrogram is, of course, that it is only due to the radiation emitted at the source region of radio bursts. Observed spectrograms of radio bursts follow the propagation of electron beams for very long distances compared to our simulations, actually comparable to the distance between the Sun and the Earth. Figure 12 shows the spectrogram of transverse waves derived from E y (t, x) for Run2 and from E z (t, y) for Run3 in the t − ω domain. Note that for Run3 similar results are found in E y (t, x), here we just take results derived from E z (t, y) for example. The PSDs are calculated in the time window t · ω pe = 90 − 800 and with the window size ∆T · ω pe = 179. • Harmonics of transverse waves, e.g., up to 6th harmonic, are generated, even with (PSD of) the 7th and 8th harmonic in a relatively weak level.
• All harmonic components of the transverse waves can be generated in the same sources region at short electron time scale t · ω pe ≤ 500.
• The fundamental and harmonic transverse waves are generated earlier than the third and higher harmonics of plasma waves, e.g., the fundamental and harmonic waves are produced at t · ω pe = 90, while the third to 6th harmonic waves occur after t · ω pe = 500.
Note that Figure  ECWs at frequency of ω = nΩ ce /γ d for Run3. PSD estimated by log 10 |E i (t, ω)/B 0 | 2 (i = y, z), respectively, is displayed in the t − ω plane, and the time window size is ∆T · ω pe = 179. The black dashed lines indicate frequency locations by ω = nω loc and ω = nΩ ce /γ d , separately. The white dashed line corresponds to ω = ω loc .
shows that higher harmonic emission as well as stronger fundamental and first harmonic emission only occur after t · ω pe ≈ 400 − 500. This agrees with the exponential growth of transverse electromagnetic wave energy (seen in Figure 3(b)), as well as with the growth of longitudinal electrostatic wave energy (see in Figure 3(a)) at this same time period. Note that longitudinal electrostatic wave energy includes mainly (forward and backward) propagating Langmuir waves as well as ion-acoustic waves.
The growth of the latter in form of ion density fluctuations can be clearly seen also at the same time (see in Figure 5(a)). All those findings are also in agreement with our linear theory calculations (based on the instantaneous distribution functions) that predict unstable exponentially growing waves in the ion-acoustic branch only after this time period (see Figure 4). In summary, those results indicate that there is a correlation between ion-acoustic waves and generation of electromagnetic wave by the plasma emission mechanism. Note that this is not necessarily a causal relation. It just provides evidence in favor of this process.
An additional evidence for the relation between ion-acoustic waves S and generation of electromagnetic emission is the wave-wave matching (beating) conditions which are a manifestation of the conservation of momentum (for the wavenumber) and energy (for the frequency). Here we focused on one possible pathway of wave-wave interaction involving ion-acoustic waves S and electrostatic Langmuir waves L. First we provide simulation evidence for the following standard interaction: here F represents the fundamental electromagnetic emission at frequency ω = ω pe . This can be interpreted as a decay process of L waves in the form L → S + F (discussed in the introduction), but somehow different because the driver are ion-acoustic waves. So it is not exactly the same process reported before in the literature (see, e.g., Thurgood & Tsiklauri 2015).
The decay process Eq. (21) implies the following beat conditions for frequency and wavenumber: In general Eq. (22) involves all components of the wavenumber vector (see Figure 7 for the angular pattern of wave emission) but for the sake of simplicity we will limit ourselves to just the parallel direction (i.e. x direction).
For ion-acoustic waves S, the wavenumber k S can be obtained from Figure 5 (b).
The region with significant (much higher than the surrounding thermal noise) spectral power along the dispersion relation curve of ion-acoustic waves is roughly confined to k · d e ∈ [0.4, 4] after t · ω pe ≈ 400. An analysis of the integrated 1D spectra k as a function of time (not shown here) shows that the spectral power is of course not constant in time: it broadens from lower to larger wavenumbers as the time evolves. There is a broad spectral peak that also shifts in time, but it is roughly located between k · d e ∈ [0.4, 1.5]. This roughly agrees with the most unstable wavenumbers at the beginning of the exponentially growing phase of the instability (see Figure 4(b2)). Therefore, the wavenumber of ion-acoustic waves is k S · d e ∈ [0.4, 1.5], and the corresponding frequency can be similarly obtained from Figure 5 (b) or be quantitatively estimated by Eq. (17).
For the Langmuir waves L, the wavenumber k L can be obtained from Figure 6(b1) and There is a clear spectral power peak at the intersection of ω = ω pe with the electromagnetic mode (green dashed curve in Figure 6(b2), located at k · d e = 0.48. However, there is still significant power in the wavenumber range k · d e ∈ [0, 1] that comes from the projection to k ⊥ = 0 due to spectral power emitted at oblique angles (see Figure 7(b)).
The spectral power outside the electromagnetic mode but around ω = ω pe , which is not described by linear theory, is a manifestation of its non-linear nature.
Note the wavenumber and frequency resolution for the spectral power of ion density fluctuations (shown in Figure 5 (b)) is ∆k = 0.044d −1 e and ∆ω = 0.0044ω loc (with 4096 samplings in time) separately, while for the power spectral density of electric fields (shown in Figure 6) is ∆k = 0.044d −1 e and ∆ω = 0.035ω loc (with 512 samplings in time) separately.
This way, one possible interaction leading to the fundamental emission as per Eq. (22) and Eq. (23) is: (ω L /ω loc = 1.02 ± 0.035) − (ω S /ω loc = 0.022 ± 0.004) → (ω F /ω loc = 0.98 ± 0.035) Note that this leads to a spectral peak in the electromagnetic wave dispersion relation These somehow symmetric spectral peaks of fundamental emission at k F d e = ±0.48, namely, positive and negative k estimated by Eqs. (26) and (24), resemble the results obtained by Ganse et al. (2012). They are, however, different than the strongly peaked with only positive k fundamental emission observed in the simulations by Thurgood & Tsiklauri (2015). One reason for those discrepancies may lie on the very different parameter set and distribution functions used for these simulations.
Another important remark is that this wave-wave interaction may in principle always occur in a plasma where those L and S waves are excited, even without a beam, like in our Run1. Those modes are just linear plasma waves, but with very low amplitude and, in the case of ion-acoustic waves, highly damped. Therefore, any kind of fundamental emission is negligible in this case. This absence of emission occurs also in the initial stage of our simulation (tω pe < 400), since ion-acoustic waves are damped and only grow later (see Figure 5). The wave-wave interactions efficiently occur after this period (tω pe > 400), as evidenced by the increased power at the fundamental mode in the spectrogram of Similar features for multiple harmonic ECWs due to ECMIs are displayed in Figure 12(b). However, all harmonics of ECWs are simultaneously formed in the source region, i.e., t · ω pe ≤ 90. Higher harmonics of ECWs are generated earlier than their counterparts of plasma emission caused by twostreaming EVDFs. For example, the third to higher harmonic of ECWs can be already generated within t · ω pe = 0 − 500, while their plasma emission counterparts have not yet formed. Different from the ECWs caused by ECMIs, Figure 12(a) shows that the 4th to 6th harmonic of transverse plasma waves (caused by electron/ion streaming instabilities) reach their maximum in t · ω pe = 520 − 700 and are damped significantly after t · ω pe = 700. Because multiple harmonic plasma emission is attributed to a series of multi-stage wave-wave interaction process as discussed before, the energy density of higher harmonic plasma emission comes from the lower harmonics. As a result, the energy density in higher harmonic plasma emission is lower than that in the lower harmonics.

SUMMARY AND CONCLUSIONS
By virtue of fully-kinetic PIC code simulations, we investigated the generation of unstably growing plasma oscillations and the excitation of radio waves due to non-thermal electron velocity space distribution functions (EVDFs) generated by 3D kinetic magnetic reconnection. The reconnection generated non-thermal EVDFs have previously been found in numerical experiments of PIC simulations and by in situ observations.
A fully self-consistent kinetic simulation of radio emissions by magnetic reconnection cannot be achieved due to the still insufficiently low temporal and spatial resolution of 3D kinetic magnetic reconnection simulations. Instead, we first carried out 3D PIC code simulations of kinetic reconnection with a relatively sparse output in time, which is still computationally expensive, to generate the characteristic non-thermal EVDFs (Yao et al. 2022). Using them, we have carried out 2.5D PIC code simulations of beam-plasma interaction, allowing a sufficiently high frequency and wavenumber resolution in the Fourier domain to derive the properties of resulting radio waves. As already explained in (Yao et al. 2022), the magnetic reconnection simulation does not show significant signatures of streaming instabilities, mainly due to the fact that the possible streaming instabilities are relatively weak in comparison with other macroscopic instabilities and plasma flows in magnetic reconnection.
Nevertheless, numerical simulations of a beam-plasma system with higher spectral resolution allows analysis of the possible radio emission caused by streaming instabilities.
In this way, we bridged the gap between the non-thermal EVDFs generated in magnetic reconnection and the formation of radio waves caused by the characteristic non-thermal EVDFs offering source of free energy. Therefore, remote observation of the radio waves, e.g., during solar flares, can be used as a tool for remote diagnostics of magnetic reconnection in astrophysical plasmas.
We first investigated the consequences of multiple harmonic plasma waves caused by two-streaming EVDFs formed along the separatrix and in the diffusion region of 3D kinetic magnetic reconnec-tion, then we investigated the consequences of multiple harmonic electron cyclotron waves caused by perpendicular crescent-shaped EVDFs generated in the diffusion region of 3D kinetic magnetic reconnection (Yao et al. 2022).
The two-streaming EVDFs (Run2) can offer sources of free energy to cause electron/ion streaming instabilities, which can generate ion acoustic waves. The relative electron/ion streaming can also generate Langmuir waves. Those EVDFs are found in the separatrices and diffusion region of kinetic reconnection. Multiple harmonic transverse plasma waves at the frequency of multiples of the local electron plasma, i.e., ω = nω loc , due to non-linear wave-wave interaction processes are generated.
The excitation of multiple harmonic plasma waves is sensitively dependent on the field-aligned drift speed of the electron beam and on the beam-to-background electron number density ratio n bm /n bg .
For small field-aligned electron drift speed, the beam-to-background electron number density ratio n bm /n bg must be large enough in order to excite harmonic(s) of plasma emission. Thus, observations of harmonic emission would allow to constraint the values of both field-aligned drift speed and beamto-background electron number density ratio of the electron beam.
We then investigated perpendicular crescent-shaped EVDFs (Run3) generated near and in the diffusion of reconnection. Such kinds of EVDFs are able to offer sources of free energy to cause electron cyclotron maser instabilities, which can generate multiple harmonic electromagnetic electron cyclotron waves due to a wave-particle interaction. After being amplified by wave-particle interactions, the latter can be isotropically emitted in all directions. Remote observations of electromagnetic electron cyclotron waves can be used for the remote diagnostics of the magnetic field as well as of the ratio of the electron cyclotron frequency to the electron plasma frequency in reconnection regions.
In both cases, the EVDFs of fast electron beams contain a two-streaming EVDF in the field-aligned direction. This results in the formation of multiple harmonic plasma emission.
The higher the harmonic number of the plasma emissions, the closer the distribution of radiation intensity approaches the field-aligned direction. Our findings have a significant practical application on remote diagnostics of reconnection by observations of harmonics of transverse emission, namely, due to the fact that the radiation intensity of harmonic of plasma emission heavily depends on the observation angle and the relative intensity variation in different locations.
Based on the angular pattern of radiation intensity of multiple harmonic plasma emission, our results can also serve as a theoretical testing ground for verifying different processes of the formation of each harmonic plasma emissions (see Appendix A). Harmonics of electrostatic Langmuir waves are excited through a coalescence process L + L n−1 → L n (Yi et al. 2007;Rhee et al. 2009;Yao et al. 2021). Via a process mediated by ion-acoustic waves S, harmonics of transverse emission at a frequency of multiples of the plasma frequency can be generated by non-linear wave-wave interaction processes as follows: • Fundamental transverse emissions F can be generated by a decay process of electrostatic Langmuir waves L and ion-acoustic waves S, i.e., L → S + F (Robinson & Cairns 1994;Reid & Ratcliffe 2014;Melrose 2008).
• Harmonic transverse emissions H can be generated by a coalescence process of electrostatic Langmuir waves L and backward scatted Langmuir waves L via L + L → H (Willes et al. 1996;Yoon 2006;Rhee et al. 2009).
Based on the above discussion, our findings can be used for the remote radio diagnostics of magnetic reconnection by which non-thermal EVDFs are generated and subsequently cause radio waves, in particular in the course of solar flares. Note that radio waves and emission generated directly in the site where magnetic reconnection takes place were observed. For example, observations of fundamental and harmonic plasma emission near the null point of reconnection in solar flares (Chen et al. 2018).
We did not consider propagation effects of the electron beam in the ambient plasma, like in the solar corona or in the solar wind, on the properties of radio waves. They should be taken into account for the sake of a correct interpretation of the observed radio signals. We just investigated the initial stage of the emission process.
Another small caveat is that we are only considering the EVDFs from magnetic reconnection as initial conditions for the beam-plasma simulations, in addition to similar parameters. We do not use other plasma conditions resulting from the magnetic reconnection evolution, e.g., density gradients, partially for simplicity but also because those effects have already been investigated (Yao et al. 2021).
In this same sense, in our beam-plasma simulations we do not initially impose any of the additional waves caused by the reconnection process. As it is well-known, not only gradients in density or bulk flow generated by reconnection can cause waves, but also Hall currents which lead to streaming instabilities, temperature-anisotropy driven instabilities, etc. Those plasma micro-instabilities can finally generate Langmuir waves, whistler waves, kinetic Alfvén waves, etc (Fujimoto et al. 2011;Treumann & Baumjohann 2013;Khotyaintsev et al. 2019). Some of those waves are caused by the EVDFs, but they can also perturb and modify the EVDFs themselves, in particular high-frequency waves with phase speed near its tail via resonant interactions. Both the waves and the distribution function must obey the Vlasov-Maxwell system of equations, which is coupled, and therefore there is always an interplay between electromagnetic fields and distribution functions. One of the most important effects of waves on EVDFs is their apparent perpendicular heating and consequent anisotropic heating, which is different from "real" heating in that it is reversible (Verscharen & Marsch 2011). Since waves and EVDFs are intrinsically linked, it is hard to disentangle the individual effects of the waves caused by the EVDFs from the waves that perturb and contribute to the formation of the EVDFs.
So we decide in the present work to analyze the EVDFs and waves caused by them, not the waves that could actually lead to the formation of EVDFs in reconnection.

ACKNOWLEDGEMENTS
We gratefully acknowledge the developers of the ACRONYM code, the Verein zur Förderung kinetischer Plasmasimulationen e.V. and the financial support by the German Science Foundation (DFG), projects MU-4255/1-1, BU 777/15-1, and BU 777-17-1. We also acknowledge the Sino-German collaboration made possible thanks to the projects with National Natural Science Foundation of China (NSFC) grant 11761131007 and 12003073 and with the Internaltional Scholarship Program of the MOST of China G2021166002L. We also gratefully acknowledge the possibility of using the computing resources of the Max Planck Computing and Data Facility (MPCDF, formerly known as RZG) at Garching and of the Max-Planck-Institute for Solar System Research at Göttingen as well as of the Technical University Berlin, Germany. We thank the referee for their useful remarks that have helped us to improve the original version of this paper.

DECLARATION OF INTERESTS
The authors report no conflicts of interest.

APPENDIX
A. RADIATION PATTERN OF MULTIPLE HARMONIC PLASMA EMISSION Rhee et al. (2009) proposed that for each harmonic of transverse emission the associated wavewave interaction processes can be verified by the angular pattern of its radiation intensity. Based on previous studies of radiation intensity pattern (e.g., see Rhee et al. 2009;Yoon 2006), we extended to talk about angular pattern of radiation intensity of multiple harmonics of plasma emission.
For the fundamental transverse emission F , it is generated by a decay process of electrostatic Langmuir waves L and ion-acoustic waves S, i.e., L → S + F . It is intuitively understood that the angular pattern of fundamental emission should be 90 • different from that of the electrostatic Langmuir waves L.
An theoretical analysis of the amplification of F emission due to L → S + F can be estimated by the transverse wave kinetic equation (Yoon 2005;Yoon et al. 2005;Yoon 2006) as follows: here I L ( k) and I F ( k) are the spectral intensities of electrostatic Langmuir waves L and fundamental emission F separately. Assume isotropic electrostatic Langmuir wave intensities, i.e., I L ( k) ≡ const, then the radiation angular pattern for the fundamental emission F is determined by the transition probability term, i.e., k × k 2 ∼ sin 2 θ. Here θ is the wave phase angle in the k − k ⊥ plane. As a result, the radiation intensity of F emission mainly occurs in angle range about 90 • in the top-half plane and 270 • in the bottom half-plane, as that shown in Figure 7(b). Note that the radiation pattern of the electrostatic langmuir waves is, however, strongly angle dependent, i.e., at angles 0 • and 180 • (see Figure 7(a)). Nevertheless, the isotropic assumption is still valid because radiation of the electrostatic Langmuir waves is very strong elsewhere (Rhee et al. 2009).
For the harmonic emission H, as we mentioned above, it is generated by a coalescence process of beam-generated Langmuir waves L and backward scattered Langmuir waves L , i.e., L + L → H (Willes et al. 1996;Yoon 2006). The amplification of H radiation associated with this coalescence process can be evaluated by the transverse wave kinetic equation (Yoon 2005;Yoon et al. 2005;Yoon 2006) here I H ( k) is the spectral intensities of harmonic emission H. For isotropic electrostatic Langmuir wave intensity, the transition probability term, i.e., k × k 2 |k 2 − 2 k · k | 2 ∼ sin 2 θ cos 2 θ, determines the angular pattern for radiation intensity of the harmonic emission H. The harmonic emission H shows the quadrupole pattern of radiation intensity, which mainly occurs in each quadrant in the k − k ⊥ plane, i.e., in angle range about θ = 45 • , 135 • , 225 • and 315 • , see Figure 7 (c).
For the third and higher harmonic of transverse emission, Rhee et al. (2009) pointed out that they are generated by a coalescence process of electrostatic Langmuir waves L and another transverse wave mode T . An candidate model satisfies this requirement is L + H n−1 → H n , which means the nth harmonic emission H n is generated by a coalescence process of electrostatic Langmuir waves L and transverse emission H n−1 (Zlotnik 1978;Cairns 1988).
We represent the wave vector of the nth transverse emission H n involved in the coalescence process L + H n−1 → H n as k n → (k n , θ n ) and electrostatic Langmuir waves as k L → (k L , θ L ). Here the wave propagation angle θ n is defined as the angle between the direction of local magnetic field and that of wavevector k n .
They satisfy the beat condition k L + k n−1 = k n . Thus the wave vector for L is k L = k n − k n−1 and the magnitude of vector k n − k n−1 lies within the range k n − k n−1 ≤ | k n − k n−1 | ≤ k n + k n−1 (A3) where k n = √ n 2 − 1 · ω pe /c is the magnitude of vector k n , it can be solved from dispersion relation of transverse wave mode for transverse emission H n (Melrose 1986).