Modeling of the Farley-Buneman instability in the E-region ionosphere: a new hybrid approach

A novel approach to nonlinear simulations of the Farley-Buneman (FB) instability in the E-region ionosphere is developed. The mathematical model includes a fluid description of electrons and a simplified kinetic description of ions based on a kinetic equation with the Bhatnagar-GrossCrook (BGK) collision term. This hybrid model takes into account all major factors crucial for development and nonlinear stabilization of the instability (collisional drag forces, ion inertia and Landau damping, dominant electron nonlinearity, etc.). At the same time, these simulations are free of noises caused by the finite number of particles and may require less computer resources than particle-in-cell (PIC) or hybrid – semi-fluid semi-PIC – simulations. First results of 2-D simulations are presented which agree reasonably well with those of previous 2-D PIC simulations. One of the potentially useful applications of the novel computational approach is modeling of the FB instability not far from its threshold.


Introduction
The Farley-Buneman (FB) instability is a low-frequency plasma instability driven by a sufficiently strong quasistationary electric field E 0 perpendicular to the geomagnetic field B 0 . This instability occurs in the weakly ionized E-region ionosphere where electrons are magnetized, while ions are unmagnetized due to frequent collisions with neutral particles. As a result, the velocity distribution of electrons is shifted relative to that of ions by the drift velocity V ≈cE 0 ×B 0 /B 2 0 (that is why the FB instability is also named the modified two-stream instability). In the equato-Correspondence to: Y. S. Dimant (dimant@bu.edu) rial and high-latitude electrojets, if the ambient electric field E 0 ≡|E 0 | exceeds a threshold value, E Thr (10 -20) mV/m, the average flow motion of electrons with respect to ions becomes supersonic. While average motion of frequently colliding ions is predominantly diffusive, their relatively small inertia is of importance. It results in excitation of field-aligned density perturbations, k k ⊥ , coupled to lowfrequency electrostatic fluctuations (here k and k ⊥ are the components of characteristic turbulence wavevectors k parallel and perpendicular to B 0 , respectively). Average amplitudes of density fluctuations usually do not exceed several percent (e.g. Pfaff et al., 1987Pfaff et al., , 1997. Typical wavelengths of density fluctuations, detected mostly through narrowband type 1 radar echoes (Cohen and Bowles, 1967;Balsley and Farley, 1971;Crochet et al., 1979;Kudeki et al., 1987;Ravindran and Reddy, 1993), are usually in a meterscale range (for review on observations at both equatorial and auroral latitudes see Fejer and Kelley, 1980;Kelley, 1989;Haldoupis, 1989;Sahr and Fejer, 1996). In addition to the equatorial and high-latitude electrojets, the FB instability can also develop at midlatitude sporadic-E layers where local electrostatic fields can exceed the FB threshold field (Schlegel and Haldoupis, 1994;Haldoupis et al., 1996;Haldoupis et al., 1997).
Starting from the two pioneering papers by Farley (1963) and Buneman (1963), the linear theory of the FB instability has been developed for many years (Lee et al., 1971;Schmidt and Gary, 1973;Ossakow et al., 1975;Fejer et al., 1984;Dimant and Sudan, 1995a,b,c;Kissack et al., 1995;Kissack et al., 1997;St.-Maurice and Kissack, 2000;Drexler et al., 2002;Dimant and Oppenheim, 2004;Kagan and St.-Maurice, 2004;Drexler and St. Maurice, 2005;Kissack et al., 2008a,b). This theory, however, has limited applications. It may give the threshold conditions for the instability, describe the initial stage immediately after the instability onset, and show some tendencies, such as field-aligned nature of irregularities. However, the linear theory per se Published by Copernicus Publications on behalf of the European Geosciences Union. 2854 D. V. Kovalev et al.: Modeling of Farley-Buneman instability in E region cannot describe the process of nonlinear saturation of the instability and provide amplitude and spectral characteristics of developed turbulence. To explain observations and make quantitative predictions, one needs a nonlinear theory. There are many important questions that must be answered by such a theory. Among those: What are the major factors that lead to nonlinear saturation of the instability? Depending on the ambient electric field E 0 and ionospheric parameters, what is the average level of density and electric field fluctuations, what are the wavelengths, phase velocities, and other spectral characteristics of the most pronounced waves in the bulk of turbulence?
Although nonlinear theory of the FB instability has been also under development for quite a long time (Skadron and Weinstock, 1969;Weinstock and Sleeper, 1972;Sudan et al., 1973;Lee et al., 1974;Rogister and Jamin, 1975;Sudan, 1983a,b;Robinson, 1986;St.-Maurice, 1993, 1995;Smolyakov et al., 2001;Bahcivan and Hysell, 2006), it is still far from completion. A mode-coupling anomaloustransport theory has been proposed by Sudan (1983a,b, for earlier nonlinear theories see references therein). Sudan's theory, following the resonance broadening approach by Dupree (1968), results in a concept of anomalous collision frequency (for review see Hamza and St.-Maurice, 1993;Robinson, 1994;Hamza and St.-Maurice, 1995). Note that some nonlinear theories employ interplay between the twostream and gradient-drift instabilities (Sudan, 1983b;Sudan et al., 1973;Keskinen, 1981). Being capable to explain some observed features of FB turbulence, resonance broadening theories of the FB instability, however, encounter certain conceptual problems and are not quite satisfactory from the theoretical viewpoint (for discussion see, e.g. Hamza and St.-Maurice, 1995).
The FB instability has been also extensively studied numerically. Nonlinear simulations of the FB instability are based on the two-fluid (Newman and Ott, 1981), particlein-cell (PIC) (Machida and Goertz, 1988;Schlegel and Thiemann, 1994;Janhunen, 1994;Oppenheim and Dimant, 2004;Oppenheim et al., 2008), or hybrid Dyrud et al., 2006) -semifluid, semi-PIC -codes. Note that for the FB instability, the kinetic effect of ion Landau damping is crucial because it gives the short-wavelength restriction on linearly unstable spectral domain. (According to the conventional two-fluid theory without ion Landau damping, the growth rate of linearly unstable Fourier harmonics would increase infinitely with the wavenumber, γ k ∝k 2 .) Therefore, two-fluid equations must include an artificial viscosity to model important kinetic effects, which is in the general case not satisfactory.
Simulations using PIC codes represent a fully kinetic treatment which automatically includes all physics of the plasma processes. However, PIC simulations have restrictions associated with the finite number of particles and the corresponding noises. In order to decrease the computational time, one usually has to aggregate the actual particles into bigger clus-ters. Besides, in order to avoid very short time steps, one often has to reduce the electron gyro-frequency and plasma frequency via the artificially increased electron mass. Even properly rescaled, such computational results not fully model the actual physical situation. Furthermore, numerical studies of wave activity have been mostly limited to 2-D cases: either in the plane parallel to the magnetic field (Machida and Goertz, 1988;Schlegel and Thiemann, 1994) or in the perpendicular plane (Newman and Ott, 1981;Janhunen, 1994;Oppenheim and Dimant, 2004;Dyrud et al., 2006;Oppenheim et al., 2008). The latter simulations are more relevant to the actual situation because they correctly take into account the dominant electron nonlinearity associated with the fluid-model term ∝δE×∇δn e , where δE and δn e are the turbulent electric field and electron density perturbations, respectively. These simulations reveal a new insight into the problem of nonlinear saturation of the FB instability. In particular, they (as well as rocket observations) show that turbulence, in addition to chaotic features, often has a pronounced ordered structure with dynamical behavior. Such 2-D simulations, however, cannot consistently model the effects of finite k which are responsible for high-latitude anomalous electron heating during magnetospheric storms or substorms (Schlegel and St.-Maurice, 1981;St.-Maurice and Laher, 1985;Providakes et al., 1988;St.-Maurice et al., 1990;Dimant and Milikh, 2003;Milikh and Dimant, 2003). While 3-D PIC simulations of the FB instability are currently under way, modeling the nonlinear stage by simpler continuous equations that sufficiently correctly take into account principal kinetic effects would be beneficial for both analytical and numerical treatment.
Our major goal is to reach a further progress in modeling nonlinear saturation of the FB instability. Fully-kinetic PIC simulations are one of the most powerful modeling tools, but they require significant computer resources and still have drawbacks described above. In this study, we are testing an alternative computational approach which by its physical implications is similar to the hybrid code by ; . In this, also hybrid, approach we consider a fluid model for electrons coupled to a continuous kinetic model for ions. The latter includes ion Landau damping and some important ion thermal factors, but does not involve discrete particles. The reason for choosing this hybrid approach is as follows. While both electrons and ions are prone to Landau damping resulting in an efficient suppression of the instability, electron Landau damping is only effective at a short-wavelength, high-frequency range where the wave growth rate is already effectively suppressed by ion Landau damping. This allows us to use for electrons the fluid model. In this paper we explore isothermal electrons, so that no electron thermal effects, as well as other kinetic effects of electrons (Dimant and Sudan, 1995a,b,c;Dimant and Sudan, 1997;Kagan and St.-Maurice, 2004;Kissack et al., 2008b,a), have not yet been included. As we explained above, ion Landau damping is of crucial importance for the FB instability, so that ions have to be treated kinetically. In our new treatment, we employ modeling technique based on continuous equations rather than on discreteparticle codes (Filbet and Sonnendrücker, 2003), so that we make use of an ion kinetic equation. Unlike electrons, heavy ions have comparable rates of collisional losses of energy and momentum, and a kinetic equation with the simple BGK collision term (Bhatnagar et al., 1954) seems to be a reasonable approximation for them. Furthermore, this equation includes ion-thermal effects which may affect significantly FB turbulence . The hybrid set of electron-fluid and ion-kinetic equations takes into account the major factors important for nonlinear saturation of the FB instability: the dominant electron nonlinearity ∝δE×∇δn e and ion Landau damping. In this paper, we restrict ourselves to the 2-D space perpendicular to B 0 , so that k =0.
The remainder of the present paper is organized as follows. Section 2 examines the physical conditions applying in the FB instability and presents the new hybrid method used in the simulations. Section 3 describes the numerical solver, whereas Sect. 4 presents the results of the performed simulations. Section 5 discusses some implications of our results, while Sect. 6 summarizes the main findings of this work, followed by an Appendix where some basics of the two-fluid linear theory are given for the purpose of completeness.

Basic conditions
The FB instability can be efficiently excited at upper D/lower E-region altitudes roughly between 80 and 120 km where electrons are strongly magnetized, while ions are essentially unmagnetized, Here ν en and ν in are the average frequencies of electronneutral and ion-neutral collisions; e,i =eB 0 /m e,i are the gyrofrequencies of electrons and ions of masses m e,i , respectively (e is the elementary charge; B 0 ≡|B 0 |; m i ≈30 amu). The collision frequencies, which are proportional to the neutral density, decrease exponentially with increasing altitude. However, throughout the entire E-region ionosphere their ratio remains nearly constant, ν en 10ν in . The most of developed FB turbulence is characterized by sufficiently low-frequency and long-wavelength waves, ω∼kv T i ν in ν en , where ω and k≈k ⊥ are characteristic wave frequency and wave number, v T i =(T i /m i ) 1/2 is the characteristic ion thermal speed, and T e,i are the temperatures (in energy units) of electrons and ions, respectively. Low-frequency and low-current plasma processes in the E-region ionosphere result in insignificant magnetic field variations. This means that these processes have an electrostatic nature and the turbulent electric field can be adequately described by an electrostatic potential , E=−∇ .

Electron fluid model
Under condition ω ν en , one can neglect both inertia and Landau damping of electrons. Assuming relatively low Eregion altitudes, we will also neglect kinetic corrections that describe electron thermal effects (Dimant and Sudan, 1995a;Kagan and St.-Maurice, 2004). For isothermal electrons, the standard continuity and momentum equations result in the diffusion-convection equation where the electron flux components are given by e ≡ n e V e = − T e m e ν en ∇ n e − en e T e ∇ δ , V e ,⊥ are the components of the electron fluid velocities, δ is the fluctuating electrostatic potential, E 0 is the ambient electric field which is practically perpendicular to the geomagnetic field B 0 , V 0 =cE 0 ×b/B 0 is the E 0 ×B drift velocity,b is the unit vector along B 0 , and the 'nabla' operators ∇ ,⊥ pertain to the coordinates parallel and perpendicular to B 0 respectively. Now we consider a purely 2-D case when all spatial variations are perpendicular to the magnetic field. Bearing in mind the characteristic time and length scales for developed FB turbulence, τ =ν −1 in and l x,y =v T i /ν in , where v T i =(T i /m i ) 1/2 , we will normalize the time and coordinates to these quantities, tν in →t; x/ l x,y , y/ l x,y →x, y, where x and y are coordinates along the V 0 and E 0 , respectively. We will also introduce a dimensionless potential, φ=e /T i , and parameters  (Farley, 1985;Dimant and Milikh, 2003). The parameter ψ ⊥ exponentially decreases as altitude increases, reaching unity at altitudes 94-97 km (e.g. Dimant and Oppenheim, 2004, Fig. 2). At the same time, the small parameter 0 remains nearly constant, 0 1.35×10 −2 . In the renormalized variables, Eqs.
(2) to (4) reduce to 1 ψ ⊥ ∂n e ∂t = ∂ 2 n e ∂x 2 + ∂ 2 n e ∂y 2 − n e ∂ 2 φ ∂x 2 + ∂ 2 φ ∂y 2 + ∂n e ∂y 2.3 Ion kinetic model As stated above, ion Landau damping is crucial for the FB instability, so that ions must be treated kinetically. Assuming for simplicity altitudes below 110 km where ions are practically fully unmagnetized, we will use a kinetic equation with the BGK collision term (Bhatnagar et al., 1954;Gross and Krook, 1956;Morse, 1964): Here v is the ion velocity, f (r, v, t) is the ion distribution function normalized according as where n i ≡n i (r, t) is the local ion density; is the ion Maxwellian function normalized to the actual density, and n 0 is the undisturbed plasma density. In the present simulations, for simplicity, we set T i =T e , where T i is the average ion temperature, although in the general case the two undisturbed temperatures may differ. Note that Eq. (8) allows for wavelike temperature perturbations associated with density fluctuations and variations of local ohmic heating. It should be noted that the BGK collisional term is a model approximation which does not follow from the rigorous Boltzmann collision operator, so its use is restricted. The BGK model is usually a reasonable approximation for lowcollisional, high-frequency regimes where the typical wave frequencies are much larger than the collision rates. However, for highly collisional, low-frequency E-region instabilities, the situation is more complicated. Indeed, in framework of the BGK collision term, collisions of charged particles with neutrals have the same characteristic rate for both momentum and energy changes. Besides, this collisional rate is not allowed to have a dependence on the particle velocity. For low-energy (E 1 eV) light electrons colliding with heavy neutrals, the characteristic rates of energy losses (mainly inelastic) and momentum changes (mainly elastic) differ by two-three orders of magnitude, and these rates depend strongly on the electron velocity (Gurevich, 1978;Schunk and Nagy, 2000). That is why the electron BGK model for low-frequency E-region processes is a rather poor approximation (for discussion, see Dimant and Sudan, 1995a,c). At the same time, ions in the E-region ionosphere collide with neutrals of about the same mass. As a result, the collisional changes of ion energy and momentum are determined roughly by the same characteristic time τ =ν −1 in which is essentially independent of the ion velocity (Schunk and Nagy, 2000). That is why the ion kinetic equation with the BGK collision rate can be considered as a reasonable approximation. Comparison of our numerical results with those for PIC ions supports this assertion.
Regardless of the spatial dimension of the problem, the ion distribution function f (v) has all three velocity components, v x,y,z . In the 2-D case, Eq. (8) involves no v z -derivatives. Integrating linear Eq. (8) over v z and passing to the dimensionless coordinates and time, we obtain where and, according to Eq. (9), we have

Electrostatic potential
For the E-region plasma processes, inequalities ω ω pi and kλ D 1 usually hold, where ω and k are the typical perturbation wave frequency and wave number, whereas ω pi (n i )=(n i e 2 / 0 m i ) 1/2 and λ D =v T i /ω pi are the ion plasma frequency and Debye length, respectively. This usually results in quasi-neutrality, n i ≈n e =n. However, we use here the more general treatment without considering the electron and ion densities equal but solving directly for the electrostatic potential via Poisson's equation, or, in terms of the above dimensionless variables,
In order to drive the FB instability, the ambient electric field should exceed at least the minimum FB threshold field obtained from the two-fluid linear theory, see Eq. (A4) in Appendix (the kinetic corrections can only increase the threshold field). Compared to the conventional two-fluid linear relations obtained under strictly quasi-neutral conditions, Eqs. (A3) and (A4) contain an additional term ∝ ν 2 in /ω 2 pi . This term arises due to small charge separation (Rosenberg and Chow, 1998); its physical nature is discussed in Appendix A. In the daytime E-region ionosphere, as well as at nighttime at altitudes above 100 km, we have ω pi ν in , so that in real ionosphere the effect of charge separation usually plays no important role. However, it should be taken into account in simulations . In our simulations, we have modeled predominantly the daytime ionosphere with n i =10 10 -10 11 m −3 (only one simulation with n i =10 9 m −3 ) and assumed ψ ⊥ <0.2 corresponding to altitudes above 98-100 km (Dimant and Oppenheim, 2004, Fig. 2). Under these conditions, according to Eq. (A5), the effect of the additional term is small.
As the initial condition, we have set a uniform fluctuation noises for electrons and ions with the electric field due to charge separation well below the FB threshold field. We have used periodic boundary conditions in both coordinates x and y. This means that the values of the particle densities and velocities are equal on the opposite sides of the simulation box. To solve these equations, we have used finite difference methods. In the 2-D computational box with the sizes L x,y , we have used a mesh with homogeneous grid sizes.
The electron equation is a nonlinear convection-diffusion equation. Diffusion part is solved with the help of a secondorder accuracy numerical scheme. Algorithm for the convection part is based on interpolation over characteristics. In our case we use Lagrange interpolation over 5 points in each direction.
Kinetic Eq. (8) has a larger dimension than the two previously described equations because in addition to the spatial variables it involves also a 2-D velocity space. We approximate the latter by a uniform grid which covers a finite domain restricted in each dimension by a maximum speed V max , −V max <v x,y <V max . We have found empirically that to adequately model the FB instability, we should choose V max ≥6v T i . We solve kinetic Eq. (10) by splitting the entire kinetic equation into two sub-equations, At each time step we solve for convection in the coordinate space with initial conditions from the previous time step and then at the same time step we solve for convection in the velocity space with the solution of convection in the coordinate space as an initial condition. A similar technique was described in Sonnendrücker et al. (1999) and Filbet and Sonnendrücker (2003) for Vlasov's equation, which differs from our ion kinetic equation only by the right-hand side. The fact that the two sub-equations form a set of 2-D convection equations allows us to use fast and efficient methods for obtaining of solution of each equation. In the simulator we use interpolation over characteristics essentially in the same way as for electron convection. We solve the full set of Eqs. (7), (10), and (13) using two different times steps, h, and h mod ≡h/N, where N is a large integer number, N=10-30. At each long time step h, we solve ion kinetic Eq. (10). At each short time step h mod , we solve alternatively electron-density Eq. (7) and Poisson's Eq. (13). We use two different time steps because solving 4-D ion kinetic Eq. (10) requires a much longer computer time than solving 2-D Eqs. (7) and (13) (more than 90% of the total computer time). That is why we need to solve Eq. (10) as seldom as possible, in other words, with the maximum time step allowed by the Courant-Friedrichs-Lewy (CFL) condition. The CFL condition for electrons imposes a much harder restriction than that for ions due to a much larger electron mobility. Trying to solve Eq. (10) as seldom as possible, we will use the maximum time step. Because redistribution of electrons induces changes in the electric field, at each modified time step we also need to solve Eq. (13).

Results of simulations
We have modeled the evolution of the 2-D FB instability using various sets of parameters. Emphasize that we have performed our simulations using the real electron mass, actual ionospheric parameters, and realistic values of the perpendicular DC electric field. Here we show the results of several runs characterized by different sets of parameters shown in the corresponding tables. In each run we have reached the nonlinear saturation of the FB instability. We verified this fact by monitoring the temporal behavior of the root-meansquare (rms) values of the turbulent electric fields and density irregularities.

Run 1
In this run, the DC electric field was more than twice the FB instability threshold field, Eq. (A4). Table 1 shows major Figures Fig. 1. Electron density for 3000, 6000, 9000, and 20000 time steps for Run 1 (see Table 1).  Table 1).
parameters. Figure 1 shows the evolution of the electron density during the development of the FB instability. In this and similar 2-D diagrams, the vertical axis is along the direction of E 0 , while the horizontal axis is along the E 0 ×B 0 -drift direction.
A heuristically expected characteristic evolution time of the FB instability, τ =ν −1 in ≈2.66×10 −4 s, in this run equals roughly 500 time steps. Simulations (especially run 5) show, however, that the time of nonlinear saturation of the FB instability proves to be up to two orders of magnitude longer than τ . A characteristic wavelength, determined roughly by λ char =2π/k char =2πC s ν −1 in with the ion acoustic speed C s 400 m/s, is λ char 0.7 m, i.e., about thirty times less than the box size. We note, however, that these characteristic scales are only for crude estimates, the actual preferred turbulent wavelengths may differ by a factor of order unity or more. We should also bear in mind that the limited sizes of the simulation box, along with the periodic boundary conditions, impose discrete-spectrum restrictions on the wavelengths: the allowed wavelengths in each direction can only equal the corresponding box size divided by an integer number, λ x,y ≡2π/|k x,y |=L x,y /N x,y , where N x,y =0, 1, 2, 3.... This limits possible flow angles for tilted waves in turbulence, θ≡ arctan(k y /k x ), to only discrete values, θ=± arctan(N y /N x ), e.g. θ=0, ± arctan(1/3), ± arctan(7/26) and the like.
The linear (exponential) growth of the instability occurs at a rather short time ∼6τ . Figure 1 shows after that during several τ the instability continues growing slowly, although not in a linear way. During this time, typical density fluctuations are less than 5%. They represent quasi-monochromatic waves with the preferred wavevector oriented practically parallel to the E 0 ×B 0 -drift direction. The preferred wavelength, λ=L x /6, where L x is the box size in x and y, is about 1.5 times λ char .
In the course of instability development, density fluctuations become deeper, and their structure becomes much more turbulent. This is a direct manifestation of nonlinear  Table 2).  Table 2). mode coupling. Although in this stage there is no clearly seen quasi-monochromatic waves, one can see that the typical wavelengths of the most pronounced waves become longer. The right panels of Fig. 1 show that the characteristic wavevector starts deviating from the E 0 ×B 0 -drift direction. The effect of wave-front tilting has been observed in previous simulations with PIC ions (e.g. Janhunen, 1994;Oppenheim and Dimant, 2004;Oppenheim et al., 2008). Dimant and Oppenheim (2004) attributed this effect to the additional instability driving mechanism of the ion-thermal nature. Due to the periodic boundary conditions and finite box sizes, the allowed discrete values of the tilt angle, θ tilt , are constrained by θ tilt = arctan(N y /N x ) with integer N x,y , as discussed above. The most pronounced values of these integers vary from N x =6 and N y =0 in the leftmost panel to N x =2, N y =1 in the rightmost panel. This corresponds to a discretestep transition from θ tilt =0 • to θ tilt ≈27 • . To the 6000 time step, the amplitude of density fluctuations increased to (15-25)% and then decreased slightly to a saturated level of (10-20)%. Note that this level of fluctuations is higher than it would be expected heuristically (Dimant and Milikh, 2003), δn e /n 0 ∼ /ν in 4%. A similar effect has been observed in previous 2-D fully-PIC or hybrid simulations Oppenheim and Dimant, 2004). A possible explanation is that the constraint to the two dimensions does not allow the turbulent energy to leak out in the parallel to B 0 direction where the instability is highly damped. Another important factor is that the turbulent energy cannot also spread over long wavelengths due to the finite box sizes. As a result, the instability is saturated at a higher level.
A run with the same basic parameters but twice as large box sizes shows that the major qualitative and quantitative characteristics of developed turbulence remain roughly the same. At the same time, there is a tendency to the development of longer wavelengths in the saturated stage, up to the simulation box sizes. The simulation box in this run was not big enough to allow longer wavelengths to develop. In the following run, we use increased box sizes.

Run 2
In this run, the box size in each direction is four times as large as the box size in Run 1. The major parameters for this run are listed in Table 2. Notice that in this run the time step is twice as short as that in Run 1. We also reduced in this run the plasma density by an order of magnitude for numerical purposes. Figure 2 shows the evolution of the density turbulence. We see that the development of instability in this run is quite similar to that in Run 1. For a sufficiently long time, as seen in the first three panels, a dominant wave is the quasi-monochromatic wave with the wavelength λ=L x /10≈0.28 m, which is directed along the E 0 ×B 0 drift. After the 50 000 time step, t≈0.013 s≈50 ν −1 in , the character of density fluctuations changes dramatically. These fluctuations become turbulent with no clearly seen preferred wave. The effect of the angular offset from the E 0 ×B 0 -drift direction, however, is still clearly seen. The wave tilting angle is about the same as in Run 1, with the similar, or slightly smaller, fluctuation level (about 20%). Figures 3 and 4 show the corresponding spectrum of the turbulent electric field which is always coupled to density fluctuations and behavior of the root-mean-square values of the turbulent electric field.
To find dominant values of the phase velocity a total spectral energy graph was constructed by analogy with Fig. 7 from Oppenheim et al. (2008). The idea is as follows. The dominant phase velocities correspond to the points of maxima in the Fourier spectral energy plot of the electron or ion density functions. Due to small charge separation and quasineutrality conditions, both densities are close to each other so that one can use any of them. After 3-D Fourier transformation of the density in space and time one obtains a function n(k x , k y , ω) which can be interpreted in terms of n(|k|, φ, ω) after transition to spherical coordinates. Then the values of n(|k|, φ, ω) can be averaged over φ and as a result two dimensional function n(|k|, ω) arises. The last step is similar to the averaging of radar data because radars measure waves with fixed wavelength but from different directions. The phase velocity by definition is v ph =ω/|k| so each value of n(|k|, ω) corresponds to the determined v ph and the dominant phase velocities are in line with the points of maxima of n(|k|, ω). The resulting plot of the total spectral energy for run 2 is presented in Fig. 5. The upper bound of the dominant phase velocities is approximately v max ph 700 m/s. It was computed as the maximum absolute value for the phase velocities corresponding to the 70 maximal values of spectrum of n(|k|, ω) (the difference in spectral values at this points is up to two orders of magnitude). The computed quantity is higher than the ion acoustic speed C s =[(T e +T i )/m i ] 1/2 406 m/s for isothermal electrons and ions, but below the velocity predicted by linear theory, V 0 /(1+ψ) 908 m/s.

Run 3
This run was intended for comparison of simulations with real electron mass and increased electron mass. For this purpose we chose parameters the same as for fully PIC simulations performed by Oppenheim and Dimant (2004), but left the plasma density and the driving field as in Run 2. The electron mass in Run 3 was artificially increased by a factor of 44 which affects the parameter 0 , Eq. (6), and grid spacings, see Table 3. The main advantage of using the increased electron mass in our simulations is that this allows modeling the instability in a larger simulation box keeping the same number of points. However, this requires changing the physical parameters.   Despite the changes in parameters, Fig. 6 shows that the instability evolution stages are qualitatively the same as in the previous runs. The initial quasi-monochromatic waves with the wavevectors parallel to E 0 ×B 0 gradually transform to turbulent waves with the tilted wave fronts. In this simulation, the preferred wavelength of the quasi-monochromatic waves is λ=L x /18≈0.44 m, when the characteristic wavelength is λ char ≈1.4 m. This wavelength is twice larger than in the run 2. The effect of wave-front tilting exists and can be seen in the density plots or spectrum graphs (Fig. 7). The density fluctuations are fairly high, ∼30−40 %. It should be noted that simulations with real electron mass give us density fluctuations lower than 25%. Figure 8 shows that after instability saturation the rms turbulent electric field is about twice the driving electric field E 0 , what is a little bit lower than in the run 2 and can be attributed to the larger box size. This figure also shows that the saturated turbulent field is much more stable than that in the simulations with the real electron (see, for example, Fig. 4).

Run 4
This run was also intended to directly compare our simulation with full PIC simulation from Oppenheim and Dimant (2004). Table 4 shows the main simulation parameters that coincide with those of Oppenheim and Dimant (2004). The difference with Run 3 is in the plasma density, driving field and box size. Note that we were able to perform simulations with the twice as large time step compared to that of Oppenheim and . Density and spectrum graphs (Figs. 9 and 10) are in agreement with the simulation by Oppenheim and Dimant (2004). It is possible to see the same wave tilting and formation of large waves after saturation with the angle about 30 • . The rms of the turbulent electric filed (Fig. 11) is in rough agreement with the results by Oppenheim and Dimant (2004). The average values are approximately the same, but there is a difference in the growth time of instability. Our growth time is about 12 ms while in Oppenheim and Dimant (2004) it was about 35-40 ms. This difference can be attributed to the electron thermal effects absent in our simulations but included in the PIC code and also to the different initial conditions.
Thus we see that our simulations without electron thermal effects give results reasonably close to those of PIC simulations. Preliminary results of simulations with the electron thermal effects included show no significant difference but show somewhat smaller rms turbulent electric field, density fluctuations, and noticeably longer time of the instability growth and saturation. Figure 12 shows contours of the ion distribution function averaged over the entire box during the instability saturation. This figure shows significant anisotropic and non-Maxwellian modifications of the ion distribution function due to the predominantly Pedersen ion response to the total electric field. The reason for the anisotropic response is that ions whose mass nearly equals that for the colliding neutral particles have comparable rates of ion-neutral collisional changes for both momentum and energy. As a result, unlike light electrons, for heavy ions collisional angular scattering in the velocity space does not lead to an effective isotropization of the ion distribution function. Furthermore, we see that the largest distortions of the ion distribution function take place at suprathermal energies. This is partially due to the kinetic effect of Landau damping which causes waves to yield their energy to resonantly interacting ions. We should bear in mind, however, that highly collisional waves excited by the low-frequency   Table 3). 31 Fig. 6. Electron density (from left to right, top to bottom) for 12 000, 26 000, 40 000, 100 000 time steps for Run 3 (see Table 3). range where the wavelengths are comparable to the ionneutral mean free path . This makes ion Landau damping for such waves less pronounced than that for plasma waves in the high-frequency, weakly collisional regime.

Run 5
Basic parameters of this run are the same as for Run 2, except for the difference in the grid sizes and, more importantly, in the value of the driving electric field, E 0 . As can be seen from Table 5, the value of E 0 is the half that in Run 2, thus being closer to the threshold value. Note that PIC simulations with the driving electric field close to the instability threshold are hardly achievable due to numerical noise comparable with density fluctuations.
According to Eq. (A4) in Appendix, for ψ ⊥ =0.1 and other parameters listed in our Table 5, the FB instability threshold field for the high latitude electrojet, B 0 5×10 4 nT, is  Fig. 7. Spectra of E 2 corresponding to Fig. 6 for 12000, 26000, 40000, 100000 time steps (from left to right, top to bottom). The maximum values of spectral modes are shown for each time step. Fig. 8. Rms turbulent electric field for Run 3. Numbers form 1 to 4 correspond to the density plots at Figure 6 for 12000, 26000, 40000, 100000 time steps, respectively. Green line corresponds to the E0 value.    8. Rms turbulent electric field for Run 3. Numbers form 1 to 4 correspond to the density plots at Figure 6 for 12000, 26000, 40000, 100000 time steps, respectively. Green line corresponds to the E0 value. E min Thr 0.022 V/m. In this run (E 0 =0.029V/m), we have (E 0 −E min Thr )/E 0 0.31, so that the excess above the minimum threshold field is not significant. In principle, this hybrid model with the fully continuous equations allows setting the driving field as close to the threshold field as possible. However, the closer the driving field is to the instability threshold, the smaller is the instability growth rate and the longer time is required to reach the instability saturation. It is important, nevertheless, that this approach allows to explore the near-threshold case without "plunging" into inherent noises caused by the finite number of PIC particles in the corresponding codes. We are planning to carefully study the FB instability dynamics near the instability threshold after we incorporate into the electron fluid module the additional energy balance equation that includes temperature perturbations (Dimant and Sudan, 1995a;Dimant and Sudan, 1997;Kagan and St.-Maurice, 2004;Kissack et al., 2008a).
Instability development stages in this run are qualitatively the same as those in the previous simulations. Initially, Fig. 9. Electron density (from left to right, top to bottom) for 2000 and 40600 time steps for Run 4 (see Table 4). The maximum values of spectral modes are shown are shown for each time step.  Table 4).  Table 4). The maximum values of spectral modes are shown are shown for each time step. Fig. 11. Rms turbulent electric field for Run 4. Numbers 1,2 correspond to the density plots at Figure 9 for quasi-monochromatic waves develop. Their length is twice as that in Run 2 and is half the heuristically expected value. During the saturation stage, the wavelengths increase, but due to restrictions on the box size their growth is limited (Fig. 13).
There are some differences in the process of instability growth. First, the total instability growth time decreases due to proximity to the instability threshold. This time is approximately 0.15 s (Fig. 14) as compared to 0.03 s in Run 2 (Fig. 4). Second, the rms turbulent electric field is on average lower than the driving electric field after saturation, while in the simulations with large enough driving field the rms of turbulent electric field was noticeably larger than the driving field. This field, however, remains larger than that expected heuristically (Dimant and Milikh, 2003). Third, the density fluctuations in the simulation are less than 10% and are closer to the heuristically expected values 4%. In future simulations, we are planning to explore a near-threshold case with the driving field much closer to the FB instability threshold field. In accord with an aforementioned remark, this may require a much longer simulation time, so that we will use highly-parallelized supercomputers.

Discussion
This paper presents first results of novel 2-D hybrid simulations which are based on continuous equations, rather than on those using discrete particles. The major objective of this  The maximum values of spectral modes are shown are shown for each time step. paper is to demonstrate feasibility and relevance of the new simulation approach. To this end, we have compared our results with results of previous PIC simulations. The comparison have shown reasonable qualitative and quantitative agreement.
This agreement has an important implication for linear and nonlinear theories of the FB and other E-region instabilities. In this paper, we have modeled ions by a kinetic equation with the simplified BGK collision term, Eq. (8). Such model have been used for analytical treatment of the FB instability since the pioneer paper by Farley (1963) and other earlier papers (Lee et al., 1971;Ossakow et al., 1975). The BGK collision model with small modifications (Morse, 1964) provides accurate transfer of the momentum and energy between the colliding particles, but this model does not follow from the rigorous kinetic theory and its applications are limited (Stubbe, 1987(Stubbe, , 1989. As we mention in Sect. 2.3, the BGK model is hardly applicable for electronneutral collisions (Dimant and Sudan, 1995a,c), but it may reasonably approximate ion-neutral collisions in the E-region ionosphere where the colliding particles have approximately equal masses. This assertion has never been verified by computer simulations, bur our results support it. In addition, in all our simulations the nonlinearly saturated stage clearly shows dominant waves with the wavevectors tilted with respect to the E 0 ×B 0 -direction. Because our simulations employ the isothermal model for electrons, see Sect. 2.2, this tilting can only be attributed to the ion thermal-driving mechanism . This suggests that the simplified BGK ion kinetic model employed here includes this mechanism automatically and hence can be successfully used in future theoretical efforts for modeling of not only the FB instability, but the ion-thermal instability (Kagan and Kelley, 2000;Dimant and Oppenheim, 2004) as well. In these simulations, we have modeled the FB instability at a high-latitude electrojet corresponding to B 0 =5×10 4 nT. These results allow a simple scaling to other locations, e.g. to the the equatorial electrojet with the twice as small magnetic field. The difference in the geomagnetic field results in the proportional change of the instability threshold field, see Eq. (A4). The latter also depends upon the altitude via the major altitude-dependent parameter ψ ⊥ ≡ e i /(ν en ν in )∝B 2 0 . Relative intensities of density and electric field fluctuations depend upon the ratio E 0 /E Thr . Typical turbulence temporal and spatial scales are determined by the same ratio and should remain invariant if expressed in terms of the characteristic parameters ν −1 e,i and C s /ν i , where ν i = i √ ψ ⊥ / 0 , while ν e = e 0 √ ψ ⊥ . The gyrofrequencies depend only upon the latitude and are practically altitude-independent within the lower ionosphere, while the parameter 0 , Eq. (6), is essentially invariant throughout the entire E-region ionosphere due to approximate constancy of ν en /ν in 10 ( Kelley, 1989).
Our results might also be applied to modeling of the FB instability at midlatitude sporadic-E layers (Schlegel and Haldoupis, 1994;Haldoupis et al., 1996;Haldoupis et al., 1997), but our current simulations do not include spatial gradients of the undisturbed ionosphere and fields. Gradients of the electric field and plasma density can play an important role in sporadic-E layers. Modeling effects of spatial gradients of the background plasma requires significant modifications of the initial and boundary conditions and, hence, of the general approach to numerical solution. This will be done in future.
Run 4 (blue). Dashed green curves show the isotropic Maxwellian function. The values of the ion distribution function at the adjacent contours (from the center to periphery) decrease by a factor of 1.7.  Table 5).   Table 5).

Conclusions
We present first results of a novel hybrid approach for 2-D simulations of the Farley-Buneman instability in the E-region ionosphere. Unlike the previous hybrid approach based on PIC technique for ions, this technique is based fully on continuous equations: fluid equations for electron density and a kinetic equation for ions with the BGK collision term. The advantage of this kinetic equation is that it includes the crucial effect of ion Landau damping while avoiding noises associated with the finite number of randomly moving particles in PIC methods. Fluid description of electrons allows modeling the real electron mass. The novel hybrid technique can be more suitable than PIC for modeling the FB instability near its threshold.
The 2-D mathematical model includes a nonlinear convection-diffusion equation for electron density, the BGK ion kinetic equation, and Poisson's equation for electrostatic potential. Our simulator can perform numerical computations of the FB instability for different ionospheric conditions. For reasonably chosen parameters, it can be implemented on a PC. However, the developed simulator is optimized for runs on computers with multiprocessor architecture. The first numerical simulations of the FB instability have shown the following major effects: nonlinear saturation of the instability, increasing wavelength in the quasi-steady saturation state, and deviation of the dominating wave vector from the direction of the E 0 ×B 0 -drift velocity of electrons. These results are in good qualitative and quantitative agreement with previous results of fully PIC or hybrid, fluid and PIC, simulations. These first results demonstrate that the new simulation method can be successfully employed in spite of its current deficiencies (simplifying assumptions of the underlying models, as well as relatively small spatial box sizes and a limited 3-D domain of ion velocities). For a more realistic and accurate description of the FB instability, we are planning to extend the new simulation technique to nonisothermal electrons, arbitrarily magnetized ions, and fully 3-D turbulence with bigger and denser simulation grids. In order to successfully implement such improvements, we will employ highly parallelized supercomputers.

Two-fluid linear theory
The FB instability can only be generated if the ambient DC electric field exceeds the minimum threshold field obtained from the two-fluid linear theory. In this Appendix we outline relevant results of this theory.
Standard linear theory implies small density and electrostatic-potential perturbations ∝ exp[i(k·r−ω(k)t)]. Here k is the wavevector andω(k)≡ω(k)+iγ (k), where real ω(k) is the linear wave frequency and γ (k) is the linear growth (γ >0) or damping (γ <0) rate. In the two-fluid model (e.g. Dimant and Oppenheim, 2004, Eq. 1) with ion inertia but no electron inertia, to the zero-order accuracy density and potential perturbations are coupled by The fluid expressions are usually valid for sufficiently longwavelength waves, kC s ν in , where C s =[(T e +T i )/m i ] 1/2 is the ion-acoustic speed for isothermal electrons and ions, so that |γ | ω. Neglecting temperature perturbations and the ion Pedersen velocity but taking into account charge separation between electrons and ions, Eq. (13), one obtains and where V 0 is the E 0 ×B 0 drift velocity of strongly magnetized electrons. The minimum threshold electric field, E min Thr , is determined by equating γ to zero and choosing the optimum wave direction, k V 0 =E 0 ×b/B 0 . As a result, one obtains For ν 2 in ω 2 pi , the altitude dependence of the FB threshold electric field is shown in Dimant and Oppenheim (2004, Fig. 5). The additional, stabilizing, term ∝ ν 2 in /ω 2 pi in Eqs. (A3) and (A4) is non-conventional. It arises due to small charge separation (Rosenberg and Chow, 1998), as explained below. For the E region where ν en 10ν in , we can express this term via the major altitude-dependent parameter ψ ⊥ as If ν in <ω pi then the additional term increases the FB threshold field, whereas if ν in ≥ω pi then it precludes generation of the FB instability for any driving electric field. The term −ν 2 in /ω 2 pi becomes important in the nighttime ionosphere below 100 km where n i 10 9 m −3 , while ψ 0.1 at high latitudes and ψ 0.4 at the magnetic equator (see Dimant and Oppenheim, 2004, Fig. 2). At daytime, this term is negligible for the entire E region. The effect of small charge separation on the FB threshold described by Eq. (A4) should always be taken into account when choosing plasma parameters for numerical simulations. This is the major effect of plasma density on the FB instability.
The physical origin of the additional stabilizing term can be outlined as follows. The small relative charge separation, δn/n 0 ∼(kλ D ) 2 , becomes even smaller as the wavelength grows (i.e. as k decreases). According to Eqs. (A3) and (A4), however, its stabilizing effect is equally significant for all wavelengths. The reason for this is as follows. The two-fluid expressions are usually valid for sufficiently longwavelength waves, kC s ν in , when the terms that describe FB driving and wave dissipation in the linear dispersion relation are second-order small terms, ∝ω 2 ∝k 2 , Eq. (A3), as compared to the first-order terms, ∝ω∝k, that determine the phase-velocity relation, Eq. (A2) (see the corresponding discussion in Dimant and Oppenheim, 2004). Unlike strongly magnetized electrons moving with the E×B 0 drift velocity, weakly magnetized ions are essentially attached to neutrals. The FB driving caused by small ion inertia is described by the first term ∝ω 2 in the RHS of Eq. (A3), while the wave dissipation caused by ambipolar diffusion is described by the last term k 2 C 2 s . To better understand these and other terms it is convenient to pass to the frame of reference moving with the wave phase velocity, V ph =ωk/k (Dimant and Sudan, 1995c;Dimant and Sudan, 1997). In this frame, ion velocity nearly equals −V ph with additional wave perturbations that eventually result in the FB driving term. Small charge separation between electrons and ions results in a slightly excessive flux of ions whose divergence is ∝k 2 V ph , where V ph ≡|V ph |∼V 0 . The contribution of this excessive divergence to the wave excitation/dissipation balance has the same k-dependence as the FB driving term, but the opposite sign. The ratio of the two terms proves to be −ν 2 in /ω 2 pi , as it appears in Eq. (A3).