Multi-instrumental analysis of the day-to-day variability of equatorial plasma bubbles

This paper presents a multi-instrument observational analysis of the equatorial plasma bubbles (EPBs) variation over the American sector during a geomagnetically quiet time period of 07–10 December 2019. The day-to-day variability of EPBs and their underlying drivers are investigated through coordinately utilizing the Global-scale Observations of Limb and Disk (GOLD) ultraviolet images, the Ionospheric Connection Explorer (ICON) in-situ and remote sensing data, the global navigation satellite system (GNSS) total electron content (TEC) observations, as well as ionosonde measurements. The main results are as follows: 1) The postsunset EPBs’ intensity exhibited a large day-to-day variation in the same UT intervals, which was fairly noticeable in the evening of December 07, yet considerably suppressed on December 08 and 09, and then dramatically revived and enhanced on December 10. 2) The postsunset linear Rayleigh-Taylor instability growth rate exhibited a different variation pattern. It had a relatively modest peak value on December 07 and 08, yet a larger peak value on December 09 and 10. There was a 2-h time lag of the growth rate peak time in the evening of December 09 from other nights. This analysis did not show an exact one-to-one relationship between the peak growth rate and the observed EPBs intensity. 3) The EPBs’ day-to-day variation has a better agreement with that of traveling ionospheric disturbances and atmospheric gravity waves signatures, which exhibited relatively strong wavelike perturbations preceding/accompanying the observed EPBs on December 07 and 10 yet relatively weak fluctuations on December 08 and 09. These coordinate observations indicate that the initial wavelike seeding perturbations associated with AGWs, together with the catalyzing factor of the instability growth rate, collectively played important roles to modulate the day-to-day variation of EPBs. A strong seeding perturbation could effectively compensate for a moderate strength of Rayleigh-Taylor instability growth rate and therefore their combined effect could facilitate EPB development. Lacking proper seeding perturbations would make it a more inefficient process for the development of EPBs, especially with a delayed peak value of Rayleigh-Taylor instability growth rate.

It is generally accepted that EPBs are normally developed via the Rayleigh-Taylor (R-T) instability at the postsunset bottomside F-region with a steep vertical density gradient after the diminishing of the E-region. The generalized R-T instability growth rate can be expressed as a function of several flux-tube integrated variables (Sultan, 1996;Equation 26): where ∑ F P and ∑ E P are the flux-tube integrated Pedersen conductivity in the F-region and E-region, respectively; V P is the vertical E×B plasma drift due to zonal electric field E at the geomagnetic equator; U P L is the neutral wind perpendicular to B in the magnetic meridian plane weighted by Pedersen conductivity; g e is the gravity acceleration (positive upward); v eff is the flux-tube integrated effective ion-neutral collision frequency; K F is the altitudinal gradient of electron density in the F-region (K F = 1/N e (∂N e /∂h)); R T is the flux-tube integrated recombination rate. For more details, please refer to Sultan (1996).
The mathematical description of R-T instability has been largely used to help understand the longitudinal and seasonal variation of EPBs (e.g., Wu, 2015;Shinagawa et al., 2018). However, understanding and predicting EPBs' short-term variation, especially its complicated day-to-day variability, has always been a major challenge for the ionospheric community. In particular, EPBs can show unusual development on certain nights but suppression on other nights even under quiet geomagnetic conditions (e.g., Carter et al., 2014;Yamamoto et al., 2018;Abdu, 2019). Considering that the R-T instability growth rate is simultaneously influenced by various ionospheric/thermospheric parameters, the primary controlling drivers of the enigmatic day-to-day variability of EPBs are the following three factors.

Prereversal enhancement
The postsunet enhancement of eastward zonal electric field and upward plasma drift due to increased zonal wind and conductivity gradient, known as PRE (Farley et al., 1986;Heelis, 2004;Eccles et al., 2015), provides a favorable condition to facilitate the growth of the R-T instability by causing the postsunset rise (PSSR) of the equatorial F-layer to higher altitudes with lower ionneutral collision frequency (e.g., Fejer et al., 1999;Sarudin et al., 2021). Along this line, some studies suggested that strong vertical plasma drift exceeding certain thresholds seems to conduct a systematic control on the development of EPBs. For example, Basu et al. (1996) and Anderson et al. (2004) indicated that the equatorial spread-F and scintillation could be generated when the postsunset upward E×B drift exceeds 20 m/s around solar minimum. In comparison, some other studies observed that strong equatorial plasma irregularities could occur near solar maximum when the postsunset upward E×B drift exceeds 50 m/s (e.g., Fejer et al., 1999;Whalen, 2001). Smith et al. (2016) suggested that the necessary PRE peak value that preceding equatorial spread-F development varies between 5 and 30 m/s across different seasons and solar activities. Moreover, a few statistical studies using satellite in-situ measurements have reported that the occurrence of EPBs is approximately proportional to the PRE magnitude, especially from a seasonal/longitudinal point of view (Su et al., 2008;Kil et al., 2009), and the irregularity occurrence probability becomes greater than 80% when the PRE is larger than 40 m/s (Huang and Hairston, 2015). Despite that the climatological correlation is strong, the PRE/PSSR does not always exhibit a clear day-today causal relationship with EPBs occurrence but shows large uncertainties (Abdu et al., 1983;Fukao et al., 2006). Moreover, it was found that local EPBs were not always effectively generated under large PRE/PSSR conditions Maruyama, 2006, 2007), but could be sometimes triggered even without large upward plasma drift Smith et al., 2016). Thus, it seems that PRE itself is not sufficient to explain the quiet-time day-to-day occurrence/suppression of EPBs, and other geophysical factors such as instability seeding sources should also play an indispensable role.

Ionospheric large-scale wave structures associated with atmospheric gravity waves
Considering that the R-T instability originally plays a catalyzing role to boost the ionospheric perturbation structures, it would be an inefficient process that may take thousands of seconds if irregularities were to develop from merely noise-like background fluctuations (Huang and Kelley, 1996). Thus, initial seeding perturbations have also been proposed as a necessary prerequisite for EPBs development to compensate for the otherwise modest strength of R-T instability, since the stronger magnitude of the initial perturbation, the less R-T instability growth rate is needed to amplify non-linearly the fluctuations into plasma bubbles (Retterer and Roddy, 2014). Many observational and theoretical studies have thus indicated that the development of EPBs could be closely related to a seeding precursor of large-scale wavelike electron density fluctuations and/or polarization electric field perturbations at the bottomside equatorial F-layer, which are potentially related to the presence of upward propagating AGWs from the lower atmosphere (e.g., Kelley et al., 1981;Singh et al., 1997;Tsunoda, 2005Tsunoda, , 2015Abdu et al., 2009;Krall et al., 2013;Huba and Liu, 2020). It is thought that AGWs with sufficient energy and vertical length can propagate into the ionosphere before being entirely dissipated, or that the primary waves may break into secondary waves which then subsequently reach ionospheric heights to cause traveling ionospheric disturbances (TIDs) (Vadas, 2007;Yizengaw and Groves, 2018). In particular, the occurrence of EPBs clusters are usually distributed quasi-periodically in longitude with interbubble distances of several hundred kilometers that are generally in agreement with the zonal wavelength of LSWS/TIDs (e.g., Röttger, 1973;Tsunoda and White, 1981;Takahashi et al., 2009;Makela et al., 2010;Taori et al., 2011;Aa et al., 2020a;Das et al., 2020). More importantly, it has been found that LSWS at the bottomside F-layer can provide not only the above-mentioned initial seeding density perturbations but also considerable upwelling of 50-100 km through large undulation of equatorial F-layer that directly contribute to an enhanced R-T growth rate (a.k.a., upwelling growth) (Tsunoda, 2010;Tsunoda et al., 2010). Thus, EPB patches are more likely triggered near each LSWS crest due to the large upwelling growth therein with elevated bottomside density gradient region, forming quasi-periodic distribution structures that are observed in ground-based and space-borne instruments (e.g., Takahashi et al., 2018;Eastes et al., 2020). Some studies suggested that the LSWS-related upwelling is comparable or even outweighs the PSSR to control EPBs development when PRE strength is weak (Tsunoda et al., 2018;Chou et al., 2020).

Neutral wind perturbation
Besides the predominant catalyzing effect of PRE and the seeding effect of AGWs as above-mentioned, EPBs' day-to-day variability might be partially complicated by the thermospheric neutral wind variation. For instance, the eastward thermospheric wind and the E layer conductivity longitudinal gradients after sunset could jointly control the PRE intensity and thereby influence the growth rate of R-T instability (Kudeki et al., 2007;Heelis et al., 2012;Abdu, 2019), although this might be more appropriate to be categorized into the PRE-related effect. On the other hand, the transequatorial wind tends to slightly decrease (increase) the low-latitude Pedersen conductivity in the upwind (downwind) side, thereby playing a destabilization (stabilization) role on the generalized R-T growth rate . The net effect is to increase the field-line integrated Pedersen conductivity and thereby suppressing the non-linear growth rate of the R-T instability (e.g., Maruyama, 1988;Mendillo et al., 2001;Krall et al., 2009). However, it should be noted that the trans-equatorial wind from summer to winter hemisphere is essentially a seasonal pattern, which are less likely subject to large day-to-day variability except when geomagnetic storms cause considerable equatorward neutral wind surge via Joule and auroral heating. During geomagnetic quiet time, the background neutral wind might be modulated by AGWs to cause ionospheric plasma density perturbation through wind components parallel and perpendicular to the geomagnetic field as follows: 1) The magnetic meridional wind perturbation (parallel components) due to AGWs move the plasma up and down along the field line via the neutral-drag effect, producing ionospheric height oscillation and density modulation though this process is not so effective around the equatorial region with small dip angles. 2) The zonal and vertical wind perturbations (perpendicular components) due to AGWs tend to create polarization electric fields Krall et al., 2013;Zhang et al., 2021) and modulate plasma density via E×B drift, which provides important seeding and undulation effects to facilitate EPB development as previously mentioned (e.g., Retterer and Roddy, 2014;Yokoyama et al., 2019). Thus, the short-term variability of EPBs associated with neutral wind perturbation could be intrinsically related to thermospheric waves.
Although significant progress on EPBs' mechanism has been obtained through those pioneering studies, our current knowledge is still limited regarding the relative importance of the concurrent or separate presence of these intertwined drivers in causing EPBs' large day-to-day variability, especially at a geomagnetically quiet time. Thus, in the present paper, we conducted a detailed event analysis of EPBs' day-to-day variability over the American sector during a geomagnetically quiet period of 07-11 December 2019, through coordinately utilizing multi-instrumental ground-based and spaceborne observations, including the Global-scale Observations of Limb and Disk (GOLD) measurements, the Ionospheric Connection Explorer (ICON) data, the global navigation satellite system (GNSS) total electron content (TEC) observation, as well as ionosonde measurements. In particular, we found that the EPBs' activity exhibited a considerably large day-to-day variation around the same postsunset time periods through the selected period. We also conducted an in-depth analysis of potential drivers by calculating the R-T instability growth rate and examining the possible background seeding perturbations, which motivated this study as a good opportunity to advance the current understanding of EPBs' enigmatic day-to-day variability.

Instruments and data description
The GOLD ultraviolet spectrometers observe Earth's airglow emissions between 134 and 160 nm through the disk, limb, and stellar occultation measurements from a geostationary orbit at 47.5°W longitude, which has an unparalleled merit of imaging the American region with an unchanged field-of-view for extended time periods (Eastes et al., 2017(Eastes et al., , 2019. In this study, we use the GOLD nighttime disk images of OI 135.6 nm emission with ∼100 km (longitude) by 50 km (latitude) resolution, which provides continuous time-evolving maps in the early evening hours that can unambiguously specify the spatial-temporal variation of low- latitude ionospheric structures in the Atlantic/American sector, especially the equatorial ionization anomaly (EIA) and EPBs (e.g., Aa et al., 2020a;Eastes et al., 2020;Karan et al., 2020).
ICON is a low-Earth orbit satellite for ionospheric and thermospheric measurements that flies at an altitude of 575 km with an inclination angle of 27°, which is equipped with four instruments: a Michelson interferometer for Global High-Resolution Thermospheric Imaging (MIGHTI) that measures the thermospheric winds and temperatures; an ion velocity meter (IVM) that provides in-situ measurements of ionospheric plasma drift velocity and number density; and two ultraviolet imagers (FUV and EUV) that measure airglow emission to derive ionospheric and thermospheric density and composition (Harding et al., 2017;Heelis et al., 2017;Immel et al., 2018). In this study, we will use the IVM ion density/drift and MIGHTI neutral wind data to investigate the possible connection between EPBs and background ionosphere and thermosphere conditions.
Ground-based GNSS TEC data are derived at Massachusetts Institute of Technology's Haystack Observatory using 6000 + worldwide receivers. The gridded vertical TEC data are provided to the community through the Madrigal distributed data system with a spatial resolution 1°(longitude) by 1°(latitude) and a time cadence of 5 min, with the data quality being examined to remove bad data points and outliers (Rideout and Coster, 2006;Vierinen et al., 2016). The sensitivity of the TEC data is a few percent of the TEC unit (1 TEC unit = 10 16 electrons/m 2 ), which is sufficiently accurate to represent EPBs. Moreover, we here also use the detrended TEC (dTEC) to examine the wavelike TID features as a proxy of AGWs, which is computed by removing the background trend for all satellite-receiver line-of-sight TEC pairs via a Savitzky-Golay low-pass filter (Savitzky and Golay, 1964) algorithm with a 30-min sliding window (Zhang et al., 2017). Besides GNSS TEC, we also use ionosonde measurements from equatorial stations of Sao Luis (2.6°S, 44.2°W) and Fortaleza (3.9°S, 39.4°W) to explore the local irregularity features and/or background ionospheric conditions. In particular, the ionograms are used to check for spread-F echoes, the iso-frequency contours of ionospheric true heights are used to check possible wave activities, and the vertical drift measurements are used to calculate localized Rayleigh-Taylor growth rate. These ionosonde data are automatically derived, which may have some limitations on their accuracy but will not considerably impact the qualitative analysis.

Figures 1A, B
show the temporal variation of interplanetary magnetic field (IMF) Bz, the longitudinally symmetric index (SYM-H), and planetary K-index (Kp) during 07-10 December 2019, respectively. Solar activity was at a very low level during this period with the F10.7 solar radio flux being around 69 SFU (1 SFU = 10 -22 W/m 2 /Hz). The geomagnetic activity was at a quiet condition during this period: the IMF Bz and SYM-H were mainly confined within ±5 and ±10 nT respectively with merely minor fluctuations; most of the Kp indices were ≤2 except for 00-03 UT on December 09 with Kp reaching 3-. This relatively quiet geomagnetic condition suggests that the possibility of significant magnetospheric driving forces causing large ionospheric day-to-day variability is less likely under such a circumstance.
Nevertheless, the low-latitude ionospheric morphology in the American sector, especially the postsunset EPBs activity, exhibited considerable day-to-day variation, even though there was no hint of a geomagnetic storm or substorm onset within a few hours before local dusk. Figures 1C-F show the GOLD nighttime ultraviolet images of OI 135.6 nm radiance at 23:55 UT on December 07-10, respectively, overlapping with ICON orbital path (red line) at close to the same time. Figure 1G-N show the corresponding longitudinal variation of ICON IVM in-situ plasma density and vertical drift measurements. On December 07, we can see from the GOLD image ( Figure 1C) the signatures of EIA crests as two bright lowlatitude zonal bands that were produced by enhanced oxygen ion emission therein. There were also noticeable meridional dark streaks embedded in the equatorial trough and cut through the EIA crests, manifesting typical EPBs structure with low density and reduced emission. The signature of EPBs can also be captured by ICON-IVM in-situ ion density measurements ( Figure 1G) as noticeable plasma bite-outs, which were associated with the strong pre-reversal enhancement (PRE) of large equatorial upward plasma drift of 100 m/s as shown in Figure 1K. For the following two nights, however, the EPBs feature was largely diminished in GOLD images around the same UT interval on December 08 ( Figure 1D) and was almost indiscernible on December 09 (Figure 1E). Similarly, the ICON-IVM measurements during these two nights showed relatively low background ion density conditions with no clear plasma bite-outs (Figures 1H, I), although the plasma vertical drift still showed some spontaneous large values of 100 m/s on December 08 ( Figure 1L). In contrast, strong EPBs re-appeared the next evening on December 10 that were simultaneously captured by GOLD image (Figure 1F) and ICON-IVM in-situ measurements (Figure 1J), associated with large upward plasma drift ( Figure 1N).
In conclusion, the coordinated GOLD and ICON measurements collectively showed that there was large day-to-day variability of EPB activity with considerable suppression on certain nights but intensification on other nights within consecutive 4 days under geomagnetic quiet conditions. Moreover, such a large day-to-day variability of EPBs intensity can also be deduced from local ionosonde measurements. Figure 2 displays a series of ionograms at the equatorial Sao Luis station between 23:30-24:00 UT on 07-10 December 2019, respectively. As can be seen, on December 07 (Figures 2A-D), there were strong diffuse echoes of range-type spread-F traces that spread across the whole F-region, which are typical characteristics suggesting the presence of large-scale ionospheric plasma irregularities that developed via the R-T instability (Abdu, 2001). At the same UT interval on December 08 (Figures 2E-H), the spread-F signatures were still considerable, but the diffuse echoes were less predominant than the previous night and were mainly confined within the bottomside F-region below 350 km virtual heights. Furthermore, on December 09 (Figure 2I-L), there were almost no spread-F features around midnight yet merely some sporadic and limited structures. Nevertheless, on December 10 ( Figure 2M-P), significant spread-F signatures were revived in the same UT interval across the whole F-region, which was much stronger than those from the previous two nights but was comparable to that of December 07. These localized ionogram results are generally in agreement with the above-mentioned GOLD and ICON measurements.
Besides ground-based ionosonde measurements and spaceborne GOLD/ICON observations, strong day-to-day variation of EPBs can also be observed in GNSS TEC maps. Figures 3A-H show the two-dimensional TEC maps over the South American area at 00:00 and 01:00 UT on 08-11 December 2019, respectively. Considering GOLD has no nighttime measurements after 00:25 UT, we here include the TEC map at a slightly later interval of 01 UT to fill the data gap and illustrate the bubble evolution more clearly. As can be seen, significant EPBs occurred on December 08 (Figures 3A, E) as depletion streaks with reduced TEC values that are marked by dashed lines perpendicular to the geomagnetic equator. In contrast, EPBs features were hardly identified on December 09 (Figures 3B, F) and 10 (Figures 3C, G), yet revived on December 11 (Figures 3D, H) at the same UT intervals. Compared with the surrounding area, the amplitude of the TEC depletion within these EPBs was approximately 5-8 TECU. Moreover, to better identify and trace the temporal evolution of EPBs, Figure 3I-L display the corresponding TEC keogram as a function of time and longitude along −10°latitude between 20 and 04 UT during 07-11 December. As encircled by the black ovals, pronounced postsunset EPBs occurred on December 07-08 ( Figure 3I) and 10-11 (Figure 3L), manifesting as parallel comb-like depletion streaks that persisted from 23 to 00 UT to at least 03-04 UT. The inter-bubble distances were estimated to be ∼400-1000 km, which are consistent with typical GOLD observations (e.g., Aa et al., 2020a;Huba and Liu, 2020;Karan et al., 2020). Nevertheless, the EPB intensity on December 08-09 and 09-10 are much weaker and less organized than the first and last days though there were some vague yet blurred streaks.

Variation of R-T instability growth rate
To discuss the potential drivers of such large day-to-day variability of EPBs' intensity, we first examine the R-T instability growth rate variation among the above-mentioned time period of December 07-10, 2019. It is known that Eq. 1 and its derivatives have been widely used in numerical simulations to estimate the flux-tube integrated R-T instability growth rate (e.g., Sultan, 1996;Wu, 2015). However, in order to maximize the usage of realistic observations and avoid the disadvantage of lack of measurements along the flux tube, we here adopt Equation 25 in Sultan (1996) to calculate the local R-T instability growth rate at bottomside F-region around the geomagnetic equator: where γ is the R-T instability growth rate; E and B are the zonal electric field and geomagnetic field magnitude, respectively; g is the gravitational acceleration; v in is the ion-neutral collision frequency; n 0 is the background F-region electron density; z is the altitude. This equation has been proved as an effective approximation to estimate the local linear growth rate of R-T instability near the geomagnetic equator, though its value might be slightly larger than those using flux-tube integrated quantities (e.g., Kelley, 1989;Otsuka, 2018;Das et al., 2021). Note that the meridional neutral wind term is not included in this simplified equation, partially because the nearly parallel-to-B wind in the magnetic meridional plane does not effectively change the conductivity around the equatorial area due to small dip angles therein. The recombination damping term is also excluded since the recombination damp is suggested to be not quite effective (Huba et al., 1996). Other flux-tube terms, such as the integrated Pedersen conductivity in the F-region and E-region, are less likely subject to significant day-to-day variability and thus will not be discussed in the current study. We here adopted a similar method as indicated by Das et al. (2021) and Kelley (1989) to calculate these parameters. Specifically, the velocity factor (E/B) is replaced by the vertical plasma drift, that is, given by Doppler drift mode measurements of Sao Luis digisonde; the inverse vertical gradient scale length of electron density, 1/n 0 (∂n 0 /∂z), is calculated from the bottomside electron density profile of Sao Luis digisonde. The ionneutral collision frequency v in is derived from Kelley (1989) as follows: where n n and n i refer to neutral and ion density. The variable A denotes the mean neutral molecular mass in atomic mass units. These parameters are calculated using NRLMSISE-00 model (Picone et al., 2002). For more details about the mathematical description, readers may refer to Kelley (1989). Figures 4A, B show Sao Luis ionosonde measurements of electron density profiles and F2-layer peak height (hmF2) as well as F-layer vertical plasma drift with error bars during December 07-10, 2019. Figure 4C displays the corresponding temporal variation of the local linear growth rate of R-T instability derived using Eq. 2 based on those Sao Luis ionosonde measurements. In general, the R-T instability growth rate is larger around nighttime, with peak values typically appearing in the postsunset hours around 19-22 LT (22-01 UT) due to the equatorial PRE effect of the uplifted ionosphere as well as the large density gradient in the bottomside F-region after local sunset and decaying of E-region, as previously interpreted in the Introduction section. Besides this typical diurnal variation pattern, one important question is whether the R-T instability growth rate also experienced similar day-to-day variability in the postsunset evening hours as that of observed EPBs activity among these days. We here compare the magnitude and timing when γ reached peak value in each evening among these 4 days: which are 15 × 10 −4 around 22:30 UT (December 07), 18 × 10 −4 around 22:30 UT (December 08), 25 × 10 −4 around 00:50 + 1 UT (December 09), and 22 × 10 −4 around 22:50 UT (December 10), respectively. As can be seen, the day-to-day variation of the peak magnitude of the R-T instability growth rate does not have a strict one-toone relationship with observed EPBs intensity. For instance, the postsunset peak intensity of γ in the middle 2 days are comparable or even slightly larger than those of the first and fourth days, yet the EPB activity was much more suppressed in the middle 2 days as shown in Figures 1-3. Moreover, the peak value of γ was the smallest on December 07 among those days, while the EPBs intensity on that day is considerably stronger than on December 08 and 09. It is worth noting that, on December 09-10, γ reached its peak value at a later time around 00:50 UT about 2 h later than the other evenings, which might partially explain the inhibited EPBs on that day due to this timing difference. However, it should also be noted that the peak intensity of γ (∼ 25 × 10 −4 ) at this time was the highest in the same UT interval among those days, while the TEC keogram in Figure 3G showed that the nighttime EPB activity even after 01 UT on December 10 was still weaker compared with other nights. Therefore, such a complicated day-to-day variability of EPBs can hardly be explained if solely considering the day-to-day variation of R-T instability growth rate. Some other parameters, such as the seeding factor of TID/AGWs, might also play an indispensable role in causing the observed day-to-day variability of EPBs in this event.

TIDs/AGWs activities
We next discuss the background ionospheric conditions in terms of wave structures to investigate the possible seeding effect due to the presence of AGWs on these days. Figure 5 shows four detrended TEC keogram as a function of time and longitude along −10 ∼−20°latitude around the observed EPBs region during 07-10 December 2019, respectively. During the evening of December 07 ( Figure 5A) and December 10 ( Figure 5D), there were noticeable wavelike fluctuations of medium-scale TID features that propagated eastward during 23-05 UT between −70 ∼−40°longitudes. These TID structures are sometimes considered a proxy that may represent ionospheric signatures of AGWs. The zonal propagating velocity of TID wavefronts is estimated to be around 80-100 m/s, which is generally consistent with the background EPB drifting speed derived from ultraviolet measurements given by previous studies (e.g., Immel et al., 2003;Park et al., 2007;Karan et al., 2020). Moreover, the wavelength of TIDs was estimated to be around a few hundred kilometers, which is generally smaller than the interbubble distances. These TIDs were sometimes embedded in the bubbles and were much more evenly distributed across latitudes than bubbles. These TIDs represent medium-scale wave structures that are different from the large-scale bubbles themselves but were present concurrently with the bubbles. In contrary to these two nights with strong TIDs, the TID intensities in the same UT interval were weak during the evening on December 08 ( Figure 5B) and were significantly suppressed on December 09 ( Figure 5C) as marked by the red oval. This lacking of ionospheric wave structures is in agreement with the observations of inhibited EPBs during the same nights.
To better check the AGW activities, Figure 6 shows the temporal variation of F-layer true heights at different plasma frequencies  at Fortaleza ionosonde between 21 and 03 UT (18-24 LT) during December 07-10. Here we use Fortaleza as a substitute since there are some data gaps in the Sao Luis results. On December 07 ( Figure 6A) and December 10 ( Figure 6D), there were identifiable downward phase propagation trends that were marked with dashed lines suggesting the possible presence of AGWs in the upper atmosphere during these two nights (Abdu et al., 2009), while such a downward phase propagation pattern was not so noticeable on December 08 and 09 (Figures 6B, C). This is consistent with those of detrended TEC results. Thus, this coincidence of strong (weak) TID/AGW characteristics together with enhanced (suppressed) EPB activities on the same night collectively illustrate a potential connection between the background wave structures and the development of EPBs. As previously described in the Introduction section, these ionospheric wave structures could not only provide important initial seed perturbations but also large undulation and upwelling of the equatorial F-layer to destabilize the density gradient. Therefore, the stronger (weaker) magnitude the background perturbations have, the less (more) R-T instability growth rate is needed to amplify the fluctuations into the non-linear regime, and the plasma bubbles will be more efficient (inefficient) to be developed (Huang and Kelley, 1996;Retterer and Roddy, 2014). For instance, in the evening of December 09-10, the TID/AGW activity was much weaker than the other nights, and the R-T instability growth rate reached a peak value more than 2 hours later than the other nights. Thus, the combination of these insufficient seeding perturbations and delayed catalyzing factors collectively caused the significant inhibition of EPBs at that night.
Last but not least, we briefly examine the thermospheric neutral wind variation during the selected period of December 07-10, 2019. Figure 7 shows the ICON-MIGHTI observation tracks and corresponding local time, zonal wind, and meridional wind profiles for consecutive orbits focusing on the American sector between 21 and 24 UT that was slightly prior to the observed EPBs on December 07 (Figures 7A-F), December 08 (Figures 7G-L), December 09 ( Figure 7M-R), and December 10 (Figures 7S-X), 2019, respectively. Although the zonal and meridional winds showed generally similar longitudinal distribution patterns with comparable amplitudes among these 4 days, there were still some discernible day-to-day differences. In particular, the horizontal winds on December 07 exhibited large wavelike fluctuations of 100 m/s that were associated with a downward phase propagation trend, especially in Figures 7C, E as marked by the parallel dotted lines. Such strong neutral wind perturbations with a downward phase propagation trend are considered to be signatures of traveling atmospheric disturbances (TADs) due to the modulation of upward propagating AGWs from the lower atmosphere, though the latitude-longitude variations were also mixed with such altitudinal variations. In comparison, the wind profiles on other days also have some moderate oscillations though not as strong as that on December 07. Recall from Figure 4 that the peak value of the linear R-T instability growth rate on December 07 was the smallest for the same postsunset period among those consecutive 4 days, while the postsunset EPB intensity on December 07 was much stronger than the following 2 days as shown in Figures 1-3. This could be possibly related to considerable initial seeding perturbations on December 07 caused by those strong neutral wind fluctuations possibly due to AGWs, which can create polarization electric fields to modulate plasma density via E×B drift and/or directly cause plasma perturbation along the geomagnetic field via ion-neutral collision (e.g., Tsunoda et al., 2010;Krall et al., 2013;Zhang et al., 2021). Such important seeding perturbations could be a necessary prerequisite to effectively compensate for the modest strength of the R-T instability growth rate and thus facilitate the development of EPBs on December 07, otherwise, it would be a more inefficient process if irregularities were to develop from merely noise-like small background fluctuations even with relatively large R-T instability growth rate (Retterer and Roddy, 2014), as was likely the cases for December 08 and 09. It is worth noting that the observed wind oscillations are not exactly around the equatorial area due to MIGHTI's north-looking remote sensing geometry, plus that there were some nighttime data gaps in MIGHTI. However, it is still reasonable to deduce a similar conclusion from the TID/AGWs results shown by the detrended TEC and ionosonde measurements in Figure 6 and Figure 7. A future theoretical simulation is needed to further quantify the neutral wind contribution, which is beyond the scope of the current study.

Summary and conclusion
In this study, we conducted a multi-instrument analysis of the day-to-day variability of EPBs over the American sector and potential drivers under a geomagnetically quiet period on 07-10 December 2019. We coordinately utilized GOLD ultraviolet nighttime disk measurements, ICON IVM in-situ plasma and MIGHTI remote sensing wind measurements, ground-based GNSS TEC observations, as well as ionosonde datasets. The main findings and conclusions are as follows: 1) The intensity of postsunset EPBs exhibited a large day-to-day variation in the same UT intervals around 23-03 UT (20-24 LT) during December 07-10. In particular, the EPBs and spread-F irregularity signatures were quite noticeable on the evening of December 07, yet considerably suppressed on December 08, and almost completely inhibited on December 09, then dramatically revived and enhanced on December 10. 2) The day-to-day variation of the local R-T instability growth rate γ), derived using some actual observations via a simple approach, does not have an exact one-to-one relationship with the observed EPBs intensity. Specifically, the peak magnitude and timing of postsunset R-T instability growth rate among those consecutive four nights of December 07-10 was 15 × 10 −4 /s at 22:30 UT, 18 × 10 −4 /s at 22:30 UT, 25 × 10 −4 /s at 00:50 UT, and 22 × 10 −4 /s at 22:50 UT, respectively. Compared with other evenings, there was a 2-h delay on December 09 when γ reached its peak value, which might partially explain the inhibited EPBs on that night due to this timing difference.
3) The day-to-day variation of TID/TADs feature has a much better agreement with EPBs' activity: there were relatively strong wavelike perturbations preceded and/or accompanying the considerable EPBs on December 07 and 10, whereas relatively weak wave fluctuations occurred on December 08 and 09 corresponding to EPB inhibition despite of large γ values. These TID/TADs structures are a proxy of atmospheric waves in the upper atmosphere, which played an important seeding role that could effectively compensate for the modest strength of R-T instability growth rate in facilitating the development of EPBs (e.g., December 07). To the contrary, the instability growth would be a more inefficient process if EPBs were to develop from merely small fluctuations of background thermosphere and ionosphere conditions even with a relatively large yet delayed peak value of R-T instability growth rate (e.g., December 09).
Our new results indicate that certain seeding factors due to wave activity in the ionosphere/thermosphere can make an important contribution to EPBs' day-to-day variability during geomagnetically quiet time. A future direction is to conduct a longer-term analysis of EPBs' variation and to analyze the contribution of drivers from above and below.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author contributions
EA was responsible for the scientific analysis of the observational results, preparing the manuscript, and organizing team efforts. S-RZ was responsible for differential TEC derivation. AC was responsible for GNSS data management. PE contributed substantially to manuscript development. WR was responsible for software development of GNSS data processing and daily GNSS data processing. All members contributed to science discussion and manuscript development.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.