Articles | Volume 41, issue 2
Regular paper
13 Jul 2023
Regular paper |  | 13 Jul 2023

Fluid models capturing Farley–Buneman instabilities

Enrique L. Rojas, Keaton J. Burns, and David L. Hysell

It is generally accepted that modeling Farley–Buneman instabilities requires resolving ion Landau damping to reproduce experimentally observed features. Particle-in-cell (PIC) simulations have been able to reproduce most of these but at a computational cost that severely affects their scalability. This limitation hinders the study of non-local phenomena that require three dimensions or coupling with larger-scale processes. We argue that a form of the five-moment fluid system can recreate several qualitative aspects of Farley–Buneman dynamics such as density and phase speed saturation, wave turning, and heating. Unexpectedly, these features are still reproduced even without using artificial viscosity to capture Landau damping. Comparing the proposed fluid models and a PIC implementation shows good qualitative agreement.

1 Introduction

Magnetized Hall-drifting electrons in the E-region ionosphere induce polarization drifts on the unmagnetized ions which tend to overshoot electrostatic equilibrium and accumulate in the crests of the local density irregularities faster than diffusion opposes them (Sahr and Fejer1996). This mechanism, which results in the amplitude enhancement of local perturbations, is the Farley–Buneman instability. This phenomenon has been shown to modify the mean state of the ionosphere in various ways as well as the magnetosphere morphology by modifying the local conductivity through anomalous heating and nonlinear currents (Wiltberger et al.2017).

Linear fluid theory of Farley–Buneman instabilities predicts some aspects of the dynamics reasonably well. Furthermore, linear kinetic theory shows that ion Landau damping effectively suppresses the growth of smaller wavelengths (Schmidt and Gary1973), which motivates the necessity of resolving kinetic ion effects. Nevertheless, the linear theory fails to explain several features observed in the experimental data obtained with rockets and coherent backscatter radars (Oppenheim et al.1996; Sahr and Fejer1996). Although particle-in-cell (PIC) simulations have been able to model many aspects of the nonlinear physics of Farley–Buneman instabilities (Oppenheim et al.2008; Oppenheim and Dimant2013; Young et al.2020), their application to non-local scales has been very challenging due to the computational cost. This limitation has motivated the exploration of more cost-effective approaches like hybrid and fluid models, which often require much less computational resources because they do not resolve the velocity distribution of the plasma. For instance, Newman and Ott (1981) and Hassan et al. (2015) proposed a fully fluid dynamical system that models Landau damping with a viscosity term on the momentum equation that damps large wavenumbers. Because of the isothermal approximation, these simulations could not capture wave turning and other thermal effects (Dimant and Oppenheim2004).

This work describes a numerical framework based on the five-moment fluid model to simulate Farley–Buneman instabilities and assesses its capability for capturing nonlinear features reproduced previously with PIC simulations. First, we used the viscosity proposed by Hassan et al. (2015), including an energy equation for both species. The simulation parameters were similar to the ones used by Oppenheim et al. (2008) to compare our results with the PIC estimates. Then, we show that most characteristic nonlinear features can be reproduced even after removing the viscosity term from the five-moment system. These results are obtained even though the standard linear theory predicts that smaller structures will grow faster. Furthermore, we argue that the nonlinear signatures obtained without the viscosity are more similar to the correspondent PIC estimates. This last result suggests that the proposed fluid framework may dramatically increase the scalability of Farley–Buneman simulations.

2 A numerical framework for Farley–Buneman instabilities based on fluid equations

The Farley–Buneman instability is an electrostatic process with the dominant dynamics mostly restricted to a 2D plane perpendicular to B. Both the magnetized electrons and unmagnetized ions collide predominantly with neutral particles (Rojas and Hysell2021). Therefore, assuming both species are locally Maxwellian, the following five-moment fluid system should be capable of capturing most of the important physics (Schunk and Nagy2009):


As usual, ns, vs, Ts, ms, and qs correspond to the density, velocity, temperature, mass, and charge of species s, respectively. The collision frequency between species s and neutral particles is represented by νsn, and μsn is the reduced mass of species s and the dominant neutral species. The frame of reference is moving with the neutral particles at a temperature Tn. On the right side of Eq. (3), the term on the left corresponds to collisional heating, while the one on the right captures collisional cooling. The term δsn captures the fraction of energy lost by particle s when it collides with a neutral (Dimant and Oppenheim2004). The electrostatic field was calculated by solving Poisson's Eq. (4).

The R(Ti,vi) term on the right-hand side of Eq. (2) is a general regularization operator, and its purpose is to dampen the growth of larger wave numbers. Here, Ti and vi are the ion temperature and velocity, respectively. Hassan et al. (2015) proposed a regularization term based on the ion viscosity operator but using the ion–neutral collision frequency instead of the Coulomb collision frequency. This term is used to damp large k modes and is only necessary for ion scales (Rojas et al.2016). If we denote the stress tensor by Π, this regularizing viscosity has the form R(Ti,vi)=Π, where

(5) Π = - n i T i ν in v i + ( v i ) T - 2 3 v i I .

We used the operator proposed by Hassan et al. (2015) to model R(Ti,vi) because it was successful in capturing several features of Farley–Buneman irregularities. Furthermore, the accuracy of this proxy will be assessed not by the quantitative estimates of the simulation but by whether it improves the resemblance to PIC simulations.

We chose a spectral solver to solve the five-moment system. Spectral methods are well-known to have outstanding accuracy and to scale very efficiently when periodic boundaries are applicable, and no shocks or discontinuities are expected (Hesthaven et al.2007). These criteria are satisfied in the case of Farley–Buneman irregularities. We build the numerical solver using the Dedalus computational framework for solving general partial differential equations using sparse spectral methods (Burns et al.2020). Our solver uses Fourier spatial discretizations with implicit integration of linear terms. The nonlinear terms are integrated explicitly with 3/2 padding for dealiasing. This solver was comprehensively described by Burns et al. (2020).

Figure 1Electron density perturbation snapshot for both models. Here, δnmax indicates the maximum perturbation for that time instant. The color bar maps the corresponding fraction of δnmax. The top and bottom rows correspond to the regularized and unregularized models, respectively.


3 Simulation setup and results

Some further simplifications can be applied to the five-moment system Eqs. (1)–(4). We omitted the gyro motion term from the ion momentum equation for the ions and used μinmi/2 and δin=1. For the electrons, μenme, no regularization is included, and δen=3.5×104 (Dimant and Oppenheim2004). Furthermore, we used the simulation parameters shown in Table 1. These parameters are the same as the ones used by Oppenheim et al. (2008) for their baseline simulation, except for the grid and box sizes, for which we used 8 times fewer grid points and half the box size, respectively. Notice that the electron mass me and the electron–neutral collision frequency νen have been artificially increased and reduced, respectively. Increasing the electron mass allows the use of larger time steps, but νen has to be reduced to maintain the same magnetization levels. The background electric and magnetic fields are defined in the y^ and z^ axes, respectively.

Table 1Simulation parameters

Download Print Version | Download XLSX

The simulation box size is several times smaller than the largest ones used for recent purely kinetic (Oppenheim et al.2008) and hybrid (Young et al.2017, 2019) Farley–Buneman simulations. In this work, we will use small simulation boxes because we assume that most of the nonlinear features of the PIC simulations are independent of the dimensions of the simulation plane. Even though an energy cascade to larger wavelengths has been reported and it seems to be an issue of insufficiently large simulation boxes, we will not address this point in this work as we are interested in the growth rate of smaller scales. Furthermore, our claim about an improvement in scalability is based on the assumption that fluid models of few moments scale better than PIC kinetic simulations for the same accuracy when kinetic effects are not dominant.

We implemented two models. We will refer to the first as “regularized” and the second as “unregularized”. Both models will solve the continuity, momentum, energy, and Poisson equations, as shown in Eqs. (1–4). The only difference is that the regularized model includes the regularization operator described in Eq. (5) and the unregularized does not.

Figure 1 shows the electron density perturbation ne/n0-1 for both the regularized and the unregularized systems at representative times. Several features are common to both models. At linear regime I, we see dominant wave modes clearly. In the mixing regime II, most wave growth is aligned close to the E×B direction, and perpendicular secondary waves start to form. After saturation, in the turbulent regime III, we see a stable evolution of the electron density perturbation at around 20 % of the background density and density structures similar to the ones obtained with PIC simulations. Moreover, we see a slight turning of the waves in the direction consistent with linear theory, and it is assumed to be present in the nonlinear regime (Dimant and Oppenheim2004). On the other hand, the dominant wavelengths for the unregularized system are smaller (≈1.5 m, similar to the PIC simulation) than for the regularized case (≈2.5 m). This difference is consistent with the idea that the regularization term not only damps larger wave numbers but affects the dynamics of all the wave modes.

Figure 2Time series for the electron temperature (blue dashed), ion temperature (blue), electron density averaged perturbation (red), background forcing field (black dashed), and root mean square of the electric field (black) for the regularized (a) and unregularized (b) models.


Although linear fluid theory predicts that smaller wavelengths will grow faster and destabilize the system without some form of regularization, we can see that the system not only remains stable but can capture several aspects of the expected nonlinear dynamics. This “self-regularization” may be a combination of several factors of numeric and physical origin. Even though spectral methods are well-known for having a minimal numerical diffusion compared to other approaches, this may play a minor role in damping the larger wave numbers. Other physical mechanisms which are captured in this model but not in some versions of the standard linear theory and which play a role in dampening smaller wavelengths are the stabilizing effect of a weak non-quasi-neutrality (Dimant and Oppenheim2011) and electron inertia (Hassan et al.2015). An investigation of the extent of each of these mechanisms is beyond the scope of this paper and its topic of future work. Although the present simulations include the stabilizing effect of thermal dynamics, we have seen similar behavior using an isothermal system, namely, the presence of a dominant mode with no evidence of a growth rate proportional to k2.

Standard metrics to diagnose the simulation, such as the root mean square (rms) of electron density, electric field perturbation, and the average heating, are presented in Fig. 2. The background electric field was raised from a value below the instability threshold at t=0 to the value shown in Table 1 in the first 2ms of the simulation. Both simulations show an increase in temperature due to Pedersen heating in the linear regime indicated as region I. All metrics reach saturation in the mixing regime characterized by region II and then stabilize in region III. The time series of both the regularized and, to a lesser extent, the unregularized simulations present an overshoot just before saturation. This behavior has been documented in hybrid (Rojas and Hysell2021) and fluid (Hassan et al.2015) simulations, but it does not seem to be present in PIC simulations.

Figure 3Time series of electron density perturbation spectral amplitude for specific wave numbers for both the regularized (a) and unregularized (b) systems.


In both simulations, the electron density and the perturbation electric rms field are larger than the corresponding PIC values. The fact that this difference is more considerable for the unregularized case (around twice the PIC metrics) suggests that the origin might be a lack of a proper damping mechanism. Although Oppenheim et al. (2008) also recovered a perturbation electric rms field larger than the background, the difference was smaller. They attributed this excess to factors like the truncation of smaller wavelengths (finite box size).

For instance, the rms field in physical space is the same as in Fourier space, so decreasing the amplitude of larger wave numbers would reduce the total rms. Moreover, the average ion temperature evolves similarly to kinetic simulations in both cases. Nevertheless, the mean electron heating is substantially lower, probably related to the simple temperature-independent collisional heating model used in both simulations. Even though there are several quantitative differences between these and the PIC results, it is interesting to notice that in the unregularized case, the saturation onset time is much closer to the corresponding time in the kinetic simulation (≈60 ms).

Figure 4Electron density perturbation spectra versus wave numbers for the unregularized system. The red dots indicate the maximum amplitude for a particular k, and the green line represents the best fit of kCs(fit) using anomalous γe,i as fitting parameters. Vd0 is the linear estimate of the phase speed using state parameters at the initialization.


Figure 3 illustrates the time series for the electron density perturbation for different wavelengths. Even though the regularized system shows strong damping for smaller wavelengths, both approaches seem to oscillate around similar amplitudes after saturation. This similarity may suggest that capturing the correct damping physics will affect the saturation onset time and the dominant wavelength value. Moreover, the amplitude of the larger wavelengths is comparatively larger for the regularized case.

The spectral properties of region III are summarized in Fig. 4 for the unregularized simulation. The spectral signatures of the regularized case are very similar but with smoother contours. The dominant wave modes propagate predominantly at phase speeds slightly above Cs. Heat capacity ratios representative of region III for both species can be estimated by fitting the expression for the ion-acoustic speed

(6) C s (fit) = γ e T e III + γ i T i III m i ,

to the most representative wave modes indicated by the red dots in Fig. 4. Here, TsIII indicates the average temperature of species s over region III. The ion's heat capacity ratio that produced the best fit with the simulated spectral peaks was γi≈2.62 in both models, which was smaller than the one reported in the corresponding PIC simulations by 14 %. On the other hand, the corresponding electron ratios were γe(r)2.32 and γe(n)2.29 for the regularized and unregularized cases, respectively. The obtained ratios for the electrons were larger than those obtained from PIC simulations by approximately 30 %. Nevertheless, this discrepancy was expected considering the oversimplification of the constant heating and cooling rates used for both species. This simplification is especially limiting for the electron thermal evolution, considering that electron heating associated with Farley–Buneman irregularities has been observed multiple times (Bahcivan2007).

Figure 5Spectra of electron density perturbations for |k|=3m-1 with respect to phase speed and flow angle α for the unregularized system in region III. The spectra were normalized with the maximum power for each flow angle. Vd follows the definition of Fig. 4.


Figure 6Simulated (ds,w(sim)), extended (ds,w(ext)), and fitted (ds,w(fit)) Doppler shift and spectral width calculated for regime III of the unregularized system. The dotted line indicates the local ion-acoustic speed.


The normalized spectra for |k|=3 m−1 are plotted in Fig. 5 with respect to the flow angle and the phase speed. Notice that the dominant modes have narrower widths and lie between the ion-acoustic and the convection speed. Both the Doppler shifts and the spectral widths can be calculated for each flow angle from these profiles. Considering that the system was assumed to be periodic, we can estimate the Doppler shifts and widths for the complete 360 by rotating the simulated ones (ds,w(sim)) appropriately. The values estimated by this rotation are labeled as ds,w(ext). These spectra were fitted to a modified version of a Doppler convection model used in Rojas et al. (2018):


This empirical model relates Doppler shifts and widths to the flow angles θ. The local convection velocity Vd, wave turning angle θ0, a, and b are fitting parameters. The red and blue lines in Fig. 6 correspond to the values of the model after estimating the parameters using a nonlinear optimization routine. We see that this empirical model agrees closely with the simulated values. Moreover, the fitted parameters are consistent with the experimental measurements by Nielsen and Schlegel (1985).

4 Conclusions

The results obtained with the proposed fluid models are consistent with the ones from PIC simulations in reproducing qualitative aspects of the Farley–Buneman instabilities. Moreover, using an unregularized five-moment fluid system, we could still reproduce most of the qualitative aspects of the diagnostics obtained with PIC simulations despite the predictions of standard linear theory. To our knowledge, this is the first time a fully fluid model is able to achieve this. Furthermore, the stabilization mechanisms of the proposed fluid electrostatic model seem to avoid the growth rates γk2. These results suggest that we may have to reconsider the necessity for capturing ion Landau damping accurately.

Even though the simulation box sizes used in this work are small compared to the dimensions used in recent PIC simulations, the potential of scalability relies on the fact that if kinetic effects are not dominant, fluid simulations usually require much less computational resources than PIC implementations to achieve similar accuracy.

Several interesting questions were raised: what nonlinear damping processes modulate the dominant wavelengths and onset saturation time? Would it be possible to build a Landau fluid proxy that captures enough of the physics we have access to by coherent backscatter radars? Would it be possible to modify the proposed fluid system to improve its scalability? Could these results be limited to spectral solvers? We will try to answer these questions in future studies.

Furthermore, we think this has significant implications for the plasma and space physics communities because it may open the door to other researchers having fluid plasma models to explore this topic.

Code availability

The simulations presented in this paper were created using the tools provided in the open-source Python library Dedalus (, Dedalus Project2023).

Data availability

No data sets were used in this article.

Author contributions

ELR proposed the mathematical approach for the simulation, did the numerical experiments, and wrote the paper. KJB contributed with the implementation and testing of the numerical model. DLH helped design the diagnostics to characterize the irregularities. KJB and DLH read the paper and provided feedback.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Financial support

This research has been supported by the Division of Atmospheric and Geospace Sciences (grant nos. 1634014 and 1818216).

Review statement

This paper was edited by Keisuke Hosokawa and reviewed by Matthew Young and Ehab Hassan.


Bahcivan, H.: Plasma wave heating during extreme electric fields in the high-latitude E-region, Geophys. Res. Lett., 34, L15106,, 2007. a

Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., and Brown, B. P.: Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Res., 2, 023068,, 2020. a, b

Dedalus Project: Dedalus – A flexible framework for spectrally solving differential equations,, last access: 11 July 2023. a

Dimant, Y. and Oppenheim, M.: Ion thermal effects on E-region instabilities: Linear theory, J. Atmos. Sol.-Terr. Phy., 66, 1639–1654, 2004. a, b, c, d

Dimant, Y. S. and Oppenheim, M.: Magnetosphere-ionosphere coupling through E region turbulence: 2. Anomalous conductivities and frictional heating, J. Geophys. Res., 116, A09304,, 2011. a

Hassan, E., Horton, W., Smolyakov, A., Hatch, D., and Litt, S.: Multiscale equatorial electrojet turbulence: Baseline 2D model, J. Geophys. Res., 120, 1460–1477, 2015. a, b, c, d, e, f

Hesthaven, J., Gottlieb, S., and Gottlieb, D.: Spectral methods for time–dependent problems, vol. 21, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press,, 2007. a

Newman, A. and Ott, E.: Nonlinear simulations of type 1 irregularities in the equatorial electrojet, J. Geophys. Res., 86, 6879–6891, 1981. a

Nielsen, E. and Schlegel, K.: Coherent radar Doppler measurements and their relationship to the ionospheric electron drift velocity, J. Geophys. Res., 90, 3498–3504, 1985. a

Oppenheim, M., Otani, N., and Ronchi, C.: Saturation of the Farley–Buneman instability via nonlinear electron ExB drifts, J. Geophys. Res., 101, 17273–17286, 1996. a

Oppenheim, M. M., Dimant, Y., and Dyrud, L. P.: Large-scale simulations of 2-D fully kinetic Farley-Buneman turbulence, Ann. Geophys., 26, 543–553,, 2008. a, b, c, d, e

Oppenheim, M. and Dimant, Y.: Kinetic simulations of 3D Farley–Buneman turbulence and anomalous electron heating, J. Geophys. Res., 118, 1306–1318, 2013. a

Rojas, E. and Hysell, D. L.: Hybrid Plasma Simulations of Farley-Buneman Instabilities in the Auroral E-Region, J. Geophys. Res., 126, e2020JA028379,, 2021.  a, b

Rojas, E., Young, M., and Hysell, D.: Phase speed saturation of Farley–Buneman waves due to stochastic, self–induced fluctuations in the background flow, J. Geophys. Res., 121, 5785–5793, 2016. a

Rojas, E., Hysell, D., and Munk, J.: Assessing ionospheric convection estimates from coherent scatter from the radio aurora, Radio Sci., 53, 1481–1491, 2018. a

Sahr, J. and Fejer, B.: Auroral electrojet plasma irregularity theory and experiment: A critical review of present understanding and future directions, J. Geophys. Res., 101, 26893–26909, 1996. a, b

Schmidt, M. and Gary, S.: Density gradients and the Farley–Buneman instability, J. Geophys. Res., 78, 8261–8265, 1973. a

Schunk, R. and Nagy, A.: Ionospheres: physics, plasma physics, and chemistry, 2nd Edn., Cambridge Atmospheric and Space Science Series, Cambridge, Cambridge University Press,, 2009. a

Wiltberger, M., Merkin, V., Zhang, B., Toffoletto, F., Oppenheim, M., Wang, W., Lyon, J., Liu, J., Dimant, Y., and Sitnov, M.: Effects of electrojet turbulence on a magnetosphere-ionosphere simulation of a geomagnetic storm, J. Geophys. Res., 122, 5008–5027, 2017. a

Young, M. A., Oppenheim, M., and Dimant, Y.: Hybrid simulations of coupled Farley–Buneman/gradient drift instabilities in the equatorial E-region ionosphere, J. Geophys. Res., 122, 5768–5781, 2017. a

Young, M. A., Oppenheim, M., and Dimant, Y.: Simulations of Secondary Farley–Buneman Instability Driven by a Kilometer-Scale Primary Wave: Anomalous Transport and Formation of Flat-Topped Electric Fields, J. Geophys. Res., 124, 734–748, 2019. a

Young, M. A., Oppenheim, M., and Dimant, Y.: The Farley-Buneman Spectrum in 2-D and 3-D Particle-in-Cell Simulations, J. Geophys. Res., 125, e2019JA027326,, 2020. a

Short summary
The standard linear fluid theory of Farley and Buneman predicts that kinetic physics are required to avoid the artificial growth of smaller structures. We explore the possibility of simulating the Farley–Buneman instability using, for the first time, a fully fluid five-moment model. This is the first time a fully fluid model has been used to simulate the Farley–Buneman instability. The results obtained with both models are qualitatively consistent with the ones from kinetic simulations.