High-Frequency Wave Generation in Magnetotail Reconnection: Nonlinear Harmonics of Upper Hybrid Waves.

MMS3 spacecraft passed the vicinity of the electron diffusion region of magnetotail reconnection on 3 July 2017, observing discrepancies between perpendicular electron bulk velocities and E → × B → drift, and agyrotropic electron crescent distributions. Analyzing linear wave dispersions, Burch et al. (2019, https://doi.org/10.1029/2019GL082471) showed the electron crescent generates high-frequency waves. We investigate harmonics of upper-hybrid (UH) waves using both observation and particle-in-cell (PIC) simulation, and the generation of electromagnetic radiation from PIC simulation. Harmonics of UH are linearly polarized and propagate along the perpendicular direction to the ambient magnetic field. Compared with two-dimensional PIC simulation and nonlinear kinetic theory, we show that the nonlinear beam-plasma interaction between the agyrotropic electrons and the core electrons generates harmonics of UH. Moreover, PIC simulation shows that agyrotropic electron beam can lead to electromagnetic (EM) radiation at the plasma frequency and harmonics. ©2019. The Authors.

Introduction

Magnetic reconnection is a fundamental plasma process that occurs ubiquitously in space, astrophysical, and laboratory plasmas. It converts magnetic field energy to kinetic energy of charged particles by dissipative electric fields. The electron diffusion region (EDR) is a key region of the magnetic reconnection process in which dissipative electric fields exist. In this region, electrons are demagnetized, so they do not follow the convection of the magnetic field line. The kinetic physics of electrons plays crucial roles in the EDR (Birn et al., 2001; Vasyliunas, 1975), and waves and electron structures are strongly associated with the reconnection process (Fujimoto et al., 2011; Hesse et al., 2011).

Investigations about electron‐scale phenomena of the EDR require high‐resolution data. The National Aeronautics and Space Administration's Magnetospheric Multiscale (MMS) mission enables observations of multiscale phenomena in the reconnection site (Burch, Moore, et al., 2016), and the four MMS spacecraft with small separations and high data rate provide detailed information about electron‐scale phenomena in or near the EDR. Various waves are observed in the vicinity of the EDR, such as upper‐hybrid (UH) or Langmuir like waves (Burch, Webster, et al., 2018; Graham et al., 2017), whistler waves (Burch, Webster, et al., 2018; Cao et al., 2017; Graham et al., 2016; Huang et al., 2017; Yoo et al., 2018), their nonlinear activities (Wilder et al., 2016, 2017), and highly localized standing whistler waves (Burch, Ergun, et al., 2018). The agyrotropic electron distribution or the crescent‐shaped electron distribution identified by Burch, Torbert, et al. (2016) is one of the most important characteristics of the EDR. Their presence was expected by previous studies (Aunai et al., 2013; Hesse et al., 2014; Scudder & Daughton, 2008), and they can generate reconnection electric field via off‐diagonal electron pressure tensor terms (Hesse et al., 2014; Hesse & Winske, 1994; Lyons & Pridmore‐Brown, 1990). Comparing spacecraft observations and analytic theory, it is known that exclusions by the cutoff energy affects the electron crescent distribution (Egedal et al., 2016; Rager et al., 2018). Moreover, it is shown that the agyrotropic electrons can be the source of free energy for beam‐plasma interactions near the EDR (Burch et al., 2019; Graham et al., 2017). Interacting with the core population, the agyrotropic electrons generate electrostatic UH waves propagating in a perpendicular direction with respect to the ambient magnetic field.

In this letter, we investigate the generation of nonlinear harmonics of UH waves by crescent electrons (agyrotropic electron beam) near an EDR of a tail reconnection event on 3 July 2017, when MMS3 spacecraft was located at [−17.6, 3.1, 2.2] in units of Earth radius (Re) in the geocentric solar ecliptic (GSE) coordinates. We use burst mode data from the fluxgate magnetometer (Russell et al., 2016), the electric field double probes (Ergun et al., 2016; Lindqvist et al., 2016), the fast plasma investigation (Pollock et al., 2016), and the search‐coil magnetometer (SCM; Le Contel et al., 2016). We compare the observation data with a two‐dimensional electromagnetic particle‐in‐cell (PIC) simulation and nonlinear theory of beam‐plasma interactions.

MMS Observation

Figure 1 presents data observed by MMS3 in the GSE coordinate. The total magnetic field is minimum at 05:26:50:730 UT and the field reverses from negative to positive (Figure 1a). The electron density slowly increases from 0.25 to 0.36 cm−3, so the plasma frequency increases from 4.5 to 5.5 kHz (Figure 1b). Near the minimum of the magnetic field, the electron perpendicular temperature increases up to 400 eV (Figure 1c), and there are deviations between perpendicular electron bulk velocities and drift speeds (Figures 1d–1f). Perturbations of electric field are observed before the X‐line crossing as shown in Figure 1g not the X‐line. As will be shown in this letter, they are linked to the perpendicular drift to y direction much faster than drift speeds. Fast plasma investigation observed cigar‐shaped electron distribution in the parallel direction and crescent‐shaped electron distributions in the perpendicular versus plane as shown in Figures 1l–1n.

Figure 1

The electron diffusion region crossing observed by MMS3. (a) Magnetic fields; (b) electron density; (c) electron temperatures; drift speeds and perpendicular electron bulk velocities in (d) x, (e) y, and (f) z direction; and (g) electric field during 05:26:50–52. The most intense waves with nonlinear harmonics are observed during 05:26:50:450–550. In this period, (h) electric fields, (i) spectrogram of the electric field, and (j) magnetic field. Dashed lines in (i) and (j) are harmonics of (the frequency of maximum intensity of the fundamental mode). (k) Phase difference between the fundamental mode and the second harmonic versus the power of second harmonic. Electron distributions in (l) B versus V, (m) B versus , and (n) versus plane when the harmonics are observed. Hodograms in (o) the x‐y and (p) the z‐y plane where black and blue dots are electric and magnetic field data, respectively. The red arrow indicates averaged magnetic field direction. MMS = Magnetospheric Multiscale; GSE = geocentric solar ecliptic; FPI = fast plasma investigation; DES = dual electron spectrometers.

We focus on the period of 05:26:50:450–550 when the most intense electric field perturbations are shown. The maximum amplitude reaches up to ~140 mV/m as seen in Figure 1h, and the oscillation is dominant. The oscillation is finite but small compared to , while is negligible. Figures 1o and 1p present hodograms in x‐y plane and z‐y plane. Black dots are electric field data, blue dots are magnetic field data, and the red arrow indicates the averaged magnetic field direction. The magnetic field is amplified by 25 times in these hodograms. They show that electric field oscillations are in the perpendicular plane with respect to the ambient magnetic field and that the waveform is linearly polarized, which is consistent with an UH wave as reported by Graham et al. (2017). In the spectrogram of the electric field, discrete harmonic structures are observed as shown in Figure 1i. Up to the fifth mode is identified in this figure. The lowest frequency mode (the fundamental mode) is spread over the frequency range between 0.6 and 1, and the maximum power is at 0.88 ( ). We note that the UH frequency can be replaced with the plasma frequency since the ratio between the electron plasma frequency and the electron gyrofrequency is very large (; thus, they are almost the same. Harmonics of are plotted by horizontal dashed lines in Figure 1i, and it seems that higher nth harmonics are located near . The power of second harmonic is on average −15dB of the fundamental, and the second harmonic electric field is aligned with the fundamental within ~1° difference. The phase difference between the fundamental mode and the second harmonics is calculated using their Fourier coefficients for each 3‐ms time window, and they clearly converge to −90° as the power of the second harmonic increases as shown in Figure 1k. Therefore, we conclude that these harmonics are real and not instrumental (Graham et al., 2014). Figure 1k presents the spectrogram of the magnetic field. Because of the SCM noise floor above 1 and that the SCM Nyquist frequency does not cover harmonics, only one mode consistent with the fundamental mode of the electric field is identified.

Simulation

Simulation Conditions

To investigate the generation of the wave and its harmonics, we use a two‐dimensional electromagnetic PIC code (Matsumoto & Omura, 1993) with a periodic simulation domain located in X‐Y plane. We assume homogeneous system because we are interested in local wave generation of high frequency and short wavelength compared to entire reconnection geometry. An ambient magnetic field lies on the X axis, and the agyrotropic beam propagates to the Y axis. Therefore, (X, Y, Z) axes in simulation correspond the GSE (X, Y, Z) axes, respectively. The grid sizes Dx and Dy both are 0.01/, and the time step Dt is 0.005 , where is the electron plasma frequency and c is the speed of light in vacuum. The numbers of grids are 1,024 × 1,024, and the number of time steps is 65,536. The proton‐to‐electron mass is set to the real value 1,836. On average, we use 1,080 particles of the second‐order shape per unit cell, so the total particle number is ~1.1 billion. Plasma parameters are based on the observed values, so they are chosen as follows: densities of the core electron and the crescent electron , where is total electron density, a thermal speed of the core electrons , the ion temperature is five times the electron temperature , the total electron plasma frequency , where is the electron gyrofrequency, a drift speed of the crescent electron , a thermal speed of the crescent electron , and the range of gyro‐angle of the crescent electron Δ = 70°.

As an initial condition, we uniformly allocate isotropic Maxwellian core electrons and ions and crescent‐shaped electrons. An initial velocity distribution of the crescent electrons is given by a beam speed , a thermal speed of the beam , and the gyro‐angle range Δ. The gyro‐phase angle of each crescent electron is defined in the ‐ plane while the axis is set as , and of each crescent electron is limited in the range of −Δ < . In addition, we consider that all particles are transient in order to model the particle motion near the EDR. Core electrons reside in the simulation domain during residence times , which are uniform randomly determined within a range of 0 <  for each particle. A core particle escapes after its own , and a new particle is introduced into the system with new and new velocities. For a crescent electron, before the particle enters into the system, and are uniform randomly determined in the range of −Δ <  and , respectively. Then, the electron enters with and escape when reaches . However, this particle distribution is not stationary. Gyration motions and replacements of particles change the distribution even if there is no effect of and field. To make a stationary distribution, we carry out a presimulation from to using only the ambient magnetic field (without and field calculations). Then, as shown in Figure 2, we obtain a stationary electron distribution with the core and the crescent while all particles are transient.

Figure 2

Electron distribution in the velocity phase space. f in the v‐v space where v (a) at the beginning (ω and (b) the end (ω. The crescent angle ϕ is limited by ±Δ, which are denoted by black dashed lines. (c) f as a function of v. Black solid line is the observation data. Red and blue solid lines are electron distribution of the simulation at the beginning and the end.Time history of (d) total energy; (e) kinetic energies of core and crescent electrons; (f) components of core kinetic energy in x, y, and z direction; (g) y‐directional kinetic energy of the crescent electrons; and (h) electric field energy of x and y direction in a log scale. They are normalized by initial total energy and ω.

Figures 2a and 2b represent electron distributions in the velocity phase space at and . They are 2‐D slices of the 3‐D distribution in ‐ plane where , namely, . The core is Maxwellian, and the crescent electrons are limited by (black dashed line). Not including entire reconnection geometry to focus on the local wave generation, the truncation by the cutoff energy is not implemented. But this initial crescent distribution could correspond to a local observation of drifting Maxwellian truncated by the cutoff energy, as studied by Egedal et al. (2016) and Rager et al. (2018). Figure 2c is a reduced 1‐D distribution in which the and coordinates are integrated. The black line is the MMS observation, and the red and blue lines are simulation data at the beginning and the end . To focus on the beam‐plasma interaction, we ignore the broadened high‐energy part of electron distribution for the initial electron distribution in this study. Thus, we use the same temperature for the core and the crescent beam, as mentioned above. If we consider meandering electrons without any other energization mechanisms, especially an effect of reconnection electric field, the temperature of the inflow core and the crescent would be the same. Therefore, by setting the same temperature, we can analyze the energization of electrons by the beam‐plasma interaction, separately. Details of the simulation will be discussed subsequently. We also carried out another simulation using a beam with higher thermal speed . This choice was made in order to better fit the higher velocity portion of the distribution. However, it did not produce a good agreement in the vicinity of the crescent beam speed . In this case, the beam‐plasma interaction was very weak, and so we do not present the results.

Simulation Results

In Figures 2d–2h, we present an energy history of particles and the electric fields normalized by initial total energy. As shown in Figure 2d, the total energy increases by 1.12% at the end of the simulation, and this energy comes from newly replaced particles. (We checked that the total energy changed by up to 0.1% without particle replacements.) In Figures 2e and 2f, the y‐directional kinetic energy of the crescent decreases, and the electric field energies increase after . They reach the minimum/maximum around then saturate undergoing a mild undulation. In the simulation, the linear growth rate  of the Ey field is 0.086, which matches very well during . Consequently, linear and nonlinear regimes can be separated by .

The kinetic energy of the core electrons (KE‐core) increases from 0.531 to 0.573 during the instability as shown in Figures 2g and 2h. This amount of energy increase is ~7.9% in units of the initial kinetic energy of core electrons. The core energy peaks when the crescent energy (KE‐crescent) is at its minimum. The increase of core energy is larger than the decrease of the crescent energy, which leads to the total energy increasing. The y‐directional core energy (KEy‐core) changes similarly to KE‐core while the z‐directional energy (KEz‐core) is delayed by approximately . Therefore, we could deduce that free energy supplied by crescent electrons is transferred to the core by a wave‐particle interaction mostly in the y direction. Then, due to the gyration, a portion of the KEy‐core changes to the KEz‐core.

Figure 3 shows and power spectra of (Figures 3a–3d) and (Figure 3e and 3f) fields as well as comparisons between the simulation and MMS observation (Figure 3g and 3h) obtained by a fast Fourier transform. The frequency and the wave number axes are normalized by and the electron skin depth , respectively, and the simulation unit is used for the power spectra. In the linear regime between and , as shown in Figure 3a, the beam‐mode grows along the y direction as a result of beam‐plasma interaction. Even though the energy increase is almost linear, the second harmonic grows due to the high amplitude of the fundamental mode. In the nonlinear regime between and , more harmonics are observed up to the fourth as shown in Figure 3b.

Figure 3

Power spectra of E and B fields. (a) E in the linear regime ω and (b) the nonlinear regime ω. Black solid lines are theoretical dispersion curves of harmonics, and the red dashed line is the effective beam speed. (c) E at fundamental (ω = 0.90ω and (d) second harmonics (ω = 1.79ω. (e) B at ω = 0.90ω and (f) at ω = 1.08ω where the radiation is observed. Comparisons of power spectrum of (g) electric field and (h) magnetic. A black line is the Magnetospheric Multisc (MMS) data, and a red line is the simulation result. Dashed vertical lines indicate harmonics of .

To analyze the beam‐plasma interaction, we employ nonlinear theory derived from the nonlinear spectral balance equation (Yoon, 2000; Yoon et al., 2003). Even though the theory explains the 1‐D nonlinear beam‐plasma interaction by a parallel beam or in unmagnetized plasma, it could be applicable to a perpendicular (agyrotropic) beam‐plasma interaction if the plasma is weakly magnetized. The full effect of the magnetic field will be briefly discussed in the subsequent section.

The nonlinear dispersion relation for nth harmonic mode, according to the theory by Yoon et al. (2003), is given by

Here is a normalized wave number defined by , which is used in Figure 3. The third term on the right‐hand side in the equation is the nonlinear shift term estimated from weak turbulence theory. We note that we replace and ω in the equation with an effective beam speed and an effective plasma frequency . Electrons with high∣have the effect of lowering the beam speed in the y direction so the effective beam speed is smaller than . The quantity in the simulation is 0.03c, and we plot this in Figure 3b by a red dashed line.

We note that the theory of Yoon et al. (2003) focuses on the harmonics generation near where n is integer number and near higher‐k region, so they used the Bohm‐Gross dispersion for the fundamental mode. In our case, the fundamental beam mode is generated over a wide frequency range  = 0.5~1.2 and lower‐k region because of high beam density and relatively slow beam speed in comparison to previous studies (Rhee, Woo, & Ryu, 2009; Umeda et al., 2003). However, our fundamental mode locates in the same dispersion curve with the previous studies, because the beam dispersion and the Bohm‐Gross dispersion in higher‐k region make a single dispersion curve based on the linear theory of beam‐plasma instabilities. Consequently, we use where the maximum power exists. The theoretical dispersion curves are plotted in Figure 3b by black solid lines. According to Rhee, Woo, and Ryu (2009), discrepancies between the dispersion curves and the spectral peaks arising in higher frequency and wave number might be attributed to the two‐dimensional effect of the simulation. As a consequence, we conclude that the power spectra of harmonics from the simulation are in good agreement with the theory, which predicts that the growth of harmonic modes should take place around the intersections between the beam mode line and the dispersion curves.

Figures 3c and 3d present power spectra associated with the fundamental mode at and the second harmonic at in space. The maxima of wave power are located in the (perpendicular) y direction while oblique propagations with respect to y direction have a certain amount of power. The spectrum is much broader than so that an ellipse‐shaped spectrum appears. When a drifting Maxwellian beam is used, the power spectrum can be modeled by an arc shape, which is proportional to Exp(cos (Willes et al., 1996; Ziebell et al., 2012); therefore, this ellipse‐shaped spectrum might originate from the crescent‐shaped beam.

Figure 3e shows the power spectra of field at . This magnetic field spectrum derives from oblique components of the fundamental beam mode in Figure 3c. We note the Bx (parallel B) component is dominant that agrees with the UH‐like properties. But we plot the Bz spectrum because it can show the whole spectrum structure. The fundamental long‐wavelength radiation is observed at as shown in Figure 3f in which spectra in the low‐k region is plotted. There are four intensity peaks, which indicate radiations propagating in the ± and ± directions. The most intense peaks are along the ± propagation directions, and their powers are almost the same. These radiations propagating along ± (normal directions with respect to the beam direction) can be explained by the electromagnetic decay process of the Langmuir waves (Edney & Robinson, 1999; Lee et al., 2018; Ziebell et al., 2015). The powers of radiations propagating to ± directions are weaker than ± propagations, and the generation mechanism of ± propagations is not clearly understood at this stage, but somewhat similar pattern is observed in fully numerical solution of the equations of EM weak turbulence theory (Ziebell et al., 2015). We note that higher harmonics of magnetic field spectrum from oblique beam‐mode components and from radio emissions are also observed (not shown here). Magnetic field harmonics of oblique beam‐mode components are observed at the same frequency with electric field beam mode, integer multiples of 0.90, and radio emission harmonics are observed near integer multiples of 1.08.

Figures 3g and 3h show a comparison of electric field and magnetic field spectra between PIC simulation and MMS observation. We obtain the MMS3 spectrum data, which is an averaged power spectral density between 05:26:50.486 and 05:25:50.490 UTC when the powers of fundamental and harmonic modes are at maximum, and we transform the simulation data to real data units. The black line is MMS3 observation, and the red line is PIC simulation, respectively, and vertical dashed lines indicate nth harmonics of . We could barely identify the fifth harmonic, which is too weak to be observable in Figures 1 and 3b. As can be seen in Figure 3g, each peak of the electric field spectrum locates slightly lower than , and the power decreases as the harmonic number n increases. The wave powers and frequencies of electric field peaks obtained from the simulation agree remarkably well with the observational counterparts.

Figure 3h plots the magnetic field data obtained from the simulation (red line) and MMS observation (black line). In the simulation data, the lower peak below comes from obliquely propagating components of nth harmonic beam modes, while the other peak slightly higher than is associated with the radio emission generated by a series of wave‐wave interactions (Rhee, Woo, & Ryu, 2009). The fundamental beam mode agrees between the simulation and the observation; however, other higher beam modes and radio emissions predicted from the simulation are not clearly seen in the MMS observation. Considering that the nominal frequency range of SCM (1 Hz to 6 kHz; Le Contel et al., 2016) does not cover harmonics (>8 kHz), this may come from the instrument limitation of the high‐frequency data.

Conclusion and Discussion

This letter shows the generation of harmonics near an EDR by agyrotropic crescent electrons. Comparing MMS observation, PIC simulation, and theoretical dispersion relation of wave modes, we identify harmonic modes generated by the nonlinear beam‐plasma interaction. Nonlinear beam‐plasma theory explains our study very well, which is applicable to zero guide‐field or weak guide‐field reconnection cases.

Since the theory employed in the present investigation assumes nonmagnetized or weakly magnetized plasmas, our findings may have a limited applicability for highly magnetized cases, however. When the ambient magnetic field is intense, then we could envision a situation where the correlation between the growing ES waves and particles interacting with the mode may be influenced by the gyration motion of electrons before a wave has a chance to grow substantially. Thus, when the growth rate is not high enough, for example, when the beam speed or the density is low, or when the distribution is broadened along the gyro‐phase angle, then for such a situation, the wave growth might be impeded even in the linear regime. We might make a conjecture that the ratio between the growth rate of the beam mode and the gyrofrequency would be a crucial factor for wave evolutions in the strongly magnetized situation, for example, as in the guide field reconnection. As such, even for a highly magnetized case, the present nonlinear beam‐plasma theory might still be applicable if the growth rate of the beam mode is sufficiently higher than the gyrofrequency.

The observation of nonlinear harmonics of beam mode implies that long wavelength electromagnetic radio emissions could be generated as a result of nonlinear processes based on previous beam‐plasma simulations. As shown in the PIC simulation, one probable scenario involves the harmonics of radio emissions near , which is explained by a series of wave‐wave interactions (Lee et al., 2018; Rhee, Ryu, et al., 2009; Ziebell et al., 2015). The fundamental radiation is generated by the electromagnetic decay process [L→S+1H] or induced scattering that involves the eletrons [L→e+1H], while the electrostatic decay process [L→S+L′] leads to a backscattered L′ mode, where L stands for Langmuir waves, S for ion acoustic waves, 1H for fundamental radiations, e denotes electrons, and L′ for backscattered modes. The second harmonics radiation is produced by a wave coalescence process between the Langmuir wave and the backscattered mode [L+L′→2H] (Willes et al., 1996), and, finally, higher harmonics nH are generated by the Zlotnik‐Cairns process between the Langmuir wave and one lower harmonics [L+(n‐1)H→nH] (Cairns, 1988; Zlotnik, 1978). All of these processes are recently confirmed with numerical analysis, which is compared against the simulation (Lee et al., 2018; Ziebell et al., 2015). Another possibility is the mode conversion of Langmuir waves (Hinkel‐Lipsker et al., 1989; Willes & Cairns, 2001) due to the highly inhomogeneous plasma near EDR.

Large amplitude electrostatic waves can in principle efficiently accelerate the particles via wave‐particle interaction. As shown in Figure 2f, the core electrons gain ~7.9% of their initial kinetic energy as a result of the beam‐plasma interaction; hence, UH waves could play a role in energization, diffusion, or scattering of core and crescent electrons near the EDR and affect the reconnection process. However, this process is apparently not sufficient to explain the accelerated electrons in the higher speed range of the MMS3 distribution as presented in Figure 2b. We speculate that these higher‐speed electrons may result from acceleration by the reconnection electric field near the EDR, but more detailed analyses are obviously called for. Gathering from the fact that the beam‐plasma interaction was found to be inefficient in the simulation run with higher thermal speed—which we have carried out but not shown—we conclude that more precise modeling of electron distributions including double Gaussian for the core population, acceleration by the reconnection E‐field, reentering of gyrating electrons, and multiple crescents is called for. Moreover, investigations of kinetic equilibrium in the EDR or in the current sheet geometry could provide where the wave generation and the accelerations via wave‐particle interaction are feasible. In this regard, it is pertinent to note the basis of previous studies of electron acceleration near the EDR (Bessho et al., 2017; Shuster et al., 2015) that modeling double or multiple crescents may be important for the higher speed portion of the observed distribution.

In summary, we have shown that nonlinear harmonics of UH wave occurs near the EDR in the 3 July 2017 event. They are linearly polarized and propagate in the perpendicular direction with respect to the ambient magnetic field. The results of the PIC simulation with transient electrons are consistent with MMS observation. They show that the nonlinear beam‐plasma interaction between the core and the crescent beam generates perpendicular modes at the integer multiple harmonics of the fundamental beam mode. Moreover, PIC simulation results show that agyrotropic electron beams near the EDR can lead to radio emission at the plasma frequency and harmonics.