Ion dynamics in electron beam – plasma interaction : particle-in-cell simulations

Electron beam–plasma interaction including ions is studied by particle-in-cell (PIC) simulations using a onedimensional, electrostatic code. Evidence for Langmuir wave decay is given for sufficiently energetic beams, as in previous Vlasov–Maxwell simulations. The mechanism for the generation of localized finite-amplitude ion density fluctuations is analyzed. Amplitude modulation due to interference between the beam-generated Langmuir waves causes random wave localization including strong transient spikes in field intensity which create bursty ion density structures via ponderomotive forces. More dense beams may quench the decay instability and generate low-frequency variations dominated by the wave number of the fastest growing Langmuir mode.


Introduction
Electron beam-generated Langmuir wave forms, observed for more than 30 years under various space plasma conditions, often exhibit a spectral double peak near to the ambient plasma frequency ω pe (Gurnett et al., 1993;Hospodarsky et al., 1994;Hospodarsky and Gurnett, 1995;Kellog et al. 1996;Bale et al., 1996;Kojima et al., 1997;Pottelette et al., 1999;Sigsbee et al., 2010;Graham and Cairns, 2013a).It is commonly assumed to be produced by the interference of the dominant Langmuir wave of the instability and a second Langmuir wave of comparable amplitude but slightly different frequency.Various efforts have been made to understand the origin of the second Langmuir wave.A well-known nonlinear wave interaction process which is able to produce such a scenario is the parametric decay of the most unstable Langmuir mode into a counterpropagating second Langmuir wave and a forward propagating ion acoustic wave.This process has been proposed long ago as basic mechanism for the interpretation of type III solar radio bursts (Cairns and Melrose, 1985).In this resonant three-wave interaction the frequency of the ion acoustic mode is very small against the plasma frequency which makes it very difficult to resolve separated peaks in the plasma frame.However, the resonance conditions require a significant difference in the wave numbers of the two high-frequency waves.That is, in case of a relative motion between the observer (spacecraft) and the main plasma (e.g., solar wind), the Doppler shift may enlarge the frequency difference such that high-resolution detectors are able to identify the two peaks.It turned out that, in several cases, using estimated beam velocities, the predictions of the three-wave parametric decay via ion sound waves are consistent with the observed beat frequencies (Hospodarsky and Gurnett, 1995;Kellog et al., 1996;Graham and Cairns, 2013a).
It should be mentioned, however, that a variety of other types of electron beam-driven Langmuir amplitude modulations have been observed which do not fit into the above parametric decay scheme.Among them are spectra with power at ω pe only (e.g., Sigsbee et al., 2010), a double peak in absence of ion acoustic activity (e.g., Sigsbee et al., 2010), multiple spectral peaks (e.g., Graham and Cairns, 2013a), a single peak with side-band structure (e.g., Gurnett at al., 1981;Thejappa et al., 2013), Langmuir eigenmodes localized in density wells (Ergun et al., 2008;Graham and Cairns, 2013b) and waveforms resulting from the interaction of Langmuir waves with electrostatic whistler/lower hybrid waves instead of ion sound (Stasiewicz et al., 1996).
The PIC simulations presented in this short contribution first confirm Langmuir wave decay as a transient process in the evolution of the beam-plasma system (Sect.3).Then the mechanism for the generation of localized finite-amplitude ion density variations is examined (Sect.4).After presenting simulation results in time representation, including the course of the energy partition (Sect.4), it is shown that a more dense beam may create a modified wave scenario, with low-frequency fluctuations dominated by k L , the wave number of the fastest growing Langmuir mode (instead of 2k L for decay) and a spectral double peak of the high-frequency field seen even in the plasma frame (Sect.6).

Simulation model
A one-dimensional electrostatic standard particle code is used with periodic boundary conditions.The initial condition is set up by Maxwellian distributions of background electrons and ions while the electron beam is represented by a drifting Maxwellian (V b drift velocity).In order to keep a tolerable degree of energy conservation, the electric field for the acceleration of electrons and ions is calculated by an implicit momentum method according to the philosophy of Langdon (1985), based on the current equation with n e , v e , n i , v i the particle density and bulk velocity of electrons and protons, respectively.The moments n e and v e are built with both beam and background electrons.In order to check the accuracy of the field calculation via Eq.( 1), the field is additionally calculated from the Poisson equation A comparison of both results did not give rise to corrections.Our frame of reference is that of the background particles at rest.The presence of an electron beam is then equivalent to a current.To avoid this initial current, one may compensate it by a small drift of the background electrons in direction opposite to the beam, as done in a number of simulation studies.The evolution of the instability then starts with zero electric field.Otherwise, in 1-D geometry, the initial current is the source of "beam-aligned" plasma oscillations with an arbitrary phase φ (Baumgärtel, 2013).This k = 0 mode is present from the very beginning prior to the onset of the instability, with an amplitude controlled by beam density and beam velocity.For thin beams (n b /n 0 0.01) where the frequency of the fastest growing mode is very close to ω pe , no additional signal can be expected from this oscillation in the power spectrum.This may change for denser beams (n b /n 0 0.05) which downshift the frequency of the dominant mode, resulting in a peak at ω pe and another below ω pe .We start the simulations generally with background electrons at rest; that is, we accept an initial current and the associated plasma oscillations.
In the simulations, time and frequencies are scaled according to the electron plasma frequency ω pe , distance is in units of the Debye length λ D , densities and velocities are normalized by the equilibrium density n 0 and the thermal velocity v the of the background electrons, respectively.As a consequence, the longitudinal electric field is then scaled by E 0 = m e /c v the ω pe .For nominal solar wind parameters n 0 = 4 cm −3 , T e = 10 eV the field unit corresponds to E 0 = 1.8 V m −1 .For this reference plasma, the mostly used beam velocity V b = 16v the is equivalent to V b ≈ 0.07c with c the light velocity.The plasma parameters listed in Table 1 are thought to be representative so that the simulation results may be used to some extent for comparison with observations.Table 1 contains in addition technical simulation parameters.

Langmuir wave decay
As known from former simulations, for a beam speed of the order of V b /v the ∼ 10, the longitudinal field starts growing after a period of some tens of ω −1 pe , dependent on the level of the initial field noise.In this initial stage the waves grow exponentially at the linear growth rate.As a result of a simulation run with V b /v the = 16, Fig. 1 depicts the spatial variation of field and proton density and the corresponding power spectra from an early part of the time development (t = 200 ω −1 pe ).The almost coherent field structure shown in panel (a) consists of a number of wave packets and is dominated by the fastest growing Langmuir mode, as indicated by the pronounced peak at k L λ D ≈ 0.075 in the power spectrum (panel b).The amplitude modulation is the result of linear interference of unstable waves with wave numbers in a small range round k L .The structure moves in beam direction with almost the group velocity of the dominating mode v gr ≈ 0.22 v the , accompanied by a continuous change of the envelope.Due to the random phases of the waves that have grown from noise, occasionally strong localizations of the field intensity may occur.
Specific low-frequency response is not yet seen in Fig. 1c although the noise level has considerably grown from its initial value 10 −10 .Up to the time t = 600 ω −1 pe the lowfrequency power has continued to increase, as seen in Fig. 2, with a peak at k s λ D ≈ 0.18, almost twice that of the maximum field intensity.This is a clear indication of parametric Langmuir wave decay.The peak exhibits an appreciable width suggesting that low-frequency fluctuations from a broad range of wave numbers are involved in the decay process.Since a spatial power spectrum at a fixed time does not distinguish between positive and negative wave numbers, backscattered daughter Langmuir waves are hidden in Fig. 2b and falsify the power spectrum.To make them visible Fig. 3 displays a dynamic spectrum for the time period 0 < t < 750 ω −1 pe and the wave number range −0.5 < kλ D < 0.5.The bottom panel of Fig. 3 presents a cut through the color graphic for a frequency slightly below ω pe , selected such that the mother wave attains its highest field strength.The backscatter wave intensity remains 1 order of magnitude below that of the mother wave and does not exhibit a similar sharp spectral peak.This is consistent with the less pronounced maximum of the low-frequency power spectrum (Fig. 2d) and suggests the Langmuir decay is broadband in nature with respect to the wave number.This property and the significantly smaller power in the backscattered wave make a subsequent decay to produce a multiple spectral peak very unlikely.
If one assumes that the waves involved satisfy the standard dispersion relations for Langmuir waves and ion acoustic waves with c s = [k B (T e + 3T i )/m i ] 1/2 the ion acoustic velocity, the resonance conditions for electrostatic decay the beam-driven wave number.One may extract k L λ D ≈ 0.075 from the bottom panel of Fig. 3 which deviates from the value k L λ D = 0.063 resulting from Eqs. ( 4)-( 7) for V b /v the = 16, T i /T e = 0.73.This discrepancy is not unexpected and can be attributed to the changes in the dispersion properties of Langmuir waves associated with the presence of the beam, particularly with the non-Maxwellian character of the electron distribution function due to the plateau formation.An indication of the beam-related change in the Langmuir wave dispersion is the fact that the frequency of the fastest growing mode is always below ω pe , with a downshift that increases with beam density.The choice of the frequency for the cut in Fig. 3 (ω/ω pe ≈ 0.99) reflects this shift.The wave number resonance condition can be satisfied with values from the ranges −0.2 < k L λ D < 0 and 0.1 < k s λ D < 0.3.

Generation of ion acoustic bursts
In the previous section evidence was given that Langmuir wave decay may proceed with some delay after the beam injection.In this section attention is drawn to a mechanism in which the high-frequency field affects heavy particles rather locally.This mechanism is based on the fact that the field energy may occasionally blow up to transient spiky events.The associated strong spatial gradients are the source of ponderomotive forces to which protons respond by the formation of isolated bursty regions of low-frequency electrostatic perturbations.Such events are already seen in the proton density profile in Fig. 2.
To get the spatial distribution of the field intensity in the simulation box, one has to build the time average of the field square E 2 at each grid point over a wave period, i.e., over the time t = 2π ω −1 pe .This quantity is henceforth denoted by < E 2 >.
Figure 4 presents profiles of < E 2 > and proton density for two subsequent times.Notable features seen in this figure are the spatial correlation of the field intensity peaks at the earlier time t = 300 ω −1 pe (a) and finite-amplitude proton density pulses at the later time t = 800 ω −1 pe (b).For example, at positions x/λ D ≈ 750, 1500, 2100 the field energy peaks and the ion acoustic pulses synchronize.This strongly suggests that the latter are generated by the former.The spiky field events appear to be of transient nature, rarely lasting longer than some hundreds of ω −1 pe .The low-frequency perturbations are supposed to propagate with ion sound speed (c s ≈ 0.023 v the ).That is, they do not move very much and may be considered as long-lived events on the timescale used here.Due to their different lifetimes the field peaks do in general not occur in coincidence with low-frequency bursts for longer time.
The detailed spatial shapes of proton density and bulk velocity, shown in the bottom panels, suggest the selected isolated event is a large-amplitude ion acoustic wave pulse with an extent of ≈ 25 λ D .
The signatures of low-frequency activity in the simulations are consistent with nonthermal, sporadic ion acoustic waves associated with bursty Langmuir waves in regions of type III solar radio sources, observed by Lin et al. (1986) and theoretically discussed by Robinson et al. (1993) and Cairns and Robinson (1995).

Time variations
Figure 5 depicts the temporal variation of the longitudinal field at a fixed point (x/λ D = 1200) in the simulation box (a), together with the power spectrum in panel (b).The initial growth phase within the interval 0 < t < 100 ω −1 pe is a result of the interaction of beam and background electrons only.Saturation is achieved by relaxation of the beam distribution to a plateau.Subsequently the system evolves randomly by linear interference of the amplified Langmuir waves.For example, the field drop, beginning at t ∼ 100 ω −1 pe , is not typical and missing in time profiles at other positions.A common characteristic, however, is the longer period of enhanced field in which the low-frequency response takes place.As expected, the power spectrum exhibits a pronounced peak close to ω pe .This peak comprises the frequencies of the beam-driven Langmuir waves, the backscattered daughter waves and the k = 0 mode.Owing to the very small frequency difference between them, no multiple peak can be resolved.Spectral peaks at the harmonics of ω pe have been commonly seen in previous simulations, the signal near the half-harmonic 3.5 ω pe is somewhat curious.
Unlike field and electron density the proton density varies on a much lower timescale.This is demonstrated in Fig. 6, which shows time profiles of the proton density at four equidistant positions within the simulation box.Each of them is Fourier analyzed, and the bottom panel presents the averaged low-frequency power spectrum with the frequency scaled by the proton plasma frequency.Most of the lowfrequency power is seen for frequencies less than the ion plasma frequency, confirming that the proton dynamics is controlled by the ion acoustic mechanism.
As a consequence of the charge neutrality at the slow timescale, the time-averaged electron density is forced to match the proton density which is generally satisfied in the simulation.
Figure 7 illustrates the time history of the energy balance.It shows how the total energy content of the simulation box is distributed over background electrons, beam electrons, ions and the electrostatic field during the evolution of the beam plasma system.The beam energy rapidly turns down in the initial phase and is almost equally supplied to field energy and kinetic energy of background electrons.Ions play a minor role in the energy exchange.Note the limited period of enhanced field energy which lasts until ∼ 800 ω −1 pe .The most intense fields, occasionally up to 2 E 0 , appear in the time period between 250 ω −1 pe and 500 ω −1 pe .During this period Langmuir wave decay takes place and ion acoustic pulses are generated.Subsequently, the field energy is continuously converted into thermal energy of core electrons.This period is accompanied by a loss of coherence.The wave packets seen in the early part of the evolution are disrupted, and the system approaches a phase of almost stationary energy partition with random spatial field fluctuations.A similar energy balance for an extremely low-energy beam was presented by Silin et al. (2007).The beam momentum is gradually transferred to the background electrons causing them to drift slowly in beam direction.

Dense beam effects
As mentioned in the introduction, not all of the frequency spectra of observed Langmuir amplitude time series are consistent with the decay scenario.In this section we present an alternative mechanism able to produce a spectral twin peak without decay.This mechanism, described earlier in the limit of rigid ions (Baumgärtel 2013), comes into play when a more dense beam (n b /n 0 0.05) penetrates the plasma.The frequency ω L of the fastest growing mode is then downshifted and separated from ω pe , the frequency of the k = 0 mode.Due to the higher beam density the latter becomes The energy is in units of the initial thermal energy of the background electrons.E tot remains constant within ∼ 2 %. comparable in amplitude to the beam-driven Langmuir wave.That is, two high-frequency modes with slightly different frequencies are present at an early stage, able to produces a beattype waveform.With respect to the ion dynamics, a different scenario is created.Langmuir wave decay is quenched, and the wave number of the excited low-frequency waves turns from 2k L (for decay) to k L , the wave number of the fastest growing Langmuir mode. Figure 8 illustrates such a situation showing results of a run with a beam of density n b /n 0 = 0.05 and beam velocity V b /v the = 8.The proton density clearly peaks near k L .This can be explained in terms of the spatial shape of the beat waves.If one of the two interfering waves is a k = 0 mode (as plasma oscillations), the spatial variation of the beat wave is determined by the wave number of the second mode.
Note that the spectral double peak here is the result of a linear interaction.That is, this process does not need a threshold field strength to occur and may be seen with weak fields and negligible ion response.It could thus be a candidate to explain the observation of beat-type waveforms without ion acoustic signatures (Sigsbee et al., 2010, their Fig.3c). Figure 9 illustrates the spatial behavior of field and proton density in a way similar to Fig. 1.
A direct comparison between simulations with and without the presence of beam-driven plasma oscillations, i.e., with and without an initial current, was made by Baumgärtel (2013) in the limit of immobile ions.

Discussion
The simulation replaces in a certain sense an observation under rigorously idealized conditions, such as one-dimensional geometry, a uniform region for the beam-plasma interaction and the electrostatic approach.A handicap of a simulation experiment is that spectral peaks are not always suitable for comparison with a real observation because of frequency changes due to Doppler shift.Another weakness concerns the unrealistically high field strengths produced in the simulation which exceed considerably the observed ones.For a real space plasma embedded in a magnetic field, the onedimensional geometry, the restriction to purely longitudinal waves and the assumption of both a uniform background plasma and a uniform beam may be not a good approximation.Thus the field intensities in the simulation should rather be considered as an upper limit of what might be reached under realistic conditions.
In addition, the electrostatic approximation prevents the model from describing Langmuir-like waves with complex polarization, e.g., Langmuir/z-modes in weakly magnetized plasmas, which appear with strong transverse field components for high beam velocities and small wave numbers.Recently observed waveforms with multiple spectral peaks were interpreted in terms of these waves (Malaspina et al., 2011;Graham and Cairns, 2013a;Kellog et al., 2013;Layden et al., 2013).
An advantage of simulations, on the other hand, is that, unlike one-point observations in space plasmas, spatial structures of field and density can be inspected and examined.Evidence for Langmuir wave decay came about by the analysis of spatial structures of field and proton density, the signatures of which were consistent with the wave number resonance conditions.
The presented results pertain mainly to the parameter set listed in Table 1, but a series of runs with other parameter combinations have been carried out.Since according to Eq. ( 7) the wave number of the fastest growing Langmuir mode is inversely proportional to V b , increasing (decreasing) beam velocity results in extended (reduced) characteristic spatial scales, e.g., the width of the wave packets.A less energetic beam with V b /v the = 8 and n b /n 0 = 0.0057 appears too weak for decay to proceed.It produces a broadband proton density power spectrum without decay signatures.A very fast and thin beam with 2 times the standard speed (V b /v the = 32, i.e., V b /c ≈ 0.14, n b /n 0 = 0.001) turned out to produce a field/ion density pattern similar to that of the standard case, including decay.
A crucial parameter for the onset of nonlinear processes such as decay or generation of ion acoustic perturbations is the initial ratio W beam /W back of beam energy to kinetic energy of background electrons, which controls the saturation level of the field intensity.It takes the value 1.47 for the standard parameters and 0.37 in the weak beam case where decay is missing.This suggests a threshold value on the order of 1.
The evolution of the beam plasma system appears to be insensitive to the ion-to-electron temperature ratio, as runs in the range 0.1 < T i /T e < 1 indicate.The same is true for the beam temperature.Simulation results do not change when the warm beam (v thb = v the ) is replaced by a cold one.
The longitudinal field calculated from the current Eq. ( 1) was repeatedly compared to the field resulting from the Poisson Eq. ( 2).There was always an excellent agreement.The relatively small amount of numerical effort required by the implicit PIC algorithm used here is at the cost of the energy conservation which is not exactly ensured.The total energy (field + electrons + protons) is conserved during runs until 2000 ω −1 pe typically by ∼ 2 %.The final electron velocity distribution shows a plateau with small slope but no high-energy tail.
Since both beam energy and beam momentum are continuously transferred to the background electrons, the electrostatic 1-D simulations suggest that the beam is eventually completely absorbed.

Summary
In this paper the self-consistent evolution of an electron beam-plasma system including ions is investigated by particle-in-cell (PIC) simulations using a one-dimensional electrostatic standard particle code with periodic boundary conditions.Initial conditions are set up by background electrons and protons at rest -thus allowing an initial current due to the beam.The initial growth of Langmuir waves is according to the classical electron beam instability without participation of protons.Saturation is achieved by rapid relaxation of the beam velocity distribution to a plateau, accompanied by a loss of beam energy.A wave packet structure develops as a result of linear interference of waves that have grown from noise.In case of sufficiently intense fields, nonlinear processes including ion response such as Langmuir wave decay can occur.Spatial profiles of field and proton density are inspected, and it was found that for low-density beams the wave numbers of the forward and backward propagating Langmuir waves and ion sound modes match the resonance conditions for decay.In time representation, Langmuir wave decay cannot be seen in the plasma frame owing to the very low frequency difference between beam-driven waves and backscattered waves.
With growing time the wave envelope undergoes significant variations and may become disrupted.Spiky field events can transiently occur with strong local increase of field intensity.Heavy particles respond to the associated ponderomotive forces by generation of ion acoustic bursts.Due to the different space/time variation of field energy spikes and lowfrequency bursts, there is in general no coincidence between them for longer time.
For dense beams where the frequency of the fastest growing mode is downshifted from ω pe , beam-aligned plasma oscillations (k = 0) come into play and establish early in the course of the instability the presence of two high-frequency waves of comparable amplitude and slightly different frequency, resulting in a beat wave structure.The corresponding spectral double peak is visible in the plasma frame.The most power of the low-frequency fluctuations is then concentrated near to the wave number of the dominant Langmuir wave instead of 2 times this value as for decay.
For nonlinear effects and ion response to become important in the beam-plasma interaction, the simulations suggest a ratio of beam energy to thermal energy of core electrons no less than 1.
The simulation confirms that the beam-plasma system goes through several characteristic phases on its way to a state without net energy exchange between particles and field.
Cases with repeatedly occurring decay resulting in multiple distinct spectral peaks as described by Graham and Cairns (2013a) could not be reproduced in the simulations, probably as consequence of the electrostatic limit.Waveforms of Langmuir eigenmodes in density cavities as observed and theoretically explained by Ergun et al. (2008) and Graham and Cairns (2013b) were not a target of the simulations.

Figure 1 .
Figure 1.Spatial profiles of field (a), proton density (c) and associated power spectra (b), (d) at an early time t = 200 ω −1 pe of the beam-plasma evolution.The almost coherent field structure is dominated by the wave number kλ D ≈ 0.075 of the fastest growing mode.Beam velocity V b = 16 v the , beam density n b = 0.0057 n 0 .

Figure 2 .
Figure 2. The same as Fig. 1 at a later time t = 600 ω −1 pe .The most power of the low-frequency response is concentrated near kλ D ≈ 0.18, almost twice the wave number of the fastest growing mode, suggesting Langmuir wave decay.The power spectrum (b) is falsified by backscattered waves.

Figure 3 .
Figure 3. ω − k dispersion diagram for the longitudinal field within the time interval 0 < t < 750 ω −1pe (top) and cut through the plot for the frequency ω/ω pe ≈ 0.99 for which the mother wave attains its maximum amplitude (bottom).Backscattered Langmuir waves appear for kλ D < 0. Parameters as in Fig.1.

Figure 4 .
Figure 4. Spatial profiles of field energy and proton density at times t = 300 ω −1 pe (a) and t = 800 ω −1 pe (b), illustrating the formation and relaxation of spiky field events and the associated proton response.Bottom panels: spatial variation of proton density and bulk velocity of the ion acoustic pulse at x/λ D ≈ 750 in better resolution.Parameters as in Fig. 1.

Figure 5 .
Figure 5.Time variation of the field at a fixed position x/λ D = 1200 in the box (a) and corresponding power spectrum (b).Parameters as in Fig. 1.

Figure 6 .
Figure 6.Time profiles of proton density at four positions x/λ D = 600, 1200, 1800, 2400 within the box (top) and power spectrum averaged over the four profiles; frequency normalized by the ion plasma frequency ω pi (bottom).Parameters as in Fig. 1.

Figure 7 .
Figure 7. Distribution of the energy content W tot of the simulation box over background electrons (W back ), beam electrons (W beam ), ions (W ion ) and the electrostatic field (W field ) as function of time.The energy is in units of the initial thermal energy of the background electrons.E tot remains constant within ∼ 2 %.

Figure 8 .
Figure 8. Beat-type waveform generated by a dense beam n b /n 0 = 0.05 at position x/λ D = 2400 in the box (a), and corresponding power spectrum (b) with two distinct peaks.Beam velocity V b = 8 v the .

Figure 9 .
Figure 9.The same as Fig. 1 but for a dense beam; n b /n 0 = 0.05, V b = 16 v the .The low-frequency power spectrum (d) peaks at the wave number of the mother wave.