Wave Emission of Nonthermal Electron Beams Generated by Magnetic Reconnection

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 PIC code ACRONYM (Kilian et al. 2012), 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.

Table 1. A Comparison of Parameters of the 3D Simulations of Magnetic Reconnection and 2.5D Simulations of a Beam–Plasma System

Parameter	3D Magnetic Reconnection	2.5D Beam–Plasma System
Grid points	256 × 512 × 1024	2048 × 2048
Physical size	$4\times 8\times 16\,{d}_{i}^{3}$	$14.3\times 14.3\,{d}_{i}^{2}$
Δx [λD]	1.56	1
mi /me	100	100
vth,e /c	0.1	0.07
Te [eV]	5110	2503
Ti /Te	1	1
vA /c	0.05 → 0	0.045
βe	0.08 → ∞	0.047
ωpe /Ωce	2 → ∞	2.2
ne [cm−3]	1.3 × 109 → (7.9 + 1.3) × 109	(7.9 + 3.95) × 109
Particles per cell	13 → (80+13)	100 + 50
${\rm{\Delta }}t\ [{\omega }_{{pe}}^{-1}]$	0.087	0.035
Total time $[{{\rm{\Omega }}}_{{ci}}^{-1}]$	9	3.6
Output period $[{\omega }_{{pe}}^{-1}]$	55.5	0.35

Note. 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 time step.

Download table as:  ASCIITypeset image

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 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\ {{\rm{\Omega }}}_{{ci}}^{-1}$. 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 analyze EVDFs at $t=5.25\ {{\rm{\Omega }}}_{{ci}}^{-1}$ when the reconnection rates reach 0.1, the same as those in the analysis of Yao et al. (2022).

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 nonthermal EVDFs for fast-moving electron beams (FEBs) at these locations are as follows: (A) two-stream 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 that 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 require to not only determine the regions with positive velocity gradients but also take into account any negative contribution to the instability from the rest of the distribution function.

Zoom In Zoom Out Reset image size

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 the unstable distributions in this direction. More precisely, the 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 timescales (i.e., ion timescales), 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).

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⁹ rad s⁻¹, which corresponds to an electron number density of n₀ = 7.9 × 10⁹ cm⁻³, 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 condition, we imposed the condition $c{\rm{\Delta }}t/{\rm{\Delta }}x=1/2\lt 1/\sqrt{2}$. Periodic boundary conditions are applied in both directions of the simulation box. The time step is ${\rm{\Delta }}t=0.035\ {\omega }_{{pe}}^{-1}$. 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 ₀ = B₀ 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 beam electrons v_th,bm = 0.08c. This beam thermal speed was calculated by fitting the parallel electron distribution function shown in Figure 1(a1) and determining its standard deviation. The perpendicular distribution exhibits similar values. Because of the strong non-Maxwellian features of the crescent-shaped EVDF shown in Figure (b1)–(b2), a Maxwellian fitting is not meaningful, so we choose to use the same temperature values as those determined for the point A. As for the background population, not shown in Figure 1, we used the thermal speed values resulting from the Maxwellian fitting carried out in Figure 8(b2) 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 0.05 ≤ v_the/c ≤ 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 us to reach a higher frequency and wavenumber resolution and provide a better understanding of the stability properties and the resulting radio waves. For the 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 us to obtain a frequency up to ${\omega }_{\max }\approx 3\pi {\omega }_{{pe}}$ in the frequency domain. The presented dispersion relation analysis is based on a time window that 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 ${\omega }_{\max }=0.025{\omega }_{{pe}}$, which is not enough to analyze the high-frequency waves such as Langmuir waves and ECWs caused by microinstabilities. Due to the large output period comparable to the reconnection timescale at ion cyclotron time (i.e., $1429{\rm{\Delta }}t=0.025{{\rm{\Omega }}}_{{ci}}^{-1}$), 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 microinstabilities that are studied here.

For the simulations in this study, we choose N_bg = 100 macroparticles per cell for the background plasma and N_bm = 50 macroparticles per cell for the electron beam. This represents a beam-to-background density ratio of n_bm/n_bg = 1/2 since the number of macroparticles per cell is proportional to the physical number density, provided a constant ratio of macroparticles 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 the simulations with a larger number of macroparticles 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 use 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 four-velocities v_∥ and v_⊥ differ from the ordinary three-velocity by a factor of γ, where the Lorentz factor is $\gamma =\sqrt{1+({v}_{\parallel }^{2}+{v}_{\perp }^{2})/{c}^{2}}$. The particle kinetic energy is E_k = (γ − 1)m_e c².

For the nonthermal 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 nonrelativistic, 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 the 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 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 macroparticles, 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, the 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, and 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 f( x , v ) in the phase space (x, y) × (v_∥, v_⊥). The distribution is spatially, homogeneously distributed and thus independent on the (x, y) coordinates. The background particle distribution function is expressed as follows,

$$\begin{eqnarray}&&{f}_{\mathrm{bg}}\left(x,y,{v}_{\parallel },{v}_{\perp }\right)={n}_{0,\mathrm{bg}}{f}_{\parallel }\left({v}_{\parallel };{v}_{\mathrm{the}}\right){f}_{\perp }\left({v}_{\perp };{v}_{\mathrm{the}}\right)\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{f}_{\parallel }\left({v}_{\parallel };{v}_{\mathrm{the}}\right)=\displaystyle \frac{1}{\sqrt{2\pi {v}_{\mathrm{the}}^{2}}}\exp \left(-\displaystyle \frac{{v}_{\parallel }^{2}}{2{v}_{\mathrm{the}}^{2}}\right)\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{f}_{\perp }\left({v}_{\perp };{v}_{\mathrm{the}}\right)=\displaystyle \frac{1}{{v}_{\mathrm{the}}^{2}}{v}_{\perp }\cdot \exp \left(-\displaystyle \frac{{v}_{\perp }^{2}}{2{v}_{\mathrm{the}}^{2}}\right),\end{eqnarray} \tag{ 3 }$$

where n_0,bg = N_bg(M/ΔV) = n₀ is the background number density, N_bg is number of macroparticles per cell of the background plasma, M is the ratio of physical to numerical particles, ΔV is the cell volume, and v_the is the thermal speed of the background electrons.

The distribution function of the beam electrons is expressed in the following form:

$$\begin{eqnarray}&&\begin{array}{l}{f}_{\mathrm{bm}}\left(x,y,{v}_{\parallel },{v}_{\perp }\right)\\ \quad =\,{n}_{0,\mathrm{bm}}\cdot {f}_{\parallel }({v}_{\parallel };{u}_{d\parallel },{v}_{\mathrm{th}\parallel }){f}_{\perp }({v}_{\perp };{u}_{d\perp },{v}_{\mathrm{th}\perp }),\end{array}\end{eqnarray} \tag{ 4 }$$

where n_0,bm = N_bm(M/ΔV) is the beam number density, and N_bm is number of beam macroelectrons per cell.

The two kinds of EVDFs of FEB for Run2 and Run3 are as follows:

• 1.  
  Two-stream EVDF.

  $$\begin{eqnarray}&&{f}_{\parallel }({v}_{\parallel };{u}_{d\parallel },{v}_{\mathrm{th}\parallel })=\displaystyle \frac{1}{\sqrt{2\pi {v}_{\mathrm{th}\parallel }^{2}}}\exp \left[-\displaystyle \frac{{\left({v}_{\parallel }-{u}_{d\parallel }\right)}^{2}}{2{v}_{\mathrm{th}\parallel }^{2}}\right]\end{eqnarray} \tag{ 5 }$$

  $$\begin{eqnarray}&&{f}_{\perp }({v}_{\perp };{u}_{d\perp },{v}_{\mathrm{th}\perp })=\displaystyle \frac{1}{{v}_{\mathrm{th}\perp }^{2}}{v}_{\perp }\cdot \exp \left[-\displaystyle \frac{{v}_{\perp }^{2}}{2{v}_{\mathrm{th}\perp }^{2}}\right]\end{eqnarray} \tag{ 6 }$$

  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.

  $$\begin{eqnarray}&&{f}_{\parallel }({v}_{\parallel };{u}_{d\parallel },{v}_{\mathrm{th}\parallel })=\displaystyle \frac{1}{\sqrt{2\pi {v}_{\mathrm{th}\parallel }^{2}}}\exp \left[-\displaystyle \frac{{\left({v}_{\parallel }-{u}_{d\parallel }\right)}^{2}}{2{v}_{\mathrm{th}\parallel }^{2}}\right]\end{eqnarray} \tag{ 7 }$$

  $$\begin{eqnarray}&&\begin{array}{c}{f}_{\perp }({v}_{\perp };{u}_{d\perp },{v}_{{\rm{t}}{\rm{h}}\perp })=\displaystyle \frac{1}{N}{v}_{\perp }\cdot \\ \,\times \exp \left[-\displaystyle \frac{{\left({v}_{\perp }-{u}_{\perp 0}\right)}^{2}}{2{v}_{{\rm{t}}{\rm{h}}\perp }^{2}}-\displaystyle \frac{{\left({\phi }-{{\phi }}_{0}\right)}^{2}}{2{{\phi }}_{{\rm{t}}{\rm{h}}}^{2}}\right]\end{array}\end{eqnarray} \tag{ 8 }$$

  $$\begin{eqnarray}&&\begin{array}{l}N=\sqrt{\pi }\beta \cdot \mathrm{erf}(\pi /\beta )\\ \quad \cdot \left[\displaystyle \frac{\sqrt{\pi }}{2}\alpha {v}_{0}\left[\mathrm{erf}({v}_{0}/\alpha )+1\right]+\displaystyle \frac{{\alpha }^{2}}{2}{e}^{-{v}_{0}^{2}/{\alpha }^{2}}\right]\end{array}\end{eqnarray} \tag{ 9 }$$

  where $\phi =\arctan ({v}_{\perp 2}/{v}_{\perp 1})$ is the polar angle with the two orthogonal components of the perpendicular speed v_⊥1 and v_⊥2, and ${v}_{\perp }=\sqrt{{v}_{\perp 1}^{2}+{v}_{\perp 2}^{2}}$. Here ϕ_th determines the angular width spread of VDF about the angle ϕ₀ in the perpendicular velocity space v_⊥1–v_⊥2. The quantity N is the normalization factor and $\alpha =\sqrt{2}{v}_{\mathrm{th}\perp },\beta =\sqrt{2}{\phi }_{\mathrm{th}}$. In our simulations, ϕ₀ = 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 above; thus v_∥ = v_x . For the components of perpendicular velocity, we simply assign v_⊥1 = v_y and v_⊥2 = v_z .

Both beam electrons, with the two-stream and perpendicular crescent-shaped EVDFs, are initialized by using Equation (4) (Figure 2). Note that both of these nonthermal EVDFs for FEBs have the two-stream EVDF along the parallel direction in velocity space (see Figure 2(a1), (b1) in v_∥–v_⊥ plane) because FEBs inherently have a positive parallel drift speed (bulk flow kinetic speed); namely, u_d∥ > 0.

Zoom In Zoom Out Reset image size

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:

• 1.  
  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 . Those visible waves include the following: Langmuir waves, ordinary (O mode), extraordinary (X-mode) waves, whistler waves, and electrostatic Bernstein waves (EBWs). Note that they have small amplitudes in comparison with the waves generated by instabilities. Waves subject to kinetic damping, such as Langmuir waves, are only visible for large wavelengths, where their damping is weak.
• 2.  
  Run2. We use an electron beam featuring a stream-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, the 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.
• 3.  
  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 ECMIs and generate multiple-harmonic ECWs 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 , which is the drift speed of the electron beam.

Here we briefly describe the evolution of both instabilities to understand the evolution of the beam–plasma system and the resulting radiation properties from magnetic reconnection EVDFs. More detailed stability analyses have been already carried out in the past. See the references in the Introduction and Yao et al. (2021), Zhou et al. (2020). We focus, in particular, on the instabilities driven by the parallel drift speed, since they are the most relevant for the generation of electromagnetic escaping waves.

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 $\mathrm{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.

Zoom In Zoom Out Reset image size

The Maxwellian electron beam of Run2 causes longitudinal electric field oscillations $\mathrm{ln}| \max ({E}_{x}){| }^{2}$, which exponentially grow between $t\approx 400{\omega }_{{pe}}^{-1}$ and $t\approx 600{\omega }_{{pe}}^{-1}$, with an 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. Figure 4(a) shows the time evolution of 1D EVDFs along the parallel direction (i.e., x direction) for Run2. The ion VDF (see the solid gray 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.

Zoom In Zoom Out Reset image size

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 the 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\approx 400{\omega }_{{pe}}^{-1}$. The electric field oscillations occur at the plasma frequency modified by thermal effects, which can be identified as Langmuir waves.

After $t\approx 400{\omega }_{{pe}}^{-1}$ 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\approx 525.0{\omega }_{{pe}}^{-1}$. 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\approx 612.5{\omega }_{{pe}}^{-1}$ 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

$$\begin{eqnarray}&&{f}_{i}({v}_{x})=\displaystyle \frac{{n}_{i}}{\sqrt{2\pi }{v}_{{thi}}}\exp \left(\displaystyle \frac{{\left({v}_{x}-{u}_{d,i}\right)}^{2}}{2{v}_{{thi}}^{2}}\right),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}\begin{array}{rcl}{f}_{e}({v}_{x}) & = & \displaystyle \frac{{n}_{{ec}}}{\sqrt{2\pi }{v}_{\mathrm{th},{ec}}}\exp \left(\displaystyle \frac{{\left({v}_{x}-{u}_{d,{ec}}\right)}^{2}}{2{v}_{\mathrm{th},{ec}}^{2}}\right)\\ & & +\displaystyle \frac{{n}_{{eb}}}{\sqrt{2\pi }{v}_{\mathrm{th},{eb}}}\exp \left(\displaystyle \frac{{\left({v}_{x}-{u}_{d,{eb}}\right)}^{2}}{2{v}_{\mathrm{th},{eb}}^{2}}\right),\end{array}\end{eqnarray} \tag{ 11 }$$

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 yields the following (Gary 1993; Che et al. 2010):

$$\begin{eqnarray}&&1+\sum _{s}\displaystyle \frac{{\omega }_{{ps}}^{2}}{{k}^{2}{v}_{\mathrm{th},s}^{2}}\left[1+{\zeta }_{s}Z({\zeta }_{s})\right]=0\end{eqnarray} \tag{ 12 }$$

where s = i, ec, or eb. Z(ζ_s ) is the plasma dispersion function,

$$\begin{eqnarray}&&Z({\zeta }_{s})=\displaystyle \frac{1}{\sqrt{\pi }}{\int }_{-\infty }^{\infty }\displaystyle \frac{\exp (-{x}^{2})}{x-{\zeta }_{s}}{dx},\qquad {\mathfrak{Im}}({\zeta }_{s})\gt 0\end{eqnarray} \tag{ 13 }$$

with arguments:

$$\begin{eqnarray}&&{\zeta }_{s}=\displaystyle \frac{\omega -{{ku}}_{d,s}}{\sqrt{2}{{kv}}_{\mathrm{th},s}}.\end{eqnarray} \tag{ 14 }$$

Here ω_ps is the plasma frequency, v_th,s is the thermal speed, and u_d,s is the drift speed of the sth species, respectively. In this subsection, we normalize our results to the electron plasma frequency ω_pe , electron thermal speed ${v}_{\mathrm{the}}=\sqrt{{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 Equation (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 Equation (12) by means of a multidimensional root-finding algorithm. This method can find all infinite roots of Equation (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 Equation (10) and the EVDF Equation (11). Note that the ion VDF does not practically vary during the evolution of the system. The solutions of Equation (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{\omega }_{{pe}}^{-1}$ and $t=525.0{\omega }_{{pe}}^{-1}$ (see dashed curves shown in Figure 4(a)), lead to a positive growth rate γ. The EVDF at $t=525.0{\omega }_{{pe}}^{-1}$ 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{\omega }_{{pe}}^{-1}$. The fitting parameters for this EVDF at $t=525.0{\omega }_{{pe}}^{-1}$ 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 Equation (12) with those parameters (see sky-blue solid curve in Figure 4(b2)) shows a positive growth rate γ > 0 in the wavenumber range about ${k}_{\parallel }\in [0,6{d}_{e}^{-1}]$ 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 Equation (12) for the linear theory analysis within this time period.

Figure 4(b1) shows the real frequency versus the wavenumber of the modes shown in Figure 4(b2). The unstable branches at $t=437.5{\omega }_{{pe}}^{-1}$ (solid light-green curve) and $t=525.0{\omega }_{{pe}}^{-1}$ (solid sky-blue 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{\omega }_{{pe}}^{-1}$). 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 is defined as follows (Baumjohann & Treumann 1997; Swanson 2003):

$$\begin{eqnarray}&&{c}_{s}^{2}=\displaystyle \frac{{k}_{B}{T}_{e}+3{k}_{B}{T}_{i}}{{m}_{i}},\end{eqnarray} \tag{ 15 }$$

where 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\approx 400{\omega }_{{pe}}^{-1}$, with the frequency and wavenumber matching the linear dispersion relation of ion-acoustic waves. Note that the 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, the 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\approx 400{\omega }_{{pe}}^{-1}$) that the strong non-Maxwellian EVDF, which is hotter than the initial one, with a heavy skewed tail and interactions with the ion EVDF, can allow the excitation of ion-acoustic waves with an amplitude that exponentially grows.

Zoom In Zoom Out Reset image size

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 that for E_x , indicating that those waves are only indirectly caused by the relative electron-ion drift speed.

Zoom In Zoom Out Reset image size

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\approx 300{\omega }_{{pe}}^{-1}$ and ends near $t\approx 400{\omega }_{{pe}}^{-1}$ 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\approx 400{\omega }_{{pe}}^{-1}$) as the longitudinal component. This implies that the positive perpendicular velocity gradients generate unstable waves from the very beginning due to the ECMI, as shown below.

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. Integrating the power spectral density (PSD) of E_x over the wavenumber k_∥ yields the power spectra with respect to frequency, namely

$$\begin{eqnarray}&&{ \mathcal P }(\omega )=\int | {E}_{x}\left({k}_{\parallel },\omega \right)/{B}_{0}{| }^{2}{{dk}}_{\parallel }.\end{eqnarray} \tag{ 16 }$$

The local plasma frequency ω_loc can then be numerically determined by the corresponding frequencies of those maxima of the ${ \mathcal P }(\omega )$ (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-stream 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 timescale 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 ${\rm{\Delta }}{n}_{i}={n}_{i}-\overline{{n}_{i}}$ within the time interval t · ω_pe = 0 − 800 in the x–t plane. Figure 5(b) shows the 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):

$$\begin{eqnarray}\begin{array}{rcl}{\omega }^{2} & = & \displaystyle \frac{{k}_{B}{T}_{e}}{{m}_{i}}\displaystyle \frac{{k}^{2}}{1+{k}^{2}{\lambda }_{{\rm{D}}}^{2}}+3\displaystyle \frac{{k}_{B}{T}_{i}}{{m}_{i}}{k}^{2}\\ & = & \displaystyle \frac{{c}_{s}^{2}{k}^{2}}{1+3{T}_{i}/{T}_{e}}\left(\displaystyle \frac{1}{1+{k}^{2}{\lambda }_{D}^{2}}+3\displaystyle \frac{{T}_{i}}{{T}_{e}}\right).\end{array}\end{eqnarray} \tag{ 17 }$$

Here c_s is the ion-sound speed defined in Equation (15). Note that this expression contains short-wavelength corrections (contained in the denominator $1+{k}^{2}{\lambda }_{{\rm{D}}}^{2}$). 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, the spectral power near this branch implies the existence of low-frequency ion-acoustic waves S.

Figure 6 shows a comparison of power spectral densities between Run1 and Run2. For Run1, which works as our control case for a thermal plasma without an energetic electron beam, we only see the electrostatic Langmuir mode in the E_x component, while the R, L, and whistler modes are in the E_y component. The electrostatic Langmuir waves are fitted by the standard Bohm–Gross dispersion relation, i.e., $\omega =\sqrt{{\omega }_{{pe}}^{2}+3{v}_{\mathrm{the}}^{2}{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:

$$\begin{eqnarray}&&{{ \mathcal P }}_{i}\left(k\right)=\int \displaystyle \frac{1}{\sqrt{2\pi \sigma }}\exp \left[-\displaystyle \frac{{\left(\omega (k)-{\omega }_{i}\right)}^{2}}{2{\sigma }^{2}}\right]| {E}_{i}(k,\omega ){| }^{2}d\omega \end{eqnarray} \tag{ 18 }$$

where σ = 0.02ω_pe is the frequency width spread along the theoretical dispersion curve.

Figure 6(c1) shows this integrated power spectra ${{ \mathcal P }}_{i}\left(k\right)$ 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⁻¹ at k_∥ = 0 to 10⁻⁷ 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⁻⁷ at k_∥ · d_e = 12. As the harmonic number increases, the magnitude of power spectrum of each harmonic electrostatic Langmuir wave decreases significantly, roughly by two orders of magnitude. Similar phenomena are found in the transverse emission at the frequency of multiple ω_loc. Figure 6(c2) shows the power spectra ${ \mathcal P }({k}_{\parallel })$ 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⁻¹¹ 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.

Figure 7 shows the PSD of E_x for electrostatic Langmuir waves and of E_y for the fundamental and up to the fifth harmonic of transverse waves in the k_∥–k_⊥ plane. The PSD of the electrostatic Langmuir waves and these transverse harmonics are anisotropic, and their radiation intensities are significantly angle dependent. An analytical analysis of the radiation intensity pattern due to the most probable wave–wave interactions for multiple-harmonic transverse emission is discussed in the Appendix, in which we solved those equations for the radiation intensity in dependence on the angle. 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. For example, the radiation intensity of the electrostatic Langmuir waves L and backward, scattered Langmuir waves $L^{\prime}$ (propagating in the negative k_∥ direction) are mainly distributed in a range about the axis of symmetry θ = 0° in the right half-plane k_∥ > 0 and θ = 180° in the left half-plane k_∥ < 0 with a broad angle width spread (Figure 7(a)).

Zoom In Zoom Out Reset image size

For the fundamental transverse mode, the maximum radiation intensity takes place about the angle θ = 90° in the top half-plane and θ = 270° in the bottom half-plane with broad width spread in angle (see Figure 7(b)). For the harmonic transverse emission H (see Figure 7(c)), the maximum wave intensity mainly occurs at each quadrant in the k_∥–k_⊥ plane, e.g., around angles θ = 45°, 135°, 225°, 315°, respectively. The wave radiation intensity of the third transverse emission H₃ (see Figure 7(d)) is observed in an angle range about angles θ = 20°, 160°, 200°, 340° and with a broad spread width in angle about Δθ ∼ 15°. For the fourth transverse emission (see Figure 7(e)), the wave radiation intensity mainly occurs in the ranges θ = 10° ± 5°, 170° ± 5°, 190° ± 5°, 350° ± 5°, 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 range θ = 7° ± 3°, 173° ± 3°, 187° ± 3°, 353° ± 3°, respectively.

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 the field-aligned drift speed of the nonthermal beam u_d∥. We thus carried out numerical experiments to test the dependence of the excitation of plasma emission on the background-to-beam number density ratio r = 2–8 and the parallel drift speed of the nonthermal beam u_d∥ = 0.1–0.5c. Note that the parallel drift speed of the nonthermal electron beams generated in the antiparallel magnetic reconnection simulation presented above is within the range u_d∥ = [ − 0.2c, 0.2c], while the background-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 nonthermal 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.

Zoom In Zoom Out Reset image size

As mentioned above, two-stream 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-stream 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.

In the perpendicular direction k_⊥, the EBWs and ECMI-generated electromagnetic ECWs coexist.

The dispersion relations of Bernstein modes are solved from the following equation (Bernstein 1958):

$$\begin{eqnarray}&&1-\displaystyle \frac{2{\omega }_{\mathrm{loc}}^{2}}{{\lambda }_{e}}{e}^{-{\lambda }_{e}}\sum _{n=1}^{\infty }\displaystyle \frac{{n}^{2}{I}_{n}({\lambda }_{e})}{{\omega }^{2}-{n}^{2}{{\rm{\Omega }}}_{{ce}}^{2}}=0\end{eqnarray} \tag{ 19 }$$

where I_n (λ) is the modified Bessel function of the first kind with argument ${\lambda }_{e}={k}_{\perp }^{2}{v}_{\mathrm{the}}^{2}/{{\rm{\Omega }}}_{{ce}}^{2}$, 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 frequency exceeds the local plasma frequency of the ambient plasma and thus can escape out (Aschwanden 2005).

The electron beam with nonthermal EVDFs offering perpendicular sources of free energy can drive ECMIs, and generate electromagnetic ECWs at the relativistic electron gyrofrequency and its harmonics (Melrose 1986), namely,

$$\begin{eqnarray}&&\omega =\displaystyle \frac{n{{\rm{\Omega }}}_{{ce}}}{{\gamma }_{d}}+{u}_{d\parallel }{k}_{\parallel }\end{eqnarray} \tag{ 20 }$$

where the Lorentz factor ${\gamma }_{d}=\sqrt{1+\left({u}_{d\parallel }^{2}+{u}_{d\perp }^{2}\right)/{c}^{2}}$ with u_d∥ and u_d⊥ as the drifting speed of the electron beam. For Run3, γ_d ≈ 1.05.

Figure 9(a1)–(a2), (b1)–(b2) shows the PSD, i.e., ${\mathrm{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 the 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.

Zoom In Zoom Out Reset image size

In order to quantitatively estimate the PSD of harmonics of ECWs, Equation (18) is applied to calculate the power spectra for each harmonic of the gyroresonance related for Run1 and ECMI-generated ECWs for Run3 along their dispersion relation curves separately. Figure 9(c1) shows the integrated power spectra of each harmonic of the ECMI-generated ECWs for Run3, while ECWs are absent in PSD of E_y for Run1. For Run3, the magnitude of power spectra of ECWs is about 10⁻⁷–10⁻², and 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, the interaction between Z mode with the harmonic ECWs and X mode with the third harmonic ECWs enhances the corresponding power spectra of ECWs. The power spectra of the fifth ECWs are 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.

Figure 9(c2) shows the power spectra of ECWs derived from PSD of E_z . Power spectra of harmonics of the gyroresonance-related ECWs for Run1 are generally by an order of 10⁻⁴ less than those of ECMI-generated ECWs for Run3. In particular, in the long-wavelength regime, the power spectra of the second and higher harmonic are higher than those of the fundamental ECWs due to the wave–wave interaction between Z and X modes with these ECWs for Run3. While in the short-wavelength regime, the power spectra of the harmonic and higher harmonics of ECWs are slightly lower but generally comparable to those of the fundamental ECWs for Run3.

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 the radiation on the top half and bottom half of the 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.

Zoom In Zoom Out Reset image size

Figure 6(b1)–(b2) shows the dispersion relation of transverse waves due to the electron/ion streaming of Run2. Similarly, Figure 9 shows 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 the radio emission.

Figure 11 shows the 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.

Zoom In Zoom Out Reset image size

Another common diagnostic for electromagnetic emission in the radio wave range from the Sun is spectrograms (frequency–time diagrams), in particular, for SRBs. 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.

Zoom In Zoom Out Reset image size

Figure 12(a) shows significant features in the spectrogram of multiple-harmonic plasma emission due to electron/ion streaming instabilities. They are summarized as follows:

• 1.  
  Harmonics of transverse waves, e.g., up to sixth harmonic, are generated, even with (PSD of) the seventh and eighth harmonic in a relatively weak level.
• 2.  
  All harmonic components of the transverse waves can be generated in the same sources region at short electron timescale t · ω_pe ≤ 500.
• 3.  
  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 sixth harmonic waves occur after t · ω_pe = 500.

Note that Figure 12, in combination with our previous discussions, also provides evidence for the role of ion-acoustic waves in the generation of electromagnetic transverse emission. Indeed, this figure 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 the 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 the ion-acoustic waves and the generation of electromagnetic waves by the plasma emission mechanism. Note that this is not necessarily a causal relation. It just provides evidence in favor of this process.

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 the simulation evidence for the following standard interaction:

$$\begin{eqnarray}&&L-S\to F\end{eqnarray} \tag{ 21 }$$

where 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 is 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 in Equation (21) implies the following beat conditions for frequency and wavenumber:

$$\begin{eqnarray}&&{k}_{L}-{k}_{S}={k}_{F}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{\omega }_{L}-{\omega }_{S}={\omega }_{F}.\end{eqnarray} \tag{ 23 }$$

In general Equation (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 the 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 Equation (17).

For the Langmuir waves L, the wavenumber k_L can be obtained from Figure 6(b1) and(c1). The red line in Figure 6(c1) shows that most of the power is concentrated between k_∥ · d_e = 0 and k_∥ · d_e = 0.75 (both of which represent a local maxima), although strong power extends toward significant higher wavenumbers (up to approximately k_∥ · d_e = 3). Note that the spectral power is directionally asymmetric: Langmuir waves with negative k are a bit stronger and extend to larger negative k than Langmuir waves with positive k. The corresponding frequency ω_L can be obtained from Figure 6(b1), or quantitatively estimated by the standard Bohm–Gross dispersion relation mentioned above.

For the fundamental emission, k_F can be obtained from Figure 6(b2) and (c2). 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 nonlinear nature.

Note the wavenumber and frequency resolution for the spectral power of ion density fluctuations (shown in Figure 5(b)) is ${\rm{\Delta }}k=0.044{d}_{e}^{-1}$ and Δω = 0.0044ω_loc (with 4096 samplings in time) separately, while for the PSD of electric fields (shown in Figure 6) ${\rm{\Delta }}k=0.044{d}_{e}^{-1}$ and Δω = 0.035ω_loc (with 512 samplings in time) separately.

This way, one possible interaction leading to the fundamental emission as per Equations (22) and Equations (23) is as follows:

$$\begin{eqnarray}&&\begin{array}{l}({k}_{L}{d}_{e}=0.96\pm 0.04)-({k}_{S}{d}_{e}=1.45\pm 0.04)\\ \quad \to \,({k}_{F}{d}_{e}=-0.48\pm 0.04)\end{array}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&\begin{array}{l}({\omega }_{L}/{\omega }_{\mathrm{loc}}=1.02\pm 0.035)-({\omega }_{S}/{\omega }_{\mathrm{loc}}=0.022\pm 0.004)\\ \quad \to \,({\omega }_{F}/{\omega }_{\mathrm{loc}}=0.98\pm 0.035).\end{array}\end{eqnarray} \tag{ 25 }$$

Note that this leads to a spectral peak in the electromagnetic wave dispersion relation curve with antiparallel propagation (negative k). The ambiguity in this equation is mainly contained in k_S (i.e., one of many possible values), since the spectral power in ion-acoustic waves is broadly distributed (while the Langmuir and fundamental emission peaks are well defined). Frequencies for each wave candidate are ω_L = (1.02 ± 0.035)ω_loc for Langmuir waves, ω_S = (0.022 ± 0.004)ω_loc for ion-acoustic waves, and ω_F =(0.98 ± 0.035)ω_loc for fundamental emission, respectively (see Equation (25)). Note the errors due to the spectral resolution are relatively large, which could explain the slight difference between the calculated ω_F and ω_loc.

The complementary spectral peak (in the dispersion relation curve of fundamental emission, ω = ω_loc) with parallel propagation (i.e., positive k_∥) can be obtained from counter-propagating Langmuir and ion-acoustic waves (often referred in the literature as a coalescence process $L^{\prime} +S\to F$), which have spectral peak values with similar wavenumber magnitudes but negative sign as per Figures 5(b) and 6(b1), namely

$$\begin{eqnarray}&&\begin{array}{l}({k}_{L}{d}_{e}=-0.96\pm 0.04)-({k}_{S}{d}_{e}=-1.45\pm 0.04)\\ \quad \to \,({k}_{F}{d}_{e}=0.48\pm 0.04).\end{array}\end{eqnarray} \tag{ 26 }$$

These somehow symmetric spectral peaks of fundamental emission at k_F d_e = ±0.48, namely, positive and negative k_∥ estimated by Equations (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.

The spectral power at the fundamental emission near k_F d_e = 0 (see Figure 6(2) and (c2)), and in general, in the region in between the peaks at the electromagnetic mode, could also be generated as a consequence of a wave–wave interaction between Langmuir and ion-acoustic waves. For example,

$$\begin{eqnarray}&&\begin{array}{l}({k}_{L}{d}_{e}=0.96\pm 0.04)-({k}_{S}{d}_{e}=0.96\pm 0.04)\\ \quad \to \,({k}_{F}{d}_{e}=0\pm 0.04)\end{array}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\begin{array}{l}({\omega }_{L}/{\omega }_{\mathrm{loc}}=1.02\pm 0.035)-({\omega }_{S}/{\omega }_{\mathrm{loc}}=0.013\pm 0.004)\\ \quad \to \,({\omega }_{F}/{\omega }_{\mathrm{loc}}=1.05\pm 0.035).\end{array}\end{eqnarray} \tag{ 28 }$$

Again, the mismatch between ω_F = (1.05 ± 0.035)ω_loc and ω_loc could be explained by the finite spectral resolution. In general, the fundamental emission in the wavenumber range k_F d_e ∈ [−0.48, 0.48] could originate primarily as a result of similar interactions with ion-acoustic waves in corresponding range k_S d_e ∈ [0.4, 1.5] within which S waves are broadband with higher power. Langmuir waves could have been also chosen at wavelengths different from k_L d_e = ±0.48, since they are broadband all the way up to k_L d_e = 0, leading to similar possible interactions. Frequencies for each wave candidate are ω_L = (1.02 ± 0.035)ω_loc for Langmuir waves, ω_S = (0.013 ± 0.004)ω_loc for ion-acoustic waves, and ω_F = (1.05 ± 0.035)ω_loc for fundamental emission, respectively (see Equation (28)).

Counter-propagating Langmuir waves can be generated by means of the process $L+S\to L^{\prime}$. For example,

$$\begin{eqnarray}&&\begin{array}{l}({k}_{L}{d}_{e}=0\pm 0.04)+({k}_{S}{d}_{e}=-0.96\pm 0.04)\\ \quad \to \,({k}_{L^{\prime} }{d}_{e}=-0.96\pm 0.04)\end{array}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&\begin{array}{l}({\omega }_{L}/{\omega }_{\mathrm{loc}}=1.02\pm 0.035)+({\omega }_{S}/{\omega }_{\mathrm{loc}}=0.013\pm 0.004)\\ \quad \to \,({\omega }_{L^{\prime} }/{\omega }_{\mathrm{loc}}=1.033\pm 0.035).\end{array}\end{eqnarray} \tag{ 30 }$$

It is also clear that a higher (and broadband) spectral power in the S waves implies the generation of counter-propagating Langmuir waves $L^{\prime}$. Equation (30) shows the frequencies for each wave candidate are ω_L = (1.02 ± 0.035)ω_loc for Langmuir waves, ω_S = (0.013 ± 0.004)ω_loc for ion-acoustic waves, and ${\omega }_{L^{\prime} }=(1.033\pm 0.035){\omega }_{\mathrm{loc}}$ for back-scattered Langmuir waves, respectively.

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 Figure 12(a).

The enhanced fundamental electromagnetic emission mode plus forward- and backward-propagating Langmuir waves after t ω_pe > 400 are the seeds for the harmonic (and higher harmonic) electromagnetic emission that appears after this time (see Figure 12(a)), as a consequence of the interaction $L+L^{\prime} \to H$. Since Langmuir waves are weak before this time, H emission is barely observed. Note that finding a combination of waves that corresponds to the spectral peaks of H near the electromagnetic mode is not straightforward, although it could be attributed to the broadband nature of Langmuir waves.

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 two-stream 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 fourth to sixth harmonics 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 the multiple-harmonic plasma emission is attributed to a series of multistage wave–wave interaction processes 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.