Modeling of the Gradient-Drift Instability in the Earth’s Equatorial Ionosphere

Abstract

The used mathematical and numerical models of the gradient-drift instability are based on mathematical and numerical models of the Rayleigh-Taylor instability in the Earth’s equatorial ionosphere. These models are constructed as systems of electrodynamically matched multidimensional nonlinear equations described equatorial F-layer of the Earth’s equatorial ionosphere. These instabilities can generate Earth’s plasma irregularities with space-and-time scales that are characteristic of equatorial F-spread. First, the space-and-time model of the developed equatorial plasma bubble is obtained by a two-dimensional numerical simulation of the Rayleigh-Taylor instability. Second, the growth rates of the plasma gradient-drift instability are obtained from the dispersion equation based on physical parameters of the developed equatorial plasma bubble. The results of numerical experiments confirm the possibility of generating the gradient-drift instability of Earth’s equatorial ionospheric plasma.

1. Introduction

Depleted concentration zones in the Earth’s equatorial ionosphere F-layer (named equatorial plasma bubbles, EPBs) are very irregular and nonstationary 3D complex structures. It is of great importance that EPBs have a very scattered range of 3D objects from huge to very little [1, 2, 3, 4, 5]. Under favorable conditions in the Earth’s equatorial ionospheric plasma for the Rayleigh-Taylor instability (RTI), medium-sized initial irregularities develop in a hyperbolic growth into huge EPBs [6, 7, 8]. It should be particularly noted that in this ionospheric process, small irregularities of small sizes do not begin to develop immediately. These small-sized irregularities can develop exclusively on the basis of sufficiently developed EPB and under favorable conditions for a variety of gradient instabilities, for example, the gradient-drift instability.

EPB groups are constantly observed from the outside of the Earth’s ionosphere. EPB in these groups differs greatly in size, as well as in the degree of development. Such measurements are also carried out within the framework of the Defense Meteorological Satellite Program (DMSP) [9]. Artificial Earth satellites DMSP provides cloud cover imagery from polar orbits that are Sun-synchronous at a nominal altitude of 830 km (Figure 1). The DMSP satellite performs valuable EPB observations, although the satellites on their trajectory are over the Earth’s equator for a short time.

Figure 1.

DPMS orbits. (https://en.wikipedia.org/wiki/File:Operational_polar_satellites.jpg).

In this work, the problem of the real feasibility of the occurrence and further deepening of secondary small-sized plasma irregularities fueled by the gradient-drift instability (GDI) [10, 11] at the EPB boundaries is investigated. The existence of such a cascade mechanism is made possible because:

• Huge differences in ionic density at neighboring longitudes occur at the boundaries of deep EPBs.

• Huge differences in ionic density at neighboring heights occur at the boundaries of deep EPBs.

It is difficult to study two instabilities of different scales at once: RTI and GDI. Therefore, the following cascade mechanism was used in this work:

• The plasma density as the basis of the ongoing processes is calculated on a two-dimensional simulation, which is satisfactory due to the condition that the EPB is strongly elongated along the Earth’s magnetic field [6, 12].

• The real possibility of the emergence and further deepening of GDI is determined by the growth rate (GR) of GDI.

2. 3D mathematic model of the Earth’s ionosphere

2.1 3D basic eqs

1. The atmosphere of Earth is the layer of gases, and the ionosphere is the ionized part of the upper atmosphere. The gas dynamics can be described in the most detail in the following three ways [12]:

   • The most detailed known description of gas motion consists of the precise determination of the instantaneous positions and velocities of all the molecules that make up the gas.

   • In practice, a molecular distribution function is introduced that satisfies the integrodifferential Boltzmann equation, by which this function is fundamentally determined by its given initial value. However, the description of the distribution function is too detailed.

   • Therefore, in reality, the equations of thermodynamics are used for calculations. In these equations, the parameters of the state are the velocity of the macroscopic mass flow, density, and temperature. In the derivation of these equations, methods dating back to Grad [13] are used.

   Note that in the equilibrium state, there is no difference between these three approximations.

   Plasma, formed at ionospheric heights, is commonly called thermal, thereby emphasizing that the distributions of electrons and ions are the Maxwell-Boltzmann distribution and can be characterized by their temperatures

   The mean free path of thermal plasma particles is usually significantly less than the characteristic spatial scales of changes in its parameters. Therefore, the thermal plasma in the ionosphere is considered as a continuous medium and is used for its modeling by the equations of magnetic hydrodynamics.

   The mathematical base of the three-dimensional numerical ionospheric model of the RTI, denoted by MI3, is the system of equations of plasma gas dynamics containing the following scalar equations for (N: 1) ion types and electrons [14, 15, 16]

2. N scalar equations of ions and electrons continuity under the assumption of the quasi-neutrality of the ionospheric plasma (MI3/1):

   $$\frac{\partial n_j}{\partial t}+\nabla \left(n_j{\mathbf{V}}_j\right)=Q_j-L_j{n}_j;$$E1

3. N scalar equations of ions and electrons motion (MI3/2):

   $$\frac{\partial {\mathbf{V}}_j}{\partial t}+{\mathbf{V}}_j\nabla {\mathbf{V}}_j=\frac{e_j}{m_j}\left(\mathbf{E}+{\mathbf{V}}_j\times \mathbf{B}\right)-\frac{\nabla p_j}{n_j{m}_j}+\mathbf{g}-{\nu }_{\mathit{jn}}\left({\mathbf{V}}_j-{\mathbf{V}}_n\right)-\sum _{j\ne k}{\nu }_{\mathit{jk}}\left({\mathbf{V}}_j-{\mathbf{V}}_k\right);$$E2

4. N scalar equations of ions and electrons thermal conductivity (MI3/3):

   $$\frac{3}{2}{n}_{j}k\left(\frac{{\mathit{dT}}_j}{\mathit{dt}}+\left({\mathbf{V}}_j\nabla \right)T_j\right)+p_j\nabla {\mathbf{V}}_j+\nabla q_j=G_j-P_j;$$E3

5. Three scalar equations of electric field potentiality (MI3/4):

   $$\nabla \times \mathbf{E}=0;$$E4

6. Three scalar equations of electric current continuity (MI3/5):

$$\nabla J\equiv \nabla \sum e_j{n}_j{\mathbf{V}}_j=0,$$E5

where j is ion and electron type, (N – 1) ions; $Q_j$ and $L_j$ are ion and electron addition and loss rates, respectively; ${\mathbf{V}}_j$, $n_j$, $m_j$, $e_j$, and $p_j$ are drift velocity, density, mass, charge, and gas pressure of ions and electrons, respectively; ${\nu }_{\mathit{jn}}$ and ${\nu }_{\mathit{jk}}$ are ions and electrons collision frequency with neutrals and with each other, respectively; $T_j$ are ion and electron temperature; $G_j$ and $P_j$ are heating and cooling rates, respectively; $q_j$ is heat flux density; the subscripts indicate the ion and electron types; k is Boltzmann’s constant; $\mathbf{E}$ is electric field strength; and $J$ is electric current density.

2.2 Dipole coordinate system

1. In the study area, which is in the altitudinal range from 150 to 950 km, it is possible to make the following assumptions about the transport processes due to the strong magnetization of the ionospheric plasma:

   • Transport processes are determined by collisions along the magnetic field.

   • Transport processes are determined by the drift motion of the plasma across the field.

     Due to the strong anisotropy caused by the Earth’s magnetic field, the processes of diffusion transfer and thermal conductivity in region F of the equatorial ionosphere occur mainly along the lines of force of the geomagnetic field. With a successful choice of the coordinate system, the cumbersome equations and boundary conditions in general will be simplified.

2. Let us choose a dipole coordinate system for the reason that the Earth’s magnetic field is approximated by the dipole approximation [17]. The coordinate components ($\alpha$, $\phi$, and $\beta$) of a curved orthogonal dipolar system are defined in terms of the spherical coordinates of the central dipole as follows:

   $$\alpha =\frac{r}{{sin}^2\theta },\beta =\frac{cos\theta }{r^2},\phi =\phi ,$$E6

   where r is the distance to the center of the Earth, $\theta$ is the latitude, and $\phi$ is the longitude (Figure 2).

3. Now let ${\mathbf{e}}_{\alpha }$, ${\mathbf{e}}_{\beta }$, and ${\mathbf{e}}_{\phi }$ be the coordinate orts in the dipolar system, and ort ${\mathbf{e}}_{\beta }$ coincides with the direction of the geomagnetic field, ort ${\mathbf{e}}_{\phi }$ is directed from west to east, and ort ${\mathbf{e}}_{\alpha }$ is chosen so that it is directed upwards at the magnetic equator. Then, we have the left-handed coordinate system:

$${\mathbf{e}}_{\alpha }\times {\mathbf{e}}_{\phi }={\mathbf{e}}_{\beta },{\mathbf{e}}_{\beta }\times {\mathbf{e}}_{\alpha }={\mathbf{e}}_{\phi },{\mathbf{e}}_{\phi }\times {\mathbf{e}}_{\beta }={\mathbf{e}}_{\alpha }.$$E7

Figure 2.

Dipole coordinate system.

Lame coefficients:

$$H_{\alpha }=\frac{{sin}^3\theta }{\sqrt{\delta }},H_{\phi }=rsin\theta ,H_{\beta }=\frac{r^3}{\sqrt{\delta }},$$E8

where $\delta =1+3{cos}^2\theta =1+3{\beta }^2{r}^4$ and $H_{\beta }={\alpha }^2{H}_{\alpha }{H}_{\phi }$. We also denote the product of the Lame coefficients by $H=H_{\alpha }{H}_{\beta }{H}_{\phi }$.

2.3 Reducing the number of differential eqs

1. Simplify the system of equations by reducing the number of differential equations and the number of unknown functions, respectively [16]. For this purpose, the following conditions of the diffusion approximation are used:

   $$\frac{d{\mathbf{V}}_j}{\mathit{dt}}\equiv \frac{\partial {\mathbf{V}}_j}{\partial t}+{\mathbf{V}}_j\nabla {\mathbf{V}}_j=0.$$E9

   Then

   $${\mathbf{V}}_j={\widehat{\upsilon }}_j\left(\mathbf{E}+\frac{B}{a_j}{\mathbf{A}}_j\right),$$E10

   where ${\widehat{\upsilon }}_j=\left(\begin{array}{ccc}{\upsilon }_{\mathit{jp}}& -{\upsilon }_{\mathit{jh}}& 0\\ {\upsilon }_{\mathit{jh}}& {\upsilon }_{\mathit{jp}}& 0\\ 0& 0& {\upsilon }_{j0}\end{array}\right)$ is the tensor of integral mobility; ${\upsilon }_{j0}=\frac{a_j^{\prime }}{B}$; ${\upsilon }_{\mathit{jp}}=\frac{a_j^{\prime }}{B\left(1+{a_j^{\prime }}^2\right)}$; and ${\upsilon }_{\mathit{jh}}=-\frac{{a_j^{\prime }}^2}{B\left(1+{a_j^{\prime }}^2\right)}$ are longitudinal, Pedersen, and Hall components, respectively; $B=\left|\mathbf{B}\right|$; $a_j=e_j\frac{B}{m_j{\nu }_{\mathit{jn}}}$; $a_j^{\prime }=a_j\frac{{\nu }_{\mathit{jn}}}{{\nu }_{\mathit{jn}}}$; ${\nu }_{\mathit{jn}}={\nu }_{\mathit{jn}}+{\nu }_t$; ${\mathbf{A}}_j=-\frac{\nabla p_j}{n_j{m}_j{\nu }_{\mathit{jn}}}+\frac{\mathbf{g}}{{\nu }_{\mathit{jn}}}+{\mathbf{V}}_n$.

   Similarly, it obtains that

   $$\mathbf{j}=\widehat{\sigma }\mathbf{E}+\mathbf{A},$$E11

   where $\widehat{\sigma }=\sum e_j{n}_j{\widehat{\upsilon }}_j\equiv \left(\begin{array}{ccc}{\sigma }_p& -{\sigma }_h& 0\\ {\sigma }_h& {\sigma }_p& 0\\ 0& 0& {\sigma }_0\end{array}\right)$ is the tensor of integral conductivity across the magnetic field; $\mathbf{A}=B\sum e_j{n}_j{\widehat{\upsilon }}_j\frac{{\mathbf{A}}_j}{a_j}$.

2. The following simplification is related to the rearrangement of the electron density in the ionosphere. Estimates [18] show that the neutrality of ionospheric plasma is maintained with high accuracy. This is enough to accept the condition of neutrality in the equations of motion and continuity of densities. Therefore

   $$\sum e_j{n}_j=0,$$E12

   moreover, the process of establishing a quasi-neutral state does not need to be investigated.

   In addition, the magnetic field is assumed to be constant in time

   $$\mathbf{B}=\text{const}$$E13

   and dipole. This is possible because the variations of the geomagnetic field created by electric currents in the irregularity and in the ionosphere as a whole are negligible compared to the geomagnetic field [18].

3. Due to the condition of electrostatics

   $$\nabla \times \mathbf{E}=0,$$E14

   the electric field is potentially (MI3/6):

   $$\mathbf{E}=‐\nabla \mathbf{\Phi },$$E15

   $\mathbf{\Phi }$ is the electrical potential.

   The dipole model geomagnetic field will not change during calculations, which makes it possible to go from 3D MI3/5 to 2D MI3/7 by integrating the electric current continuity equation:

   $${\nabla }_{\perp }\left(\widehat{\sigma }{\nabla }_{\perp }\mathrm{\Phi }\right)={\nabla }_{\perp }\mathbf{A},$$E16

   ${\nabla }_{\perp }$ acts perpendicular to the force lines.

   Since the electric current density, as well as the velocities of charged particles, can be calculated based on algebraic expressions, they are also not included in the final system.

   The result is 2 N three-dimensional nonstationary nonlinear equations with 2 N unknowns $n_2,\dots ,n_N,T_1,\dots ,T_N,\mathrm{\Phi }$:

   $$\left\{\begin{array}{l}\frac{\partial n_i}{\partial t}+\nabla \left(n_i{\mathbf{V}}_i\right)=Q_i-L_i,i=\overline{2,N},\\ \frac{3}{2}{n}_{j}k\frac{{\mathit{dT}}_j}{\mathit{dt}}+p_j\nabla {\mathbf{V}}_j+\nabla q_j=G_j-P_j,\\ j=\overline{1,N},\\ {\nabla }_{\perp }\left(\widehat{\sigma }{\nabla }_{\perp }\mathrm{\Phi }\right)={\nabla }_{\perp }\mathbf{A}.\end{array}\right.$$E17

2.4 3D final eqs

1. The resulting system of nonlinear three-dimensional strongly coupled multicomponent hydrodynamic equations with nonlocal properties allows for further simplification. The greatest difficulties to solve are the three-dimensional inhomogeneous asymmetric elliptic equation for the electric field potential obtained from the electric current continuity equation. Therefore, this equation will be simplified.

   First, let us take advantage of the fact that there is a dedicated direction in the ionospheric plasma—the direction of the geomagnetic field. The plasma along it at the heights of the equatorial region F is a highly conductive medium due to a large amount of longitudinal conductivity ${\sigma }_0$. Therefore, let the geomagnetic field lines be equipotential, then the electric potential will be independent of the direction of the geomagnetic field

   $$\mathrm{\Phi }=\mathrm{\Phi }\left(\alpha ,\phi \right).$$E18

2. Finally, the dimension of the three-dimensional equation for the potential is reduced by integrating it along the lines of force. At the same time, the integration region should take into account the parameters of the E-region: its lines of force should begin under the E-region in the southern hemisphere, pass through the geomagnetic equator, and end under the E-region in the northern hemisphere (see Figure 2).

   Then, at the ends of the power lines, natural boundary conditions can be set, consisting of the impermeability of the lower boundary of the E-region to the electric current. As a result, the integral of the derivative in the direction of the geomagnetic field turns to zero, and only the following integral remains:

   $$\underset{\mathrm{E}+\mathrm{F}}{\int }H{\nabla }_{\perp }\mathbf{j}\mathit{d\beta }=0,$$E19

   where E and F are E- and F-regions.

3. So, a three-dimensional model MI3 suitable for numerical solution is constructed

   $$\left\{\begin{array}{l}\frac{\partial n_i}{\partial t}+\nabla \left(n_i{\mathbf{V}}_i\right)=Q_i-L_i,i=\overline{1,N},\\ \frac{3}{2}{n}_{j}k\frac{{\mathit{dT}}_j}{\mathit{dt}}+p_j\nabla {\mathbf{V}}_j+\nabla q_j=G_j-P_j,j=e,i,\\ \frac{\partial \left(\underset{\mathrm{E}+\mathrm{F}}{\int }{\sigma }_p\frac{H}{H_{\alpha }^2}\mathit{d\beta }\cdot \frac{\mathrm{\partial \Phi }}{\partial \alpha }\right)}{\partial \alpha }+\frac{\partial \left(\underset{\mathrm{E}+\mathrm{F}}{\int }{\sigma }_p\frac{H}{H_{\phi }^2}\mathit{d\beta }\cdot \frac{\mathrm{\partial \Phi }}{\partial \phi }\right)}{\partial \phi }+\\ \frac{\partial \underset{\mathrm{E}+\mathrm{F}}{\int }{\sigma }_h{H}_{\beta }\mathit{d\beta }}{\partial \alpha }\cdot \frac{\mathrm{\partial \Phi }}{\partial \phi }-\frac{\partial \underset{\mathrm{E}+\mathrm{F}}{\int }{\sigma }_h{H}_{\beta }\mathit{d\beta }}{\partial \phi }\cdot \frac{\mathrm{\partial \Phi }}{\partial \alpha }=\underset{\mathrm{E}+\mathrm{F}}{\int }{\nabla }_{\perp }\mathbf{A}\mathit{d\beta }.\end{array}\right.$$E20

   The well-known global empirical model MSIS was required to determine the characteristics of the uncharged component of the Earth’s ionosphere [19, 20].

   Before starting calculations of the development of irregularities, the model was established by running a background simulation of ionospheric plasma until it reached a periodic solution with a period of a day.

   This 3D MI3 model has been successfully used (e.g., [21]).

3. 2D mathematic models of the Earth’s ionosphere

1. The following simplification is related to the following property of the majority of sufficiently developed medium-scale equatorial irregularities obtained from observational data and numerical three-dimensional experiments on the previous three-dimensional model. Irregularities are strongly elongated along the lines of force and remain so during their movement. Let the irregularities be two-dimensional, then their dynamics are described in the plane of the geomagnetic equator by two-dimensional equations.

   By discarding the terms with longitudinal components in the equations of the three-dimensional model and translating all variables to the coordinates $\alpha$ and $\phi$, a two-dimensional model MI2 of developed irregularities is obtained, suitable for numerical solution:

   $$\left\{\begin{array}{l}\frac{\partial n_i}{\partial t}+{\nabla }_{\perp }\left(n_i{\mathbf{V}}_i\right)=Q_i-L_i,i=\overline{1,N},\\ \frac{3}{2}{n}_{j}k\frac{{\mathit{dT}}_j}{\mathit{dt}}+p_j{\nabla }_{\perp }{\mathbf{V}}_j+{\nabla }_{\perp }{q}_j=G_j-P_j,\\ j=e,i,\\ {\nabla }_{\perp }\left(\widehat{\sigma }{\nabla }_{\perp }\mathrm{\Phi }\right)={\nabla }_{\perp }\mathbf{A}.\end{array}\right.$$E21

   This 2D MI2 model has been successfully used (e.g., [6]).

2. The last significant simplification is based on observational data and numerical experiments on previous models, which show that sufficiently developed irregularities are strongly elongated in height. But this circumstance is used only in the equation for the potential, where only the zonal gradient is taken into account:

   $$\frac{\partial j_{\phi }}{\partial \phi }=0.$$E22

   In addition, let $N=1$, considering only one positive ion, which usually does not lead to incorrect conclusions in the dynamics of irregularities [14].3. Because due to the magnetization of the plasma of the F-layer

   $$a_1=e_1\frac{B}{m_1{\nu }_{1n}},$$E23

   then

   $$\frac{\partial j_{\phi }}{\partial \phi }=\frac{\partial \left({\sigma }_p{E}_{\phi }+A_{\phi }\right)}{\partial \phi }=0,{\sigma }_p{E}_{\phi }+A_{\phi }=\text{const}\left(\alpha \right)=j_{0\phi }\left(\alpha \right),$$E24

   and

   $$E_{\phi }=\frac{j_{0\phi }\left(\alpha \right)-A_{\phi }}{{\sigma }_p},$$E25

   where $j_{0\phi }\left(\alpha \right)$ is the longitude component of the background electric current ${\mathbf{j}}_0=\mathbf{A}$, which must depend on the height only.4. In the simplest case

   $$j_{0\phi }=A_{\phi }\cong n_{10}{m}_1\frac{g}{B},$$E26

   where $n_{10}\left(\alpha \right)$ is the background altitude profile of the ion concentration. Now the zonal component of the electric field $E_{\phi }$ can be found from the algebraic expression.

   The following four equations with four unknowns $n_1,T_1,T_e,E_{\alpha }$ form an express model MI0 of developed irregularities, suitable for both numerical and analytical solutions:

   $$\left\{\begin{array}{l}\frac{\partial n_1}{\partial t}+{\nabla }_{\perp }\left(n_1{\mathbf{V}}_1\right)=Q_1-L_1,\\ \frac{3}{2}{n}_{j}k\frac{{\mathit{dT}}_j}{\mathit{dt}}+p_j{\nabla }_{\perp }{\mathbf{V}}_j+{\nabla }_{\perp }{q}_j=G_j-P_j,j=e,i,\\ \frac{\partial E_{\alpha }}{\partial \phi }=\frac{\partial E_{\phi }}{\partial \alpha }.\end{array}\right.$$E27

   This 2D MI0 model has been successfully used (e.g., [22]).

   An analytical solution of the 2D MI0 model was obtained [23].

4. Gradient-drift instability model

4.1 Mathematic model

1. The two-dimensional model MI2 is used.

   When a deep EPB is formed, very large density differences of charged plasma particles occur on the entire contour of this EPB over the Earth’s magnetic field:

   • Longitude very large density differences of charged plasma particles.

   • Altitude very large density differences of charged plasma particles.

   These huge differences in plasma density favor the emergence and intensification of drift waves growing in the hyperbolic growth regime, when the linear increment is sufficiently large. Figure 3 shows schematic drawing of the strongly deep EPB structure in the plane of the Earth‘s equator with numbers:

   • he horizontal straight line across EPB indicates the maximum of the background F-layer of the Earth’s ionosphere.

   • “1” is located at the place of the greatest drop in plasma density.

   • “2,” “3,” and “4” are in place of the lateral, upper, and lower boundaries with huge drops in plasma density, respectively.

2. To be able to work with the GDI, two equations are added. First, the following notations are introduced: $N_0$ is the unperturbed electron concentration; ${\mathrm{\Phi }}_0$ is the unperturbed potential of the electric field, made dimensionless by dividing by the induction of the magnetic field B; $V_{0\phi }$ and $V_{0\alpha }$ are the undisturbed drift speeds; and ${\nu }_r$ is the recombination rate of electrons. The perturbed values are also denoted but without index 0 [24].

   Then for the perturbed values $N$ and $\mathrm{\Phi }$ taking into account the fact that the conditions ${\nabla }_{\perp }\left(V_{\phi }\right)=0$ and ${\nabla }_{\perp }\left(V_{\alpha }\right)=0$ for the ionospheric plasma are satisfied, two equations are added:

   $$\frac{\partial N}{\partial t}+V_{0\phi }\frac{\partial N}{\partial \phi }+V_{0\alpha }\frac{\partial N}{\partial \alpha }+{\nu }_{r}N+\frac{\mathrm{\partial \Phi }}{\partial \phi }\frac{\partial N_0}{\partial \alpha }-\frac{\mathrm{\partial \Phi }}{\partial \alpha }\frac{\partial N_0}{\partial \phi }+\frac{\mathrm{\partial \Phi }}{\partial \phi }\frac{\partial N}{\partial \alpha }-\frac{\mathrm{\partial \Phi }}{\partial \alpha }\frac{\partial N}{\partial \phi }=0,$$
   $$\mathrm{\Delta \Phi }+\frac{V_{0\phi }}{N_0}\frac{\partial N}{\partial \alpha }+\frac{V_{0\alpha }}{N_0}\frac{\partial N}{\partial \phi }=0.$$E28

   The system of two differential equations is nonlinear since the first equation is nonlinear. The system of these two differential equations can be solved using only the first two harmonics of variation:

   $$\begin{array}{l}N=N_{1}exp\left(-\mathit{i\omega t}+{\mathit{ik}}_{\phi }\phi +{\mathit{ik}}_{\alpha }\alpha \right)+N_{2}exp\left(-2\mathit{i\omega t}+{\mathit{ik}}_{\phi }\phi +{\mathit{ik}}_{\alpha }\alpha \right)+\\ N_{3}exp\left(-\mathit{i\omega t}+2{\mathit{ik}}_{\phi }\phi +2{\mathit{ik}}_{\alpha }\alpha \right),\end{array}$$
   $$\begin{array}{l}\mathrm{\Phi }={\mathrm{\Phi }}_{1}exp\left(-\mathit{i\omega t}+{\mathit{ik}}_{\phi }\phi +{\mathit{ik}}_{\alpha }\alpha \right)+{\mathrm{\Phi }}_{2}exp\left(-2\mathit{i\omega t}+{\mathit{ik}}_{\phi }\phi +{\mathit{ik}}_{\alpha }\alpha \right)+\\ {\mathrm{\Phi }}_{3}exp\left(-\mathit{i\omega t}+2{\mathit{ik}}_{\phi }\phi +2{\mathit{ik}}_{\alpha }\alpha \right).\end{array}$$E29

3. In accordance with refs. [10, 11], under the condition $k_{\phi }<<k_{\alpha }$, the instability increment obtained from this system of two equations has the form

$${\gamma }_g=\frac{V_{0\phi }}{L_{\phi }}\cdot \frac{k_{\alpha }^2}{k_{\phi }^2+k_{\alpha }^2}-{\nu }_r,L_{\phi }={\left(\frac{\partial ln{N}_0}{\partial \phi }\right)}^{-1}.$$E30

Figure 3.

Schematic structure of the developed EPB.

Such a nonlinear increment of increase can be used in the case when an irregularity of the plasma density occurs at the lateral boundaries of the deep EPB. The frequencies of plasma waves in the presented mathematical simulation should not be less than the frequencies corresponding to ionospheric irregularities at the boundaries of which they exist. Also, these frequencies do not rise above the wave frequencies of the charged components of the ionosphere. In this case, we have wave sizes of 10–1000 m.

4.2 Numerical model

1. The computational simulation occupied the plane of the magnetic equator not lower than 100 km and not higher than 1700 km. In width, the two-dimensional simulation space occupied approximately 400 km. A relatively large two-dimensional simulation space serves to exclude the impact of processes at the edges on the results of simulation of the development of irregularities. Numerical simulation was carried out at medium $F_{10.7}=150$ and small $k_p=3$ [24].

   The boundary values of the electric potential are obtained from the undisturbed electric potential:

   • On the right, the model’s undisturbed electric potential was only 0.001 V/m. This value sets the movement from the bottom up with a model speed of 40 m/s.

   • At the bottom of the simulation area, the density of charged particles corresponds to the chemical equilibrium of the processes, and the heating of charged particles is equivalent to the heating of neutral plasma particles.

   • The model value of 0 is taken from above and on the sides for the plasma particle flow rates.

   The model seed EPB is given by the following equation:

   $$n_i=n_{i0}\left(1+\left(a_n-1\right)exp\left(-{\left(\frac{\phi -{\phi }_0}{l}\right)}^2-{\left(\frac{\alpha -{\alpha }_0}{l}\right)}^2\right)\right),$$E31

   where $n_{i0}$ is the background concentration of ions; $a_n$ is the ratio of the concentration at the center of the irregularity to the background value of the concentration at this height ($a_n=0.5$ for the simulation); $\left({\phi }_0,{\alpha }_0\right)$ are the model seed EPB center coordinates; and l is the model seed EPB radius ($l=5$ km for the simulation).

2. The system of differential equations of the mathematical model was approximately solved by approximating derivatives with finite difference on quasi-uniform grids, grinding to the middle of the two-dimensional simulation space to the size of 1.5 km. The solution of the equations of motion has been simplified by the fact that they are algebraic. The solution of the equation for the electric potential was difficult, because it is not just elliptical, but also:

   • It is not self-adjoint due to the presence of Hall conductivity components in the equation.

   • It significantly depends on the ion density of the ionosphere.

   Due to the complexity of solving finite-difference equations, it was necessary to develop approximate second-order equations of accuracy. On the other hand, the use of only two-dimensional simulation of processes in the ionosphere allowed the use of small finite-difference grids, which made it possible to satisfactorily calculate the little scale of some parts of the EPB.

   Numerical solution of mathematical models of Rayleigh-Taylor instability has some difficulties that not all authors manage to overcome [25]:

   • Computational finite-difference schemes for solving transfer equations in numerical experiments with RTI simulation should have satisfactory accuracy (second order).

   • The traditional erroneous way, practiced by some authors for more than half a century, is to neglect the fact that RTI develops in a mode with aggravation. Such a development of events with a poor choice of finite-difference schemes is the reason for a critical increase in the errors of approximate calculations, which is the reason for the nonphysical results of numerical experiments.

3. The finite-difference transfer equations included in the numerical model are solved as follows

   • Nonlinear 2D equations are solved according to the so-called splitting scheme, and this scheme is symmetric, as a result of which the scheme has the second order of accuracy.

   • The 1D equation is solved by a specially developed monotonic method with nonlinear correction of flows for which the minmod limiter was chosen [26, 27].

   The electric potential equation included in the numerical model is solved using the multigrid method with the specially developed W-scheme.

   According to numerous calculations by various authors, the deep EPB has the following characteristics:

   • The width of the “leg” of the mushroom-shaped EPB is approximately $L_{\phi }\approx 1\mathrm{km}$.

   • The width of the “cap” of the mushroom-shaped EPB is approximately $L_{\phi }\approx 10-100\mathrm{km}$.

A numerical experiment to simulate the development of RTI made possible that both the charged particles drift velocity and the background irregularities sizes are obtained.

4.3 Results of numerical experiments

1. The finite-difference simulation MI2 of the RTI development constructed above was the first step necessary to determine the values of the simulated deep EPB parameters of the Earth ionosphere and their differences, which, in turn, are necessary to calculate the values of the GDI development rate. Figure 4 presents the calculated contours of the ionospheric components with the following characteristics:

   • The numerical experiment was carried out on a rectangle with a width of 100 km (longitude) and a height of 800 km (from 150 to 950 km).

   • The process snapshots are shown for the times 2500, 2800, 3100, 3400, and 3700 s, respectively, since the launch of RTI by creating the model seed EPB [24].

2. Figure 5 represents the values of the following parameters, which were obtained during a computational experiment using mathematical-numerical simulation of MI2, varying in time, which is counted since the launch of RTI by creating the model seed EPB

   • Figure 5(a) represents the calculated maxima of the GDI growth increment $1000{\gamma }_g$ (1/s).

   • Figure 5(b) represents the calculated height of the maximal value of the GDI growth increment ${\gamma }_g$ (km).

3. In accordance with Figure 5, which shows the results of a numerical experiment, the calculated linear growth rate of GDI can rise to values of 0.006 1/s, or the characteristic linear development time of GDI can be 170 s. It is obvious that the obtained values of the linear increment of GDI may well be the reason for the rather rapid growth of irregularities of small size, which, in fact, are the main factor in the development of F-spread in the equatorial latitudes of the Earth at ionospheric heights.

Figure 4.

Distribution of electron concentration at different times from the beginning of initialization in the plane of the magnetic equator. The marks on the left at the boundaries of the zones are equal to the decimal logarithm of the electron concentration.

Figure 5.

The values calculated by model MI2 and depending on time from the moment of initialization of RTI: (a) maximum value of the GDI growth increment $1000{\gamma }_g$ 1/s and (b) the height of the maximum value ${\gamma }_g$ km.

5. Conclusions

Using the MI2 model for numerical calculations on a personal computer leads to the following conclusions about the nature of F-spread:

1. Density differences of charged particles of equatorial ionospheric plasma across the Earth’s magnetic field at the lateral boundaries of the cap of a deep mushroom-like EPB can most likely be the initial link of the cascade mechanism from density differences to GDI.

2. The linear growth rate of secondary GDI may well rise to values of 1 /(170 s) and higher, being the trigger for the appearance and deepening of irregularities of small size.

3. The constructed mathematical and numerical simulation is designed to simulate irregularities of tens of meters in size, so the numerical experiment showed that GDI could well provoke equatorial F-spread.