The Changing Eigenfrequency Continuum During Geomagnetic Storms: Implications for Plasma Mass Dynamics and ULF Wave Coupling

Geomagnetic storms are one of the most energetic space weather phenomena. Previous studies have shown that the eigenfrequencies of ultralow frequency (ULF) waves on closed magnetic field lines in the inner magnetosphere decrease during storm times. This change suggests either a reduction in the magnetic field strength and/or an increase in its plasma mass density distribution. We investigate the changes in local eigenfrequencies by applying a superposed multiple‐epoch analysis to cross‐phase spectra from 132 geomagnetic storms. Six ground magnetometer pairs are used to investigate variations from approximately 3 < L < 7 and across the whole dayside sector. We find that at L > 4, the eigenfrequencies decrease by as much as 50% relative to their quiet time values. Both a decrease in magnetic field strength and an increase in plasma mass density, in some locations by more than a factor of 2, are responsible for this reduction. The enhancement of the ring current and an increase in oxygen ion density could explain these observations. At L < 4, the eigenfrequencies increase due to the decrease in plasma mass density caused by plasmaspheric erosion.


Introduction
Geomagnetic storms are prolonged periods of intense solar wind-magnetosphere coupling that compress the Earth's magnetosphere, energize the magnetosphere-ionosphere system, and result in strong enhancements of the ring current, among other phenomena. These have been documented since the nineteenth century (Stewart, 1861). Coronal mass ejections and corotating interaction regions are the most common drivers of geomagnetic storms. The most intense coupling happens if the interplanetary magnetic field (IMF) is orientated southward and results in enhanced plasma convection within the magnetosphere (Akasofu et al., 1963;Gonzalez et al., 1994). During these disturbances, large quantities of energy are released from the Earth's magnetotail, often causing increased substorm activity (e.g., Hutchinson et al., 2011) and enhancing the Earth's radiation belts (e.g., Murphy et al., 2016).
The presence of a geomagnetic storm is usually monitored using ring current indices, most commonly the Dst https://wdc.kugi.kyoto-u.ac.jp/dstdir/dst2/onDstindex.html or the Sym-H (Iyemori, 1990) index. Ring current indices are determined from the perturbation of the horizontal magnetic field strength using a set of ground-based magnetometers located at midlatitudes. These locations map to inside the ring current region of space where they largely avoid other current systems. A storm usually has three observable phases in the ring current indices (see Figure 1). During the initial phase, the ring current index increases due to the enhancements in the magnetopause currents as the magnetosphere is compressed, which increases the magnetic field strength on the ground. During the main phase of the storm, the ring current intensifies, The red section is the initial phase of the storm, the blue section is the main phase, and the green section is the recovery phase. An additional 12 hr of data has been shown in gray at either end of the storm interval.
which results in a sharp drop in the ring current index due to the weakened magnetic field. The recovery phase usually lasts several days as the ring current decays. Hence, a sharp decrease in the Sym-H index is indicative of a geomagnetic storm (e.g., Akasofu et al., 1963;Chapman & Ferraro, 1930).
Magnetic flux tubes, or field lines, can support standing Alfvén waves. The eigenfrequencies on any given flux tube depend upon the flux tube length, wave polarization, and the Alfvén speed v A , which in turn depends on the magnetic field strength, B, and plasma mass density, , of the flux tube (Cummings et al., 1969). The dependencies on v A are given by where 0 is the permeability of free space. The eigenfrequencies of any given flux tube thus depend upon the relationship between these three parameters: Eigenfrequencies will increase if the magnetic field strength increases, the plasma mass density decreases, or the flux tube decreases in length.

Eigenfrequencies During Storm Times
Several studies have observed that the eigenfrequencies of closed magnetic field lines decrease as the ring current index decreases (becomes more negative). Wild et al. (2005) predicted this by using the time of flight technique to calculate eigenfrequencies in a T96 magnetic field model (Tsyganenko, 1996) parameterized by the Dst index.  used the same technique but with the more realistic plasma mass density model of Sandhu et al. (2017) based on 12 years of Cluster data (L > 5.9) that could be parameterized by the Dst index and found the same result. Most recently, Wharton, Wright, Yeoman, James, et al. (2019) performed a statistical analysis of eigenfrequencies measured by the cross-phase technique applied to ground magnetometers and found that the eigenfrequencies decreased when the Sym-H index was reduced for L > 4.
Other studies have investigated single geomagnetic storm events. Kale et al. (2009), Takasaki, Kawano, et al. (2006), and Chi et al. (2005) all looked at the Halloween storm event in October 2003, and each found that the eigenfrequency decreased after a sudden storm commencement. Du et al. (2005) also found an eigenfrequency decrease for a storm in July 2000, Lee et al. (2007) for the March 1991 storm, andRae et al. (2019) for the March 2013 storm. Warner and Orr (1979), Engebretson and Cahill Jr. (1981), and Takahashi et al. (2002) are among others to have observed eigenfrequencies decrease with increasing geomagnetic activity at the latitudes investigated in this study.
10.1029/2019JA027648 number density decreased. A reduced plasma mass density would increase the eigenfrequencies. However,  used this model and the T96 magnetic field model to show that the eigenfrequencies decreased for low Dst values. This was because the reduction in the magnetic field strength created by an enhanced ring current was dominant in determining the eigenfrequencies over the reduced plasma mass density. Several other studies have shown that the plasma mass density decreases during enhanced geomagnetic activity (e.g., Corpo et al., 2018;Dent et al., 2006;Denton et al., 2006;Menk et al., 2014).
However, it has also been shown that the plasma mass density increases in response to increased geomagnetic activity (Takahashi et al., 2002(Takahashi et al., , 2010Takahashi, Denton, et al., 2006), and Min et al. (2013) showed that the plasma mass density barely changes. Chi et al. (2005), Takasaki, Kawano, et al. (2006), and Kale et al. (2009) also determined that the plasma mass density increased after the sudden storm commencement. Equation 1 shows that if the plasma mass density increases, then changes in both the plasma mass density and magnetic field strength are acting in cooperation to decrease the eigenfrequencies by reducing the Alfvén speed. In this scenario, just considering the change in magnetic field strength would not fully account for the decrease in eigenfrequency.

The Objectives of This Study
Studies based on correlating eigenfrequencies with a ring current index (e.g., Sandhu et al., 2017;Wharton, Wright, Yeoman, James, et al., 2019) do not fully capture the eigenfrequency evolution throughout a geomagnetic storm's three phases. Eigenfrequency data taken during the main and recovery phases will be combined when binning by only the value of a ring current index. For example, a bin containing all samples taken when the Dst index was −80 nT would include samples corresponding to both the main and recovery phases of a storm. It is important to consider the phases separately as each storm phase is associated with different magnetospheric processes. The statistics will also be dominated by data in the recovery phase, as this is typically the longest phase. Hence, the results will not typically describe what is happening in the main phase of the storm.
Other studies have looked at how the eigenfrequencies vary for specific storms (e.g., Lee et al., 2007;Rae et al., 2019) to explore the time evolution. However, these studies do not provide information on whether the storms are representative of the statistically average storm, especially as these studies focussed on more extreme events. Any storm related trends in the data will also be contaminated by diurnal effects. For example, if observations beginning in the morning observe an eigenfrequency decrease throughout the day as a storm is progressing, it is difficult to determine how much of this decrease is due to the storm and how much is due to the expected diurnal variation for a single event. The eigenfrequencies normally decrease when moving into the afternoon sector (e.g., Chappell et al., 1971;Chi et al., 2013;Sandhu, Yeoman, James, et al., 2018;Wharton, Wright, Yeoman, James, et al., 2019). Data from multiple storms are required to disentangle these two factors influencing the change in eigenfrequency.
A statistical study of a large number of storms enables us to determine the average trends in eigenfrequencies during storms and to reduce the diurnal variation. One hundred and thirty-two storms are identified using the algorithm of Walach and Grocott (2019). The cross-phase technique of Wharton et al. (2018), and Wharton, Wright, Yeoman, and Reimer (2019) is used to measure the eigenfrequencies, and the superposed multiple-epoch analysis method employed by Hutchinson et al. (2011) is used to find the average variation for all the storms. By multiple epoch, we mean that we have analyzed the initial, main, and recovery phases of the storm independently. This is the first reported study to employ a superposed epoch analysis to investigate the eigenfrequencies in cross-phase spectra. The superposed epoch analysis also puts all of the storms onto a common time grid so we can see how the average eigenfrequency varies relative to the storm phase and it prevents data from different storm phases being mixed together. By only including data in the analysis from a given magnetic local time (MLT) sector, the diurnal variation can be removed from the average. Different MLT sectors can then be analyzed separately. A chain of magnetometer stations allows us to investigate changes in eigenfrequency for different latitudes as well as MLTs.
Previous studies have provided a valuable understanding of how mass density changes. However, due to varying instrument coverage between studies, differing techniques, and limitations of the methodology as discussed earlier, the typical storm time variations are inconclusive with reports of both enhanced mass density and depletions of mass density during geomagnetic storms. We use our data to determine whether the plasma mass density increases or decreases for the average geomagnetic storm by solving the wave equation of Singer et al. (1981) in conjunction with the TS05 magnetic field model (Tsyganenko & Sitnov, 2005). It is found that whether the plasma mass density increases or decreases depends on the L-shell observed.

Geomagnetic Storm Intervals
To study the average changes in eigenfrequencies during geomagnetic storms, we have used a storm list created using the method of Walach and Grocott (2019). Their algorithm provides times that mark the four boundaries associated with the three phases mentioned in Section 1: initial, main, and recovery.
This algorithm searches the Sym-H index for minima less than −80 nT. This is the same criteria employed by Hutchinson et al. (2011). Some weak storms may have a Sym-H minimum greater than −80 nT, but it was deemed best to remove these to ensure that the minimum was due to a geomagnetic storm. This minimum defines the time at the end of the main phase. An example of this procedure is illustrated in Figure 1. The beginning of the main phase is defined as the most recent time before the Sym-H minimum crossed the quiet time level of −15 nT. Walach and Grocott (2019) found that this was a more reliable definition than using the preceeding maximum Sym-H value. To define the beginning of the initial phase, the maximum value of Sym-H before the beginning of the main phase is identified. The most recent time before this when Sym-H is −15 nT is the beginning of the initial phase. Note that this definition of the start of the initial phase is different to that of Hutchinson et al. (2011) and Murphy et al. (2018), who use a definition based on dynamic pressure enhancements. The end of the recovery phase is defined as the next time after the end of the main phase when the Sym-H index returns to −15 nT. This algorithm is also discussed in Walach and Grocott (2019), where the list is used to investigate geomagnetic storms with the Super Dual Auroral Radar Network (Chisham et al., 2007;Greenwald et al., 1995) from 2010 to 2016.
This algorithm was applied to Sym-H data from 2002 to the end of 2018. After visually checking for false detections, this gave a list of 132 storms. The dates and Sym-H minima associated with each storm are given in the supporting information.

Magnetometer Data
This study used 10-s resolution magnetometer data from the International Monitor for Auroral and Geomagnetic Effects (IMAGE) array (Luhr, 1994). North-south component data from the stations listed in Table  1 were used in the analysis, which corresponds to the toroidal mode if we assume a 90 • rotation of the wave ellipse by the ionosphere (Hughes & Southwood, 1976). Midpoint dipolar L-shell values are given to help visualize the latitudinal coverage of the pairs used in this study.

Cross-Phase Superposed Epoch Analysis
To investigate how the eigenfrequencies of magnetic field lines change during geomagnetic storms, a different approach was used to that of Wharton, Wright, Yeoman, James, et al. (2019). They simply binned measurements according to their Sym-H value. However, the Sym-H index only reaches storm levels, here defined as Sym-H < −80 nT, for a very small fraction of the total time available. Data from both the main and recovery phases will also be combined together, and no information about how the eigenfrequencies The eigenfrequencies were determined using the cross-phase technique of Waters et al. (1991), which requires two latitudinally and closely spaced ground-based magnetometers. Two are required to detect the phase change with latitude that occurs at the resonant frequency of the midpoint of the magnetometers. Several papers have automated this technique (Berube et al., 2003;Chi et al., 2013;Sandhu, Yeoman, James, et al., 2018;Wharton et al., 2018).
Further to this, Wharton, Wright, Yeoman, and Reimer (2019) developed a Lomb-Scargle (LS) cross-phase technique that could process unevenly spaced data and use a higher frequency resolution because the frequency grid is independent of the properties of the data used. The chosen resolution was 4 times that achievable with a discrete fourier transform as this value was adopted in previous studies which used LS (e.g., Bland et al., 2014;Wharton, Wright, Yeoman, James, et al., 2019;Wharton, Wright, Yeoman, & Reimer, 2019). The dynamic cross-phase spectrum uses a 40-min sliding window, which gives a frequency resolution of 0.417 mHz without oversampling for the IMAGE magnetometer resolution. Changing the oversampling factor to 4 reduces the frequency resolution to 0.104 mHz, allowing us to be more precise in determining the frequencies of waves. The LS technique is not intrinsically more accurate than using the discrete Fourier transform, but it does allow us to create a denser frequency axis. This was the motivation for using it in this paper.
The superposed mulitple-epoch analysis method used by Hutchinson et al. (2011) was then applied to the derived cross-phase spectra. This method treats the three phases of geomagnetic storms separately by calculating the mean duration of each of the three storm phases (initial, main, and recovery). A superposed epoch analysis is then applied to each storm phase. This method is better than applying a superposed epoch analysis to the whole storm interval because the durations of each phase are not always in the same proportion to each other for different storms. Treating the three phases independently accounts for this. These mean durations were 33.16, 11.76, and 53.10 hr for the initial, main, and recovery phases, respectively. This enabled us to create a common time grid to which the three phases of each storm could be normalized to so we could observe average changes in each of the three storm phases.
For each storm, the LS cross-phase spectrum was calculated. Figure 2 shows an example of how the rescaling was performed on one of the storm intervals. Figure 2a contains the filtered magnetometer data (800-s boxcar filter (Wharton et al., 2018)) for two magnetometer stations. Figure 2b shows the LS cross-phase spectrum for this storm. The initial, main, and recovery phases have been shaded red, blue, and green, respectively. An additional 12 hr of data has been included on each end, shaded in gray, to represent the quiet time periods and will be used to compare variations to quiet time conditions.
Before the rescaling and averaging process could be performed, there is a complication concerned with MLT that had to be considered. The full duration of a geomagnetic storm is usually a few days. This means that the diurnal variation of the eigenfrequencies must be taken into account. For example, eigenfrequencies are always higher in the morning sector than the afternoon (e.g., Takahashi et al., 2016;Wharton et al., 2018). This is because the plasma mass density is higher in the afternoon sector relative to the morning due to flux tubes refilling from the ionosphere throughout the day. If this is not accounted for, changes in the eigenfrequencies that are diurnal and storm related will be mixed together. To compensate for this, a MLT sector is first defined (e.g., 10-14). We arbitrarily chose 4 hr as a compromise between minimizing diurnal variability and not excluding too much data. Only data taken when the magnetometer pair were within that MLT sector were included in the superposed epoch analysis described above. The analysis was then performed independently for other MLT sectors. This minimized contamination of the spectra by diurnal variation so that only the storm-related variations were observed.
To perform the superposed epoch analysis, each frequency in the periodogram was treated as its own time series. Each "frequency time series" was interpolated and then rescaled to the common time axis for each phase of the storm. This gave the rescaled LS cross-phase spectrum in Figure 2c. This process was repeated for all storms so that the spectra could then be superposed and averaged using the common time grid. In this example, data from all 24 hr of MLT have been included to more clearly show the rescaling process. For the forthcoming analysis, only data taken from 4-hr MLT sectors are included. The equivalent plot here would contain considerable white space due to the excluded data.
The superposed multiple-epoch analysis method used by Hutchinson et al. (2011) was also applied to variables from the OMNI data set, which describes the upstream solar wind and contains geomagnetic indices. This is shown for the Sym-H index in Figure 3a, where the three phases of the storm have the same color (red, blue, and green). The mean, median, and upper and lower quartiles of the Sym-H index have been overplotted. Figure 3 gives an example of the superposed epoch analysis applied to the 132 storms for the Pello-Muonio magnetometer pair in the 10-14 MLT sector. Figure 3b shows the number of storms that contributed data to each time on the common time axis. This number varies because the storms were not equally distributed in MLT. On average, 22 storms were expected at any given time due to the size of the MLT sector. Figure 3c shows the results of the superposed epoch analysis on the LS cross-phase spectra taken in the 10-14 MLT sector for all 132 storms. The LS cross-phase spectra for all of the storms have been circularly averaged in each frequency-common time bin because the cross phase is a circularly distributed variable (Mardia & Jupp, 2000). A variation in the fundamental eigenfrequency is visible. The trend is least clear in the initial phase but begins at ∼8 mHz at the beginning of the phase and reaches ∼4 mHz at the start of the main phase. Other bands can also be seen at higher frequencies, which complicates the interpretation of the harmonic numbers. The eigenfrequency reaches a minimum at ∼4 mHz in the main phase, where it remains steady before increasing again throughout the recovery phase back to 8 mHz. Only the fundamental was observable in the recovery phase. Red shading indicates initial phase, blue is main phase, and green is recovery phase. Gray is 12 hr before and after the storm.
Identifying the fundamental eigenfrequency in the initial phase is ambiguous and complicated by the possible presence of a higher harmonic. This could be the third harmonic, based on the arguments in Wharton et al. (2018), and it also decreases in value toward the main phase from ∼15 to ∼5 mHz. It could also be the fundamental eigenfrequency increasing around 24 hr then decreasing again. The selection algorithm in Section 4.1 tries to address this ambiguity with a manual process.
The same analysis was performed for the other station pairs at the latitudes given in Table 1. Although there were some differences in whether the third or fundamental harmonic were being observed during the initial phase, a consistent pattern of decreasing eigenfrequencies was observed at all but the lowest latitudes. This decrease in eigenfrequency with increasingly negative Sym-H is in agreement with observations by  and Wharton, Wright, Yeoman, James,, et al. (2019). Other latitudes and MLTs are presented in the next section, and the results will be summarized in Figures 6 and 7, which separate out the storm time and diurnal variation.

Controlling Factors for the Eigenfrequency
Changes in the eigenfrequency values depend on the magnetic field structure and plasma mass density, as discussed in Section 1. Calculating the average changes in these parameters allows an understanding of why the eigenfrequencies vary throughout the storm cycle. To do this, a range of appropriate frequencies centered on the eigenfrequency must be extracted from the data and expressed as a function of the common time axis to enable calculation. This is difficult using the data in the format displayed in Figure 3. Section 4.1 describes this extraction process. In Section 4.2, we parameterize the TS05 magnetic field model (Tsyganenko & Sitnov, 2005) and calculate the plasma mass density. In Section 4.3, we repeat this for all station pairs in Table 1 and for a range of MLTs. This enables the derivation of eigenfrequency maps and, by accounting for magnetic field changes according to the TS05 model, estimations of equatorial plasma mass density. The results provide a fuller understanding of how the structure of the dayside magnetosphere changes during geomagnetic storms, at least within the region covered by the magnetometers in this study.

Isolating the Fundamental Eigenfrequency
The technique to extract the fundamental eigenfrequency profile is based on that of Wharton et al. (2018) with an additional step, which is illustrated in Figure 4. Figures 4a and 4b are reproduced directly from  (Figure 4c). To do this for any given value, it is averaged with values within ±2 mHz and ±2 hr, giving a range or box size of 4 mHz by 4 hr. At the edges of the spectrum, these ranges will be capped by the spectrum boundaries (0 or 50 mHz for frequency and the start and end times), shrinking these ranges to half their size for the boundary values. Hence, the standard error will be larger at the plot boundaries. The standard deviation is also calculated and stored for the purposes of the t test explained below.
Next, the mean and standard deviation of the cross phase for each common time are calculated. Only values less than the mean minus one standard deviation are retained. The result of this is illustrated in Figure 4d. It can be seen that a narrow frequency band where the fundamental frequency lies dominates the remaining data. Then a t test is applied to remove the nonresonant features. The statistic is the smoothed cross-phase value divided by its standard deviation and multiplied by minus 1. Only cross-phase values with a t statistic greater than one are kept. Nonresonant signatures randomly had a highly negative cross-phase value, but their large standard deviation can be used to remove them. Only the stable values then remain, which are presumed to represent real eigenfrequencies. This process is also discussed in Berube et al. (2003), the algorithm upon which the Wharton et al. (2018) algorithm was originally based. These are shown in Figure  4e.
We then use a process that searches through the spectrum and isolates the remaining independent features and assigns them a number. This process searches up each column in Figure 4d, beginning at the bottom left, until it finds a data point. It then searches this column and adjacent columns for data in adjacent cells and groups this data together. This group is a feature, and these data are then discounted from the rest of the search. Once all points that are grouped together are found, the search up the columns continues looking for the start of the next feature and so forth. The final list of identified features is labeled with the numbers in the yellow boxes in Figure 4e that are automatically attached to their respective features. "F" stands for feature and the number corresponds to the order in which they were identified. Hence, lower numbers are to the left where the search began. The user can then manually select which "features" they think contribute to the fundamental eigenfrequency. This is a subjective process to a certain extent. In some cases, the labels overlap because two features begin in very similar locations on the plot. In such circumstances, trial and error is unfortunately required to select the desired features. For the example in Figure 4e, the features 1, 2, 10, 12, 15, 21, 22, and 23 were selected by the user. These features are then fitted with a smooth line to yield an eigenfrequency profile, which is the thick red line in Figure 4e.
The band identified in Figure 4e is spread over a range of frequencies. An uncertainty on the eigenfrequency is found by calculating the mean frequency width of the selected features that composed the fundamental eigenfrequency. This range is illustrated by the two thin red lines in Figure 4e and represents an estimate in the uncertainty in measuring the fundamental frequency using this method.
This method was considered more reliable than employing a complicated automated algorithm. The primary issue is that there are often multiple harmonics, which could contaminate the fundamental, so automatic identification becomes difficult. This method was a way of resolving this problem. It is noted that advanced machine learning methods could be applied in future but such methods were not needed for the small number of cases in this study.

Determining the Plasma Mass Density
Once the fundamental frequency has been determined as a function of the common time as described in Section 4.1, the next step is to describe the magnetic field line. Once achieved, the magnetohydrodynamic (MHD) wave equation of Singer et al. (1981), Equation 2, can be solved to determine the plasma mass density.
We chose the TS05 magnetic field model (Tsyganenko & Sitnov, 2005) to model the magnetic field lines, which was designed to include the variability introduced by geomagnetic storms and specifically includes a ring current contribution. Huang et al. (2008) have shown that the performance of this model is better than previous empirical magnetic field models. The parameters for this model are the Dst index (we have used the Sym-H index, which can be considered as a high resolution Dst index (Wanliss & Showalter, 2006)), solar wind dynamic pressure and velocity, IMF x and y components and the dipole tilt angle.
To parameterize the TS05 model across the common time axis, the superposed multiple-epoch analysis method described in Section 3 was applied to the Sym-H, pressure, IMF B z , and solar wind velocity x component data. This process is illustrated in Figure 5 for the Pello-Muonio pair in the 10-14 MLT sector. The The mean values of these four parameters at each common time were used for parameterizing the TS05 model. The solar wind velocity y and z components, dipole tilt angle, and IMF B y data were all set to 0, as this was their expected average values. Over long time scales, the solar wind flow is expected to be radial and the dipole tilt angle will average to 0. Therefore, the IMF B y component should not affect the storm process either. This was confirmed from applying the superposed epoch analysis method to these variables to check. The superposed epoch analysis cross spectrum from Figure 3b is reproduced in Figure 5e with the eigenfrequency profile from Figure 4e plotted over it.
The TS05 model was evaluated for every tenth point available on the common time axis. Figure 5f shows the maximum radial distance of any point along the field line, assumed to represent the equatorial crossing 10.1029/2019JA027648 point. As the initial phase begins, the field line moves earthward until the main phase, where it then moves outward beyond its original position. The field lines recover toward their original position at the end of the recovery phase.
Plotted in Figure 5g is the magnetic field strength at the maximum radial distance of the field line. Note that in the equatorial plane, the magnetic field strength is primarily dominated by the B z component. This is observed to increase throughout the initial phase as the magnetopause currents are enhanced and the closed magnetic field line strength increases, then drops throughout the main phase as the enhanced ring current weakens the background magnetic field. It then returns to normal during the recovery phase. These results show that the TS05 model is capturing the expected aspects of a geomagnetic storm.
The MHD wave equation, derived by Singer et al. (1981), has been used to calculate the equatorial plasma mass density. The justification for using this equation is given below. Solving this equation requires knowledge of the magnetic field structure and the eigenfrequency. We have used the TS05 magnetic field model and the eigenfrequencies measured from our data. In Equation 2, is the field line displacement, h is a scaling factor describing the separation of the field lines along the field, s is the field-aligned coordinate, B 0 is the ambient magnetic field strength, 0 is the magnetic permeability of free space, is the angular eigenfrequency, and is the plasma mass density. The separation h is calculated in the toroidal direction using a closely spaced field line. Singer et al. (1981) state that Equation 2 can be applied to any arbitrary geometry. The derivation of this equation implicitly assumes that orthogonal field-aligned coordinates can be used. However, Salat and Tataronis (2000) have shown that such fields only exist in special cases, such as a dipolar field. In general, fields that can support field-aligned orthogonal coordinates do not exist and this includes the TS05 fields. Using orthogonal coordinates requires B · (∇ × B) = 0 or that there is no shear in the magnetic field. An absence of shear requires that the magnetic field is 0 in one of the three coordinate directions. In a dipole field, this would be the azimuthal direction. For a more realistic magnetospheric field line, there will be variation in the magnetic field in all three field-aligned coordinates due to shear, and thus, an orthogonal field-aligned coordinate system is not possible. Therefore, solving the wave equation in a nonorthogonal coordinate system such as that of Rankin et al. (2006) would be more accurate, despite the more complicated formulation of the wave equation. In this study though, we only consider magnetic field lines in the dayside magnetosphere away from the magnetopause, which can be considered to be quasi-dipolar. Hence, any error from using Equation 2 is expected to be minimal, and many other authors have used this equation when tackling this problem (e.g., Berube et al., 2006;Waters et al., 1996;Wharton et al., 2018). Therefore, we have applied this equation to our problem.
Equation 2 is solved using the Runge-Kutta method described in Wharton et al. (2018). The plasma mass density is assumed to have the power law form for its distribution along the field line given by This is a common assumption made by many authors (e.g., Angerami & Thomas, 1964;Carpenter & Smith, 1964;Cummings et al., 1969;Denton et al., 2006;Menk et al., 1999Menk et al., , 2004Orr & Matthew, 1971;Wharton et al., 2018). In Equation 3, is the plasma mass density and r is the radial position along the field line. The subscript eq represents the equatorial value. The exponent determines the rate at which the plasma mass density changes with radial distance. A positive exponent causes the plasma mass density to decrease radially, so that it reaches a minimum along a field line in the equatorial plane. An exponent of 3 has been used, and it is assumed this is a reasonable value based on other studies and has been used by many other authors (e.g., Poulter et al., 1984;Menk et al., 1999Menk et al., , 2004. This distribution contains a minimum at the equator and maxima at the base of the field lines in the ionosphere. The size of this exponent has been shown to make little difference to the calculated equatorial mass densities (e.g., Menk & Waters, 2013;Orr & Matthew, 1971).
The equatorial plasma mass density profiles calculated with the mean fundamental frequency and the upper and lower bounds are shown in Figure 5h. The lower bound eigenfrequency corresponds to the highest plasma density and vice versa. In this example, the plasma density remains stable throughout the initial phase, which could suggest that as the magnetosphere is compressed, that plasma loss processes increase. The plasma mass density then increases to its highest value at the beginning of the main phase. It then decreases back to quiet time levels throughout the recovery phase. Small differences in the eigenfrequency correspond to large differences in the calculated plasma mass densities, in agreement with previous studies (e.g., Takahashi & Denton, 2007).

Mapping the Plasma Mass Density
The process of estimating the equatorial mass density described in Section 4.2 was repeated for each station pair in Table 1 for three MLT sectors, 6-10, 10-14, and 14-18 MLT. Only the dayside was considered because the cross-phase technique is ineffective at night (e.g., Wharton, Wright, Yeoman, James, et al., 2019). There equatorial magnetic field data in the middle, and equatorial plasma mass density data on the right. Plots of the form of Figure 3 show that there are a sufficient number of storms in each MLT sector at all times to do this analysis.
For each MLT-radial bin, the quiet time average value of either the eigenfrequency, magnetic field strength, or plasma mass density is calculated from the 12-hr period before the initial phase and the 12-hr period after the recovery phase. These data are shaded in gray in Figures 1-5. This value is then used to normalize the respective data in that bin so we can see how the values change relative to their quiet time average. This approach made it much easier to discern the general trends as opposed to observing individual plots such as Figure 5. Each color scale is centered on 1. Hence, blue indicates a reduction in that value, and red an increase relative to quiet times.
An equivalent plot with absolute eigenfrequency, equatorial magnetic field, and plasma mass density values is shown in Figure 7 in the same format as Figure 6. This shows that the variation in eigenfrequency with L-shell and MLT is much greater than the storm time variation at any given position, making such variations difficult to discern. Hence, we focus our discussion on the normalized data in Figure 6. Here, an L value refers to the radial distance in the equatorial plane to the field line in Earth radii. Figure 6a shows the beginning of the initial phase, when values are near the quiet time average. Figure  6b shows the end of the initial phase, where the eigenfrequencies have decreased in the 14-18 MLT sector and the magnetic field strength has increased for the outer L-shells. There has also been an increase in the plasma mass density, mostly clearly in the afternoon sector. Figure 6s shows the end of the main phase. The eigenfrequencies have decreased by as much as 50% across most L-shells and MLTs. At the innermost L-shells, the eigenfrequencies are actually greater than the local quiet time average. The magnetic field has weakened, especially in the 6-10 and 14-18 MLT sectors. The plasma mass density has also increased across all sectors, except at the innermost L-shells where the opposite is observed. Figure 6d shows partway through the recovery phase. The eigenfrequencies have risen back toward mean levels again but at the innermost L-shells, they are still higher than during quiet times. The magnetic field is still weaker than during quiet times. The plasma mass density has decreased again and remains above average only at the outermost L-shells. It is still below average at the inner distances. Finally, Figure 6e shows the end of the recovery phase, where for the most part, all three parameters have returned back to their quiet time averages.

Discussion
Based upon this analysis, Figure 6 demonstrates that there is a clear pattern in the behavior of measured toroidal eigenfrequencies during an average geomagnetic storm. The Alfvén speed, and hence a given field line eigenfrequency, is proportional to the magnetic field strength and inversely proportional to the square root of the plasma mass density (Equation 1). Hence, understanding both the field line structure and the plasma mass density distribution is necessary to determine why any field line eigenfrequency changes with time.
Section 5.1 discusses the changes seen at L > 4 in the context of changing magnetic field strength and plasma mass densities. Section 5.2 then discusses changes observed for L < 4, which had different variations in eigenfrequency and plasma mass density. Figure 7 shows that during all phases of the storm, the eigenfrequencies are consistently greater in the morning sector than the afternoon sector, which is the same across all observed L-shells. With asymmetries in magnetic field configuration as described by the TS05 model removed, the asymmetry persists and consequently must be caused by plasma mass density asymmetry. This asymmetry in MLT is well documented (e.g., Poulter et al., 1984;Sandhu, Yeoman, James, et al., 2018;Takahashi & McPherron, 1984;Takahashi et al., 2016;Wharton et al., 2018;Wharton, Wright, Yeoman, James, et al., 2019), and we confirm that this asymmetry persists throughout all phases of a geomagnetic storm. The higher plasma mass densities in the afternoon lower the Alfvén speed in this sector. This extra plasma could be due to the existence of a plasmaspheric bulge, plasmaspheric plume, or increased refilling of ionospheric flux tubes. The average ion mass is also greater in the afternoon sector due to convection of heavy ions from the nightside plasma sheet (Sandhu et al., 2016).

Causes of the Variation in the Eigenfrequencies: L > 4
We now consider the storm time variations that occur for L > 4. Figure 6b shows the change that occurs during the initial phase. For L > 4, the eigenfrequencies first decrease in the afternoon sector, which corresponds with plasma mass density increases of approximately 100% in some cells. This change is also most prevalent at L < 7. At radial distances greater than this, the magnetic field has increased due to the enhanced magnetopause currents, which would increase the eigenfrequencies.
Changes during the main phase are shown in Figure 6c. Here, the eigenfrequencies have decreased by as much as 50% across all MLT sectors. This corresponds with a magnetic field decrease, most prevalent away from the noon sector where magnetospheric compression is strongest, and plasma mass density increases across all MLT sectors. The decrease in eigenfrequency during the main phase when the ring current index is low agrees with the observations by Chi et al. (2005), Takasaki, Kawano, et al. (2006), Kale et al. (2009), Sandhu, Yeoman, and, Wharton, Wright, Yeoman, James, et al. (2019), and Rae et al. (2019). Figures 6d and 6e show that the eigenfrequencies recover back to their quiet time levels as the recovery phase progresses. The return to quiet time values takes longest in the afternoon sector. This may be due to the existence of a lasting plasmaspheric plume.
The storm time variation in the equatorial magnetic field strength is as expected, due to the increase of the ring current magnetic field. However, this study demonstrates that the role of the plasma mass densities in determining the storm time evolution of the eigenfrequencies is greater than previously thought by  and Rae et al. (2019). Previous studies have differed about whether the plasma mass density increases or decreases during storms. Sandhu et al. (2017) found using an empirical model that the plasma mass density decreased during intervals of low Dst index, despite the average ion mass increasing. This was because the electron number density decreased due to increased magnetospheric convection. Several other studies have also put forward evidence that the plasma mass density decreases during enhanced geomagnetic activity (e.g., Dent et al., 2006;Denton et al., 2006;Corpo et al., 2018;Menk et al., 2014). However, there are other studies that have found the plasma mass density increases (e.g., Chi et al., 2005;Kale et al., 2009;Takahashi et al., 2002Takahashi et al., , 2010Takasaki, Kawano, et al., 2006;Takahashi, Denton, et al., 2006), which agrees with our findings. This increase was attributed to an increasing average ion mass (e.g., Kronberg et al., 2014). We show that at distances of L > 4, the plasma mass density doubles during the main phase of the storm across the whole dayside MLT sector and a wide range of L-shells. This shows that the plasma mass density and the magnetic field change in cooperation instead of opposition to lower the eigenfrequencies.
Where does this extra plasma mass originate from during the geomagnetic storm? The density of the hot O + population has been shown to increase during the main phase for L < 7 (Yue et al., 2019), and the partial pressure contribution from O + can be as great as that of H + (e.g., Nosé et al., 2011;Mitani et al., 2019). A rise in the density of O + in the ring current could explain our observations of enhanced plasma mass density with lower magnetic field. The cold O + population is also thought to increase due to higher outflows on the nightside during storms (e.g., Gkioulidou et al., 2019;Kistler et al., 2016). However, it is not possible to determine from our data if these populations are responsible for the plasma mass density changes, and if so, which populations dominate. Investigating these changes is beyond the scope of this study.
From a theoretical viewpoint, the lower eigenfrequencies during storm times would allow fast mode waves of given frequency to couple to field lines deeper in the magnetosphere, where they would produce standing Alfvén waves through the field line resonance mechanism (e.g., Southwood, 1974). Figure 7 also shows that these locations would be deeper in the afternoon than the morning due to the magnetospheric asymmetry about noon. Single event case studies such as Lee et al. (2007) and Rae et al. (2019) show that decreased eigenfrequencies during storm times coincide with enhanced ultralow frequency (ULF) wave power across large regions in the inner magnetosphere. Although this study shows that lower eigenfrequencies at L > 4 are a consistent feature of geomagnetic storms, it does not provide measurements of ULF wave power. The decreased eigenfrequencies are expected to contribute to the increased ULF wave power due to the power law power distribution of ULF waves (e.g., Rae et al., 2012), but other factors may also play important roles including the initial power of the fast mode waves, their frequency spectrum, their azimuthal wave numbers, or changing ionospheric boundary conditions. Establishing the importance of and isolating the effect of decreased eigenfrequencies on ULF wave propagation will be the focus of a future study.

Causes of the Variation in the Eigenfrequencies: L < 4
By including the Tartu-Nurmijärvi pair, changes within the plasmasphere have also been observed. Figure  6 shows that the eigenfrequencies and plasma mass densities had a different behavior on field lines where L < 4 than L > 4. Instead of decreasing, the eigenfrequencies increase, and the plasma mass densities decrease instead of increasing. The timing of this change is also delayed until the main phase, whereas changes for L > 4 begin during the initial phase. Relative to quiet time, the eigenfrequencies have increased by more than 50% partway through the recovery phase. The magnetic field change is much weaker for L < 4 compared to L > 4, suggesting that the reduced plasma mass density is almost entirely responsible for the increase in Journal of Geophysical Research: Space Physics 10.1029/2019JA027648 eigenfrequencies. It also suggests that the effect of the ring current magnetic field is weaker for these inner shells. Figure 6d shows a decrease in plasma mass density of more than 50%.
During a geomagnetic storm, the plasmapause boundary is eroded away by the increased convection of the Dungey cycle and a plasmaspheric plume forms (e.g., Chappell et al., 1971;Menk et al., 2014). Regions that previously were within the plasmapause are now outside and have a much lower plasma mass density. This would increase the eigenfrequencies (Dent et al., 2006;Kale et al., 2007;Menk et al., 2004). It is difficult to say where the plasmapause is in Figure 6 or Figure 7. Figure 8 shows the plasma mass density data from Figure 7 as radial profiles. Each column represents one of the three MLT sectors, and each row a point along the storm time axis shown at the bottom, the same as in Figures 6 and 7. The blue line shows the plasma mass density profile at that time, and the dashed red line is from the time/row before to help show how the plasma mass density has changed. If there is a consistent location for the plasmapause during all storms, it would be visible as a classic "knee" in this plot but because these data show the average values of 132 storms, there is likely to be considerable variation in the plasmapause position between different storms. This has the effect of blurring out the plasmapause position, and hence, the "knee" is not seen in Figure 8. However, we still expect the plasmapause to retreat in this region due to enhanced magnetospheric convection. In Figures  8c and 8d, for the innermost station pair, the plasma mass density decreases as expected from plasmapause erosion. The O'Brien and Moldwin (2003) plasmapause model shows that the plasmapause distance at noon will decrease to L < 3 for a −100 nT storm, which supports our observation of decreasing plasma mass density for L ∼ 4.
There may still be an increase in the average ion mass due to higher densities of O + or an increase in the mass density of the energetic population, as suggested in Section 5.1. However, the observations suggest that any increase in the average ion mass is inferior to the cold plasma loss of mass from plasmapause erosion observed for L < 4. Hence, we see the opposite behavior in the eigenfrequencies for the innermost L-shells.

Conclusions
This study determined how the eigenfrequency continuum changed during geomagnetic storms and how this was controlled by the underlying magnetic field strength and plasma mass density. We analyzed the eigenfrequencies of magnetic field lines during 132 storms using a superposed multiple-epoch analysis on cross-phase spectra derived from ground-based magnetometer stations in the range 3.15 ≤ L ≤ 6.42. This region corresponds to a key region during geomagnetic storms because of its proximity to the ring current and radiation belts. This has enabled us to determine the average time evolution of the fundamental eigenfrequency during each of the three phases of a geomagnetic storm. We then solved the Alfvén wave equation of Singer et al. (1981) using the TS05 magnetic field model of Tsyganenko and Sitnov (2005) to determine how the plasma mass density changes directly for the first time.
It was found that the fundamental eigenfrequency decreased in all dayside MLT sectors during the main phase for L > 4 and this corresponded to a weaker magnetic field and an increased plasma mass density, with a clear MLT dependence on the reduction. The changes in both of these dependencies would decrease the eigenfrequency. We suggest that the increase in plasma mass density is caused by an increase in Oxygen ion density.
For L < 4, the eigenfrequency increased relative to the quiet time value. This was accompanied by a decrease in the plasma mass density, which is thought to be caused by plasmapause erosion during the main phase of a storm. However, the magnetic field strength still decreased, suggesting that the plasma mass density has a stronger influence over the eigenfrequency than the magnetic field at these locations.

Data Availability Statement
The data are available at the website (http://space.fmi.fi/image/beta/). The OMNI solar wind and geomagnetic indices data are publicly available from the NASA Space Physics Data Facility, Goddard Space Flight Center (http://omniweb.gsfc.nasa.gov/ow.html). This research used the SPECTRE High Performance Com-