Equatorial plasma bubbles (EPBs) are large-scale plasma density depletions in the equatorial ionospheric F region, typically forming post-sunset due to the development of the Rayleigh–Taylor instability (e.g., Kelley, 2009; Woodman, 2009). This phenomenon was named EPB because the lower density region grows nonlinearly and penetrates through into the top of the F region. These bubbles can severely disrupt radio wave propagation by inducing scintillation in amplitude and phase, which affects communication and navigation systems that rely on ionospheric propagation. The concept of EPB was proposed by Woodman and LaHoz (1976) based on radar observations and supported by numerical simulations on a magnetic equatorial plane (Scannapieco and Ossakow, 1976). There have been a number of simulation studies of EPBs since the first outcome reported by Scannapieco and Ossakow (1976). The historical review of the numerical simulation studies of EPBs was presented by Yokoyama (2017).

Despite their critical importance, predicting the day-to-day variability of EPB occurrence remains a significant challenge due to the complex interplay of contributing factors. Several studies have addressed the day-to-day variability of the occurrence of EPBs (e.g., Abdu et al., 2009; Carter et al., 2014; Aa et al., 2023), but it was quite difficult to determine a key factor that controls their occurrence. From the modeling approach, the EPB occurrence characteristics were investigated by using a global atmosphere–ionosphere coupled model (Wu, 2015; Shinagawa et al., 2018; Pedatella et al., 2024). The linear growth rate of the Rayleigh–Taylor instability estimated from the simulated parameters shows reasonable seasonal and longitudinal patterns and strong day-to-day variability. Shinagawa et al. (2018) attributed the day-to-day variability to the forcing from the lower atmosphere. The ionospheric altitude variation driven from above (solar and geomagnetic activities) and below (atmospheric activities) makes the occurrence conditions of EPBs more complicated.

This paper aims to address this challenge by utilizing the 3D high-resolution bubble (HIRB) model, which provides a detailed framework for simulating EPB evolution under a range of ionospheric conditions and thereby improves our understanding of their behavior and predictability (Yokoyama et al., 2014; Yokoyama et al., 2015; Yokoyama et al., 2019). The spectral characteristics of the irregularities inside EPBs have been studied using the HIRB model (Rino et al., 2018b; a; Rino et al., 2023), and a comparison with radar observations has been conducted (Tulasi Ram et al., 2017; Tulasi Ram et al., 2020). In this study, we concentrate on the impact of the ionospheric E-region on the generation of EPBs. It has been known that the E-region conductivity contributes to the flux-tube-integrated linear growth rate of the Rayleigh–Taylor instability because the equatorial F region is coupled with the off-equatorial E region by the magnetic flux tube. To the best of our knowledge, however, such contribution of the E-region conductivity has not been carefully studied. Understanding the importance of E-region conductivity will help in understanding the day-to-day variability of EPB occurrence and the prediction of EPB occurrences in the future.

The high-resolution bubble (HIRB) model developed by Yokoyama et al. (2014) is used in this study. It incorporates an advanced 3D numerical simulation framework to accurately replicate the growth and dynamics of EPBs in the equatorial ionosphere. The governing equations in the model are the continuity (Equation 1) and momentum (Equations 2, 3) equations for O + (F region) and NO + (E region), and electrons, and the divergence-free current condition (Equation 4), which are written as: ∂ N j ∂ t + ∇ ⋅ N j V j = S j . (1) e E + V j × B + M j g − ∇ N j k B T N j + M j ν j n U − V j = 0 (2) − e E + V e × B + M e g − ∇ N e k B T N e + M e ν e n U − V e = 0 (3) ∇ ⋅ J = ∇ ⋅ e ∑ j N j V j − N e V e = 0 . (4) where j stands for each ion species, N j , e is the ion/electron density with quasi-neutrality condition ( ∑ j N j = N e ) , V j, e is the ion/electron velocity, S j represents the chemical terms, e is an electron charge, E = E 0 − ∇ ϕ is the electric field, E 0 is the background electric field, ϕ is the electrostatic polarization potential, B is the dipole magnetic field, M j , e is the ion/electron mass, g is the gravitational acceleration, k B is the Boltzmann constant, T = T j = T e is the ion/electron temperature (isothermal condition), ν j n , e n is the ion/electron collision frequency with neutrals, U is the neutral wind velocity, and J is the total current density. Background parameters are obtained from NRLMSISE-00 and IRI-2007: F10.7 is 150, local time is 2000, the day of the year is 83, and the longitude is 135 ° .

The simulation setting in this study is basically the same as those conducted in Yokoyama et al. (2014), except for the plasma density ( NO + ) in the E region. Six different initial conditions were set by increasing the production rate of NO + ions in Equation 1 by factors of 2, 3, 10, 20, and 30. This modification is applied to all latitudes so that the magnetic field lines with any apex altitudes over the dip equator penetrate the E region at the corresponding latitudes. Increasing NO + has negligible impacts on the collision frequency and flux tube electron content gradient. The uniform eastward neutral wind of 120 m s − 1 is applied in the F region, and the background electric field was set to be zero for simplicity.

Figure 1 shows plasma density profiles at the beginning of the simulations. Six solid lines indicate NO + densities for six different simulation conditions, and a dotted line indicates the common O + density for all cases. The difference of the N O + density in the E region between the highest and the lowest cases is less than one order. Then, the initial sinusoidal perturbation resembling a large-scale wave structure (e.g., Tsunoda and White, 1981) is applied by raising the density profile perpendicular to B in the same way as Yokoyama et al. (2014). Figure 2 shows plasma density distribution on a magnetic equatorial plane at T = 3600 s after the beginning of the simulation for the six cases described above. Results at T = 7200 s are shown in Figure 3. It is clearly seen that larger NO + density in the E region weakens the growth of EPBs. The initial seedings in the top three cases eventually turned into structured EPBs in the top of the F region, while the seeding stayed at the bottom of the F region in the bottom three cases. Although the difference of the N O + density in the E region between the highest and the lowest cases is less than one order, it has a strong impact on the growth of EPBs. Our simulations reveal critical insights into the dynamics of EPB growth and its sensitivity to E-region plasma density or conductivity.

The local linear growth rate of Rayleigh–Taylor instability ( γ L ) was given as Equation 5 (e.g., Kelley, 2009) γ L = E z B − g ν i n 1 N ∂ N ∂ z . (5)

This formula does not have an E region contribution to the growth rate. In the equatorial and low-latitude ionosphere, the equatorial F region is coupled with the off-equatorial E region along the magnetic flux tubes. The flux-tube-integrated linear growth rate of Rayleigh–Taylor instability ( γ FT ) was derived by Sultan (1996) and given as Equation 6 γ FT = Σ P F Σ P E + Σ P F E z B L 3 − g e ν eff F K F , (6) where Σ P F and Σ P E are the flux-tube-integrated Pedersen conductivities in the F region and the E region, respectively, L is the McIlwain L -parameter, g e is the downward gravity acceleration, ν eff F is the flux-tube-integrated effective ion-neutral collision frequency weighted by the electron density, and K F is the vertical gradient of flux-tube-integrated electron content in the F region. The recombination rate that would appear in this formula is ignored for simplicity.

Figure 4 shows the local and flux-tube-integrated linear growth rate of Rayleigh–Taylor instability at the initial stage for the six simulation cases. The maximum value of the flux-tube-integrated growth rate for the six cases was 1.093 × 1 0 − 3 , 0.980 × 1 0 − 3 , 0.905 × 1 0 − 3 , 0.668 × 1 0 − 3 , 0.539 × 1 0 − 3 , 0.471 × 1 0 − 3 in descending order. The difference in the growth rate comes only from the factor Σ P F / ( Σ P E + Σ P F ) and stays within approximately a factor of 2 among them. Needless to say, the local growth rates are exactly the same in all cases.

Even minor changes in the linear growth rate could lead to significant differences in EPB growth after a few hours. This finding is particularly important in real applications, where variability in E-region conductivity due to factors such as geomagnetic activity or lower atmosphere phenomena can lead to significant changes in EPB behavior. Furthermore, the impact of E-region conductivity on the temporal characteristics of EPBs may suggest that real-time measurements could be valuable for improving EPB forecasting. By integrating E-region conductivity data into predictive models, it may be possible to enhance the accuracy of forecasts and provide more reliable warnings for communication and navigation systems affected by EPBs.

Our simulation results, unfortunately, emphasize the difficulty of forecasting EPBs based on the growth rate estimation, even though we have access to multiple real-time observations. First, we need to obtain the flux-tube-integrated growth rate, which means ionospheric parameters along the magnetic flux tube, such as E-region plasma density at the off-equatorial regions. This information may only be available at limited longitude sectors where sufficient instruments have been installed. Second, even if sufficient observations are available to estimate the flux-tube-integrated growth rate, the threshold of the growth rate by which the evolution of EPBs should be judged is difficult to define. As shown in this study, a moderate variation of the growth rate becomes a significant difference in EPB growth into the top of the ionosphere. This is a major factor that has made forecasting EPB generation quite difficult for several decades.

This study advances our understanding of equatorial plasma bubble (EPB) dynamics by employing the 3D high-resolution bubble (HIRB) model to simulate and analyze EPB growth under various ionospheric E-region conditions. A key finding is the significant impact of E-region conductivity on EPB development, even when linear growth rates of the Rayleigh–Taylor instability (RTI) show moderate variation. Increased E-region conductivity leads to weaker EPB growth. This underscores the importance of considering E-region conductivity as a crucial factor in EPB forecasting models. Integrating real-time E-region conductivity measurements into forecasting models could further enhance their accuracy and reliability, offering better predictions and mitigation strategies for communication and navigation systems affected by EPBs. However, our results highlight that traditional linear growth rate analyses alone may not fully capture the complexities of EPB behavior and suggest the difficulty of predicting EPB generation in advance. Overall, this research contributes valuable insights into the intricate relationship between ionospheric parameters and EPB formation, emphasizing the need for a holistic approach to EPB modeling. Future work should focus on refining these models and incorporating additional factors to improve forecasting capabilities and better understand the nuances of EPB behavior.