Occurrence characteristic of plasma bubble was studied using ground-based GPS receiver networks. The occurrence rate of plasma bubble derived from the global GPS network has higher spatial and temporal resolution than that derived from the other observational techniques because of its wide coverage of the observation. The physical characteristics of plasma bubble occurrence were studied in detail with this novel data set. Twenty-three GPS receivers around the dip equator were used to reveal the occurrence from 2000 to 2006. Characteristics of the monthly occurrence rates were different among the regions. Although it was found that sunset time lag effect plays an important role for the monthly variation, two asymmetries which could not be explained with the sunset time lag scenario were found: (1) asymmetry between two solstices and (2) asymmetry between two equinoxes. The seasonal variation of the F-region conductivity integrated along the geomagnetic field line could partially explain the solstice asymmetry. Semiannual occurrence rates from 2000 to 2006 were used to study the year-to-year variation from the high solar activity period to the low solar activity period. The dependency of the occurrence on the solar activity was different among the regions. Occurrence rates against the latitude/altitude were investigated in the Asian region in 2004. It was found that the occurrence was high and constant for a station whose height on the dip equator (HODE) was less than 700 km. They began to decrease when HODE was higher than 700 km and was almost zero where it was higher than 900 km.
 Several physical mechanisms have been proposed to explain the seasonal-longitudinal variation. Plasma bubble is widely believed to be generated by the Rayleigh-Taylar instability. Uplift of the F-layer after the sunset by the eastward electric field, which is known as the prereversal enhancement, makes favorable condition for the Rayleigh-Tayler instability. Tsunoda  showed that the occurrence of plasma bubble was maximum when sunsets in the conjugate E-regions is simultaneous. Simultaneous sunset is expected to intensify prereversal enhancement. The angle between the geomagnetic declination and sunset terminator line, α, was introduced as a proxy of the simultaneous sunset in the conjugate points [Burke et al., 2004]. On the other hand, Rottger  and McClure et al.  emphasized that the initial perturbations by atmospheric gravity waves play an important role to generate plasma bubble. Recently, Kil et al. [2004a, 2004b] suggested that the seasonal longitudinal variations of plasma bubbles could be determined by the plasma density in the nighttime F-region. Among many models on seasonal longitudinal variations, Tsunoda model is recognized as the most plausible model to explain the seasonal-longitudinal variation of plasma bubble occurrence. While this model explains some features of the observed global distributions of plasma bubbles, seasonal-longitudinal occurrence variability has not been understood well. This discrepancy is studied in detail in this study using a large observational data set.
 One of the reasons that the occurrence characteristics have not been fully understood is difficulty of observations for wide area. In order to understand the occurrence distributions of plasma bubbles, global and continuous observational data is required. Statistical data were obtained from ground-based observations, such as ionosondes, radio scintillations, and backscatter radars. Global coverage is, however, difficult to obtain from ground-based observations because the field-of-views and the number of its observational locations has been limited. In situ measurements by satellite were also carried out by Hinotori satellite, and DMSP satellites. Although satellite measurements can cover all longitudes, observational areas are fixed in the geographic latitude and the satellite altitudes. For plasma bubble detection it is necessary for a satellite to fly within ±15 degrees of the geomagnetic latitude after the sunset for plasma bubble detection.
 New observational technique of plasma bubble is ground-based GPS receivers data. Observation with GPS receivers covers large area since the number of the receivers is large. GPS receivers data covers all longitudinal sectors continuously. They measure total electron content (TEC) continuously with little lack of data. It is a strong tool to reveal the occurrence of plasma bubbles. In this study, we used GPS data to detect plasma bubbles. The occurrence of plasma bubble was investigated globally and continuously. The observations of ground-based GPS receiver are described in section 2. In section 3, the occurrence distributions of plasma bubble detected by GPS are presented. Temporal variations of occurrence rates were investigated from 2000 to 2006, and each month among five longitudinal areas. Dependence of the occurrence rates on latitude was also investigated.
2.1. Data Set
 A global observation of the equatorial ionosphere is necessary to clarify the global activity of plasma bubble. GPS receiver networks such as International GNSS Service (IGS), a GPS network by Japan Agency for Marine-Earth Science and Technology (JAMSTEC), and Scripps Orbit Permanent Array Center (SOPAC) can observe wide longitudinal area. The number of GPS receivers of these networks was about 1300 in 2004, all over the world. Twenty-three receivers out of them around the magnetic dip equator were used for the detection of plasma bubble in this study: 17 receivers of IGS network, three of JAMSTEC, two of SOPAC, and one of Solar-Terrestrial Environment Laboratory, Nagoya University. The station codes, their networks, and the geographic longitudes and latitudes are listed in Table 1. Distribution of the receivers is shown in Figure 1. Star marks and filled circles represent the locations of these GPS receivers. The solid curve indicates the magnetic dip equator.
The station code, the network which provides the data, the geographic longitude and latitude, and HODE are listed.
 One of the advantages of ionospheric observation with GPS is its wide coverage. Occurrence characteristics of plasma bubble are known to be different among longitudinal sectors. The difference was interpreted to be caused by the geometry of the geomagnetic field line and the other ionospheric parameters [e.g., Maruyama and Matsuura, 1984]. We studied the occurrence of plasma bubble in five longitudinal regions. These regions are defined based on the geomagnetic declination angle and longitude, Central Pacific (CP; 160°E–120°W), Eastern Pacific (EP; 120°W–60°W), Atlantic (AT; 60°W–20°W), African (AF; 20°W–60°E) and Asian (AS; 60°E–160°E) regions as shown in Figure 1. More than two GPS receivers were used in each region in this study. The GPS receivers for each region are listed in Table 2. Global occurrence of plasma bubble was studied with receivers represented by star marks for 7 years from 2000 to 2006.
Table 2. Five Regions and Their GPS Stations Whose Data Is Used in Figures 5–6a
Averaged HODE (km)
The region was classified according to the geographic longitude and the geomagnetic declination (see text). HODE for each station is averaged and shown on the third row.
IISC, SAMP, NTUS, BAKO, TVST, GUAM, CNMR
THTI, ASPA, FALE
AREQ, ARE2, RIOP, GLPS, GALA
 Latitudinal dependence of plasma bubble was also studied using a meridional array of GPS receivers in the Southeast Asian region, which are represented by filled circles in Figure 1 around 100°E of longitude. As a proxy of the geographic latitude, “height on the dip equator (HODE)” is introduced. The schematic illustration of HODE is presented in Figure 2. HODE is defined as an altitude of the geomagnetic field line on the magnetic dip equator which passes 400 km altitude above a GPS receiver site. HODE for each receiver is listed on the third row of Table 1.
2.2. TEC Variation and Rate of TEC Index
 TEC was derived in every 30 s from GPS data of the GPS receivers. At the time of plasma bubble occurrence, GPS-TEC fluctuates rapidly. Examples of the diurnal TEC variation with and without plasma bubble are shown in Figures 3a and 3b. The variations were observed with the receiver at Bukittinggi (BUKT), which is located at 0.2°S and 100.3°E in the geographic coordinates, on 22 and 23 March 2004, respectively. The horizontal axis represents the local time for 1 day. The vertical axis shows slant TEC in TEC Unit (TEC; 1TECU = 1016el/m2). The TEC data were derived for all ray paths between GPS satellites and the BUKT GPS receiver whose elevation angle larger than 45°. The instrumental biases of each transmitter and the receiver are contained in the TEC values in the top of Figures 3a and 3b. U shapes in TEC variations were due to the change of the ray path length inside of the ionosphere. In Figures 3a and 3b, daily variations were seen, that is, high TEC in the daytime and low in the nighttime. Intense TEC depletions and fluctuations were observed soon after the sunset in Figure 3b. These were caused by plasma bubble that is characterized by large depletion of the ionosphere electron density.
 The severe electron density depletions inside of plasma bubble cause small-scale fluctuations of TEC. The small-scale fluctuations of TEC were utilized to detect plasma bubble using ground-based TEC in this study. The TEC depletions of plasma bubble are difficult to be distinguished from the other observed slant-TEC variations because the timescale of TEC depletions of plasma bubble is an order of 1 h that is comparable to the timescale of slant-TEC variations caused by the change of the ray path length inside of the ionosphere. Furthermore, the amplitudes of them are comparable. The small-scale fluctuation of TEC with 5 min period has a large amplitude inside of plasma bubble [Beach and Kintner, 1999], while they have little amplitudes outside of plasma bubble. Therefore the small-scale fluctuation of TEC is a good proxy to measure the activity of plasma bubble using the ground-based GPS-TEC data. The time derivative of TEC is often used to measure small-scale ionospheric fluctuations. Time differential TEC is called rate of TEC change (ROT) [Pi et al., 1997], that is, R(t) = S(t + Δt) − S(t), where R(t) and S(t) are ROT and slant TEC at time = t, respectively. The sampling interval is Δt = 0.5 (min). The value of S(t) contains instrumental biases of each transmitter and receiver. Since the instrumental biases are almost constant during the sampling interval of 30 s, subtracting S(t) from S(t + Δt) can cancel the biases. It is also necessary to remove the effect of the radio wave ray path length change inside of the ionosphere from R(t). The ionosphere was assumed as a thin shell, and corrected ROT is approximately given as R′(t) = cosθ(t)(S(t + Δt) − S(t)), where θ(t) is a zenith angle of the line-of-sight from the receiver to the satellite.
 To identify small-scale fluctuations, the standard deviation of ROT was used. This value is called rate of TEC change index (ROTI) [Pi et al., 1997]. It is often used to investigate ionospheric fluctuations [e.g., Beach and Kintner, 1999; Bhattacharyya et al., 2000]. The time window of the standard deviation of ROTI reflects the scale of the phenomenon. Five minute window is widely used in studies for plasma bubbles and high-latitude irregularities [Beach and Kintner, 1999]. We calculated ROTI with 5 and 30 min windows and found no significant difference between them. Therefore we used 5 min window ROTI in this study in order to make it easy to compare with the previous studies using ROTI. The spatial scale of the phenomenon which is detected by 5 min ROTI corresponds to the displacement of the pierce point of the GPS radio wave at the ionospheric height for 5 min. When the ionospheric height is assumed at 400 km, velocity of the pierce point of the GPS radio wave at the ionospheric altitude is approximately 70 m/s around the zenith, therefore the scale of the phenomenon is ≃20 km. Size of plasma bubble is 100 km ≃ 400 km in zonal direction [e.g., Fukao et al., 2006; Tsunoda, 2005]. ROTI detects the substructures inside plasma bubbles.
 Field of view of the GPS observation whose zenith angle was less than 45° was ≃ 800 km scale. Radio wave signals from more than four GPS satellites can be received within the field of view at most of times. The maximum value of ROTIs derived from each GPS satellite signal was used as the ROTI value for the 5 min period. Considering the scale size of plasma bubble and FOV of GPS receiver, the maximum values is better proxy than the average value.
 Eastward zonal velocity of plasma bubble is reported as 100 m/s ≃ 200 m/s [Fukao et al., 2006; Lin et al., 2005]. Therefore, plasma bubble could be observed for at least 30 min by a GPS receiver. In order to detect plasma bubble, one 30 min ROTI was derived from six values of 5 min ROTI. Five minute ROTI has also large values when errors in TEC measurement happen. The largest ROTI among the six values was excluded to remove the errors. The average of the remaining 5 min ROTIs was defined as a 30 min ROTI. In this paper, this 30 min ROTI is referred as ROTI and used for the analysis.
 In the bottom of Figures 3a and 3b, ROTI on 22 and 23 March 2004 were plotted, respectively. The error bars of ROTI show the standard deviation of five values of 5 min ROTI. The vertical dashed lines show the local noon time, sunset time, and midnight time at 400 km altitude over the receiver site. It is seen that the value of ROTI was low and constant during the day time in both days. ROTI was almost constant on 22 March 2004, even in the nighttime. While on 23 March 2004, it was large after the sunset time. The ROTI value for the satellite number (Pseudo Random Noise: PRN) 30 satellite was plotted with the thick line. In the top of TEC plot, data of PRN 30 was also plotted with a thick line. The large value of ROTI corresponds to the sharp depletion of TEC. Good correspondences between plasma bubbles detected in TEC data and ROTI were seen for the other satellites. ROTI was used to identify the occurrence of plasma bubble using the ground-based GPS-TEC data.
2.3. Detection of Plasma Bubbles
 Plasma bubble is a phenomenon which appears after the sunset. It causes TEC fluctuations and high ROTI value. The difference of the daytime ROTI and the evening ROTI were used for detection of the occurrence of plasma bubbles in this study because noise level of ROTI is different for each GPS receivers. Three steps were conducted to identify plasma bubble for 1 d at one station, as follows:
 1. Rev, evening ROTI, was derived with ROTIs from sunset to midnight. This time period is defined as the period when the solar zenith angle at 400 km altitude over the station was from 90° to maximum. The average of 30 min ROTIs during this period was defined as Rev. In the case of 22 and 23 March 2004 at BUKT which area shown in Figures 3a and 3b, Rev were 0.07 and 0.68 TECU/min, respectively. Using ROTIs from sunset to midnight as Rev, avoids the effect of a loss-of-lock of radio wave according to plasma bubble. When plasma bubble occurs and the ionosphere is disturbed severely, a loss-of-lock of GPS radio waves occurs, and GPS measurement becomes impossible for 3 h at longest. This data gap has been one of the problems of TEC data when ionosphere is disturbed [Aarons, 1993; Kelley et al., 2002]. Using 6 h data from sunset to midnight makes it possible to detect plasma bubbles even if data lacks occurs even for a few hours.
 2. Daytime ROTI, Rday was the average value of ROTI during 3 h from the noon. In the case of Figures 3a and 3b, Rday were 0.04 and 0.05 TECU/min, respectively.
 3. The differential value, Rev − Rday, was used as an index of plasma bubble activity for one day at one station. In the case of Figures 3a and 3b, this value was 0.03 and 0.63 TECU/min, respectively. Figure 4 (top) presents Rev − Rday at BUKT from 1 March 2004 to 28 February 2005, though there were data lacks in February 2005. There were maxima of Rev − Rday in March and September. It is known that plasma bubbles occur most often in March and September in the Asian region [e.g., Aarons, 1993; Burke et al., 2004]. Distribution of Rev − Rday from March 2004 to February 2005 is shown as a histogram in Figure 4 (bottom). The histogram shows that most of Rev − Rday were between −0.05 TECU/min and +0.05 TECU/min. The number of days when ∣Rev − Rday∣ was less than 0.05 TECU/min was 267 d out of 332 d. Variation of Rev − Rday within ±0.05 TECU/min could arise from difference of GPS satellite trajectories and temporal variations of small fluctuations of electron density which were not related to plasma bubble. Larger values of Rev − Rday than 0.05 TECU/min could be caused by plasma bubble. In order to identify plasma bubble, a threshold was set to Rev − Rday = 0.075 TECU/min. The threshold 0.075 TECU/min was determined by adding a margin of 0.025 TECU/min to 0.05 TECU/min to exclude the TEC variations that were not generated by plasma bubble. The number of days when Rev − Rday was less than 0.075 TECU/min was 259 d, while the number of days when Rev − Rday was more then that was 73 d.
 Plasma bubbles are known to be affected by geomagnetic disturbance [Fejer and Scherliess, 1997; Scherliess and Fejer, 1997]. In order to investigate occurrence rates under the geomagnetically quiet condition, the data in the geomagnetical active period was not used. Kp index was used to determine the geomagnetic active days. Two Kp index values from 1500 to 2100 in local time were averaged for the Kpsunset. The analysis of plasma bubble detection was done when Kpsunset ≤ 5.
3.1. Monthly Variation
 Monthly occurrence rates of plasma bubble were derived from 2000 to 2006 with the procedure described in the previous section. Five regions are defined according to their geographic longitude and the geomagnetic declination angle, δ. It has been reported that the occurrence of plasma bubble was affected by the geomagnetic declination angle [Tsunoda, 1985] The five regions and their range of δ are as follows: Asian (−2° < δ < 2°; AS), African (−1° < δ < 4°; AF), Atlantic (−20° < δ < 5°; AT), Eastern Pacific (−2° < δ < 5°; EP), and Central Pacific (10° < δ < 12°; CP). The GPS receivers in each region, are listed in Tables 1 and 2. HODE is averaged in each region and listed in the third column of Table 2. Monthly occurrence rates are shown in histograms in Figure 5. The ratio of days when plasma bubble appeared to the observation days in a month at one receiver was averaged in each region. The averaged value was defined as a monthly occurrence rate in the region. The scale is shown on the left-hand side axis of Figure 5. The numbers of observation days, that is, the number of days when the data were available in each month was averaged in each region and represented with lines in Figure 5. The scale of observation days are shown on the right-hand side axis.
 Two maxima were seen in annual occurrence variations in the African, Asian, and Atlantic regions. The maxima of occurrence rates in the Asian region were in March/April and September/October. Tsunoda  explained that the seasonal-longitudinal variation of plasma bubble is controlled by the geometry between the geomagnetic field line and the terminator line. In this study, the sunset time lag between the geomagnetic conjugate points was introduced to represent the geometry of the geomagnetic field line and the terminator line. Sunset time lag is shown with solid curves on Figure 5 (bottom). The geomagnetic field line where the sunset time lag was calculated is one which passed 400 km altitude on the dip equator, where longitude is the same as the GPS receiver site. The difference between the sunset time at the foot points of the geomagnetic field line at 100 km altitude in the both hemisphere was defined as sunset time lag. Solid curves in Figure 5 were averaged sunset time lags for all stations listed in Table 2. The time lag was absolute value in unit of hour with the scale of the left-hand side axis. Occurrence rate of plasma bubble is expected to be high during the period when the sunset time lag is small, in other words, when the sunsets of the geomagnetic conjugate points are synchronized. We defined this day as “magnetic equinox” in this paper. The magnetic equinoxes in the Asian region is around the end of March and September, and the occurrence rate of plasma bubble was maximum in March/April and September/October in the Asian region. Although the maxima were around the magnetic equinoxes in the Asian region, some features of the annual variation of plasma bubble occurrence were found to be inconsistent with the sunset time lag scenario. Major discrepancies between the annual variation of the occurrence, and the sunset time lag scenario are following two asymmetries: (1) asymmetry in occurrence between two solstices and (2) asymmetry in occurrence between two equinoxes. In the Asian region, sunset time lags are the similar between around two solstices as shown on Figure 5. On the other hand, occurrence rate of plasma bubble was larger around June solstice than around December solstice from 2000 to 2003. Asymmetry 2 was seen in 2005 and 2006 in the Asian region. Occurrence rates in March equinox (ME) were larger than in September equinox (SE). These asymmetries were seen in the other regions than the Asian region. Asymmetry 1 was clearly seen in the African region from 2000 to 2006, that is, occurrence rates were larger in June solstice than in the December solstice though sunset time lags had no difference between two solstices. In the Eastern Pacific and Atlantic regions, the occurrence rates in the December solstice were larger than those around the June solstice. The asymmetry in the Atlantic region could be explained with sunset time lag scenario because the sunset time lag around the December solstice was smaller than that around the June solstice. In the Eastern Pacific region, however, the sunset time lag around the December solstice was larger than that around the June solstice. The other mechanism than the sunset time lag effect are seems to contribute to the occurrence of plasma bubble in the African, Asian, and Atlantic regions. Asymmetry 2 was clearly seen in the Central Pacific region. In the Central Pacific region had no maximum around the ME. The annual variation without ME maximum repeated for years in this region. Even the regions where there were two maxima, asymmetry between the occurrence rate between in the ME and the SE were seen. In the African region, the occurrence rate in the SE was larger than those in the ME from 2002 to 2005. In the Eastern Pacific region, the occurrence rates in the ME were larger than those in the SE in 2005 and 2006.
3.2. Year-to-Year Variation
 Occurrence rate for 7 years from 2000 to 2006 was derived for the five regions. Figure 6 shows the semiannual occurrence rates in each region during these 7 years. The GPS receivers in each region were listed on the Table 2. Occurrence rate was defined as a ratio of the number of days when plasma bubble was detected to the number of observation days in a half year from January to June and from July to December. The ratio was plotted as semiannual occurrence rate of plasma bubble with the scale on the left-hand side axis. The triangle, diamond, square, asterisk, and cross marks represent Atlantic (AT), Eastern Pacific (EP), Central Pacific (CP), Asian (AS), and African (AF) regions, respectively. In order to compare the occurrence rates with the solar activity, F10.7 was plotted in Figure 6. F10.7 was represented with filled circles with the scale on the right-hand side axis.
 Variation of semiannual occurrence rate for seven years was different among the regions. In the Asian region, the semiannual occurrence rate was maximum in early 2002, that was 70%. It decreased year by year to be about 10% in 2006. On the other hand, in the Atlantic region, the occurrence rate was between 20% and 50%. It did not decrease linearly while F10.7 decreased. In the African region, the maximum 70% in the early 2002 decreased to 40% in 2003. It did not decrease after 2003 when F10.7 was below 140. In the Eastern Pacific region, the occurrence rate was between 20% and 40% from 2000 to 2005 except for early 2002. It decreased to be 10% in 2006. The occurrence rate in the Central Pacific region were 30% at most, which was smaller than the other regions, and decreased as the solar active decreased.
 The low occurrence rates in the Central Pacific region could be explained by the location of the GPS station in this region. A proxy of latitude of GPS stations, HODE was listed in Table 2. The averaged HODE in the Central Pacific region was 983 km, while that of the other regions were low than 700 km. The effect of HODE on occurrence rate is studied in the next section.
3.3. Latitudinal Variation
 In addition to temporal variations, it is necessary to study the latitudinal effect on the occurrence rate of plasma bubble. Occurrence rates of plasma bubble are expected to be higher around the geomagnetic equator than the off-equatorial region.
 Latitudinal dependence of plasma bubble was investigated using meridional array of receivers in the Southeast Asia region. The GPS receivers used for this study were Phuket (PHKT; HODE = 400 km), Sampali (SAMP; HODE = 454 km), Nanyang Technological University (NTUS; HODE = 515 km), BUKT(HODE = 580 km), Chiangmai (CHMI; HODE = 662 km), Padang(PDNG; HODE = 612 km), and Bakosurtanal (BAKO; HODE = 952 km). The locations are shown in Figure 1 with filled circles. Figure 7 shows occurrence rates of these stations in 2004. Occurrence rates were plotted as a function of HODE, whose scale is on the horizontal axis. Names of stations are shown on the top side of the figure. Filled circles and star marks indicate occurrence rate averaged for 2 months; March–April and September–October, respectively. Open square represents occurrence rate for the whole year of 2004.
 Occurrence rates were almost constant below 600 km of HODE. They were 30–35% and 35–45% in March–April and September–October, respectively. Occurrence rates were smaller above 600 km of HODE than that below 600 km of HODE. They were lower by one order at BAKO whose HODE is 952 km than those of the other stations.
4.1. Solstice Asymmetry and Equinox Asymmetry
 In order to illustrate the monthly variation of plasma bubble, monthly occurrence rates in 2003 and 2004 were presented in Figure 8. During these 2 years, year-to-year variations were smaller than in any other years (see section 3). The receivers used were the same as those used in Figure 5 and listed in Table 1. The histograms showed the averaged occurrence rates in 2003 and 2004. Two dots on each month were the occurrence rates in 2003 and 2004. The dotted curves showed the sunset time lag between the geomagnetic conjugate points (see Figure 3).
 Some of previous studies have interpreted that monthly variation of plasma bubble is caused by the seasonal change of the sunset time lag between the geomagnetically conjugate E-regions [e.g., Tsunoda, 1985]. Some characteristics of the monthly occurrence shown in Figure 8 agree with the sunset time lag scenario while several characteristics do not agree. Major discrepancies were two asymmetries: (1) asymmetry between solstices and (2) asymmetry between equinoxes.
 Asymmetry between solstices was clearly seen in the African region where occurrence rates in the June solstice season were larger than those in the December solstice season. This feature was also seen in the Asian region from 2000 to 2002 when the solar activity was high. In the Eastern Pacific and Atlantic regions, occurrence rates in the December solstice season were larger than those in the June one. The difference in occurrence rates in the Atlantic regions could be explained by the sunset time lag scenario because the sunset time lag was larger in the June solstice than in the December one. In the Eastern Pacific region, however, the higher occurrence rate in the June solstice than the December could not be explained by the sunset time lag scenario.
 Seasonal-longitudinal variation of plasma bubble has been studied with scintillation data between ground-based receivers and radio beacon satellites. Aarons  compiled scintillation data and studied the longitudinal variation of plasma bubble occurrence. Figure 6a in the work of Aarons  shows occurrence rates of scintillation with 137 MHz radio signal from 1969 to 1970 at Huancayo in Peru, which is located at (12°S, 75°E). According to the figure, the occurrence rate was more than 80% in December while it was 0% in June. Occurrence of scintillation at Natal in Brazil at 6°S and 35°E was shown in Figure 8a of Aarons . The occurrence of scintillation with 254 MHz was more than 60% in December while it was 0% in June. Higher occurrence rate in June than in December was found at Huancayo and Natal, which are in the Atlantic and Eastern Pacific regions, respectively. These results are consistent with the plasma bubble occurrence rate derived with GPS-TEC data.
 Higher occurrence rate in June than in December in the Asian region also were reported in the previous scintillation studies [Paulson, 1981; Aarons, 1993]. Paulson  studied scintillation occurrence with 250 MHz radio wave at Guam(14°S, 143°E) in the Asian region. The occurrence rate was more than 10% in June while it was 0% in December during the solar maximum periods. Scintillation occurrence was investigated in Manila (14°N 121°E), in the Asian region by Aarons . He used 250 MHz radio wave from satellites and found that occurrence rate was about 20% and 0% in June and December in 1980, respectively.
 On the other hand, the asymmetry between two solstices in the African region which we found with the GPS-TEC observation was not seen in the scintillation observation from 1968 to 1970 at Accra in Ghana (5°S, 0°E). Koster  investigated satellite scintillations at Accra and found that occurrence rate derived from beacon satellite scintillation at Accra was smaller in June than that in December. This discrepancy could be attributed by the difference of measurement techniques and/or observational periods and/or the locations.
 An explanation for asymmetry between solstices was not given to the solstice asymmetry. One of physical parameters candidate that makes the asymmetry is seasonal variation of the flux tube integrated conductivities in the F-regions, ΣFP. Plasma bubble is widely believed to be generated by the growth of electron density fluctuations by the Rayleigh-Taylor (R-T) instability. The growth of the generalized R-T instability is controlled by the flux tube integrated conductivities in the E- and F-regions, ΣEP and ΣFP, and other parameters. The growth rate of the generalized R-T instability, γ, given by Zalesak et al.  was
where VP is the vertically upward component of plasma drift by the E × B/B2 drift due to the eastward electric field, UnP is the vertically downward component of the neutral wind velocity g represents the gravity acceleration, νin is the ion-neutral collision frequency, L is the scale length of the vertical gradient of the flux-tube integrated plasma density in the F-region on the magnetic equator, and R is the recombination rate that is integrated along the flux tube.
 Favorable conditions for the R-T instability are established soon after the sunset. The conductivity in the F-region is higher than in the E-region after the sunset. Under this condition, the F-region dynamo due to the eastward neutral wind controls the electro dynamics of the ionosphere. As a result, the eastward electric field enhances around the sunset terminator. The enhanced eastward electric field makes the E × B drift speed large. This enhancement of the eastward electric field, that is, prereversal enhancement, makes the collisions frequency, νin, and the recombination rate, R, small [Kelly, 1989]. It makes the growth rate large. The effect of prereversal enhancement is important. The F-region dynamo is high when the field line integrated Pedersen conductivity in the F-region is large. Solid curve in Figure 8 represents ΣFP in each region with the scale on the right-hand side axis. The conductivity of each point was derived from the conductivity model (http://swdcwww.kugi.kyoto-u.ac.jp/ionocond/index-j.html). The IRI 90 and the CIRA 72 model were used for ionized and neutral atmosphere models, respectively. The conductivity was integrated along each geomagnetic field lines which were traced from the points at 300 km and 400 km altitudes of the dip equator with the IGRF model. The field lines integrated conductivity was averaged in each region. In the African and Asian regions, ΣFP had maximum in June and minimum in December. In these regions, occurrence rate of plasma bubble were higher in the June solstice than in the December one. On the other hand, in the Eastern Pacific and Atlantic regions, ΣFP has maximum in December and minimum in June. In the Eastern Pacific region, occurrence rate of plasma bubble was higher in the December solstice than in the June solstice. This asymmetry between ΣPF for two solstices in these regions could cause the seasonal variation of plasma bubble occurrence rates.
 The equinoctial asymmetry was clearly seen in the Central Pacific region. In the Central Pacific region, occurrence rate in SE was much larger then that in ME. Large occurrence rates in SE than ME were found in the African region from 2002 to 2005 and in the Atlantic region in 2004 and 2005. In the Asian and Eastern Pacific regions, in 2005 and 2006, the occurrence rates in ME were larger than in SE. The asymmetry in the Central Pacific region was also reported in the previous studies [Aarons, 1993; Burke et al., 2004]. Larger occurrence rate in SE than in ME in Atlantic region and larger occurrence rate in ME than in SE in the Asian and Eastern Pacific regions agreed also with the previous studies [Aarons, 1993; Burke et al., 2004]. In the African region, however, scintillation data in the work of Koster , shows larger occurrence rate in ME than in SE. The difference of observational techniques or periods might make the discrepancy between scintillation data at Accra and ground-based GPS data.
 ΣFP cannot explain the equinoctial asymmetry since ΣFP have little asymmetry between two equinoxes. One candidate that may cause this asymmetry is the neutral wind in the F-region at the sunset time. The neutral wind is a source of the F-region dynamo. Effect of neutral wind in the F-region is not revealed in detail since measurement of neutral wind is difficult. Further investigation is necessary.
4.2. Year-to-Year Variation
 Characteristics of semiannual occurrence rates of plasma bubble were different among the five regions as shown in Figure 6. The region where the occurrence rate decreased the most was the Asian region, whose occurrence rate decreased from 70% in 2002 to 10% in 2006. Occurrence rate decreased in every region other than the Atlantic region. It decreases from 70% to 30%, from 50% to 10%, and from 30% to 10% in the African, Eastern Pacific, and Central Pacific regions, respectively. In the Atlantic region, it was between 20% and 40% and did not decrease with decreasing the solar flux.
 Previous studies suggested that yearly occurrence rate of plasma bubble tended to correlate with yearly averaged F10.7 [Huang et al., 2002; Sahai et al., 2000]. In the work of Huang et al. , probability of plasma bubble occurrence detected by DMSP satellite, PDMSP, was derived as a function of F10.7 in each longitudinal bin. Using data on Table 2 of Huang et al. , PDMSP in the Asian, Central Pacific, and Atlantic regions are estimated as PDMSP = −8.8 + 0.12F10.7, PDMSP = −6.8 + 0.1F10.7, and PDMSP = −12.6 + 0.22F10.7, respectively. Using F10.7 = 220 in 2001 and F10.7 = 70 in 2006, PDMSP was estimated for 2001 and 2006 and shown in the second row of Table 3. The annual occurrence rates derived from GPS data, PGPS, shown in Figure 6, are shown in the first row. The years 2001 and 2006 were regarded as high solar activity (HSA) and low solar activity (LSA) periods, respectively. In every three region, occurrence rate measured with GPS data was larger than PDMSP. The difference PGPS > PDMSP can be interpreted to be caused by the difference of observational altitude, and sensitivity of the measurement. Occurrence rate derived by the DMSP satellite observation was a ratio of a number of passes when plasma density was depleted to a number of all passes whose magnetic latitude were within ±30°. It was impossible to detect plasma bubbles lower than the orbital altitude, about 840 km. The ground-based GPS data can detect plasma bubble at any altitude above the GPS station. Therefore occurrence rates with GPS data tend to be larger than those derived from the satellite observation.
Table 3. Occurrence Probability of Plasma Bubble Using GPS-ROTI Data (PGPS), DMSP Data (PDMSP), and Airglow Data (PAG)a
 Another study on occurrence rate of plasma bubble and its dependence on the solar cycle was developed using airglow images [Sahai et al., 2000; Makela et al., 2004; Kim et al., 2002]. Plasma bubble was detected with all-sky OI 630nm images in Cachoeira Paulista (2.7°S, 45.0°W) in the work of Sahai et al. . Their definition of HSA and LSA were 1988–1991 and 1994–1998, respectively. Figure 2 of Sahai et al.  gives the monthly occurrence rates in HSA and in LSA in the Atlantic region. The annual occurrence rates, PAG, were 26% and 19% for HSA and LSA, respectively, which are shown in the third row of Table 3.
 Airglow observations were also conducted at Haleakala, Hawai [Makela et al., 2004; Kim et al., 2002]. Makela et al.  derived monthly probabilities of nights with plasma bubble for 2002 and 2003 using Cornell All-Sky Imager (CASI). Annual occurrence rates were derived from Figure 9 of Makela et al.  and shown in the third row of Table 3b as HSA data. Analysis of data for LSA was done by Kim et al. . Monthly occurrence rates were presented on Table 1 of Makela et al. . As PAG for LSA, the occurrence rates of 1966 and 1967 were averaged. It was 17% as shown on the third column on Table 3 for LSA. In the Central Pacific region, PGPS > PAG > PDMSP was found in the both HSA and LSA. In the Atlantic region, on the other hand, PDMSP ≃ PGPS > PAG was found in HSA while it turned to be PGPS > PAG > PDMSP in LSA. In the both region, PGPS tended to be closer to PAG than to PDMSP. PDMSP was smaller than PGPS and PAG especially in LSA. PGPS in LSA was more than one fourth of PGPS in HSA in every region. On the other hand, PDMSP decreased more than one order through HSA to LSA. Huang et al.  discussed that the growth rate for bottom side irregularities was fairly constant over the solar cycle while for the topside it was much larger in the solar maximum than in the solar minimum. Therefore, occurrence rates derived by the DMSP satellite, whose observation altitude was on the topside could be more sensitive to the solar activity than those derived by the airglow data whose observation altitude was on the bottom side. It can explain the larger dependence of PDMSP on F10.7 than PGPS.
 The difference of apex altitudes on the solar activity agrees with study of solar cycle dependence of plasma bubble with a Jicamarca unattended long term investigations of the ionosphere and atmosphere (JULIA) radar [Hysell and Burcham, 2002]. JULIA, which is a 50 MHz coherent scatter radar located at 12°S and 77°W in the geographic coordinate, was used to observe 3-m-scale field-aligned irregularities (FAI) of plasma bubble. The altitude that JULIA observed FAI was from 95 km to 900 km on the dip equator. Statistical study with JULIA FAI data from 1996 to 2000 exhibited that the altitude that plasma bubble appeared were different among solar activity under geomagnetically quiet conditions. According to Figure 2 in their paper, apex altitudes of plasma bubbles were less than 600 km, around 700 km, and higher than 800 km for low, moderate, and high solar activity conditions, respectively. Low, moderate, and high solar activity conditions were defined as F10.7 ≤ 80, 80 < F10.7 ≤ 180, and 180 < F10.7, respectively. This difference of apex altitudes on the solar activity was interpreted understood climatological changes in the behavior of the electric field [Hysell and Burcham, 2002; Fejer et al., 1999].
4.3. Latitudinal Variation
 Occurrence rates were plotted against HODE of the GPS station in Figure 7 in order to study the latitudinal variation. The occurrence rates were almost constant where HODE was between 400 km and 700 km. Above 700 km of HODE, occurrence rates decreased. According to studies of in situ measurements by satellite in the solar maximum, the occurrence rate was constant where the apex height was less than 900 km. It decreased from 900 km of the apex height and become almost zero at 1400 km [Singh et al., 1997a, 1997b; Watanabe and Oya, 1986]. Huang et al.  reported that the apex height where plasma bubble reaches depended on the solar activity. Dependence of the apex altitude on the solar activity was also observed by JULIA VHF radar [Hysell and Burcham, 2002]. The altitude of 700 km, where occurrence rates of plasma bubble derived by GPS data began to decrease was lower than the apex height where occurrence rate derived with satellite and radar data began to decrease. One of the reasons why the altitude dependence was different could be the solar activity effect on the altitude where plasma bubble reaches, though it is not understood well. Another reason is a difference between the field of view of a GPS station and in situ measurements. Since the elevation mask is set to 45° (see section 2), the field of view of a GPS station is a circle with 400 km of radius over the station. This field of view corresponds to ±150 km of HODE. That is, HODE for one station has ambiguity of ±150 km. Because the field of view includes the area where HODE is larger than that of the station, the occurrence rate can be smaller than that of in situ measurement.
 Making use of advantage of ground-based GPS receivers data, occurrence characteristics of plasma bubble were studied. Monthly, semiannual occurrence rates for all longitudinal regions were derived. Characteristics of the monthly occurrence rates were different among the regions. Although it is believed that sunset time lag effect plays an important role for the monthly variation, there were some characteristics that could not be explained with the sunset time lag scenario. One of them was the difference of occurrence rates between two solstice seasons. In the Africa and Asian regions, occurrence rates in the June solstice season were larger than those of December. In the Eastern Pacific and Atlantic region, occurrence rates in the December solstice season were larger than those of June. It was found that ΣFP could explain the solstice asymmetry. In addition to the solstice asymmetry, asymmetry between two equinoxes was found. In the Central Pacific region, the occurrence rate in the September equinox was much larger than that of the March equinox. Not only the Central Pacific region but other regions had the equinoxes asymmetry. Semiannual occurrence rates from 2000 to 2006 showed year-to-year variation. The dependence on the solar activity was different among the regions. The occurrence rates in any regions except for the Atlantic region had a good correlation with F10.7. In the Asian region, where the best correlation was found, the occurrence rates decreased to 30% from 70%. Occurrence rate against latitude was investigated in the Asian region in 2004. It began to decrease when the HODE was higher than 700 km and was almost zero where the HODE was higher than 900 km.
 There have been many studies on occurrence of plasma bubble. The occurrence rates were different among the data sets. One of data sets for statistical occurrence study is data of in situ observation, which covers all longitudes. However, observations with in situ measurements were largely affected by the satellite orbit and altitude. Although ground-based observation has a limit in the location, ground-based GPS receivers data has a lot of advantage because the number of receiver is large. In this paper, plasma bubble occurrence was studied with GPS data for the first time. This will lead to an unprecedented view of occurrence characteristic of plasma bubble, which will help to understand the process of plasma bubble generation.
 We thank International GNSS Service, Japan Agency for Marine-Earth Science and Technology, Scripps Orbit Permanent Array Center, and Solar-Terrestrial Environment Laboratory, Nagoya University for providing the GPS receiver data.
 Amitava Bhattacharjee thanks J. Vincent Eccles and Archana Bhattacharyya for their assistance in evaluating this paper.