1 Introduction

In the past over 70 years, considerable efforts have been taken to describe and explain the equatorial plasma bubbles (EPBs) phenomena from both an experimental and a theoret- ical viewpoints. It is now generally believed that the physical mechanism responsible for the growth of EPBs is Rayleigh-Taylor (R-T) instability^{[1, 2, 3]}. The R-T instability is excited in the bottomside F region and evolves into a plasma bubble that penetrates the F peak to the topside F region. More recently, with the continuing installation of, and our increased dependence on, global navigation satellite systems, such as the American Global Positioning System (GP(S)), the European Galileo, the Chinese COMPA(S)S/Beidou, the Russian Global Navigation Satel- lite System, and the Japanese Quasi-Zenith Satellite System constellations, the importance of understanding and predicting the occurrence of EPBs has significantly increased^[4]. EPBs oc- currence has complex temporal and spatial variations. To better forecast the generation and development of EPBs, it is necessary to study the R-T instability in the equatorial ionosphere.

Hydrodynamic instability is one of the important topics in the fields such as fluid mechan- ics, plasma dynamics, and atmospheric dynamics. The linear instability theory based on the normal mode method has much more achievements in the hydrodynamic instability research. Rayleigh^[5] deduced the linear instability necessary condition for the parallel shear flow of invis- cid fluid, which is the well-known Rayleigh criterion. Fjortoft^[6] extended the Rayleigh criterion, but the sufficient condition has not been derived. Howard^[7] presented the semicircle theorem, that is, all eigenvalues of the instability models must locate in a given semicircle. In the recent years, researchers all around the world have made significant contributions to hydrodynamic instability studies, including linear, weakly nonlinear, and nonlinear hydrodynamic instability studies^{[8, 9, 10, 11, 12, 13, 14]}.

In the present study, the R-T stability and instability sufficient conditions in the equatorial F region are deduced based on the electron/ions continuity equations, momentum equations, and the current continuity equation^{[15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]}. First, small perturbation linear equations that are vari- able coefficient partial differential equations are obtained. Second, these equations are changed to the eigenvalue equations by the normal mode method, and the sufficient conditions of R-T stability and instability are derived by the variational approach. Last, an approximate numer- ical method of solving systematic eigenvalues is introduced, and a modified numerical model of the development of EPBs is used to validate the sufficient conditions.

2 Basic theory

Using a geomagnetic west, vertically up and geomagnetically north Cartesian coordinate system (x, y, z). The geomagnetic field B is in the z-direction. In the F region, O⁺ is the dom- inant ion, and other ion species are neglected. The electrodynamics of the nighttime equatorial ionosphere can be expressed by the ion/electron continuity equations, momentum equations, and the current continuity equation as follows^{[26, 27]}:

where n is the electron or ion number density, e is the electronic charge, the subscripts e and i represent the electron and ion species, respectively, υ_R is the recombination rate, υ_αn is the electron or ion-neutral collision frequency (the subscript α stands for electron e or ion i), E is the electric field, B is the magnetic field, g is the acceleration due to gravity, W is the neutral wind, and q_α, m_α, and Vα are the electron or ion charge, mass, and velocity, respectively.Ω_α = |B|q_α/m_α is the gyro-frequency, b = B/|B| is a unit vector of the magnetic field, and |B| is the magnetic field strength. In what follows, any motions or variations parallel to the magnetic field are ignored^[15]. E region short circuiting effects are ignored for two reasons. At night, the E region conductivity is reduced, and also when the E region recombination is taken into account, the morphological evolution of F region density perturbations is only slightly altered by the E region short circuiting^{[28, 29]}. Inertial terms in (2) are dropped, and this approximation is certainly valid on the bottomside and up to 450 km or so^[2]. Then, the steady velocity of each species is given by where K_α = E/|B| + g/Ω_α + υ_αnW/ α. Equation (4) is rather unwieldy and needs to be simplified. In the F region, υ_en/Ω_e ≈ 0 and υ_in/ Ω_i ≪ 1. For electrons, the term g/ Ω_e e can be neglected^[15]. Then, the ion and electron velocities are simplified to

Equations (1) and (2) are now perturbed and linearized by using n = n₀+n₁, E = E₀+E₁ = -∇ψ₀ - ∇ψ₁, and '0 is the zeroth-order electric potential. Supposing that υ_in, n₀ and υR are the functions of the altitude (y), and neglecting any second-order perturbation terms, small perturbation linear equations can be readily obtained

where g′ = |B|(g - υ_in²w_x/ Ωi - υ_inw_y). Denote $\hat \nu $_in = ( $\frac{{{\partial _{\nu in}}}}{{{\partial _y}}}$/υ_in, $\hat g$ = (E_0x - g′/υ_in)/n₀, $\hat n$₀ = ( $\frac{{{\partial _{{n_0}}}}}{{{\partial _y}}}$/n₀, and b₀ = n₀/|B|. (7) can be rewritten as

3 Sufficient conditions of R-T instability and stability

To study the electrostatic instability of (8) in the presence of a vertical zeroth-order den- sity gradient, we can write the perturbed electric potential and the plasma density as n₁ = N(y) exp[iα(x - γ t)] and φ₁ = Φ(y) exp[iα(x - γ t)], respectively, where N(y) and Φ(y) are the eigenfunctions, α is the wave number, γ= γ_r + iγ_i is the eigenvalue, and γ_r and γ_i are the real and imaginary parts, respectively. Substituting n₁ and φ₁ into (8) yields

Using the local approximation α≫($\hat n$₀+$\hat \nu $_in)^{[20, 30, 31]}, the term ($\hat n$₀+$\hat \nu $_in)Φ′ in (9) is neglected. We have

Substitution of (10) into (9) leads to

Equation (11) is similar to the Orr-Sommerfeld (O-S) equation in fluid dynamics at the first glance, but it has some differences in the first brackets. The differences can make the necessary condition derived in the fluid dynamics while the sufficient conditions are derived here. (11) is multiplied by $\bar \Phi $ ($\bar \Phi $ is the complex conjugate function of Φ), and the resulting expression is integrated with respect to y over the interval (h, H) with the boundary condition Φ(h) = Φ(H) = 0, which leads to

As it holds that we can obtain ) where Re represents the real part of a complex. From the following relation: it is obtained that Taking the real part on both sides of (12) yields Define where $\hat g$ = (E_0x-g′/υ_in)/n₀, g′ = |B|(g-υ_in²w_x/ i-υ_inw_y), and $\hat n$₀ = ( $\frac{{{\alpha _{{n_0}}}}}{{{\alpha _y}}}$/n₀. The expression can be rewritten as

3.1 Sufficient condition of R-T instability

Let

If γ_i is positive, the ionosphere is R-T unstable, that is, the sufficient condition of the equatorial R-T instability is λ_min > 1.

3.2 Sufficient condition of R-T stability

Similarly, let max

If γ_i is negative, the ionosphere is R-T stable. Therefore, the sufficient condition of the equatorial R-T stability is λ_max < 1. Equations (20) and (21) are equivalent to the following variational equation:

We introduce the functional J[Φ] as

To get the functional extreme values, the first-order variation of J[Φ] is computed as follows:

Then, the functional extreme value problems can be changed to the eigenvalue problems of the second-order differential equation, considering that δ$\bar \Phi $ is arbitrary, i.e.,

The eigenvalues reflect the systematic and comprehensive R-T instability and stability of the whole equatorial ionosphere altitude, rather than the R-T instability and stability of the specific ionosphere altitude or parameter. For eigenvalue problems, many methods and theories have been put forward. Here, we use an approximate method to estimate the systematic eigenvalues.

4 Approximate numerical calculation method

Choose the basis function ψ_j (y)=$\sqrt {\frac{2}{{H - h}}} $sin $\frac{{j\pi (y - h)}}{{H - h}}$ in the L₂(h, H) space, which satisfies the following orthogonal relations:

Let Φ(y) =$\sum\limits_{j = 1}^N {{c_j}{\varphi _j}(y)} $, which satisfies the boundary conditions, where cj are the unknowns. If we denote the matrix C = (c₁, c₂, · · · , c_N)^T, the following balance can be obtained: where B₀ = (b_jk)_N×N is an N × N matrix, and b_jk= ∫_h^H υ_R(α²ψ_jψ_k + ψ′_jψ′_k)dy. Similarly, where A = (a_jk)_N×N is an N × N matrix, a_jk = ∫_h^H-G(y)ψ_jψ_kdy.

Then, (22) can be simplified to the following equations:

Introduce J[C] as

To get the minimum or maximum values, the first-order variation of J[C] shall be equal to zero, i.e.,

Since δC is arbitrary, AC = λB₀C is derived. Substituting AC = λB₀C into (29) leads to

If the eigenvalues of AC = λB₀C are λ₁≤λ₂≤· · ·≤λ_N, λ₁ is the minimum eigenvalue, and λ_N is the maximum eigenvalue. The ionosphere is R-T unstable when λ₁ > 1, and the ionosphere is R-T stable when λ_N < 1.

5 Validation and discussion

To validate the sufficient conditions, some numerical experiments are performed in the follow- ing. In the numerical calculations, thermospheric quantities are obtained from the NRLM(S)I(S)E- 00 model^[32], ionospheric parameters are from the 2007 International Reference Ionosphere (IRI- 2007) model^[33]. The initial ambient electron density profile is superimposed on each contour figure for references, and the altitude ranges from 250 km to 450 km with the electron density peak altitude (hmF₂) at ∼375 km. To make G(y) have the same sign, we take h as 250 km and H as 375 km. For simplicity, the magnetic field strength is regarded as uniform (3×10^-5 T) in the simulation, and the horizontal electric field strength is supposed to be altitudinal uniformity in the initial background ionosphere.

Simulation 1 We just change the horizontal electric field strength, and other parameters are kept the same. The disturbance wave number is 1/km (the wavelength is ∼6 km), and the neutral wind velocity is zero. The minimum and maximum systematic eigenvalues versus the horizontal electric field strength are shown in Fig. 1. The unstable area derived from the sufficient condition is marked in white, and the stable area is in dark. The unstable area is located at the large eastward electric field strength regime. When the eastward electric field is present, the plasma can be lifted under the action of E × B, which steepens the electron alti- tudinal gradient further, and is beneficial to the development of R-T instability. The westward electric field can suppress the electron altitudinal gradient below the peak altitude, which is not conducive to the R-T instability.

Ossakow et al.^[15] have set up the basic equations of the development of EPBs in the equa- torial ionosphere, including the effects of neutral wind and electric field. Huba et al.^[26] showed that the R-T instability is not damped by recombination in the F region, and there is no source term υ_Rn₀ to maintain the electron density at the value n₀ in the nighttime. Therefore, the models of EPBs development can be expressed as follows:

Among many previous works^{[15, 29, 31]}, superimposed on the background-density at t = 0 was a cosine-like perturbation of the form

where l is the initial perturbation wavelength, and x = 0 is the center of the mesh x. That is, when the bottomside electron density has a 5% amplitude initial perturbation, whether can it initiate R-T instability to generate the plasma bubbles or not? The two cases are carried out under different horizontal electric field strength conditions. One is in the unstable area, while the other is in the stable area, as shown in Fig. 1.

Case 1 Figure 2 shows the isodensity contours at different moments with l = 6 km initial perturbation, and the east-west extent is 15 km. The horizontal electric field strength is -2.0× 10-3 V/m. The initial disturbance amplitude is 5%. The white dashed line is the ambient electron density profile at the time t without the initial disturbance. Due to the recombination effects, the ambient electron density decreases gradually, especially in the bottomside F region. At 3 000 s (see Fig. 2(a)), the plasma bubble is quite obvious, and the maximum electron depletion is ∼80%. At 6 000 s, the plasma bubble has penetrated the peak altitude with an innermost depletion contour of ∼100% (see Fig. 2(d)). We also note that the other electron depletions in the wings are formed. The given ionosphere is R-T unstable, which is in accord with the conclusions derived in Fig. 1.

With the 0.5% initial perturbation, the bubble has been clearly formed in the bottomside F region, and the maximum electron depletion is ∼60% at 3 000 s (see Fig. 3(a)). The bubble has penetrated the F peak at 6 000 s (see Fig. 3(d)). When the ionosphere is R-T unstable, even the smaller amplitude initial perturbation (such as 0.5%) can trigger the ionospheric R-T instability and form EPBs.

Case 2 When the horizontal electric field strength is E_0x = 4 × 10^-4 V/m, the ionosphere is R-T stable according to the sufficient condition. Figure 4 exhibits contour plots of electron density at 300 s and 3 000 s with l = 6 km initial perturbation. At 300 s, a cosine-like perturbation can be clearly found on the bottomside electron density with the 5% (see Fig. 4(a)) and 10% (see Fig. 4(c)) amplitude initial disturbances. It is not obvious that the electron depletion and enhancement at 3 000 s are under 5% initial perturbation conditions. The initial perturbations gradually become weakened, and the plasma bubble is difficult to grow and reach the topside F region. With the 10% initial perturbation (see Fig. 4(d)), no obvious electron depletion and enhancement exists. The initial perturbation gradually tends to be stabilized and is difficult to develop. When the ionosphere is R-T stable, even larger amplitude (such as 10%) initial perturbations cannot initiate the ionospheric R-T instability.

The EPBs can be initiated or suppressed when the systematic eigenvalues are in the range of λ_min <1 and λ_max >1 (marked in light dark in Fig. 1), that is, the ionosphere may be R-T unstable or stable in this eigenvalue extent, the corresponding figures of the development of EPBs are not shown here.

Simulation 2 We just change the zonal and vertical wind velocities, respectively, and other parameters are kept the same. The electric field strength is E_0x = -6 × 10^-4 V/m, and the wave number is 1/km. The minimum and maximum systematic eigenvalues are plotted versus the zonal wind velocity (a) and vertical wind velocity (b) in Fig. 5. The minimum systematic eigenvalues increase with the eastward wind velocity. According to the sufficient condition, the eastward wind is conducive to the development of R-T instability. In the evening sector, the horizontal pressure gradient across the terminator produces large thermospheric eastward zonal neutral winds. The equatorial zonal neutral winds interact with Earth’s magnetic field resulting in a vertical polarized electric field. At nighttime, E region plasma densities rapidly decay, the E region dynamo becomes negligible, and the F region dynamo dominates the nighttime plasma drifts^[2]. Under such conditions, the neutral and ionized portions of the atmosphere at F region altitudes move eastward at the same speed. Additionally, the eastward neutral wind at the equator may explain the phenomena that the upper portions of bubbles will take on a westward tilt, while the lower portions will tilt eastward, giving rise to the ‘fishtails’ and ‘Cs’ observed by coherent backscatter radar measurements of field-aligned small-scale irregularities^[34]. Cha- pagain et al.^[35] showed that the nighttime and night-to-night variations in the zonal neutral winds, equatorial plasma bubble velocities, and plasma drifts were well correlated. Kudeki et al.^[36] presented experimental evidence and modeling results which indicated that the eastward thermospheric wind was the primary controlling factor of equatorial spread-F (E(S)F) initiation in the post-sunset ionosphere.

Sekar and Raghavarao^[37] with the linear theory, and Raghavarao et al.^[38] with the nonlinear numerical modeling, suggested that vertical downward (upward) winds have the potential to cause (inhibit) the ESF/bubble phenomenon. Sekar et al.^[39] hypothesized that the variability in the occurrence characteristic of ESF was in connection with the possible day-to-day variability of vertical winds. Raghavarao et al.^[40] showed evidence for the presence of downward winds collocated with irregularities in electric fields and plasma densities as revealed by a unique combination of highly accurate measurements with instruments onboard DE-2 satellite. From Fig. 5(b), it can be seen that the minimum systematic eigenvalues increase with the downward wind velocity, that is, the downward wind is beneficial to the development of R-T instability. The conclusion derived from the sufficient condition is consistent with the experimental and modeling results mentioned above.

6 Conclusion

Sufficient conditions of the R-T stability and instability in the equatorial ionosphere are preliminarily derived as follows: if the minimum systematic eigenvalue is greater than one, the ionosphere is R-T unstable, while if the maximum systematic eigenvalue is less than one, the ionosphere is R-T stable. To validate the sufficient conditions, we design some numerical experiments under different ionospheric parameters conditions. When the minimum systematic eigenvalue is greater than one, the plasma bubble arising from the initial perturbation can develop rapidly and penetrate the peak altitude, indicating that the ionosphere is R-T unstable. While the maximum systematic eigenvalue is less than one, the initial disturbance is difficult to grow and tends to be stabilized gradually, indicating that the ionosphere is R-T stable. The numerically experimental results show that our sufficient conditions are reasonable.

Sufficient conditions are deduced under the assumption of α≫($\hat n$₀+$\hat \nu $_in). In general, the scales of $\hat n$₀ and $\hat \nu $_in are both ∼10^-2/km, and therefore the sufficient conditions are suitable to the initial perturbation of small wavelength. However, observations show that the large scale initial perturbation (∼10² km) can also cause the occurrence of EPBs^[41]. How to take out the precondition α ≫ ($\hat n$₀ + $\hat \nu $_in) is the future work. Moreover, when the systematic eigenvalues are in the range of λ_min < 1 and λ_max > 1, the ionosphere can be R-T stable or unstable, which needs further study.