Back to the home page
Previous simulations of the instability have not been able to answer these questions definitively. In the simulations of Newman and Ott [1981], a numerical error was made in the potential field solver, calling into question the validity of the results. The simulations of Machida and Goertz [1988] never explicitly addressed questions associated with the saturation level. They were instead interested in electron heating dynamics, a topic we do not consider in detail here. They modeled the plane parallel to magnetic field, whereas we model the perpendicular plane. The simulations of Janhunen [1994] employed particle electrons as well as ions. The system was therefore required to be quite small in physical dimensions, and thus was not able to provide the larger view of how the saturated state interacts with the longer wavelengths.
The electrostatic potential was found by solving div(J) = 0 on a grid. The total current density, J, consisted of an ion current, accumulated from the ion particles onto the grid, and the electron current, described by the inertialess electron momentum equation. The resulting equation is formally the same as Poisson's equation, with the exception of two terms: one which takes the form grad(n).grad(phi), and one which is of magnitude of kappa*grad(n) x grad(phi), where kappa is the the ratio of the electron cyclotron frequency to the electron collision frequency. This second term represents the divergence of the nonlinear electron Hall current. The factor kappa is huge in the topside equatorial electrojet--typically of order 100. This makes the nonlinear electron Hall current term the dominant nonlinear term in the field equation.
When the waves in the simulation system are of small amplitude, this term can be small compared to the linear terms. In this case, an iterative method works, in which the two nonlinear terms are evaluated using the potential from the previous iteration and added into the source term. The resulting equation is then formally Poisson's equation and is easily solved for the next potential iterate using, for example, a pseudospectral method.
For relatively modest wave amplitudes, ((delta n)/n_0 of order 0.01), however, this method no longer converges because of the large size of kappa. The two nonlinear terms must then be grouped together with the del^2(phi) term. The resulting linear operator on phi varies from nearly symmetric when wave activity is low to nearly antisymmetric during moderate wave activity, making the equation a difficult one to solve. After trying several matrix inversion packages, an effective method was found to be one described by Freund and Nachtigal [1994], termed the quasi-minimal residual (QMR) method. This package has been generally effective and acceptable in speed in our particle simulations.
Another complication encountered in the field solve comes from the fact that the field equation is singular. If phi(x,y) is a solution, then so is phi(x,y)+const. This leads to the following numerical pitfall: it is easy to construct a finite difference scheme which does not reflect the singular quality of the differential operator. The corresponding finite difference operator would typically be "nearly" singular and thus could potentially generate very inaccurate solutions. Alternatively, the matrix representing the finite difference operator could be reduced by one row and column, thereby reducing the complete set of finite difference equations down to a non-singular set. The problem in this case is that the eliminated equation may not be "close to" a linear combination of the remaining equations, implying that the grid point to which the removed equation corresponded may not be treated properly. These were the sorts of issues that were not dealt with correctly in Newman and Ott [1981].
These problems can be bypassed by making sure that the finite difference equations are singular in the same sense as the original differential equation. We employed a finite volume method, which preserves the mathematical structure of the two derivatives involved in the differential equation. We can show rigorously that the resulting difference equations have the following necessary properties: (1) if phi(x,y) is a solution to the finite difference equations, then so is phi(x,y)+const., and (2) the sum any N-1 of the N finite difference equations yields the remaining equation. A consistent method may thus be obtained by removing one row and column from the finite volume set of difference equations.
All the simulations we describe here are available as movies on the World Wide Web. The address is http://buteo.colorado.edu/~meers/fbi-movies.html.
At the beginning of the standard simulation, the rms density perturbations are observed to grow linearly with time with a growth rate consistent with linear theory. The instability then saturates with an rms density perturbation of about 5%.
A number of interesting features appear in the saturated wave phase. Perhaps the most prominent of these is the appearance of strong secondary waves which travel perpendicularly to primary wave propagation direction. The waves take the form of regions of density enhancements and depletions propagating along the wavefronts of the primary wave. The simulation movies show quite clearly that these disturbances travel upward along the primary wave density peaks and downward in the valleys.
We also observe modifications in the primary waves themselves. Portions of the primary wavefronts tend to turn their propagation direction away from horizontal (i.e., along the direction of primary electron drift), progagating instead either somewhat upward or downward. We will refer to this phenomenon as "wave turning." Further, this tendency is not symmetric; the waves more often redirect their propagation upwards than downwards. This asymmetry may account for radar observations, which see more power when looking eastward than westward in the equatorial electrojet.
The waves also propagate more slowly once saturation occurs. The phase velocity lies significantly below the phase velocity predicted by linear theory but above the sound speed.
The conflict among all these effects tends to tear the primary wavefronts apart, thereby effecting an end to linear growth. Torn segments of the primary wavefronts are also observed to reconnect with neighboring wavefronts. We are at present examining the interaction among the various effects in an attempt to provide a detailed explanation of how saturation occurs.
While no clear explanation yet exists as to what causes saturation, we can say with certainty that, within the limitations of our simulation, saturation must be occurring through the action of the nonlinear electron Hall drift. When the corresponding term in the potential equation, kappa*grad(n) x grad(phi), is removed from the simulation (along with the much smaller term grad(n).grad(phi)), saturation does not take place until the wave grows to amplitudes ten times that seen when the terms are included. Further, the behavior once saturation occurs is completely different. Neither turning of the propagation direction nor breakup of the primary wavefronts occur. The waves maintain an essentially monochromatic profile with no change in phase speed as the waves saturate. Examination of the ion phase space suggests that saturation is now occurring through the substantial trapping of ions by the wave.
Other properties of the wave saturation process were uncovered by varying one or more of the parameters of the above simulation. When the driving electric field was lowered to its typical value in the topside equatorial electrojet (15 mV/m, corresponding to a horizontal drift of 600 m/s), a very large number of ion simulation particles (3.6 million) running on the Cray machines at LANL was necessary to lower particle noise sufficiently to see clear wave saturation. The saturation rms density variation was then found to be 3% to 4%. This case was also characterized by a more modest tendency of the wave to turn away from its primary propagation direction.
We were also able to able to check the wave power as a function of elevation angle for these realistic physical parameters. At a wavelength of 3 meters, the power falls at the rate of 0.25 db per degree of elevation angle. This compares favorably with 3-meter radar observations, where a fall-off rate of 0.3 db / degree was measured [Ierkic et al, 1980].
When the driving electric field was increased further, the simulations showed a much faster linear growth rate followed by a more dynamic saturated state. The secondary waves were so strong that distinguishing them from the primary wave became difficult. The disparity between the primary wave phase velocity and linearly-predicted phase velocity also became unmistakable: with the difference between the linear phase velocity and the sound speed large in this high drift case, it was clear that the observed waves were propagating well below the the linear phase velocity but above the sound speed. These strongly driven cases were also characterized by a significant coupling of energy to long-wavelength modes. These modes were observed to travel with phase velocities significantly below the sound speed.
The effect of varying the electron temperature was observed to be similar to that of varying the driving electric field. The reason is simple: decreasing the electron temperature decreases the sound speed. Since the linear growth rate of the instability increases as the difference between the squares of electron drift speed and the sound speed, it follows that a lower electron temperature should lead to a more vigorously unstable state. This is exactly what is observed. Lower electron temperatures produce more wave turning and more turbulence in the saturated state, just as occurs in the higher electron drift case. One difference, of course, is that the waves travel more slowly, since both the linear phase velocity and sound speed are smaller when the electron tempeature is low. We do not expect that our simulations are telling the complete story in these vigorously unstable cases. In the actual plasma, these instabilities probably lead to significant electron heating, modifying the sound speed. This may be the situation, for example, in the auroral electrojet.
A two-fluid code was also developed for this project; however, its results are currently unreliable except for purposes of comparison to the particle code. The two-fluid code should produce the same results as the particle code, as the results of the particle code suggest that the only significant ion kinetic effect present is occurring at short wavelengths, where ion Landau damping is suppressing the instability. This effect was mimicked in the fluid code by imposing artificial damping at those wavelengths. Despite the inclusion of this damping, however, the two-fluid code developed numerical instabilities at the scale of the grid spacing, rendering the code unreliable. Before this instability developed, though, we were able to see many of the same features displayed by the particle code, including wave turning and the development of secondary waves.
M. Oppenheim and N. Otani, Spectral characteristics of the Farley-Buneman instability: simulations vs. observations, J. Geophys. Res., in review (1996).
M. Oppenheim, Nonlinear simulations and theory of the Farley-Buneman instability in the equatorial E-region, J. Atmos. Terr. Phys., in review (1996).
M. Oppenheim, Nonlinear simulations and theory of the Farley-Buneman, Instability in the E-region ionosphere, Ph. D. Thesis, Cornell University, (1995).
M. Oppenheim, N. Otani and C. Ronchi, Nonlinear simulations and theory of the Farley-Buneman instability in the equatorial E-region ionosphere, Proc. of the Ninth Inter. Symp. on Equatorial Aeronomy, A3-1, (1995).
M. Oppenheim, N. Otani and C. Ronchi, Hybrid simulations of the saturated Farley-Buneman instability in the ionosphere, Geophys. Res. Lett. 22, 353 (1995).
Farley, D. T. A plasma instability resulting in field-aligned irregularities in the ionsophere, J. Geophys. Res. 68, 6083 (1963).
Freund, R. W. and N. M. Nachtigal, An implementation of the QMR method based on coupled two-term recurrences, SIAM Journal on Scientific Computing 15, (1994).
Ierkic, H. M., B. G. Fejer, and D. T. Farley, The dependence on zenith angle of the strength of 3-meter equatorial electrojet irregularities, Geophys. Res. Lett. 7, 497 (1980).
Janhunen, P., Perpendicular particle simulation of the E region Farley-Buneman instability, J. Geophys. Res. 99, 11461 (1994).
Newman, A. L., and E. Ott, Nonlinear simulations of type I rregularities in the equatorial electrojet, J. Geophys. Res. 86, 6879 (1981).