Turbulent plasmas in space, laboratory experiments, and astrophysical domains can often be described by weak turbulence theory, which can be characterized as a broad spectrum of incoherent interacting waves. We investigate a fundamental nonlinear kinetic mechanism of weak turbulence that can explain the generation of whistler waves in homogeneous plasmas by nonlinear scattering of short wavelength electrostatic lowerhybrid (LH) waves. Two particleincell (PIC) simulations with different mass ratios in two dimensions (2D) were performed using a ring ion velocity distribution to excite broadband LH waves. The wave modes evolve in frequency, and wavenumber space such that the LH waves are converted to whistler waves. The simulations show the formation of quasimodes, which are lowfrequency density perturbations driven by the ponderomotive force due to the beating of LH and whistler waves. These lowfrequency oscillations are damped due to resonant phase matching with thermal plasma particles. By comparing the phase and thermal speeds, we confirm the nonlinear scattering mechanism and its role in the 2D evolution of the ring ion instability. Although the nonlinear scattering is theoretically slower in 2D than in 3D due to the absence of the vector nonlinearity, these simulations show that quasimodes are an important diagnostic for nonlinear landau damping in PIC simulations that has not been utilized in the past. The nonlinear scattering mechanism described here plays an important role in the generation of whistler waves in active experiments, which will be experimentally studied in the upcoming Space Measurement of a Rocket Release Turbulence experiment.
I. INTRODUCTION
Whistler and lowerhybrid (LH) waves can be excited through a variety of physical mechanisms in both space and laboratory plasmas. The most wellstudied mechanisms involve linear instabilities. For example, temperature anisotropies^{1} or heat flux instabilities^{2} can drive whistler waves linearly unstable and heavy ion beams^{3,4} or sheared electron flows^{5} can drive lowerhybrid waves unstable. More recently, observations in magnetic reconnection events indicate the generation of whistlers by mechanisms such as electron beams^{6–8} and temperature anisotropy by LH drift instability (LHDI), whereby the LHDI creates changes in the electrons leading to temperature anisotropies that generate the whistler waves.^{9} In turbulent plasmas, the excitation of whistler waves can also occur through nonlinear mechanisms such as scattering of electrostatic LH waves through particles, also known as nonlinear Landau damping or nonlinear induced scattering.^{10–12} In particular, whistler wave generation by the nonlinear induced scattering mechanism has been studied theoretically,^{12–15} experimentally,^{16} and through simulations.^{12,15,17,18}
Studying the nonlinear generation of whistlers in turbulent plasmas is one of the main objectives of the Space Measurement of a Rocket Release Turbulence (SMART) experiment, which is based on previous active experiments.^{19–21} The SMART experiment^{4,14,15,22} will inject approximately 1.5 kg of Barium (Ba) in Earth's ionosphere at an altitude of $ \u2248 500$ km. The Earth's magnetic field at this altitude is approximately $ 4 \xd7 10 \u2212 5$ T, and the plasma beta parameter^{14} $ \beta e \u2250 2 \mu 0 n e k B T e / B 0 2 \u2243 10 \u2212 5 \u2013 10 \u2212 4$. The neutral Ba particles will expand, ultimately forming ion ring velocity distributions as they are photoionized. The ion rings will be unstable to LH waves,^{4} and these will lead to nonlinear scattered whistler waves.
While there are investigations of LH turbulence with particleincell (PIC) simulations, there has been a lack of studies with both particle and field data analysis that can corroborate the generation of whistler waves by induced scattering. For example, in the studies of Winske and Daughton,^{17,18} a ring of heavy ions in velocity space linearly generates lowerhybrid waves and whistler waves were seen later in the simulation. However, in those simulation studies, a detailed analysis of the role of nonlinear scattering by particles was not considered. In particular, analysis of the electron and ion density fluctuations in kspace as a function of time reveals the presence of field aligned density fluctuations, which are the physical manifestation of socalled “quasimodes” or virtual modes, which are the signature of nonlinear induced scattering. We believe that the observations of these quasimodes are an underutilized tool in identifying nonlinear induced scattering in PIC simulations, and in this manuscript, we use carefully designed 2D simulations to show this.
The purpose of this investigation is the study of the physics of the generation of whistler waves in homogeneous weakly turbulent plasmas. We generate the LH waves through an energetic ion ring velocity distribution and study the impact of nonlinear scattering on the evolution of the instability. We show that the instability is saturated by the scattering of LH to whistler waves. We present detailed analysis for the evolution of the LH and whistler (W) waves. We also study the electron and ion velocity distribution functions (VDFs) and show their evolution in time. The data analysis shows agreement with the whistlers being generated by the nonlinear induced scattering mechanism. Specifically, we identify the quasimodes and the individual constituent waves involved in the nonlinear induced scattering.
In Sec. I A, we review the basic concepts of weak plasma turbulence in order to clarify terminology and to highlight different perspectives. In Sec. II, we describe the methods used to simulate such a process. In Sec. III, we present the results of two case studies. In Sec. IV, we analyzed the results obtained. We then discuss our results in Sec. V, and we conclude with Sec. VI.
A. Brief review of weak plasma turbulence theory and nonlinear induced scattering
Nonlinear induced scattering (NLIS) is a class of nonlinear wave–particle interactions in a weakly turbulent plasma. NLIS is a process that allows transfer of energy from one unstable mode into another (potentially) stable mode and the particles that satisfy the nonlinear Landau resonance condition. For NLIS of LH waves to whistler waves, most of the energy is transferred to the whistler wave, with only a small portion going to the particles. Thus, the energy flows from the linearly unstable mode to a stable mode, effectively saturating the system.^{23}
There is another important ratio^{11} and that is $ \epsilon 2 = \tau a c / \tau int$, which is the ratio of the autocorrelation time $ \tau a c = [ k \Delta ( \omega / k ) ] \u2212 1$ of the waves and the interaction time of particles and waves τ_{int}. Here, $ \Delta \omega $ and $ \Delta k$ are the spectral widths of the incoherent spectrum of the random waves considered. The τ_{int} can be taken to be, for example, the trapping time $ t t r = ( m / ekE ) = \omega b \u2212 1$, where ω_{b} is the bounce frequency. The τ_{ac} is the duration in which a particle “feels” the force by a wave with wavenumber $ k + \Delta k$. Therefore, if $ k \Delta ( \omega / k ) \u226b \omega b$, the wave is not able to capture (trap) the particle, and so $ \epsilon 2 \u226a 1$. In the weak turbulence limit, both $ \epsilon 1 < 1$ and $ \epsilon 2 < 1$. Notice this means that the turbulence can be characterized by a broad spectrum of lowamplitude waves. All other states are characterized as the strong turbulence regime ( $ \epsilon 1 \u2273 1$) where there is a set of large stochastic nonlinear modes. The other case is the largeamplitude single wave regime ( $ \epsilon 2 \u2273 1$) where the analysis used is nonperturbative (see, e.g., Ref. 25).
In this study, we try to simulate the conditions of weakly turbulent plasmas and thereby test the set of equations abovementioned, and also show that whistler waves can be nonlinearly generated in low beta plasmas by this scattering process.
II. PIC SIMULATIONS OF NONLINEAR SCATTERING: SIMULATION SETUP
We use the TristanMP ParticleinCell code.^{34} The TristanMP is a massively parallel (MP) code that has been widely used in simulations of whistlers, e.g., whistler heat flux instability^{35} and magnetosphere whistlermode chorus waves^{36} and shocks.^{37} For this investigation of LH to whistler wave scattering, we have carried out electromagnetic simulations with two spatial dimensions (2D), but we retain the three velocity and field components. We use periodic boundary conditions to simulate an infinite system. The different iontoelectron mass ratios are as follows: $ m i / m e = 100$ and $ m i / m e = 1024$. We load a cold ion ring distribution perpendicular to the plane of the simulation box that should excite LH waves.^{4} The simulation box has dimensions in x–y, with N_{x} = 1944 and N_{y} = 648, number of cells in each direction, for the low mass ratio case. The length ( $ \Delta x$) of each cell is $ \Delta x = d e / 20$, where $ d e \u2261 c / \omega p e$ is the electron skin depth (c is the speed of light). The background magnetic field $ B 0 = B 0 x$ is uniform. The number of particles in each cell (ppc) is 1024; this is $ \u2248 1.29$ × 10^{9} particles for the low mass ratio case.
We load a cold Maxwellian background plasma with ions (i), electrons (e), and a cold energetic ion ring (r), with $ n e = n i + n r$ the total electron density. The time step in both simulations is $ \Delta t = 0.0225 \omega p e \u2212 1$, where ω_{pe} is the electron plasma frequency. The ion ring density is $ n r = 0.3 n i$ and the mass of the ring ions is m_{r} = m_{i}, so $ \alpha = 0.3$. The ratio $ \omega p e / \Omega e = 1.7$ is used throughout.
III. RESULTS
A. Low mass ratio case, $ m i / m e = 100$
For the case where $ m i / m e = 100$, we use the following background plasma temperatures: $ T e = T i = 10$ eV, which implies $ \beta e = \beta i = 1.1 \xd7 10 \u2212 4$ and $ \theta i = 6.26 \xd7 10 \u2212 4 c$, and $ \theta e = 6.26 \xd7 10 \u2212 3 c$. The thermal speed of the ring ion is $ \theta r = 1.98 \xd7 10 \u2212 4 c$, with a mean ring speed of $ v R = 22 \xb7 \theta i$. The speed of the ring is chosen to be higher than in the actual SMART experiment^{4} for computational purpose only. The initial ion and electron distributions are shown in Fig. 1.
As shown in Soto et al.,^{4} the ion ring is unstable to LH waves. Figure 2(a) shows the time history evolution of the average of the electric field energies, which we calculate by summing the square of the field's components over the entire simulation box and dividing it by the dimensions of the box, e.g., $ E x 2 = 1 / ( N x N y ) \u2211 i , j N x , N y  E x [ i , j ]  2$. We see a rapid exponential growth for the E_{y} component, following it by a decay and then a secondary, but smaller, growth phase. The $ E x = E \u2225$ component shows growth and decay but at a much lower level. The E_{z} component, which is purely electromagnetic, does not grow during the linear phase of the instability but as we can see the growth starts near region II, indicating that this is a consequence of nonlinear scattering. Panel (b) shows the time series of the squared magnetic field components.
The nature of E_{y}, as an LH wave, can be established by plotting the fast Fourier transform (FFT) of it in $ k \u2212$space (as shown in Fig. 3). There, panel (a) shows a snapshot of the E_{y} field at the time $ \omega L H t = 23$, marked as region II in Fig. 2(a). The white dashed line in Fig. 3 panels (a) and (c) was obtained analytically by using the double resonance condition,^{4} i.e., by equating $ \omega = k \u22a5 v R = \omega L H 1 + ( m i / m e ) ( k \u2225 / k \u22a5 ) 2$, where $ \omega L H \u2243 0.05 \omega p e$ is the LH frequency oscillation. (Henceforth, we normalize our frequencies to the LH oscillation for the plots only, unless otherwise indicated. Also, the wave numbers are henceforth normalized to $ k \xaf = k d e$, but we drop the bar in k). Panel (b) shows the E_{y} field but in real space at the same time as in (a), notice that the wave is highly oblique. Figure 3—panels (c) and (d)—shows E_{y} now at the later time marked as region III in Fig. 2(a), i.e., $ \omega L H t \u2243 60$. Note the wave becomes more perpendicular (with respect to the $ B 0 x$) at these later times. It is clear from the figures that a broadband spectrum of LH waves has been excited by the ion ring.
On the other hand, returning to the Bfields in Fig. 2(b), we see a rapid exponential growth followed by a secondary peak in B_{x}, the transverse electromagnetic components B_{y} and B_{z} grow quasimonochromatic afterward. These components are the whistler modes. In order to verify this, we plot the FFTs of B_{y}, in ω vs k format, in Fig. 4. The panels (1)–(3) correspond to a window time of $ \Delta \omega L H t = 20$ centered at the regions marked in Fig. 2(b), i.e., at times $ \omega L H t \u2243 50 \xb1 10 , 100 \xb1 10 ,$ and 176 ± 10 (1, 2, and 3, respectively). Note, the whistler dispersion relation is indicated with the white dashed line. Also, panel (4) shows a zoomin region at the maximum whistler power. We see the broadband whistler waves evolving toward maximum power at frequencies just below the LH oscillation, i.e., $ \omega \u2192 \omega w < \omega L H$, in agreement with theory.^{23}
The evolution of the ion ring distribution is shown in Fig. 5. We see a feature of the 2D geometry, whereby the ring loses energy to the wave asymmetrically Fig. 5(b). The second peak in the electric field energy E^{2} is due to the ions rotating back to the plane of the simulation and becoming, again, unstable to the LH waves; this has also been discussed elsewhere.^{17} For completeness, the weak turbulence limit is valid in this simulation as the ratio of the wave energies to the plasma electron energy is: $ W / n e k T e < 1$ (as shown in Fig. 5).
The evolution of the electronVDF as a function of $ v \u2225$ is shown in Fig. 6. There, the vertical dashed lines show the plasma thermal velocities and the ionring speed magnitudes. We see that as the wave grows the eVDF slowly broadens and by $ \omega L H t \u2273 46$ a flattening in the VDF tails is formed. This flattening at suprathermal speeds ( $ v \u2225 \u2273 0.05 c \u226b \theta e$.) is not due to the nonlinear (NL) scattering but due to the linear Landau damping of the LH wave itself (Sec. IV B), and it is because of the low mass ratio used, which allows for larger $ k \u2225$ to satisfy the linear Landau resonance condition at this late times. The NL scattering with the electrons happens at much lower speeds, i.e., at thermal speeds or less, as we will show further below.
IV. ANALYSIS
In the linear exponential regime, the system is governed by the generation of the LH waves as we discussed in the linear theory section. The LH waves are excited by the ring at the double resonance condition. Figure 7 shows the fields E_{ky}, B_{ky}, and E_{kx} (FFTs) and the electron density perturbations n_{ek} at $ \omega L H t \u2248 12$, i.e., at the region marked as I in Fig. 2(a). First, a broadband LH turbulence is generated by the ring. The LH turbulence extends from $ \omega = \omega L H \u2243 0.05 \omega p e$ (at $ k \u2225 = 0$) to $ \omega \u2243 8 \omega L H \u2248 0.4 \omega p e$ (at $ k \u2225 \u2243 15 , \u2009 k \u22a5 \u2243 20$), indicating broadband turbulence: $ \delta \omega = ( \omega max \u2212 \omega min ) \u2273 0.3 \omega p e \u226b \omega L H$.
Third, at about $ \omega L H t \u2273 25$, we begin to see intense parallel electron density perturbations forming in the region $ 10 \u2009 \u2272 \u2009 k \u2225 \u2009 \u2272 \u2009 15$ ( $ k \u22a5 \u223c 0$). These parallel density perturbations are created by the ponderomotive force on the electrons, and they can be seen in Figs. 8 and 9.
The ponderomotive force arises from the beating between the LH waves created by the ionring and the scattered intermediate LH/W waves, and it is a defining feature of the nonlinear induced scattering. The electron nonlinearity is sensitive to the ratio of parallel phase velocity and the electron thermal velocity.^{24,38} When $ \Delta \omega / ( k \u2225 \theta e ) \u2272 1$ ( $ \Delta \omega \u2261 \omega 1 \u2212 \omega 2$ is the low beat mode's frequency), the electrons respond efficiently to the parallel ponderomotive force $ F p \u2225 = F p \xb7 x$. Because the electrons are strongly magnetized, $ k \u22a5 \rho e \u226a 1$, a fluid description can be invoked to highlight the main features of the $ F p \u2225$. In this fluid limit, the nonlinear $ F p \u2225$ is $ \u221d \u2212 m e ( v \xb7 \u2207 ) v \u2225 + q ( v \xd7 B ) \u2225$. Here and below for notational clarity, we do not write out the implied symmetric terms in the ponderomotive force (for more details, see, e.g., Ref. 38).
The beat wave frequency $ \Delta \omega $ should be reflected in the electron (ion) density perturbations, $ n e k \u2009 ( n i k )$. Figure 10—panel (a)—shows the density perturbations, n_{k}, vs time (at a fixed $ k \u2225 = 10.12 , k \u22a5 \u226a 1$). First, there is an exponential growth at $ \omega l h t \u2273 23$ and, later, a decay with a repetitive pattern as time goes on. From the density peaks, we obtain the frequency $ \Delta \omega \u2243 1.2 \xd7 10 \u2212 2$. The associated phase velocity to this lowfrequency oscillation is therefore $ v \varphi n \u2261 \Delta \omega / k \u2225 \u2243 1.2 \xd7 10 \u2212 3$ and thus $ v \varphi n \u2009 \u2272 \u2009 \theta e$, as expected from the theoretical grounds described in Sec. I A.
The value obtained in Eq. (12) is therefore in agreement with the NL density perturbations associated phase speed, since $ v \varphi b = O ( v \varphi n )$, and importantly $ v \varphi b \u2264 \theta e$, which is indicative of the NL scattering.
In Fig. 12(a), we show these beating modes corresponding to the foregoing frequencies together with the density perturbations as a function of time. We can note: (1) the B_{ky} arises later in time with respect to E_{ky}, indicating a scattered mode, and (2) that the density perturbations begin growing due to the interaction of these modes.
Finally, we track the evolution of the density perturbations at $ k \u2225 = 8 , \u2009 k \u2225 = 6$ and $ k \u2225 = 3$ [Fig. 12(b)]. Note that as time progresses, the quasimode at $ k \u2225 \u2243 8$ peaks first with the largest amplitude subsequently the quasimode with $ k \u2225 \u2243 6$ peaks with a smaller amplitude than the first, and so on. We find that the associated phase speeds are $ v \varphi n k = 8 = 1.22 \xd7 10 \u2212 3 , \u2009 v \varphi n k = 6 = 1.14 \xd7 10 \u2212 3 , \u2009 v \varphi n k = 3 = 1.08 \xd7 10 \u2212 3$, respectively. Notice that the quasimodes' phase speed remains approximately the same, $ \u223c O ( 10 \u2212 3 )$, despite the fact that the frequency decreases, and the reason is that k is also decreasing. The quasimodes' phase speeds are on the same order of magnitude as the ion acoustic speed, i.e., $ c s = ( T e + 3 T i ) / ( m e + m i ) \u2243 9 \xd7 10 \u2212 4 \u223c 10 \u2212 3$, indicating that they are ion acoustic modes damped by the electrons.
A. Nonlinear scattering rate
B. Particles
Note that the scattering (and the energy transfer) is toward lower frequencies and long wavelengths. This lower whistler frequency can also be seen in Fig. 4—panel (4), where the maximum whistler power is concentrated just below the LH frequency.
The development of the flattening at times $ \omega L H t \u2273 40$ in the parallel velocity of the eVDF tails at suprathermal speeds is mainly due to the LH waves. To see this, we take two representative values along the double resonance LH curve (white dashed line in Fig. 3), with the following corresponding k's: $ k 1 \u2225 = 1.35 , k 1 \u22a5 = 7.55$, giving the $ \omega 1 = 0.10$, and $ k 2 \u2225 = 2.61 , k 2 \u22a5 = 10.11 , \u2009 \omega 2 = 0.14$. Because the electrons are magnetized, we then have: $ v 1 \varphi = \omega 1 / k 1 \u2225 = 0.10 / 1.35 \u2248 0.075 c$ and $ v 2 \varphi = \omega 2 / k 2 \u2225 = 0.14 / 2.61 \u2248 0.053 c$, in agreement with the flattening region seen in $ v \u2225$ (Fig. 6). Note that this flattening in the eVDF is long after the saturation time t_{s}, first achieved by the nonlinear scattering, as indicated in Fig. 13.
The time energy histories of the fields and particles are shown in Fig. 14(b). There, we see on the top: the ion ring energy E_{r} and bottom: the background electrons and background ions, and in this last panel, the particle energies are normalized to their initial values. We can see that E_{r} has lost 6% of the initial energy at $ \omega L H t \u2243 30$, and approximately 10% toward the end of the simulation. The slight increase after the dip in E_{r} is due to the quasilinear evolution of the ring, and in more realistic boundary conditions, in an open system, the energy of the ring should continue decreasing.^{42}
C. Large mass ratio case, $ m i / m e = 1024$
So far, we used a small mass ratio in order to understand the detailed physics of the evolution of LH turbulence and nonlinear scattering. We now demonstrate the essential physics for a larger mass ratio. For this case, we use the same background plasma parameters and temperatures as before, i.e., $ \beta e = \beta i = 1.1 \xd7 10 \u2212 4$, and same $ v R = 22 \xb7 \theta i$. However, due to the larger ion mass, we have: $ \theta i = 1.96 \xd7 \u2212 4 c$. This makes the system evolve in a similar way. However, there are a few differences. For example, Fig. 15—panel (a) shows the evolution of the field energies (note, $ W / n e T e < 1$). There, we can see that the saturation of the first peak in the electric field energy is at the time $ \omega L H t \u2243 32$, which is at approximately the same saturation point as in the low mass ratio case: $ \omega L H t \u2243 31$. The LH waves are excited with a growth rate of $ \gamma L \u2243 0.14 \omega L H$ (the growth rate was obtained by the best fit to $  E y ( t ) $ calculated over the entire box, not shown). This is in agreement with theory and also very close to the value obtained from the low mass case. The reason for this is that the wave growth rate only depends on the density ratios, as in our case the ringion and background ion masses are equal. The FFT of E_{y} in $ k \u2212$ space [Fig. 15—panel (b)] taken during the early exponential phase ( $ \omega L H t \u2243 14$) shows that the $ k \u2225$ values are not as large as the low mass ratio case; that is, the wave starts off more perpendicular. The LH turbulence now extends from $ \omega = \omega L H \u2243 0.016 \omega p e$ (at $ k \u2225 = 0$) to $ \omega = 10.3 \omega L H \u2248 0.16 \omega p e$ (at $ k \u2225 \u2243 8 , \u2009 k \u22a5 \u2243 25$), as shown in Fig. 15(b).
The fact that the wave is more perpendicular means that the $ k \u2225$ values excited are smaller than in the low mass ratio simulation. Therefore, in order to increase the resolution and see the nonlinear scattering at short k's, we needed to increase the size of the simulation box to: N_{x} = 3744 and N_{y} = 1944. Because of the larger mass, the dimensions of the simulation box in terms of the ion skin depth are $ 6 d i \xd7 3 d i$, which is not that far from the size of the low mass simulation ( $ 9.6 d i \xd7 3.2 d i$), note $ d i = c / \omega p i$ is the ion skin depth. However, in terms of grid size the box is almost twice as large in the parallel direction and three times larger in the perpendicular direction. This made the total number of particles increase from $ \u2243 1.3$ × 10^{9} (low mass ratio case) to $ \u2243 7.5$ × 10^{9} particles.
In this second simulation, the whistler waves are also present. In Fig. 16, the nonlinearly excited whistler waves are shown. There, we can see their evolution toward low frequencies and large wavelengths.
Figure 17 shows the fields and the electron density in the same format as in Fig. 7. There, we can see in panel (a) that at $ \omega L H t \u2243 14$, we have broadband LH wave activity, created by the ring. Note that in the time history evolution of the fields, the $ B y 2$ peaks after the $ E y 2$ component, indicating scattering of $ B y 2$, as expected. At about $ \omega L H t \u2243 39$, panel (b), we begin to see the parallel electron density perturbations in the n_{ek} panel. Note this is past the peak of E^{2} and before saturation of the B^{2}. However, in fact, the n_{ek} and n_{ik} begin growing before this time as we can see in panel (c), where the densities at $ k \u2225 = 5$ are plotted vs time. These are the quasimodes. Panel (d) shows a zoomin in kspace of the n_{ik} showing the parallel density perturbations at lower kvalues. From the peak of the quasimodes, we find the following associated phase velocity $ v \varphi n = 1.96 \xd7 10 \u2212 4$, which implies $ v \varphi n \u2243 c s = 1.94 \xd7 10 \u2212 4 \u2243 \theta i$. This suggests that the quasimodes are now damped by the ions.
V. DISCUSSIONS
Our PIC simulations clearly demonstrate the excitation of LH waves (by an ion ring distribution) and importantly the subsequent excitation of whistler waves by nonlinear wave–particle scattering in a weakly turbulent plasma. From the fields' data, we were able to obtain through detailed dispersion analysis the frequencies and wave numbers that confirm the wave modes in the system. In particular, dispersion analysis in ω vs k space shows how well the data follow the whistler modes' dispersion relation. We also are able to track their evolution as the simulation progresses to ultimately cascading down in frequency to just below the ω_{LH} oscillation, in agreement with theory.^{23}
In our investigation, analysis of the parallel density perturbations reveals the presence of damped oscillations at very low frequencies $ \Delta \omega \u226a \omega L H , \omega w$. These lowfrequency oscillations, or quasimodes, were postulated in the literature analytically,^{11} but here we show how they come out of the simulation. We are able to track their evolution in time and find that in the case of the low mass ratio simulation, their associated phase space velocity is $ v \varphi n \u2243 c s \u2009 \u2272 \u2009 \theta e$. This suggests that the quasimodes are ion acoustic modes damped by the electrons, that is, the scattering is off the electrons. In contrast, in the larger mass ratio case, the analysis suggests that the quasimodes are damped by the ions, $ v \varphi n \u2243 c s \u2243 \theta i$, that is, the scattering is off the ions. In both cases, however, we find $ v \varphi n \u2243 c s$.
The growth rates obtained from our data are in good agreement with theory. As we showed, the nonlinear growth rate obtained balances with the linear growth rate, in agreement with Eq. (13). However, realistic comparisons between theory and simulations for the nonlinear scattering rate, which includes the vector nonlinearity, require a full 3D simulation.
In this study, we are limited by the iontoelectron mass ratios (which make the simulations very costly to run). However, as far as the NLIS process is concerned, it is likely that simulations with realistic mass ratios will not show much different results qualitatively speaking as we have found when we varied the mass ratio from 100 to 1024. As we found in the present study, increasing the mass ratio may require further increase in the size of the simulation box.
Another limitation of all simulations, both 2D and 3D, using periodic boundary conditions is that in reality, most instabilities occur only in some finite region of space. Lowerhybrid turbulence being electrostatic in nature will be confined to the region of instability; however, the scattered electromagnetic whistler waves have large group velocities and will quickly escape. This changes the nature of the saturation of ring instabilities and thus has profound implications for active space experiments such as SMART. This can be modeled by including a convective loss term in the wavekinetic equation for the whistler waves.^{14,26} The resulting system of equations has a predator–preylike oscillatory behavior that is controlled by the convective loss rate. This may be modeled using open boundary conditions, which we will leave to a future publication. The objective of this article was to understand and highlight in detail the physics of NLIS.
Implications. Although the nonlinear scattering in turbulent plasmas has been studied for several decades now, few PIC simulation studies have carried out the analysis in sufficient detail necessary for understanding the interrelated kinetic processes described here. For example, to the best of our knowledge, the quasimodes have never been observed (or properly diagnosed) in a PIC simulation. Previous studies have also not been able to find both the beating modes and the scattering rates in a simulated computer experiment. For example, in the simulations of Winske and Daughton,^{17,18} they concluded that the generation mechanism for the whistlers they saw was due to 3w decay, and missed the important role of the NLIS process in the saturation of the system. Of course 3w decay is always a possibility in turbulent plasmas as we discussed in Sec. I A; however, it is often the case that it is difficult to satisfy the resonant 3w matching conditions because of domain size issues. In our simulations, we did not observe the resonant 3w matching conditions. So, the 3w decay in their simulations could be due to the parameters used, specifically (1) their beta is larger than in our simulations, as they used $ \beta = 6 \xd7 10 \u2212 3$ in the first study^{17,18} and $ \beta = 6 \xd7 10 \u2212 2$ in the second study,^{18} whereas our plasma beta is $ \u223c 1 \xd7 10 \u2212 4$ (which is 60 and 600 smaller than their case, respectively). (2), They did not show the presence of quasimodes. So, without an analysis of the density perturbations in kspace, it is difficult to know if they played a role on the nonlinear whistler excitation via NLIS. Furthermore, with the large mass ratios, they used: $ m r = 4 m i$, with m_{i} = 1024, the LH wave is more perpendicular, and so, the $ k \u2225$ values are much smaller requiring a larger box to possible excite quasimodes at very short k. In upcoming work, we will compare 2D and 3D simulations using the original Winske and Daughton^{17,18} parameters and look for the signatures of nonlinear induced scattering that we have found in this study. Therefore, this simulation study paves the way to understand more complex 3D simulations of weakly turbulent plasmas.
Regarding the implications for active experiments, such as the SMART experiment, we can say that this study confirms one of the basic principles upon which SMART was developed, that is, the conversion of LH to whistler waves in low beta plasmas.
In regard to the generation of LH and whistler waves in low beta plasmas, we think that this mechanism can compete with other established mechanisms for the generation of LH waves (e.g., see Ref. 43 and references therein) and whistler waves in space and solar wind (SW) environments.^{35,44,45} For example, in the turbulent SW at close distances to the Sun, LH waves may be created in solar flares^{27} and our study implies that this could be an efficient method to generate whistler waves by nonlinear scattering through particles in a uniform plasma. The Parker Solar Probe (PSP) is now observing very close to the Sun, and many recent studies indicate the presence of LH and whistler waves detected by PSP.^{43,46}
VI. CONCLUSIONS
The main conclusions are as follows:

The nonlinear generation of whistler waves, in these simulations, was demonstrated via nonlinear induced scattering of shortwavelength LH waves through the background plasma particles in a uniform low beta turbulent plasmas.

Nonlinear plasma theory predicts the formation of quasimodes; here, we showed for the first time how the quasimodes naturally evolve in a 2D PIC simulation. We tracked their evolution and found that the quasimodes can be interpreted as damped modes by the electrons in the low mass ratio case and by ions in the large mass ratio case. Importantly, the quasimodes' phase speed is always at the ion acoustic speed in both simulations.

As the iontoelectron mass ratio is increased, the excited LH waves tend to be more perpendicular. Therefore, in order to resolve the quasimodes at very short k's, the simulation box had to be increased significantly in comparison with the low mass ratio case.

Finally, we proposed that this mechanism could have implications in the solar wind (and perhaps in other astrophysical environments) as an effective mechanism of whistler wave generation.
ACKNOWLEDGMENTS
We would like to thank Dr. Leonid Rudakov for his insights into the investigations of nonlinear induced scattering and also for helpful discussions. This work was supported by the U.S. Naval Research Laboratory (NRL) Base Program and the Defense Advanced Research Projects Agency.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
A. Rualdo SotoChavez: Conceptualization (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Writing – original draft (lead); Writing – review & editing (lead). Chris Crabtree: Formal analysis (equal); Investigation (supporting); Writing – review & editing (equal). Gurudas Ganguli: Funding acquisition (lead); Investigation (supporting); Writing – review & editing (supporting). Alex Fletcher: Investigation (supporting); Writing – review & editing (supporting).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.