Large-scale simulations of 2-D fully kinetic Farley-Buneman turbulence

Currents flowing in the Earth’s ionospheric electrojets often develop Farley-Buneman (FB) streaming instabilities and become turbulent. The resulting electron density irregularities cause these regions to readily scatter VHF and UHF radar signals. Many of the observed characteristics of these radar measurements result from the nonlinear behavior of this plasma. This paper describes a set of high-resolution, 2-D, fully kinetic simulations of electric field driven turbulence in the electrojet. These show the saturated amplitude of the waves; coupling between linearly growing modes and damped modes; the evolution of the system from dominance by shorter (1 m–5 m) to longer (10 m–200 m) wavelength modes; and the propagation of the dominant modes at phase velocities that lie below the linearly predicted phase velocity and close to but slightly above the acoustic velocity. These simulations reproduce many of the observational characteristics of type 1 waves. They provide information useful in accurately modeling FB turbulence and demonstrate the significant progress we have made in simulating the electrojet.


Introduction
In the E-region ionosphere, turbulent processes driven by strong ambient DC electric fields create plasma density irregularities responsible for type 1 radar echoes.These irregularities have been studied experimentally and theoretically for five decades.In the last decade, numerical simulations became an important tool in exploring the nonlinear behavior of E-region instabilities.However, these simulations were investigation and then presents results from a number of simulations.We then conclude with a brief discussion and summary.We plan to follow this paper with another one discussing the implications of these results for theory and observations.In the near future, we also plan to run a set of 3-D simulations and compare the results to these 2-D simulations.

Background
Plasma irregularities in the E-region ionosphere were first detected shortly after the invention of radar in the 1940s and detailed radar studies of this phenomenon began in the 1950s (Bowles, 1954).Since that time hundreds of studies, both experimental and theoretical, have been undertaken to further our understanding of the origin, evolution, and effects of E-region irregularities.These studies have been reported on in hundreds of papers and a number of books (Fejer et al., 1984;Farley, 1985;Kelley, 1989).In the last 15 years, simulations have become a useful tool for exploring the nonlinear evolution of E-region waves.For a general discussion of background material appropriate to this paper, please refer to Oppenheim and Dimant (2004).In this paper, we will limit our discussion of background to material that has appeared in the last few years.
On the experimental side, there has been a concerted effort over the last few years to simultaneously use multiple frequency radars, imaging radars, and Faraday rotation to study electrojet turbulence at the Jicamarca Radio Observatory (Hysell et al., 2007).They found structuring of the system at the 3-5 km scale size and east-west asymmetries which may be explained by thermal effects (Dimant and Oppenheim, 2004).Woodman and Chau (2002) and Hysell et al. (2007) present compelling evidence that the phase velocities of the waves go as C s cos θ where C s is an estimate of the ion acoustic speed and θ is the angle between the direction of the fastest waves and the measured wave.
A number of theories have been advanced to better explain the behavior of Farley-Buneman (FB) turbulence in the last few years.Otani and Oppenheim (2006) further developed the theory of three or four mode coupling mechanisms of saturation (Hamza and St.-Maurice, 1993;Sahr and Farley, 1995;Dimant, 2000).This paper shows how a FB system, limited to very few modes will saturate and show behavior such as density perturbation levels and phase velocities which match observations quite nicely.They also show how mode-coupled systems have constraints on phase velocities and wave directions.However, these systems never resemble the fully turbulent simulations that include diffusive and kinetic effects neglected in the mode coupling models.Hysell and Drexler (2006) further developed the non-spectral "blob" approach to modeling the saturated FB system (St.-Maurice and Hamza, 2001), showing how the blobs rotate and distort.While these models show some interesting aspects of non-linear dynamics of FB turbulence, they also exclude diffusion and kinetic effects.The net result is that the "blob" models do not resemble the simulated FB turbulence.
Finally, there have been a number of recent papers commenting on interesting aspects of FB turbulence.Haldoupis et al. (2005) looks at the role of density gradients on FB waves.Dyrud et al. (2006) shows the similarities between rocket data and simulations of FB turbulence where the electrons are modeled as an isothermal fluid.Bahcivan and Hysell (2006) discusses the physics of gradient drift wave driven FB turbulence.Fontenla (2005) suggests that a form of FB turbulence heats the Sun's chromosphere.

Simulation methods and limitations
To study the growth, evolution and saturation of E-region waves, we use a particle-in-cell simulator to follow the dynamics of a collection of ions and electrons as they respond to a driving electric field and their own self-generated electric field.To reduce the computational expense, we simulate motion within the 2-D plane perpendicular to the geomagnetic field.This includes the important nonlinear wave coupling between perpendicular modes, but precludes possible excitation of modes oblique to the geomagnetic field, B 0 (Oppenheim et al., 1996).This means that the saturated amplitudes of the resulting waves probably do not reflect the amplitude in a true 3-D system.The 2-D simulations presented allow us to explore the competing roles of 2-D mode coupling, linear growth, and thermal effects.Further, these simulations present the most physically accurate simulations of field-driven electrojet turbulence yet performed.
In Oppenheim and Dimant (2004) we discuss in detail the use and limitation of our fully kinetic 2-D PIC code.For this research we improved the parallel processing approach by domain decomposing the problem.This means that we can have different sets of processors working on different spatial regions of the problem.This allows us to reach much larger sizes and numbers of PIC particles.We call this code the Electrostatic Parallel PIC or EPPIC code.

Simulation configuration
We start with an analysis of a baseline simulation with parameters appropriate for a run at 101 km altitude in the auroral region.This run has the parameters described in Table 1.All subsequent simulations will be a variation of this baseline simulation.While these runs use auroral parameters, they apply equally well to the equatorial electrojet at an altitude approximately 4km higher.In particular, the more weakly driven runs probably show the dominant physics of the equatorial electrojet.The simulations use the neutral wind frame of reference where x points in the E 0 ×B 0 direction, ŷ along E 0 , and ẑ along B 0 .
This simulated density evolves from an initial state of white noise through a linear growth stage to a turbulently saturated final state.Figure 1 shows the density at two times in this evolution, the first shortly before saturation and the last well into saturation.During the evolution many processes occur, from the initial anisotropic heating of the ions and electrons to cascade processes which couple energy across many scales.By comparing many simulations, we will describe the dominant processes involved in saturation.While the details of how these simulations evolve depends on the parameters, qualitatively all runs show essentially the same behavior as our baseline case.
Tracking electrons as they gyrate dominates the computational cost of these simulations.To reduce this expense, we use an electron mass which is 44 times more massive than the physical electron mass, making the ion to electron mass ratio close to that found in a H + plasma.Our theoretical analysis says that one can raise the electron mass without substantially impacting this system until the point where electron Landau damping becomes important, at about 80 m e .We performed a test to demonstrate that this elevated electron mass has little effect on these simulations.We ran the baseline case with m sim e =6 m e .This run showed no appreciable changes from the baseline case, showing the same growth rate, saturation amplitudes, and dominant modes.As we reduce the electron mass, we also reduce the frequency of updating the ion velocities and positions.This allows us to increase the number of PIC ions (not changing the density) without increasing the computational cost.The increased number of ions reduces the level of particle noise, allowing  at Boston University.We will show images from this large run when appropriate.

Saturation amplitude
The average perturbation electric field amplitude, E 2 1/2 gives information on the nature of the turbulence.Figure 2a shows that the average field reaches 0.08 V/m for the 4096×4096 run.This exceeds the driving electric field of 0.05 V/m by 60%.For the 1024×1024 baseline run with the same parameters, the rms electric field reached a more modest 0.05 mV/m, only 25% above the driving field.It remained the same for the run with the more realistic ion/electron mass ratio.This high rms field is a common feature of 2-D kinetic simulations of the FB instability.The 12% maximum density perturbation seen in Fig. 2 arises from the PIC noise throughout the mesh.Since we have 16 million cells, this maximum should reach approximately 4 times the RMS particle noise level.The RMS value is a far better indicator of the actual noise.The difference between the maximum electron and ion densities results from using 8 times as many PIC particles to model ion than electron dynamics.We do this because the slow ions iterate only once for every 32 electron timesteps and, hence, we use 8 times as many ions at little computational cost.Despite this difference in maximum values, the system remains highly quasineutral and the saturated electron and ion RMS densities track each other well.
The E 2 1/2 value reflects the combination of the turbulent electric field plus that resulting from individual particle noise.In EPPIC, the ions are cycled less frequently than the electrons since the ions do not need to resolve the electron cyclotron frequency.The ions need a time step which resolves ion time scales such as the ion acoustic frequency, the ion plasma frequency, and the ion collision rate.The electric field applied to the ions is the field averaged over the intervening electron time-steps.Therefore, the ions experience a less noisy electric field than do the electrons.Because the ions time step is much slower than the electron time step, the number of PIC-particles representing the ion density can be increased without substantial computational expense.
The averaged density perturbation of the turbulent plasma can be compared to the perturbation amplitude of the electric field, theoretical estimates, and in-situ measurements by rockets.Figure 2b shows the time evolution of the peak plasma density reached as well as the averaged, (n−n 0 ) 2 /n 2 0 dx.From these simulations, we estimate that the saturated average density perturbations is ∼12% while the initial, non-turbulent perturbation is ∼2.5%.The maximum density perturbation is typically ∼40% while the predisturbance perturbation level is ∼12%.
Examining these perturbation levels for a series of simulations with different electric field drivers should reveal the relationship between this driver and the perturbation amplitudes.Figure 3 shows the average electric field amplitude.One sees that the RMS electric field exceeds E 0 by 20% to greater than 50% where the higher driving fields cause the greater excesses.The RMS density perturbations start at a modest 6% and rise to 40%.The RMS electric field does not change as a function of electron mass.
The 4096×4096 run with baseline simulation parameters shows a 50% higher saturation electric fields E 2 1/2 , and a similar maximum |E| field as the baseline case.The average ion density perturbation n for the larger run drops 50% to 8% but the peak maximum density perturbation remains similar.Unfortunately, the large run also shows evidence that it has not yet saturated at all wavelengths as max(n) maintain a small upward trend though E 2 1/2 does not.

Spectra of the saturated system
Radar data typically probes 1 or a few wavevectors of the system with sensitivity to the phase velocity of those waves.Rockets can measure many wavelengths in-situ for a limited time and at a single moving point.Both these measurements gather spectral data which we plan to compare to the simulation spectra.
Spectra from two spatial dimensions plus time such as density, n i (x, y, t), form 3-D spectra, n i (k x , k y , ω), which are challenging to display.One can show cross sections or various averages of the data.We will present a number of these as we attempt to extract the most pertinent information.
Radars typically measure electrojet waves at a single wavelength but frequently from a number of directions.For comparison purposes, we will show analogous information generated from these simulations.Figure 5 shows the spectral information from the 4096×4096 simulation for 4 wavelengths (i.e.|k| values).We can see from these spectra that the phase velocity changes roughly as a cos(θ−φ) relationship where φ is the angle of the phase shift due to the turning of the wave.This turning results from a modification of the linear FB relationship due to thermal effects, mostly ion heating (Oppenheim and Dimant, 2004).One can also see this turning in the density images of Fig. 1.One sees from the top 4 images of Fig. 5 that the largest amplitude signal will return from nearly in the E×B direction and will fall off as one measures in the direction perpendicular to this.The power falls off more rapidly for modes near from the dominant, linearly growing, wavelength (∼1.5 m) then for longer or shorter wavelengths.While the cosine relationship of phase velocity agrees with the results of Woodman and Chau (2002), the fall off in energy may not.First Quarter Last Quarter Fig. 4. Spectra showing the t 1 t 0 n i (k x , k y , t)dt for the first fourth of the simulation (t 0 =0, t 1 =t final /4) and the last fourth (t 0 =t final * 3/4, t 1 =t final ).The logarithmic scale tends to deemphasize the strongest modes which in the right image are at small values of k.
To illustrate the average behavior of these simulations, we show the time averaged spatial spectra from the first quarter of the simulation and the last quarter of the simulation in Fig. 4. One sees from these the downward tilt of the waves and how the system evolves from having predominantly short wavelength modes as one would expect from linear kinetic theory to longer waves distributed inhomogeneously but smoothly through k-space.When we look at animations of (k x , k y ) evolving in time, we see a group of waves growing with a range of wavenumbers initially centered at k=(3.4,−0.1).Then within a few milliseconds the dominant mode slides in k-space towards k=(0, 0) eventually settling on the distribution show on the right of Fig. 4.  Ion Thermal Energies Fig. 6.Electron and ion thermal energy of large baseline simulation.These are expressed as m v 2 j /K b where j is the direction unit vector, m is the species mass and K b is the Boltzmann constant.The temperatures rise from 300 K, initially because of frictional heating due to E 0 ŷ, and further due to wave heating.Our predicted electron and ion temperatures are very similar to those that develop in the simulations.Note that these simulations include v z values for each particle.These values change only because collisions with neutrals slowly make the distributions more isotropic (no E z exist).
The peak velocity is also listed on each frame of the figure.This velocity exceeds the acoustic speed of the system if one calculates that speed using the temperature of the neutrals and for isothermal electrons and ions, C sn =407 m/s.However, neither species is truly isothermal, nor do they remain at the neutral temperature of 300 K. Figure 6 shows the temperature evolution of the large baseline simulation.If we adjust for the enhanced temperature and assume the electrons are fully adiabatic while the ions are isothermal, then C s ≈500 m/s.While these assumptions are not precisely correct, a dramatically higher acoustic speed seems incorrect.For short wavelengths, the ions will not have time to fully isotropize during a wave cycle.If the ions behaved as a fully 3-D adiabatic gas then Cs≈550 m/s, while if they were a 1-D adiabatic gas then C s ≈650 m/s.In all cases, these peak velocities remain well below the velocity predicted by linear theory, V d /(1+ )=950 m/s.We will analyze the phase velocities more completely in the discussion section.
If these simulations were fully 3-D, we would expect the electron temperatures to rise more dramatically as electrons accelerate due to the small wave-driven parallel electric field.Even in a 2-D simulation with no E || , the moments of the electron distribution increases because of wave induced compression and inelastic collisions between electrons and neutrals (see Dimant and Oppenheim, 2004, for more details).
To examine phase velocity over a range of wavelengths, we can compare the wave frequency, ω, from the saturated component of the simulation to wave number, |k|. Figure 7 shows these spectra.Clearly most of the energy falls within the acoustic modes with a ∼700 m/s velocity.One can also see that there is substantial spectral broadening, particularly at low k.At high k, the acoustic speed may increase as  angles from the largest simulation.We generate these using two methods.The first mode takes the second moment of n(|k|, ω) for a given k, where ω 0 is the first moment of n(|k|, ω) and )dω is the normalization factor.Since low amplitude noise at large velocities can cause the second moment to become dramatically larger than one would obtain if the noise level were lower, we also fitted n(k, ω) to the sum of a Gaussian plus a quadratic polynomial.The standard deviation of the fitted Gaussian gives a better representation of the spectral width.One sees that the simulation generates substantial spectral broadening in all directions and at all wavelengths.The spectrum is narrowest where it is strongest, in the direction of greatest linear growth, and broadest perpendicular to that direction.Shorter wavelengths have less spectral broadening but the difference between 1.5 and 6 m modes is less than a factor of two.
Rockets generate a spectrum by comparing k to ω for the time-series data they collect.Simulations can generate similar spectral information.Figure 9 shows such a comparison generated by Fourier transforming the n i (x, y, t) array from the saturated part of the simulation into a n i (k x , k y , ω) array.We then interpolate along paths of fixed |k| to obtain a n i (|k|, ω) array.We then find the maximum and average values of these arrays and plot them to make Fig. 9. Making a direct comparison between these simulation results and the rocket spectrum will require substantial further analysis.
These simulations all show energy flowing into longest horizontal waves which the simulation supports.For the 4096×4096 simulation, this means that energy flows into waves with wavelengths approaching 160 m.This mode eventually becomes the dominant mode in the system.Figure 10 compares the growth of the longest mode to one of the dominant linearly driven modes from the baseline simulation.One can see that the longer mode develops much later than the shorter one.However, its growth rate is similar or faster than the shorter mode.For long wavelengths, the predicted linear growth rate is proportional to k 2 , so that, if linear growth were the cause of the long wavelength modes growth, we would expect this mode to grow at 1/436 the speed of the longer mode (neglecting kinetic reduction of the short wavelength growth rate).The 40 m wave (dashed line) also grows at the same rate as the 160 m wave, implying that these are not linear FB growth processes, but instead result from coupling to shorter modes.
For wavelengths much less than the size of the simulation, saturation appears to result from coupling to other modes or from turning of the wave (not a completely separate idea).
Waves which span the entire simulation model cannot saturate in the same way since they cannot couple to longer waves nor can they smoothly turn since the simulation does not support a broad spectrum of waves at such long wavelengths.These longer waves can saturate through a alternate mechanisms including coupling to shorter wavelengths or particle heating.We are not certain of the mechanism.

Discussion and conclusions
These simulations show the nonlinear behavior of saturated FB waves in 2-D.While one can modify the parameters to change the driving field and the effective altitude, these changes do not dramatically impact the qualitative behavior of the waves.This is because the waves dissipate their energy by coupling to other modes both at longer and similar wavelengths.In this section, we briefly discuss some interesting features of the saturated FB waves found in our 2-D simulations, estimate basic quantitative characteristics of these waves, and compare our estimates with these simulation results.
The most intriguing feature of our simulations is the transfer of turbulent energy from short-wavelength waves (preferred by the linear instability growth rates) to much longer waves.In the saturated state, the maximum of turbulent energy always shifts toward the longest waves allowed by the simulation box size.As we increase our box size, not only does the most energetic wave become longer but the rms turbulent electric field becomes larger.In the real electrojet, the 3-D geometry and natural gradients of the background medium may impose a preferred maximum wavelength.
What physical mechanism causes inverse cascading from short to long wavelength?For long wavelengths the linear instability growth rate scales as the wavenumber squared, γ ∝k 2 .If one waits long enough, one would expect long wavelengths to develop.However, in the simulations these waves develop much faster than predicted by the linear theory, suggesting that nonlinear mode-coupling plays a crucial role in their development through an inverse cascade mechanism.We do not yet have a model explaining this aspect of the simulations.
We ran these simulations with a range of driving electric fields, E 0 and found that the threshold electric field necessary to generate instabilities is nearly twice the minimum isothermal fluid-model value, (for the undisturbed temperatures T e0 =T i0 =T n =300 K, m i ≈30 amu, and B 0 =0.5×10 5 nT).There may be several reasons for the discrepancy.First, the driving field itself heats both electrons and ions, increasing the threshold.In our simulations, however, this initial heating was fairly small, ∼10−20 K for each species.Second, the small charge separation between electrons and ions, disregarded in quasineutral fluid theory, can also increase the threshold field by a factor Mode (20,-6) Fig. 10.n i (k, t) vs. time for two modes, the solid line shows the (1,0) mode which is at k=(2π/41, 0) m −1 , the dotted line shows the (20,−6) mode at k=(2π/2.0, 2π/6.8)m −1 and the dashed line shows the (4,0) mode. of , where ω pi is the ion plasma frequency and the critical density is given by n cr =m i ν 2 i 0 /e 2 ≈4.84×10 7 m −3 (Rosenberg and Chow, 1998).In our simulations, we had (1−n cr /n 0 ) −1/2 ≈1.05, so that this factor should not play an important role.Third, and most importantly, the non-isothermal and kinetic behavior of plasma appears to substantially elevate the instability threshold.If both electrons and ions behaved adiabatically, this would increase the threshold value by a factor of γ 1/2 e =γ 1/2 i =(5/3) 1/2 ≈1.29.In these collision-dominated, low-frequency, E-region processes, however, the situation is more complex and neither the electrons nor ions behave a simple adiabatic or isothermal fluid.This is also true in these simulations as discussed below.
Both ion and electron thermal behaviors modify the instability threshold.The simplified fluid-model analysis in Dimant and Oppenheim (2004) shows how ion thermal effects cause the predominant F-B waves to travel in a different direction from E 0 ×B 0 -drift velocity but impacts the instability threshold only a little.Electron thermal factors calculated using both kinetic and fluid models give a corrected value of electron adiabatic index γ e which will raise the threshold substantially (Dimant and Sudan, 1995a,b;Kagan and St.-Maurice, 2004).In these simulations, we modeled electron-neutral collisions using an elastic hardsphere interaction module which gives a kinetic collision frequency proportional to the electron velocity.We matched the inelastic collision and heat transfer rates by colliding the electrons with a reduced mass neutral species.Using Eqs. ( 61) to (64) from Dimant and Sudan (1995b) with α=1/2 and x=0, we obtain γ e =1+P ω,k ≈2.58.This quasi-adiabatic coefficient is valid for waves in the range of wavelengths, ω−kV 0 ≈ kV 0 /(1+ ) ν i .For shorter wavelengths, the kinetic effects of ion Landau damping increases the threshold drastically.Note that for the longer wavelengths -the ones that dominate these simulationswhere (2m e /m * n )ν i ω−kV 0 , we expect electron cooling to reduce γ e which will decrease the threshold.The combination of all listed factors can increase the threshold by a factor somewhat less that 2, though there may be some additional unidentified factors.Also, bear in mind that it is difficult to determine the exact value of the threshold electric field from simulations because small instability growth rates close to the threshold require long times for the instability to develop.Finally, fluctuation noise associated with the finite number of PIC particles can preclude the measurement of waves driven near their instability threshold, though with these massively parallel simulations, we can drive that noise down to a fraction of a percent.
Figure 6 shows that heating occurs in a two-step process for both electrons and ions.The initial heating stage results from the frictional heating caused by the driving electric field itself.Then, after several tens of milliseconds, the wave electric field causes additional heating.Equation (34) from Dimant and Milikh (2003) predicts an initial heating for ions of 12 K which agrees reasonably well with the anisotropic temperature increases of 20 K for V y and 11K for V x and V z .The anisotropy occurs because the ion Pedersen drift velocity (v Ped 95 m/s) is nearly one third of the mean thermal speed defined by v T i =(T i /m i ) 1/2 .In the highlycollisional regime, when v Ped is not small compared with v T i , the distribution function can easily deviate from a simple Maxwellian and become anisotropic.
The electron velocity distribution is much more isotropic than that of ions because both collisions of light lowenergy electrons with heavy neutrals results mainly in angular scattering and relatively slow energy changes.The energy exchange rate (δ en in Dimant and Oppenheim (2004)) has been modeled in our code by an artificially reduced neutral mass, m * , so that in our baseline simulation δ en =2m e /m * n ≈8×10 −2 (this value exceeds the actual E-region values by one to two orders of magnitude but compensates appropriately for the artificially massive electrons).Using Eq. ( 37) from Dimant and Milikh (2003), we obtain T i −T 0 25 K, which agrees well with the initial electron heating.The characteristic times of this initial heating, (δ en ν e ) −1 5 ms for electrons and ν −1 i 0.6 ms for ions also agree well with the observed times of the initial heating, as seen in Fig. 6.
The 2-D turbulent heating level can be estimated from the same formulas using the rms values of the turbulent electric field E instead of E 0 .According to Fig. 2, E reaches 84 mV/m, which yields an estimate for the final electron heating of 95 K and for ions 45 K.These estimates give a slightly lower value for electron temperature and exaggerate significantly ion heating.This may be due to oversimplification of the frictional heating model in the theory which does not account for deviations from Maxwellian distributions or the rather large ratio of the drift velocity to the thermal speed.Note that the rms turbulent electric field in these simulations was larger than one would expect from simple heuristic reasonings and larger than those observed by rockets.We hope that 3-D simulations will reduce the rms field value because the third dimension creates an important degree of freedom for instability saturation.
These heating effects may help explain the elevated value of the wave phase velocity compared to the cold-plasma ion-acoustic speed.Using the final elevated temperatures along with kinetically calculated γ e , as we did above, we can easily obtain a renormalized ion-acoustic speed of C s 700−750 m/s.Assuming that most energetic waves settle on the margin of linear stability for, V ph ≈C s , we obtain the peak values listed on Fig. 5.A more accurate analysis is needed, but even these simple estimates show that a conventional assumption of peak V ph ≈C s with the underlying reasonings do not contradict our 2-D simulations.
The principle limitation on these results arises from the simulations being only 2-D.We expect to address this and further analyze these 2-D results in future papers.

Fig. 1 .
Fig. 1.Plasma density as a function of position at two different times.The color bar to the right indicates the ion density perturbation amplitude.Note that the maximum values of the density perturbation, δn i reflect approximately four times the RMS value of δn i because of statistical variations on large meshes.

Fig. 2 .
Fig. 2. Electric field energy and density perturbation maximum and average from baseline simulation.Left image shows the time evolution of the maximum (top) line and averaged electric field squared in (V/m) 2 .Right image shows the evolution of the maximum (top lines) and RMS of the perturbed density (bottom lines).The electron densities are shown with solid lines while ion densities with dashed lines.The base-line noise level of the simulation can be assessed from 10 ms to 25 ms.
Fig. 3. Saturated RMS electric field and density perturbation levels from four simulations.Left image shows the RMS electric field while the right image shows the RMS of the perturbed density versus driving electric field.Note that the label of the left image, E represents E 2 1/2 .

Fig. 5 .
Fig. 5. Spectra comparing v ph =ω/|k| to angle with respect to the E 0 ×B 0 direction for four wavelengths.The top four figures indicate the strength of the signal with the color bar to the right of the image.The bottom four normalize the amplitude of the signal at each angle to the maximum signal strength at that angle, in order to better compare with radar data which generally show normalized spectral amplitudes.The curve drawn through the signal is the moment of the signal at each angle and this indicates the mean phase velocity.The peak value of v ph is listed on the bottom of the figure.

Fig. 7 .Fig. 8 .
Fig. 7. Power as a function of |k| and ω.This is generated by summing over all the angles from spectra like those seen in Fig. 5. Lines for C sa =500 m/s, V d /(1+ )=950 m/s and C smax =650 m/s (dashed) have been placed on the figure for reference.

Fig. 9 .
Fig.9.Power as a function of |k| and ω.This is generated by summing over all the angles from spectra like those seen in Fig.5.The "max" curve shows the n(|k|, ω) with the highest ω while the "avg" curve averages n(|k|, ω) over all ω.

Table 1 .
Parameters of Baseline Simulation.Here ≡ν e ν i / e i and N e =5.12×10 8 and N PIC i =4.1×10 9 ).This is required 1024 processors running for approximately 35 h the equivalent of over 4 years of single CPU processor time.These runs were performed on the IBM Blue Gene supercomputer www.ann-geophys.net/26/543/2008/Ann.Geophys., 26, 543-553, 2008