Plasma parameter analysis of the Langmuir decay process via Particle-in-Cell simulations

The beam-plasma mechanism, based on the Langmuir decay process, has been proposed to explain naturally enhanced ion-acoustic lines (NEIALs), which are spectral distortions in incoherent scatter radar (ISR) data frequently observed in the vicinity of auroral arcs. In this work the effect of the Langmuir decay process on the ISR spectrum is studied and compared with an analytical model for different plasma parameters by using an electrostatic parallel particle-in-cell (EPPIC) code. Simulations show that the code is working in accordance with theory for a wide range of beam and plasma values and that the features of the spectrum are sensitive to changes of those values. These results suggest that the EPPIC code might be used to build a spectrum-plasma parameter model which will allow estimation of beam and plasma parameters from observed spectra. Simulations also confirm that background electron density (ne) plays an important role in determining the maximum detectable wavenumber of the enhancement. Specifically, results demonstrate that an increase in ne makes the enhancements of the ion acoustic more likely line at large wavenumbers, a finding consistent with statistical studies showing more frequent NEIAL occurrence near solar maximum. Finally, the simulations expose some inaccuracies of the current theoretical model in quantifying the energy passed from the beam to the Langmuir waves as well as with the range of enhanced wavenumbers. These differences may be attributable to the weak Langmuir turbulent regime assumption used in the theory.


Introduction
Incoherent scatter radar (ISR) is a technique used to estimate the plasma parameters of the ionosphere (Evans, 1969).Accurate estimation of plasma density, drift, and temperature, depends on detailed knowledge of the spectral characteristics of the received signal.Under certain circumstances and in certain locations, the ISR spectrum has strongly enhanced, denoted Natural Enhancement of Ion Acoustic Lines (NEIALs) (Sedgemore-Schulthess and St Maurice, 2001).Although standard ISR fitting techniques are not applicable for these cases, these spectral anomalies are a valuable remote sensing diagnostic since they provide insight into kinetic processes relevant to small scale magnetosphereionosphere coupling.There are four major explanations of NEIALs.Two are based on two-stream instabilities (e.g.Foster et al., 1988;Rietveld et al., 1991;Wahlund et al., 1992), and two are based on wave-wave interaction (Forme, 1993;Bahcivan and Cosgrove, 2008).The latest experiments (Grydeland et al., 2004;Blixt et al., 2005;Strømme et al., 2005;Michell et al., 2009;Akbari et al., 2012;Isham et al., 2012) suggest that the Langmuir wave decay process (Forme, 1993(Forme, , 1999) ) is a plausible explanation to this phenomenon.In this explanation, the free energy is provided by an electron beam which can trigger either weak Langmuir turbulence (WLT) or strong Langmuir turbulence (SLT), depending on the energy of the beam.Recent results presented by Akbari et al. (2012) and Isham et al. (2012) show evidence of the first ISR observations of SLT.
The Langmuir decay process explanation comprises two steps.First, a down-going beam of electrons destabilizes Langmuir waves of the ionosphere through a plasma-beam instability.The beam-plasma instability can favor, depending upon the parameters of the beam and the plasma, either the Langmuir beam instability (LBI), where the Langmuir waves are enhanced ( b < 1), or the beam mode instability (BMI), where the beam mode is enhanced ( b > 1), where b is defined as (Gary, 1985;Diaz et al., 2010) b ≡ where v b , v b , n b and n e are the beam bulk velocity, the beam velocity spread, the beam density and the background electron density, respectively.Second, if the beam/plasma parameters produce the LBI, the enhanced Langmuir waves, with wavenumber k L0 and frequency ω L0 , can potentially decay into up-going Langmuir waves (k L1 ,ω L1 ) and downgoing ion acoustic waves (k s1 ,ω s1 ).To produce the decay, the electric field amplitude of the excited Langmuir wave, given by Melrose (1986) (assuming a WLT regime and v b < v b ), has to overcome the threshold given by Fejer (1979), where n i is the ion density, k B the Boltzmann constant, T i the ion temperature and C thr is a constant value that accounts for effects that can prevent the cascading, such as Landau damping.The value C thr is obtained from the expression presented by Fejer (1979) and Diaz et al. (2010), and is approximately 0.05 for the beam/plasma parameters of the simulations presented in this work.If the up going Langmuir waves gain enough energy from this process, they might also trigger a second cascading to downgoing Langmuir waves (k L2 ,ω L2 ) and up-going ion acoustic waves (k s2 ,ω s2 ).Since downward traveling Langmuir waves excite downgoing ionacoustic waves, and upward traveling Langmuir waves excite upgoing ion-acoustic waves, ion-acoustic waves traveling in both directions could be destabilized in regions of intense electron precipitation.Diaz et al. (2010) summarize the theoretical model of the Langmuir decay process, which emphasizes the possibility that the first Langmuir harmonic should be present if Langmuir decay is producing NEIALs.Due to the complexity of the beam-plasma instability, numerical simulation is the most important method to study the characteristic of ion acoustic enhancements due to a beam-plasma instability (Kasaba et al., 2001;Guio and Forme, 2006;Yi et al., 2010;Ziebell et al., 2011;Diaz et al., 2011).In recent research based on an electrostatic parallel particle-in-cell (EPPIC) code, Diaz et al. (2011) showed that the first Langmuir harmonic is present and might have enough power to be detected with current ISRs.The work of Diaz et al. (2011) represents the starting point for the present paper which presents more systematic treatment of the Langmuir decay process.The purpose of this exercise is to analyze simulation results from the EPPIC code for different beam-plasma scenarios, and to synthesize these results in a way that will facilitate future development of an estimation scheme for distorted spectrum.

Simulation description
The simulations performed in this work use an electrostatic parallel particle-in-cell (EPPIC) code (Oppenheim and Dimant, 2004).The code and the procedures needed to obtain the simulation parameters for different plasma states are described in more detail in previous works (Oppenheim and Dimant, 2004;Diaz et al., 2008Diaz et al., , 2011)).The 2-D simulation scenario of this work utilizes three species: electrons and protons (H + ) which compose the background plasma, and a weak electron beam.Representative values for these species and computational parameters are summarized in Table 1.In simulations presented below some of these parameters will be changed while leaving all others constant in an effort to characterize spectral distortions under various circumstances.
In order to simulate the behavior of the plasma when a beam is injected through it, particles leaving the volume must be discarded, while new particles with the original statistical characteristics must be continuously injected at x = 0 uniformly along y.The beam is set to travel on the positive x-direction (see axes in Fig. 3).Because the injection is in this direction, an open boundary condition is applied to the x-direction only, while a periodic boundary condition is used for the y-direction.Since the code is electrostatic the Poisson's equation is used to obtain, through finite difference approach, the potential as function of x and y for each instant of time, which is used to obtain the electric field in each point of the volume.Since the volume has open boundaries along x, the Poisson's solver uses a direct method of solution of the system of equation obtained by the finite difference method for this dimension, while a Fourier transform approach is used along y taking advantage of the periodicity of this dimension.With the electric field, the particles are advanced by using a leap frog method, obtaining a new charge density distribution for the next time.The simulation continues executing these basic steps (update particle positions, then fields) until a stop condition is reached (Birdsall and Langdon, 1985).
No constraints are imposed on the velocity distribution of the particles except at the beginning of the simulation and at the boundaries of the volume.The initial velocity distribution is Maxwellian for all species.Zero drift is used for the background plasma while the beam has a bulk velocity in the x-direction.The current code can perform analysis of the velocity distribution evolution, however, this process is very expensive if high velocity resolution is needed since it implies more outputting which is one of the most costly processes of the code.Therefore, velocity distribution analysis was not performed in this work.Table 1 also summarizes the values of the desired spectral specifications which together with the parameters of the plasma and the beam define the inputs with which the simulation is operated (Diaz et al., 2011).The total size of the box (N x x × N y y) is related to the resolution in the kspace, while spatial steps x and y are related to the maximum detectable wavenumber (k).When a beam injection is set (specifically v b ), the time step t is adjusted to be small enough to guaranty that no particle in the volume crosses more than one cell in each advancement of the time ( t < x/v b ), ensuring the stability of the simulation (Diaz et al., 2011).In order to save computational time, data is written to files only at the sample times which are given by the maximum frequency scattered by the plasma.The entry N s is the number of time steps per sample period.The total time of the simulation is related to the resolution in the ω-space -longer simulations (with more time steps, N t ) yield higher resolution in the IS spectrum ion line at low wavenumbers.N lags is related to this spectral resolution, since it is the number of points between the shoulders of the ion acoustic line at a low wavenumber.The actual "coherent" behaviour of the particles interacting with the wave modes is leveraged by grouping the particles into particle clusters, called macroparticles, which reduces the computational time and simulation noise.The number of a macroparticles per processor is denoted N p .The number of particles in each macroparticle can be obtained using the particle density, the volume of the box and the number of macroparticles per processor.
The main physical parameters were selected based on averaged measurements observed at high latitudes (Diaz et al., 2011).The beam parameters were, in general, selected to have current densities (|J|) between 100 to 300 µ m −2 , values consistent with observed current densities over the discrete aurora (Stasiewicz et al., 1997).Since spacecraft potentials are usually a few volts negative, the electrons with energies from ∼0.08 to ∼6 eV are repelled.Because of this, the thermal electron distribution at those energies has never been measured effectively in the nighttime mid-altitude auroral zone.The simulations presented in this work theoretically explore possible signatures of NEIALs produced by these low-energy types of electron beam distributions.The distribution function of the electron beam was selected to satisfy the current density and the enhanced radar wavenumbers observed by experiments.
Although the physical parameters (e.g. the mass of the species) are changed to facilitate simulations, these input parameters were selected to obtain growth rate curves similar to those computed with realistic values.The background electron density is selected to be close to the peak of the F-region since this value set (through the plasma frequency) has the maximum enhanced wavenumber, given a certain beam velocity, that can trigger a decay (k = ω pe /v b ).The selection of this large n e imposes a separation between the ion acoustic and plasma lines requiring a restrictively short sampling rate to avoid aliasing between all modes present in the IS spectrum (including the Langmuir harmonics).In order to bring the ion acoustic and plasma lines closer in frequency (smaller ω pe ), a smaller ion to electron mass ratio (actually m i /m e ∼ 2000) is used.An analysis performed by Diaz et al. (2011) shows that the optimal mass ratio to avoid artificial behavior of the decay process and to preserve spectral accuracy of the ion acoustic line is 500.Usually, the theory uses two limits for the temperature ratio, T e /T i = 1 and T e /T i 1.However, with the simulation it is possible to use more realistic ratios.For the present work, the ratio T e /T i = 2 is chosen, since lower temperature ratio produces more "flat" thermal ion acoustic spectra, facilitating identification of some enhancements (Diaz et al., 2008).Even though the geomagnetic field has no effect on the longitudinal unstable waves caused by the electron beam, a constant value of 50 000 nT is used to be in accordance with the reality.
The simulation performed with the parameters presented in Table 1 was able to trigger a Langmuir decay event (Diaz et al., 2011), a result consistent with theoretical calculations.In this paper Langmuir decay events are presented in three different manners: (1) as perturbations of the relative ion density variation, (2) as an overcoming of the decay threshold by the average electric field, and (3) as an enhancement of the ion acoustic line (IS spectrum).The relative density variation is defined as where n j (x, t) is the density at every point of the mesh at each sample time of species j , and n 0j the average density of species j .The total electric field as a function of time is calculated as where N x and N y are the total points of the mesh in xand y-directions, respectively; x y the step in x-and ydirections, respectively; l and m are indexes for points in the x-and y-directions, respectively, and t is the time index.Finally, the IS spectrum is defined as where N j (k, ω) is the relative density variation of the species j , in k and ω space, and the angle brackets represent an average over multiple independent samples of the spectrum.
The spectral convergence to the mean value of the spectrum for a plasma in thermal equilibrium can be facilitated by taking advantage of the angular independence of the spectrum over a two-dimensional plane (k x ,k y ) (Diaz et al., 2008).In the present work, the angular independence is lost due to the injected electron beam and the magnetic field.However, Diaz et al. (2011) show that the asymmetries due to the beam and the magnetic field are small for angles close to 0 • (or along x).Hereafter, all spectra obtained throughout this work use the integration of 44 spectra obtained by two simulation runs.Each run is composed of 22 angular integrated spectra obtained assuming a 3 • angular resolution between angles 15 • to −15 • and 165 • to 195 • over the (k x , k y ) plane (see Diaz et al., 2008).The spectrum figures are presented in dB (10 log 10 [S j (k, ω)]) to facilitate detection of the enhanced part of the spectrum together with the part of the spectrum where no enhancement is present.As a reference the ∼ −90 dB represents the thermal equilibrium power spectrum level (Fig. 5).At this point a clarification has to be done for the use of the term "IS spectrum" in this work since the simulated scattered spectrum is not incoherent for the cases simulated here.The term "IS spectrum" in this work refers to the power spectrum obtained from the relative electron density fluctuations and how an ideal multi-frequency wide-band IS radar spectrum would appear under the simulated scenarios.
Although most of the simulation parameters remain constant throughout this work (as outlined above), some simulation parameters, such as the sampling rate, are modified between simulations to take computational advantage offered by the particular physical parameters under investigation.Changes from the values presented in Table 1 will be highlighted where they occur.

Parameter analysis of the Langmuir decay
This section presents a systematic study of the sensitivity of the Langmuir decay process to the physical parameters of the background plasma and the beam.The conditions given by b < 1 and E 2 L > E 2 thr can be rearranged as and by using Eqs.( 1), ( 2) and (3).From Eqs. ( 4) and ( 5), it is possible to see that the critical parameters are the ratios n b /n e and v b / v b as shown in Fig. 1.All the simulations presented in this section satisfy the condition given by Eq. ( 4), analyzing the validity of the condition given by Eq. ( 5).In order to keep the factor α 0 in Eq. ( 5) constant (at 5.9, with the values presented in Table 1), v b will be kept constant through the simulations presented in this section.The ratio v b / v b is varied by altering v b .Therefore, the beam density (n b ), velocity (v b ) and background electron (ion) density (n e ) are chosen to explore the response of the decay process.

Density of the electron beam (n b )
Four different values of n b are simulated in this subsection.Besides the different densities, only the parameters N s and N t were changed from the original selection.Those parameters were changed to 16 and 34 000 respectively.The four cases are selected to produce enhancement of the Langmuir waves  over wavenumbers estimated from the growth rate equation (Diaz et al., 2010, Eq. 7) and shown in Fig. 2. Three of the cases satisfy the theoretical decay condition presented in Eq. ( 5) as shown in Table 2. Thus, those densities are expected to produce a decay from excited Langmuir waves into ion acoustic waves.The last density is selected to have E 2 L − E 2 thr less than zero (∼ −5 V 2 m −2 ), therefore no decay is expected for that density.
All simulated beams enhanced Langmuir waves which can be observable in the electron density variations.The growth of all four E 2 L (t) curves from thermal level also confirm this fact in Fig. 4. Figure 3 shows that ions are perturbed, as expected, for the three largest densities and that larger beam density produces larger and faster perturbations of the ions (see video 1 in supplementary material).The weakest electron beam, unexpectedly also generates enhanced ion acoustic waves.The simulated E 2 L (t) curves for all beam densities are able to overcome the decay threshold producing perturbations of the ions.It is important to notice that the inaccuracy in the prediction is in E 2 L (Eq. 2) and not in the threshold value (Eq. 3) since the curves for E 2 L (t) show that the lowest beam density crosses the threshold value at ∼0.12 ms, which agrees with the time when perturbations start to appear in the ion density variation (Fig. 3).Even though the spectral resolution is not the best (the resolution can be improved with longer simulations), it is possible to see the same basic characteristics in the IS spectrum as in the density variation plots.Larger beam densities produce stronger wave mode enhancements on the spectrum.First Langmuir harmonic (at ∼ +2ω pe and k ∼ +2ω pe /v b , see Yoon et al. (2003) for more detailed estimations for ω and k of the first harmonic) is also clearly produced and its enhancement is higher for larger densities (see Fig. 5).An interesting characteristic shown by the simulations is that the first Langmuir harmonic is present for negative frequencies when the Langmuir decay occurs.The possible detection of the first Langmuir harmonic might be useful to constrain the values of the beam and the electron density.When the decay is established, the values of those parameters are restricted to a much smaller region over the plane (v b / v b ,n b /n e ) as it is shown in Fig. 1.The Langmuir harmonic at negative frequency also serves as evidence that the Langmuir decay process is happening at lower wavenumbers (or radar frequencies).Figure 6 shows a zoom on the ion acoustic line.Larger beam densities clearly produce larger amplifications of the ion acoustic waves and over a wider range of wavenumbers.This may be due to a secondary cascade over a larger range of wavenumbers.
A more detailed analysis on the spectral characteristics will be conducted for one specific case.Similar analysis applies to the rest of the cases of Sect.3.1.rizes the predicted wavenumbers where the enhancements should be produced for the plasma line as well as for the ion acoustic line of a beam with parameters presented in Table 1. Figure 5b shows the incoherent spectrum for this beam.The prediction of the wavenumber of the excited positive Langmuir wave (k L0 ) is in approximate agreement with the one obtained by the simulation (Fig. 5b).The power of the enhancement of the positive Langmuir mode is ∼50 dB above the thermal case (Diaz et al., 2011).The calculated range for the excited wavenumbers of the first cascaded Langmuir wave (k L1 ∈ [10.4,31.4]m −1 ) also approximately agrees with the range obtained by the simulation.Given the low resolution of the ion acoustic line, it is difficult to say if the wavenumber of the first cascaded ion acoustic wave (k s1 ), which should enhance the shoulder at positive frequencies, exhibits the predicted behavior.The excited wavenumbers of the ion acoustic line covers most of the wavenumbers from 0 to 60 m −1 (Fig. 6b).However, the highest wavenumber enhanced of the positive Langmuir mode (k max L0 ∼38 m −1 in Fig. 2) would produce an enhancement at k max s1 ∼ 68.4 m −1 (Table 3), which seems to agree with the highest wavenumber excited in the ion acoustic mode.The first cascaded Langmuir wave (L1) seems to trigger a secondary decay.The secondary decay produces an excitement of the positive Langmuir mode (k L2 ) close to the theoretically calculated values of [2.8,22.8]m −1 .This decay should also produce an enhancement on the negative shoulder of ion acoustic mode (k s2 ) at wavenumbers close to k s2 ∈ [13.2 ,53.2] m −1 , which can be observed.Figure 6b shows a predominant enhancement over both shoulders of the ion acoustic line between the wavenumbers 0 to 60 m−1, the enhancement of low wavenumber might be an indication of a third decay process of the second enhanced Langmuir wave k L2 .

Velocity of the beam (v b )
Three beam velocities are simulated leaving the remaining parameters unchanged.Only n b , N s , and N t are changed from the original configuration to 7.5 × 10 8 m −3 , 16 and 34 000 respectively.Table 4 presents the electron beam velocities under investigation.
The three speed values are selected to produce enhancement of the Langmuir waves over wavenumbers estimated  from the growth rate equation (Diaz et al., 2010, Eq. 7), similarly to those obtained in Sect.3.1 (Fig. 2).The two highest selected velocities should produce Langmuir decay, unlike the lowest selected electron beam velocity.Figures 7a and 3b present the ion density variation for two of the three simulated velocities, showing that for the two largest velocities the decay is triggered.The lowest selected electron beam ve-locity, as expected, does not produce a decay.However, the relation E 2 L once again appears to be inexact in quantifying the amount of energy transferred from the electron beam to the LWs.On the other hand, the theoretical threshold for the magnitude squared electric field (E 2 thr ) once more shows its accuracy.Figure 7b shows that for the two largest beam velocities, E 2 L (t) cross the threshold at t ∼ 40 µs and t ∼ 90 µs, respectively, which agrees with the appearance times of the structures in the ion density variation (Figs.7a and 3b).
Figure 8a shows IS spectrum, for the beam velocity of v b = 6.75 × 10 5 m s −1 .The wave number (|k|) where the enhancement of the positive Langmuir mode starts to be evident (|k| ≈ 18 m −1 ) approximately agrees with the theoretical calucalations |k| ∼ ω pe /v b ≈ 22 m −1 ).However, because the maximum growth rate is reached at larger wavenumbers, the maximum enhancement is reached in the simulations at |k| ≈ 29 m −1 .The expected result for a slower electron beam velocity is to shift the enhancement toward larger wavenumbers.Figure 8b shows the IS spectrum for the case when v b = 4.9 × 10 5 m s −1 , where it is possible to see that the enhancement starts around |k| ≈ 24 m −1 (∼ ω pe /v b ≈ 30 m −1 , from the model).In addition to the shift, the intensity of the power is reduced (∼2 to 3 dB) along the whole range of wavenumber where the Langmuir waves are excited.
Characteristics similar to those exhibited for the positive Langmuir mode can be seen in the ion acoustic mode in Fig. 9.The enhancement covers a large range of wavenumbers for both beam velocities with a lower intensity for the lower beam velocity.Some preference can be noticed for the positive ion acoustic mode (positive shoulder), although the resolution is not good enough for more detailed analysis.
Figure 8a together with the theoretical expressions of the decay made by Diaz et al. (2010) (similar to those presented in Table 3) can be used to make a more quantitative analysis of the spectrum of one of the cases simulated in Sect.3.2.Figure 8a shows the incoherent scatter spectrum for the beam velocity v b = 6.75 × 10 5 m s −1 .As expected, the range of the excited wavenumbers for the positive Langmuir mode is narrower and shifted toward larger wavenumbers compared with the case where v b = 7.5 × 10 5 m s −1 (Fig. 5b).The prediction for the lowest wavenumber for the excited positive

Langmuir wave (k min
L0 ∼ 22 m −1 ) approximately agrees with the one obtained from the simulation (Fig. 8a).The power of the enhancement is 40 dB higher than that of the thermal case (Diaz et al., 2011).The wavenumber range of the first cascaded Langmuir wave (k L1 ∈ [14.4,33.4] m −1 ) also approximately agrees with those obtained from the simulation.Those wavenumbers are associated with the negative Langmuir mode (−ω pe ).The enhancement of the ion acoustic mode (k s1 ) appears at wavenumbers lower than those calculated theoretically (k s1 ∈ [36.4, 74.4] m −1 ), developing a maximum at wavenumbers close to 30 m −1 (Fig. 9a).This reveals a likely problem with the theoretical expressions.The first cascaded Langmuir wave (k L1 ) does not seem to generate a secondary decay, or if it does, it is too weak to be detected.If the secondary decay had occurred, an excitation of the positive Langmuir mode would be visible at wavenumbers close to k L2 (∼6.5 m −1 ), which is not present.Figure 9a shows a close up of the ion acoustic line and presents a tendency to enhance the positive ion acoustic mode (positive shoulder), as expected.The enhancement of the ion acoustic line increased the power level over 20 to 30 dB above the thermal level.

Background electron density (n e )
The velocities simulated previously in Sects.3.1 and 3.2 put most of the distortions in the wavenumber range between 8 m −1 and 60 m −1 , a wavenumber range where most of the IS radars operates.However, the strongest enhancements of the incoherent scatter spectra (the most likely to be detected) appeared far below |k| = 54 m −1 which is the wavenumber of Sondrestrom radar.For the beam presented in Table 1 (first case) and n e = 2.5 × 10 11 m −3 , the simulation showed a decay of the Langmuir waves.This case would place the approximated wavenumber of the enhanced Langmuir waves at k ∼ ω pe /v b = 19 m −1 .If we want the wavenumber to be higher, in order to detect the enhancements with ISRs of higher frequency, the beam velocity should be decreased.For instance, if the beam velocity is decreased to v b = 3.9×10 5 m s −1 , the approximated wavenumber of the enhanced Langmuir waves would be k ∼ ω pe /v b = 38 m −1 .This low velocity would not be able to trigger the decay with n b = 7.5×10 8 m −3 and n e = 2.5×10 11 m −3 (E 2 L − E 2 thr = −14 < 0).However, if n e is increased to 1×10 12 m −3 and used with the original beam (Table 1), the LBI is enhanced ( b −1 = −0.69), the decay is triggered (E 2 L −E 2 thr = 181 V 2 m −2 > 0) and the approximate wavenumber of the enhanced Langmuir waves would be placed at k ∼ ω pe /v b = 38 m −1 as desired.Therefore, this section will explore if increasing the background density (n e = 10 12 m −3 ), which determines the plasma frequency, will allow higher velocities (Table 1) to destabilize ion acoustic waves at larger wavenumbers.The physical parameters, changed from those in Table 1, are n b = 2.5 × 10 9 m −3 , n e = 10 12 m −3 , v b = 7.8 × 10 5 m s −1 , |J| = 312 µA m −2 , E 2 L − E 2 thr = 181 V 2 m −2 and b −1 = −0.7.The new background density changes the Debye length; therefore, the simulation parameters must be changed.Those parameters are summarized in Table 5.
Figure 10a shows the ion density variations for a simulation with a larger background density.This figure shows structures, which demonstrate that a beam triggers a bump on tail instability of Langmuir waves and a decay of those into ion acoustic waves.It is possible to see that structures start to appear at t ≈ 0.014 ms, which agrees with the time when the E 2 L (t) curve crosses the theoretical threshold E 2 thr = 154 V 2 m −2 (Fig. 10b).Even with the spectral resolution of this simulation ( t was reduced by half, making it more difficult to reach larger times), the spectral analysis (Fig. 11) shows, as expected, that the enhancement over the ion acoustic and positive Langmuir waves were shifted to higher wavenumbers with respect to the previous simulations due to the larger background density used.

Density variation for a low velocity beam
This section will explore the precondition b < 1 (Eq.1).In Sect.3.2 was shown that lower beam velocities places enhancement of the Langmuir waves at higher wavenumbers, which eventually might produce enhancement of the ion acoustic waves also at higher wavenumbers (or radar frequencies).However, Sect.3.2 also showed that for a density of the electron beam equal to 0.3 % of the background density, a beam velocity too low ( 5 × 10 5 m s −1 ) cannot trigger a decay.In this section v b and v b are changed to 3.9×10 5 m s −1 and 7.8×10 4 m s −1 , respectively, leaving the ratio v b / v b = 5 equal to that used in Sect.3.1.If the density of the beam is selected to be n b = 7.5 × 10 8 m −3 , this beam would be unable to trigger the decay (E 2 L − E 2 thr = −14 < 0).Increasing n b to 1.16 × 10 9 m −3 , the decay should be produced (E 2 L − E 2 thr = 0.12 > 0).Increasing n b increases the value of E 2 L −E 2 thr ensuring the decay (however it has a limit) when the precondition b − 1 < 0 is violated.This happens when n b overcomes the value of n * b ≈ 2 × 10 9 m −3 (0.8 % of n e ).For values over this beam density n * b , no decay should be produced; therefore, to verify this prediction made by theory, two values of the beam density much larger than 2×10 9 m −3 are selected.Besides the two used beam densities (Table 6), v b , v b , N s and N t are the only changed parameters, changed to 3.9 × 10 5 m s −1 , 7.8 × 10 4 m s −1 , 16 and 22 000, respectively, for these simulations with respect to the original values.Table 6 also shows the theoretical E 2 L − E 2 thr calculated with the new beam parameters, which reach extremely high values compared to the threshold needed to trigger the decay (> 0), but with values of b −1 larger than zero.As expected, none of the beam densities triggers a decay.For both n b the electron density variations are perturbed, however no perturbations are perceived on the ion variation density for both cases.Figure 12a shows that the BMI is being enhanced.Figure 12b shows no enhancement on the ion acoustic line.The curve E 2 L (t) converges to a value ∼ 3 V 2 m −2 for both densities, below the threshold to trigger the decay and far below the values theoretically obtained, confirming that Eq. ( 2) is only valid if b < 1 (see video 2 in supplementary material).
The current densities of the two latest simulations are as ∼81 µA m −2 (Sect.3.2) and ∼31 µA m −2 (lowest n b in Sect.3.1) achieved the decay.Consequently, it is possible to conclude that the current density is not a good indicator to evaluate the decay capacity of a beam.

Discussion
The simulations performed in previous sections show that, for some beams, it is possible to simultaneously observe enhancement of the plasma lines, the ion acoustic line, and the Langmuir harmonic line (e.g.Fig. 5a).Although these results agree with simultaneous enhancements of plasma line and ion acoustic line reported previously (Strømme et al., 2005;Akbari et al., 2012;Isham et al., 2012), those events are not common.On the other hand, no Langmuir harmonic has been reported to be produced naturally and simultaneously with enhancements of either plasma or ion acoustic lines.Simulations also show that NEIALs could be observed at different wavenumbers (or radar frequencies).However, there have never been simultaneous NEIALs observations reported at EISCAT UHF and VHF radars, which operate at two different frequencies but observe a common scattering volume.Cabrit et al. (1996) suggests that this is due to the difference between the UHF and VHF radar beam widths.Guio and Forme (2006) argue that the enhancement in the IS spectra can vary by orders of magnitude between two different frequencies.This difference of the enhancements can be increased by the integration process of the radar (seconds) rendering the enhancement at one of the frequencies undetectable.The simulations performed at this work do not include the possible diminishing effect that the integration process can have over the spectrum, which might be important given the ephemeral nature of the phenomenon.Therefore, the simulated enhancements may be larger than those observable by radar.
Although the integration procedure (temporal and spatial) might explain some of the lack of detections at the current ISRs; it is also possible that the simulations of this work are overestimating the electron beam, thus giving the impression that the simultaneous detection is likely.The beam was selected to be in agreement with some of the data collected for beams at the top of the ionosphere.However, low energy beams at lower altitudes have not yet been measured accurately, which necessarily led us to an arbitrary selection of the beam.Future work will attempt a closer connection of the model to beam measurements at the top of the ionosphere.Furthermore, measurements in the topside ionosphere might be extrapolated to altitudes where NEIALs are being detected by modeling the beam degradation due to collisions.
Section 4 tested the validity of the precondition b < 1.If this condition is not satisfied, the decay cannot be triggered and the decay condition E 2 L > E 2 thr cannot be used to predict the Langmuir decay process.On the other hand, in Sect.3.3 it was shown that larger background densities (n e ) can help to place the enhancement at higher wavenumbers (k ∼ ω pe /v b ) and to compensate for the increament of the v b / v b ratio of slow-cold beams by the n b /n e ratio ensuring that b stays less than one.At high altitudes the beam is expected to be fast and warm, placing the possible enhancement at low wavenumbers.While the beam is penetrating the atmosphere, it might become slower and colder.Although this is speculative, the lack of data allows the assumption.
On the other hand, the background density increases until a maximum at 300 km of altitude.The combination of these two effects could produce an increase of the enhanced wavenumbers (k ∼ ω pe /v b ) with decreasing altitude down to 300 km.For altitudes lower than 300 km the background density decreases with altitude and the beam velocity will keep decreasing.Whichever parameter decreases faster will determine the maximum wavenumber.For instance, if √ n e (ω pe ∝ √ n e ) decreases faster than v b , then the wavenumber will be lower than that obtained at 300 km (maximum n e ).If √ n e decreases at the same rate as v b , then the wavenumber will be similar to that reached at ∼300 km.Finally, if √ n e decreases slower than v b , the wavenumber will be higher than that reached at 300 km.Since the beam is becoming slower and colder, it is also possible that the parameter b overcomes one, passing the energy to the beam mode instead of the Langmuir mode, thus inhibiting the decay process.Even so, some detailed investigation of the evolution of the beam while it penetrates the atmosphere is necessary to see if Eqs. ( 5) and ( 4) are satisfied.It seems likely that the maximum enhanced wavenumber is obtained at the maximum background density n e .Statistical studies of radar data (see Buchert et al., 1999) had shown that it is more likely to observe coherent echoes in incoherent scatter radars of lower frequency or wavenumber than for higher frequency radars.The simulations presented in this work agree with this fact.Since Sondrestrom operates at a high wavenumber, it has never observed a coherent echo or NEIAL.The simulations and the previous discussion suggest that such a high wavenumber cannot be reached because the maximum values of ionospheric n e are too low to favor instabilities at such high wavenumbers.For example, the maximum background density of n e = 2.5 × 10 11 m −3 cannot favor the decay at high altitudes and high wavenumbers because low beam velocities would be needed.However, the simulations show that low velocities cannot trigger the decay of Langmuir waves because they are too weak (or warm not satisfing Eq. 5).Slower and colder beams might be reached at lower altitudes where low n e is found making the beam too strong (not satisfying Eq. 4).It is difficult to estimate the precise background density needed to actually place enhancements of the ion acoustic mode close to k ∼ 54 m −1 .However, it is possible to argue that the actual value of the background density needed to trigger detectable instabilities at high wavenumbers may be larger than the value obtained in this work, 10 12 m −3 .This is because, although the plasma frequency would be higher for realistic m e , allowing larger beam velocities for large wavenumbers, the actual mass of the ions (O + instead H + ) would increase the threshold needed to trigger the decay (Eq.3).Hence, this situation would require even larger beam velocities to produce a destabilization of the ion acoustic waves.Rietveld et al. (1996) showed that an increment of the solar activity increases the number of events of ion acoustic enhancement detected by a given incoherent scatter radar.
Apparently, there is a contradiction of this result, with the beam model explanation of NEIALs, since more solar activity might produce more energetic beams.These energetic beams might generate enhancement of ion acoustic waves at lower wavenumbers (or radar frequencies) making the detection of this enhancement unlikely with the same incoherent scatter radar.However, an increment in solar activity would also produce an increment of the background density (n e ) so that the enhancement of those more energetic beams can be detected at the same (or even larger) wavenumbers as was shown in Sect.3.3.
It is also important to notice the behavior exhibited for the ion acoustic enhancement at high wavenumbers.For most of the simulations performed in this work, the enhancement stops quite suddenly at large wavenumbers, before decreasing to an unnoticeable level.The enhancements exhibit a Gaussian noise-like shape.This kind of enhancement might lead to misestimation of the plasma parameters, especially the electron density, because this kind of enhancement is not strong enough to be noticed and it does not produce any significant change in the shape of the ISR spectrum.Guio and Forme (2006) got results where the upgoing plasma line were more enhanced than the downgoing plasma line for a downgoing beam for both weak and strong turbulent regimes.However, it is not clear that the simulations performed in this work agree with this result.Figure 11a shows an enhancement on the negative plasma line stronger than the one produced on the positive plasma line for k < 36 m −1 .Nevertheless, this feature can be attributed to the fact that after the decay the Langmuir waves traveling in opposition to the beam are enhanced at lower wavenumbers than the Langmuir waves traveling in the same direction that the beam (Diaz et al., 2010, Eq. 16).This possible disagreement has to be analyzed in more detail in future works.
The model summarized by Diaz et al. (2010) assumes a weak Langmuir turbulent regime, although simulations can be developing either weak or strong Langmuir regimes.The main goal was to explore the accuracy of the weak Langmuir model for different plasma and beam parameters.Since the beam has to be appropriate to trigger the decay (Beam Mode Instability limitation) it was hypothesized that the weak Langmuir turbulent model was accurate enough to represent the phenomenon.However, the inaccuracies in the quantification of the transferred energy from the beam to the Langmuir waves for both regimes (weak and strong) imposes the need to review the validity of Eq. (2).Problems with the accuracy of the enhanced wavenumbers (after decay) are also revealed.However, these issues might be due to the lack of resolution in the current simulations.
All of the expected characteristics of enhanced ISR spectra are present in the simulations.This implies that the simulator might be used in a systematic way to formulate an empiricalsimulated model of the ISR spectrum for different plasma and beam parameters which could be used to estimate the parameters with actual distorted ISR data.However, if the code wants to be used to obtain a simulation-empirical model of the distorted spectrum, the frequency (ω) resolution has to be improved on each line of the radar making longer simulations.Although a simulation of millisecond can actually take many hours of computation, the computing time can be significantly reduced by reducing the need of outputting data which can be done discarding the initial time and focusing in the simulated time where the decay is taking place.

Conclusion
The EPPIC code was successfully used to study natural enhancement of ion acoustic waves via injection of an electron beam into a background plasma.This work shows that the signatures of the spectrum are sensitive to the parameters of the beam and plasma, reinforcing the idea that the parameter estimation is possible for beam-distorted ISR spectra.The EPPIC code might help to build a semi-empirical spectrumplasma parameter model which might allow a parameter estimation process.Simulations also show that a multi-channel (0, ±ω pe , ±2ω pe ) and multi-wavenumber (multiple radar frequencies) IS radar system would help to constrain the estimation procedure and discriminate the driver of the ion acoustic enhancement.Other explanations of NEIALs would not produce some of the signature of the beam-plasma process.For instance, the first Langmuir harmonic would not be present if a beam is not the cause of the observed NEIAL.Therefore, design of experiments that look for Langmuir harmonics arise as a future task to do.
Simulations confirm that the main parameter driving the Langmuir enhancement-decay process is the ratio v b / v b .While the velocity of the beam, v b , is relevant to set up the wavenumber enhanced of the plasma and ion acoustic lines, the v b / v b ratio is relevant to select which wave mode is enhanced by the beam (either the Langmuir or the beam mode).
The simulations suggest that this strong dependence of the decay on the ratio of beam velocity to beam velocity spread and background density limits the maximum value of wavenumber (near the F-region peak).This might explain why radars of high frequency, such as Sondrestrom, are unable to detect NEIALs.
Simulations also show that the first Langmuir harmonic is evident for negative frequencies when the decay is triggered, which might be used by ISR of higher wavenumbers as indication of a Langmuir decay process at lower wavenumbers in case those radars could detect the first up-going Langmuir harmonic (−2ω pe ).Although NEIAL model, in general, works well, these simulations found some issues with the model.The main problems are as follows: Decay prediction of waves close to the threshold due to the inaccuracies in the value of the energy transferred from the beam to the Langmuir waves; and the wavenumber placement of the ion acoustic waves that will be excited, especially those coming from a secondary decay.Both issues might be related to the validity of the equations used for the different regimes.For instance, Eq. ( 2) uses quasilinear approach (Vedenov, 1963;Boyd and Sanderson, 2003) assuming WLT regime to calculate the energy transferred from the beam to the Langmuir waves.However, this equation shows some inaccuracies not only for the SLT regime (as expected) but also for the WLT regime tested with the simulation (e.g.weakest beam in Sect.3.1).

Fig. 1 .
Fig. 1.Area over the v b / v b -n b /n e plane for which both equations (Eqs.4 and 5) ensure LBI and Langmuir decay.The area of Langmuir decay is increased when α 0 is increased ( v b is increased).

Fig. 4 .
Fig. 4. Magnitude squared electric field versus time for simulations with parameters presented in Sect.3.1 for the four beam densities.

Fig. 7 .
Fig. 7. (a) Snapshots of the ion density variation for simulations with parameters presented in Sect.3.2 for v b = 6.75 × 10 5 m s −1 .The ion density variation associated with v b = 7.8 × 10 5 m s −1 is shown in Fig. 3b.The beam with velocity v b = 4.9 × 10 5 m s −1 does not triger the decay; therefore, the ion density variation has a noise-like appearence being omitted for this reason.(b) Simulated magnitude squared electric field (E 2 t ) versus time for simulations with parameters presented in Sect.3.2.

Fig. 10 .
Fig. 10.(a) Snapshots of the ion density variation for simulations with parameters presented in Sect.3.3 (n e = 10 12 m −3 ).(b) Magnitude squared electric field E 2 t versus time for simulation with parameters presented in Sect.3.3.

Fig. 11 .
Fig. 11.(a) Incoherent scatter spectrum and (b) ion acoustic line of the incoherent scatter spectrum for simulations with parameters presented in Sect.3.3.

Fig. 12 .
Fig. 12.(a) Incoherent scatter spectrum and (b) close up to the ion acoustic line of the IS spectrum for simulations with n b = 1.25 × 10 10 (Sect.4).

Table 1 .
Summary of the parameters used to simulate the decay of Langmuir waves.

Table 5 .
Summary of parameters changed from those presented in Table1.