Variations of High‐Latitude Geomagnetic Pulsation Frequencies: A Comparison of Time‐of‐Flight Estimates and IMAGE Magnetometer Observations

The fundamental eigenfrequencies of standing Alfvén waves on closed geomagnetic field lines are estimated for the region spanning 5.9≤L < 9.5 over all MLT (Magnetic Local Time). The T96 magnetic field model and a realistic empirical plasma mass density model are employed using the time‐of‐flight approximation, refining previous calculations that assumed a relatively simplistic mass density model. An assessment of the implications of using different mass density models in the time‐of‐flight calculations is presented. The calculated frequencies exhibit dependences on field line footprint magnetic latitude and MLT, which are attributed to both magnetic field configuration and spatial variations in mass density. In order to assess the validity of the time‐of‐flight calculated frequencies, the estimates are compared to observations of FLR (Field Line Resonance) frequencies. Using IMAGE (International Monitor for Auroral Geomagnetic Effects) ground magnetometer observations obtained between 2001 and 2012, an automated FLR identification method is developed, based on the cross‐phase technique. The average FLR frequency is determined, including variations with footprint latitude and MLT, and compared to the time‐of‐flight analysis. The results show agreement in the latitudinal and local time dependences. Furthermore, with the use of the realistic mass density model in the time‐of‐flight calculations, closer agreement with the observed FLR frequencies is obtained. The study is limited by the latitudinal coverage of the IMAGE magnetometer array, and future work will aim to extend the ground magnetometer data used to include additional magnetometer arrays.


Introduction
Local excitation of a given field line by a source of wave energy at its eigenfrequency is known to result in the resonant oscillation of the field line (Dungey, 1954a(Dungey, , 1954b. These oscillations, referred to as FLRs (Field Line Resonances), are an important mechanism for the transfer of energy and momentum in the magnetosphere (Elkington et al., 1999;Rae et al., 2007). Therefore, it is of scientific interest to have an understanding of the expected frequencies of FLRs and their spatial variations.
The FLR process is the result of a standing Alfvén wave on a closed geomagnetic field line. Consequently, the fundamental eigenfrequency for a field line is determined by the length of the field line and the local Alfvén speed at all points along the field line. By definition, the local Alfvén speed is defined by the magnetic field strength and the plasma mass density. Since these values vary with radial distance, the fundamental eigenfrequency varies across the field lines, forming the Alfvén continuum. Energy input into the Alfvén continuum at a discrete frequency will drive an FLR at the field line where the driving frequency matches the field line eigenfrequency. Large-scale spatial variations in the frequency of standing Alfvén waves are thus attributed to the magnetic field configuration and the distribution of plasma mass density. Although work by Vellante and Förster (2006), for example, has provided a detailed and accurate examination of eigenfrequencies for field lines located well within the plasmasphere, the study presented here will focus on closed field lines with larger L values located in the outer plasmasphere, plasmatrough, and near-Earth plasma sheet. Previous studies have explored estimating fundamental FLR eigenfrequencies for closed geomagnetic field lines using the time-of-flight approximation (Lee & Lysak, 1990;Warner & Orr, 1979;Wild et al., 2005), 10.1002/2017JA024434 involving the assumption of a magnetic field model and plasma mass density model. Notably, a study conducted by Wild et al. (2005) presented time-of-flight calculated frequencies for high-latitude closed field lines, using the T96 magnetic field model (Tsyganenko, 1996) and a simple radial power law mass density model. The results demonstrated that the use of the realistic T96 magnetic field model provides more representative estimates of the eigenfrequencies, compared to previous calculations using more simplistic magnetic field models (Lee & Lysak, 1990;Warner & Orr, 1979). The results account for local time variations, field configuration at high latitudes, and dependences on geomagnetic activity and solar wind conditions. However, the mass density model assumed by Wild et al. (2005) is relatively simplistic. It has been shown by Sandhu et al. (2016a) that the magnetospheric mass density distribution contains strong variations with MLT (Magnetic Local Time) and cannot be described validly with a simple radial power law mass density model. Therefore, the time-of-flight analysis of fundamental FLR eigenfrequencies presented by Wild et al. (2005) can be further refined by including a more realistic plasma mass density model to account for the mass density contribution, in addition to the realistic T96 magnetic field model. This forms the motivation of this study.
In order to assess the validity of time-of-flight calculations, a comparison to observed FLR frequencies would prove useful. A brief overview of previous studies to examine average large-scale spatial variations in FLR frequency is now presented. Among the first studies to consider the spatial distribution of FLR frequency was Obayashi and Jacobs (1958), where the latitudinal variation was assessed. Geomagnetic pulsations observed using ground-based magnetometers were analyzed, for a period spanning about 80 days and covering a geomagnetic latitude range from approximately 40 ∘ to 60 ∘ . Although a large amount of scatter was present in the data, the results indicated that a clear decrease in pulsation frequency with increasing latitude was observed. The latitude dependence was also observed in subsequent studies (Orr & Matthew, 1971;Samson et al., 1971), leading to a notable study conducted by Samson and Rostoker (1972). Samson and Rostoker (1972) presented an analysis of Pc4 and Pc5 pulsations comprising approximately 185 h of geomagnetic pulsation activity, observed using ground-based magnetometers covering latitudes from approximately 58 ∘ to 78 ∘ . A linear trend in the data was identified, such that on average, the resonant pulsation frequency linearly decreased with latitude. As well quantifying the L value dependence, some consideration of dependences on MLT was included, with results indicating the dayside FLR frequency is statistically greater than the nightside FLR frequency. Variations in FLR frequency with MLT were further examined in studies by Yumoto et al. (1983) and Poulter et al. (1984). Using ground-based observations obtained by the STARE (Scandinavian Twin Auroral Radar Experiment) radar and the Slope Point radar, Poulter et al. (1984) identified 64 pulsation events, covering an interval of 1978-1980 and a latitude range of approximately 58 ∘ to 70 ∘ . The diurnal variation in average pulsation frequency was clearly identifiable, where frequencies were observed to maximize at dusk. Subsequent studies provided support for the observed latitudinal and MLT variations in FLR frequency (Engebretson et al., 1986;Glassmeier et al., 1984;Junginger & Baumjohann, 1984;Mathie et al., 1999;Plaschke et al., 2008;Takahashi et al., , 2004Takahashi et al., , 2015, using both ground-based and in situ data sets. A study conducted by Takahashi et al. (2002) provided further insight into the spatial variations of the resonant frequencies of field lines by presenting latitude-MLT maps. Using electric field and magnetic field observations from the AMPTE/CCE (Active Magnetospheric Particle Tracer Experiment / Charge Composition Explorer), a data set comprising of 3,171 events was obtained. A decrease in frequency with L shell, and frequency values minimizing in the noon-dusk sector were the key features observed. In a similar manner, Liu et al. (2009) presented an average map of FLR frequency distribution in L-MLT space, based on observations of electric fields and magnetic fields obtained by THEMIS (Time History of Events and Macroscale Interactions during Substorms) over the interval from November 2007 to December 2008. The average frequencies of Pc4 and Pc5 toroidal waves indicate a decreased frequency with L value and peak values located at approximately dawn, in agreement with prior studies. Although these studies clearly demonstrate the variability in FLR frequencies with L value and MLT, the temporal and spatial extent of the data sets are the key limitations. Overall, a review of previous studies indicates that there is a lack of large-scale, statistical studies of the FLR frequency distribution. This introduces difficulties in providing any in-depth comparisons of the time-of-flight calculations to direct observations, and consequently, a statistical analysis of FLR frequencies observed by IMAGE (International Monitor for Auroral Geomagnetic Effects) is also conducted here.
is a function of local plasma mass density, , and magnetic field strength, B, so will vary significantly along closed field lines. Taking this into account, the fundamental mode frequency, f , of a given field line can be approximated using the time-of-flight technique (Warner & Orr, 1979): where the integral is taken over the full field line length. Equation (2) calculates the time taken for an Alfvén perturbation to traverse along a given field line. This technique requires assumptions of the magnetic field line length and field-aligned variations of the magnetic field strength and plasma mass density, which determine the Alfvén speed distribution along the field line (equation (1)).
In their study, Wild et al. (2005) assumed the T96 magnetic field model, and a simple radial power law mass density model. Estimates of fundamental FLR eigenfrequencies (calculated using equation (2)) are shown in Figure 1, illustrating the time-of-flight technique. As the mass density models used throughout this study represent the average conditions in the closed magnetosphere, the T96 magnetic field model parameters used in all time-of-flight calculations presented correspond to typical conditions at spring equinox. The inputs are a solar wind dynamic pressure of 2 nPa, IMF (Interplanetary Magnetic Field) B y and IMF B z equal to 0 nT, and a Dst index value of 0 nT. The frequencies are shown plotted at the field line positions in the T96 magnetic equatorial plane (left panel) and the field line footprint latitude (right panel), considering field lines with L values spanning 4 ≤ L < 10. The radial power law mass density model is defined to have a mass density of 100 cm −3 at L = 4, decreasing with radial distance in the equatorial plane with a r −4 dependence, as used in previous studies (Warner & Orr, 1979;Wild et al., 2005). The mass density model is based on OGO (Orbiting Geophysical Observatories) observations of H + densities, presented in a study by Chappell (1972). Local time variations in the magnetic field configuration are represented as variations in the frequency with MLT.
In order to improve upon the estimates of Wild et al. (2005), this study will use the T96 magnetic field model combined with a more realistic mass density model. The T96 magnetic field model (Tsyganenko, 1996) is parameterized by the solar wind dynamic pressure, IMF B y , IMF B z , and Dst index. The mass density model used will now be briefly described. Using Cluster observations obtained between 2001 and 2012, Sandhu et al. (2016bSandhu et al. ( , 2016a obtained empirical models describing field-aligned variations of the total electron density, based on WHISPER (Waves of High frequency and Sounder for Probing of Electron density by Relaxation) data, and the average ion mass, based on CODIF (ion Composition and Distribution Function analyser) data. The empirical models included dependences on L value and MLT (Magnetic Local Time) and cover all MLT sectors. The L value represents the equatorial point of the field line, and details concerning the determination of the L value are given in Sandhu et al. (2016b). The empirical model describing field-aligned variations in electron density, n e , is given by where R norm describes the position along a given field line and is defined as the radial distance divided by the L value of the field line. The model is parameterized by the power law index, , power law equatorial electron density, n e0 , and the equatorial peak height, a. The electron density model describes variations for field lines within 4.5 ≤ L < 9.5. It is noted that the functions are not restricted to be continuous at R norm = 0.8. Previous empirical models describing electron density variations assume a power law form along the full field line length (Cummings et al., 1969;Denton et al., 2002;Denton, Menietti, et al., 2004;Goldstein et al., 2001;Sheeley et al., 2001), whereas the key feature of the Sandhu et al. (2016a) model is the presence of an equatorial enhancement, which is described through the inclusion of the Gaussian form function. The physical contributions to the electron density distribution, and the potential cause of the equatorial enhancement, are considered in detail in Sandhu et al. (2016a). The corresponding model for field-aligned variations in average ion mass, m av , is where is the power law index and m av0 is the equatorial average ion mass. This model covers the region spanning 5.9 ≤ L < 9.5. It is noted that the CODIF data set used to determine the empirical average ion mass model covers high-latitude points along field lines and is mapped to cover the equatorial region of the field lines. By combining the empirical models for total electron density, n e , (equation (3)) and average ion mass, m av , (equation (4)), using the simple relation = n e m av (5) the total plasma mass density, , can be obtained. Under this approximation, the plasma is assumed to be quasi-neutral, and the mass of electrons is assumed to be negligible compared to the mass of ions. It is also assumed that the CODIF ion measurements are representative of the full ion population. Observations of the average ion mass taken in the usual MAG (Magnetosphere) mode operation (0.025-40 keV/q) were compared to observations for the cold population, taken by CODIF operating in the RPA (Retarding Potential Analyser) mode (0.7-25 eV/q), and the results indicate that no dependence of the average ion mass on energy is present for the energy range (0.7 eV/q-40 keV/q) and the L values considered. A detailed justification of this assumption and analysis is included in Sandhu et al. (2016b). The key advantage of the method used by Sandhu et al. (2016bSandhu et al. ( , 2016a to estimate the total plasma mass density is that the contributions of both plasma number density and ion composition are included through the combination of the two models. The empirical mass density model, which covers all MLT sectors within 5.9 ≤ L < 9.5, is assumed in the time-of-flight calculations.

10.1002/2017JA024434
The main motivation for the development of the Sandhu et al. (2016a) mass density model and the choice for its inclusion in this study is due to the use of a more representative description of field line dependences, as opposed to a simplistic power law form used in previous models (Berube et al., 2005;Maeda et al., 2009;Menk et al., 1999;Takahashi et al., 2002Takahashi et al., , 2006Takahashi et al., , 2010Takahashi et al., , 2014. Denton, Takahashi, et al. (2004) first identified that the mass density distribution can be locally peaked at the magnetic equator, decreasing off-equator and then increasing at higher latitudes, demonstrating that a simple power law form does not fully describe the mass density variations along field lines. Following models (Denton et al., , 2009Takahashi & Denton, 2007;Takahashi et al., 2004) have accounted for the peaked field line distribution similarly to Sandhu et al. (2016a); however, they are not without their limitations. These studies use measurements of multiharmonic toroidal Alfvén wave frequencies to infer the field-aligned distribution of mass density, where an inversion technique is utilized to indirectly estimate the mass density. In addition, the data sets used all have limitations in terms of the spatial coverage provided. Takahashi et al. (2004) used CRRES (Combined Release and Radiation Effects Satellite) data covering only 1200-1800 MLT, and while this data set was extended by Denton et al. (2006) to a wider range, only 11% of the additional measurements lie outside this MLT sector. In contrast, Takahashi and Denton (2007) utilized GOES (Geostationary Operational Environment Satellite) 5 magnetometer data to examine mass density variations with good MLT coverage over all sectors. However, this particular study was limited to field lines at geosynchronous orbit only, significantly restricting the L range of the resulting model. The model presented by Sandhu et al. (2016a) is based on data sets that provide higher spatial coverage, for both L and MLT, as well as comprising direct measurements as opposed to indirectly inferring mass density values from magnetic field data. Furthermore, although previous studies account for a local equatorial peak in the mass density distribution along field lines, the peak was thought to result solely from increased heavy ion concentrations toward the magnetic equator Denton et al., 2006Denton et al., , 2009Takahashi & Denton, 2007;Takahashi et al., 2004). Sandhu et al. (2016bSandhu et al. ( , 2016a demonstrate that the equatorial peak in mass density is due to both the electron density and the average ion mass contributions, where the peak is most prominent on the nightside in the premidnight MLT sector. The physical contribution to the electron density equatorial enhancement is suggested to be the convection of plasma from the nightside plasma sheet toward Earth, as discussed in Sandhu et al. (2016a). As expected, the average ion mass is locally peaked at the equator across all local times due to the result of centrifugal forces acting more effectively on heavy ions Lemaire & Gringauz, 1998;Sandhu et al., 2016b). Cold plasma dynamics are dictated by the balance of corotational and convective electric fields, which generally means the flux tubes move more slowly and hence allow for more plasma buildup on the duskside compared to the dawnside (Chappell, 1972) but still convect plasma through to all local times. Hence, the combination of these factors results in a mass density equatorial peak that has an MLT asymmetry peaked in the premidnight MLT sector, as reported by (Sandhu et al., 2016a) and accounted for in the empirical mass density model.
An alternative method of determining the fundamental eigenfrequency for a given field line is using the Singer et al. (1981) method, where the MHD wave equation can be numerically solved using a fourth-order Runge-Kutta algorithm (Waters et al., 1995(Waters et al., , 1996. In the study conducted by Wild et al. (2005) the frequencies estimated using the time-of-flight approximation (equation (2)) were directly compared to frequencies calculated using the Singer et al. (1981) method, where the T96 magnetic field model and radial power law mass density model were used for both. Although the time-of-flight technique is relatively simplistic compared to the Singer et al. (1981) method, Wild et al. (2005) showed that the frequencies were in good agreement up to field line footprint latitudes above approximately 75 ∘ , where the estimates diverged. As the empirical mass density model of Sandhu et al. (2016aSandhu et al. ( , 2016b covers field lines with footprints well below 75 ∘ , the simpler and quicker time-of-flight approximation is chosen over the Singer et al. (1981) method to estimate the fundamental FLR eigenfrequencies.

Refined Estimates of FLR Eigenfrequencies
The time-of-flight technique is applied to calculate the fundamental FLR eigenfrequencies for a range of field lines, using the T96 magnetic field model and the empirical mass density model, and the calculated frequencies are presented in Figure 2. For comparison, the frequencies have been calculated using both forms of the electron density model: the power law form (upper panels) and the power law form combined with the Gaussian function at low latitudes (lower panels). The left panels show the frequencies plotted at the field lines' position in the T96 magnetic equatorial plane, and the right panels show the frequencies plotted at the field lines' footprint position, following the same format as Figure 1. The spatial coverage of field lines shown is restricted by the L range of the mass density model used (5.9 ≤ L < 9.5). Typically, the standard error  (4)).
for the mass density in a given L-MLT bin is ∼ 5.1 amu cm −3 , which translates to an error in the time-of-flight calculated frequencies shown in Figure 2 of approximately 0.03 mHz, assuming no error in the magnetic field model.  Wild et al. (2005), is included. The gray profiles represent the time-of-flight calculated frequencies using the radial power law mass density model, corresponding to values shown in Figure 1. The green profiles correspond to the empirical mass density model, where the Gaussian function was used at low-latitude points along a given field line in the electron density model. The purple profiles correspond to the empirical mass density model, where a power law form was extrapolated and used at low-latitude points along a given field line. Many models using the power law form for plasma number density are based on spacecraft observations, which are used to fit for the equatorial electron density model parameter, n e0 , using a power law function (equation (3a)). Assuming the power law index is described by equation (3d), electron density observations obtained by a spacecraft located at high latitudes along field lines, where the power law dependence is observed, would correspond to an equatorial electron density value of n e0 under a power law distribution. In contrast, for a spacecraft The gray lines correspond to a radial power law mass density model, and the green and purple lines correspond to the empirical mass density model (a combination of the electron density model defined by equation (3) and the average ion mass model defined by equation (4), as discussed in section 2). The purple lines show calculated frequencies using the extrapolated power law form for the low-latitude electron density dependence, based on hypothetical spacecraft observations of electron density off-equator (solid) and in the equatorial plane (dashed). The green solid lines show calculated frequencies using the Gaussian function for the low-latitude electron density dependence. Each panel corresponds to field lines at different MLTs, as indicated. located in the magnetic equatorial plane, the observed equatorial electron density would be equal to n e0 + a, due to the localized enhancement. Therefore, assuming a power law form for the field-aligned electron density can produce two models with different equatorial values, due to spacecraft locations. The corresponding mass density models were used to calculate the fundamental FLR eigenfrequencies using the time-of-flight technique, as represented by the purple profiles in Figures 3 and 4. The purple solid profiles correspond to the power law form using an equatorial electron density equal to n e0 , representing results for high-latitude spacecraft observations, and the purple dashed profiles correspond to the power law form using an equatorial electron density equal to n e0 + a, representing results for equatorial spacecraft observations. Therefore, the electron density models used for the solid and dashed profiles are n e = n e0 R − norm and n e = (n e0 + a)R − norm , respectively. The parameters n e0 , a, and are defined in equation (3). The spatial variations in frequency exhibit some notable features, as well as several clear differences in estimated frequencies depending on the choice of mass density model, which will now be discussed. Note that the reduced coverage of the green and purple profiles relative to the gray profiles is due to the L range of the empirical mass density model, as only field lines with an L value within 5.9 ≤ L < 9.5 have a corresponding time-of-flight calculated frequency value. Consequently, the MLT range is more limited at higher footprint latitudes (Figure 4) because the change in L value for a given footprint latitude with MLT is larger, so field lines over all MLT do not fall within the empirical mass density model L range.

Discussion
The results of the time-of-flight calculations indicate several features and differences in profiles through using various mass density models. The aim of the time-of-flight analysis is to improve the estimates of Wild et al. (2005), where a simple radial power law mass density model was used with the realistic T96 magnetic field model. Through a comparison of the frequency profiles corresponding to the radial power law mass density model (solid gray) and the empirical mass density model (solid green), shown in Figures 3 and 4, the contribution of the mass density distribution to fundamental FLR eigenfrequencies can be assessed. Figure 3 shows the frequency variation with footprint latitude. In agreement with multiple previous FLR observations and time-of-flight calculations, the estimated frequency for the radial power law mass density model (solid gray profile) is observed to decrease with increasing footprint latitude, across all MLT sectors (see also Figure 1), which is a feature attributed to the magnetic field contribution, as field line length increases and magnetic field strength decreases for increasing footprint latitude (Engebretson et al., 1986;Glassmeier et al., 1984;Liu et al., 2009;Mathie et al., 1999;Obayashi & Jacobs, 1958;Orr & Matthew, 1971;Plaschke et al., 2008;Poulter et al., 1984;Samson & Rostoker, 1972;Samson et al., 1971;Takahashi et al., 2002Takahashi et al., , 2004Takahashi et al., , 2015Wild et al., 2005;Yumoto et al., 1983). This feature is also observed for the profile corresponding to the empirical mass density model (solid green). As the empirical mass density describes decreased values with increasing L value, this acts to increase the Alfvén speed (equation (1)) and therefore increase the fundamental FLR eigenfrequency (equation (2)  the empirical mass density model profile (solid green) in Figure 3, and the decrease in frequency with increasing footprint latitude indicates that magnetic field variations (changes in field line length and magnitude along the field line) are the dominant contribution. Figure 4 illustrates the frequency variations with MLT, and it can be seen that the profile corresponding to the use of the radial power law mass density model (solid gray) demonstrates strong dependences with MLT, with frequency values maximizing on the dayside. The diurnal variation in standing Alfvén wave frequency is a feature observed in previous studies (Junginger & Baumjohann, 1984;Liu et al., 2009;Mathie et al., 1999;Plaschke et al., 2008;Poulter et al., 1984;Takahashi et al., , 2015Yumoto et al., 1983). The radial power law mass density model profiles only includes MLT asymmetries for the magnetic field contribution and assumes the plasma mass density contribution has no MLT dependence. Therefore, the magnetic field results in higher frequencies on the dayside, where the magnetic field lines are shortest and the magnetic field strength maximizes (Mathie et al., 1999;Wild et al., 2005). By replacing the radial power law mass density model with the empirical mass density model, plasma mass density contributions to the MLT asymmetry in time-of-flight frequencies are included. The solid green profile in Figure 4 demonstrates the MLT dependences of the estimated frequencies using the empirical mass density model, which, similarly to the use of the radial power law mass density model (solid gray profile), exhibits variations with MLT. In comparison to the solid gray profile, the solid green profile has a larger amplitude for the MLT dependence and the peak frequency MLT location is shifted toward dawn with the lowest values present toward dusk. This is due to the mass density values maximizing in the evening sector, as a result of the plasmaspheric bulge and average ion mass enhancement (Sandhu et al., 2016b(Sandhu et al., , 2016a, which acts to decrease the time-of-flight calculated frequencies (Junginger & Baumjohann, 1984;Takahashi et al., 2015). Furthermore, the mass density contribution increases the MLT asymmetry of the frequency profiles, a feature observed through a comparison of the solid gray profiles and solid green profiles in Figure 4. It is noted that the MLT asymmetry in frequency is greatest for lower latitudes (upper panel), corresponding to the strongest MLT asymmetries in mass density (Sandhu et al., 2016a).
A comparison of the profiles corresponding to the empirical mass density model (green profile) and the radial power law mass density model (gray profile), in Figures 3 and 4, indicates some large deviations, such that the empirical mass density model profile consistently displays lower frequency values. This suggests that the radial power law mass density model tends to underestimate mass density values compared to the empirical mass density model (equations (2) and (1)). The reason for this could be due to the radial power law mass density model neglecting contributions of heavy ions to the total plasma mass density, whereas the empirical mass density model accounts for ion composition variations. The average ion mass model clearly indicates that it is inappropriate to assume a H + plasma for this region and that O + ions are a significant factor (Sandhu et al., 2016b), which act to increase the total plasma mass density.
Overall, Figures 3 and 4 demonstrate that the use of a realistic empirical mass density model, as opposed to a simple radial power law mass density model, refines the time-of-flight calculations. The consideration of ion composition and MLT asymmetries in the empirical mass density model provide further features in the spatial variations of fundamental FLR eigenfrequencies, resulting in decreased frequency values, enhanced MLT asymmetries, and frequency values peaking at earlier MLTs compared to the use of a radial power law mass density model. This indicates that the mass density distribution represents a significant contribution to determining the fundamental eigenfrequency of a field line, and using a simplified power law mass density model can produce large deviations in the calculated frequencies.

Contribution of Equatorial Enhancement
A key feature of the empirical mass density model is the inclusion of an observed localized enhancement in number density at the magnetic equatorial plane, neglected by previous models. As mentioned in section 2, previous studies assume that plasma number density follows a power law form along closed magnetic field lines, such that the number density maximizes at the field line ionospheric footprint and approaches a minimum at the magnetic equatorial plane. As shown by the empirical model for total electron density (equation (3)), Sandhu et al. (2016a) found that although this distribution is appropriate for high-latitude points along field lines, close to the magnetic equatorial point of a field line, a statistically significant localized enhancement in total electron density is observed. This was described by including a Gaussian form function at low latitudes (equation (3b)). Figure 2 shows time-of-flight calculations corresponding to electron density models using the power law form extrapolated to low latitudes along field lines (upper panels) and using the Gaussian function at the equatorial region (lower panels). A comparison of the estimated frequencies presents some notable differences, indicating that the equatorial enhancement in plasma number density represents a nonnegligible contribution in determining fundamental FLR eigenfrequencies. This will now be examined further, through a consideration of the latitude and MLT profiles. Figure 3 shows differences between frequencies calculated using the empirical mass density model (green profile) and the empirical power law forms (purple profiles). As previously mentioned, the solid purple profiles are intended to represent hypothetical results for high-latitude spacecraft observations, and the dashed purple profiles represent hypothetical results for equatorial spacecraft observations. The solid purple profiles, representing a power law electron density distribution with an equatorial value of n e0 , is larger than the dashed purple profiles, representing a power law electron density distribution with an equatorial value of n e0 + a. This is expected, because the dashed purple profiles correspond to larger number density values, and hence larger mass density values compared to the solid purple profiles, across the full field line length. Therefore, as shown by equations (2) and (1), the decreased Alfvén speeds result in decreased eigenfrequencies for a given field line. Figure 3 shows that the green solid profiles, corresponding to use of the empirical mass density model accounting for the equatorial enhancement in plasma number density, have frequency values above the purple dashed profiles and below the purple solid profiles. This is due to the number density values from the electron density model. Considering the corresponding electron density values used for each profile, the solid purple profiles neglect the equatorial enhancement and provide values less than the solid green profiles. On the other hand, the dashed purple profiles do not account for the decrease in electron density away from the peak in the Gaussian region, assuming a power law form across the full field line, and provide values higher than the solid green profiles. Therefore, the resulting mass density values for the solid purple profiles will be less than the solid green profiles, and the mass density values for the dashed purple profiles will be greater than the solid green profiles. The time-of-flight calculations shown in Figure 3 show that the solid purple profiles overestimate the frequency and the dashed purple profiles underestimate the frequency, through a comparison with the solid green profiles. Therefore, by assuming a power law form for the plasma number density, deviations from the corresponding observed peaked number density distribution occur. Figure 3 indicates that the deviations are negligible at lower footprint latitudes for the range considered here and increase at higher footprint latitudes. This is due to the observed equatorial enhancement in electron density becoming increasingly prominent at larger L values, attributed to the presence of the nightside plasma sheet population.
In addition to the latitude profiles (Figure 3), the MLT profiles shown in Figure 4 also demonstrate that the solid purple profiles show frequency values greater than the solid green profiles, and the dashed purple profiles show frequency values lower than the solid green profiles. The difference between the values is smallest at around noon and greatest in the nightside MLT sectors for field lines with a fixed footprint latitude. This is due to the relative magnitude of the equatorial number density enhancement being greatest for nightside MLT sector, corresponding to the presence of the plasma sheet population. Therefore, the purple profiles underestimate/overestimate frequencies most significantly for nightside MLT sectors. It can be concluded, from Figures 3 and 4, that using a power law form for number density in the low-latitude region can result in large deviations compared to including a Gaussian form function to account for the equatorial number density enhancement, particularly at high L values for nightside MLT sectors.

Statistical Analysis of IMAGE Data
In order to test the validity of the calculated frequencies, this study aims to compare the time-of-flight estimations to ground magnetometer FLR observations. The formation of standing Alfvén waves on closed geomagnetic field lines excites resonant oscillations of the field line, which can be observed in corresponding magnetic field perturbations on the ground. Toroidal standing Alfvén waves drive perturbations in the north-south direction, or x-component, at the ground, and by identifying the frequency of the perturbations, the natural frequency of the relevant field line is obtained. By identifying the fundamental FLR eigenfrequencies using ground magnetometer observations, the variations in average frequency with L and MLT are determined. The distribution of average frequencies can then be compared to the calculated time-of-flight frequencies.
The IMAGE ground magnetometer array is utilized here, providing observations from 35 stations located within 50 to 75 ∘ magnetic latitude, with a time resolution of 10 s (Lühr, 1994). The data set employed in this study comprises observations obtained between 2001 and 2012, corresponding to the same time period as the Cluster data set used for the mass density model. In order to routinely identify resonant oscillations in the observed north-south magnetic field component, B x , an automated FLR identification method was developed, which will now be detailed.

Automated FLR Identification
The cross-phase technique (Waters et al., 1991) is chosen to identify resonant toroidal oscillations. Consider field line oscillations monitored by two ground magnetometer stations, latitudinally separated and approximately located along the same longitudinal meridian, in the presence of a fast mode incoherent driving wave. The fast mode wave consists of a continuum of frequencies. For each station, the corresponding field line will resonate independently at its eigenfrequency, determined by the field line length, as well as the variation in magnetic field strength and plasma mass density along the field line (equation (2)). Due to the variation in eigenfrequency with latitude, field lines located at the northern station and southern station will resonate with slightly different frequencies. Assuming a linear variation in the fundamental eigenfrequency, evaluating the phase difference between the north-south (B x ) components from the two stations provides the cross-phase difference function. The cross-phase difference function will be maximized at the eigenfrequency corresponding to the field line located at the midpoint of the two stations considered. Therefore, using the north-south (B x ) component from two closely spaced stations, the fundamental eigenfrequency of the midpoint field line can be determined via the cross-phase technique. The cross-phase technique makes some notable assumptions. Magnetic field lines are approximated as uncoupled, such that the oscillations are purely toroidal; only damping through ionospheric dissipation is present; and the coupling between the fast mode and toroidal wave mode is identical for all field lines. Despite these simplifying assumptions, the method has proved relatively robust for identifying FLR frequencies. It remains reliable at low signal strengths and is less sensitive to amplitude variations compared to alternative techniques (Waters et al., 1995).
The cross-phase technique is implemented in this study as part of an automated identification method. The automated method is demonstrated in Figure 5, which shows a representative sample of IMAGE data obtained between 00:00 and 12:00 UT on 10 November 2006, for the Abisko (ABK) and Tromsø (TRO) ground magnetometer stations. The steps taken to identify resonances and the frequencies of resonant oscillations in the B x data are now detailed, where the relevant panels corresponding to each step of the process are as numbered on the right in Figure 5. 1. Taking the magnetic field data for two ground magnetometer stations that are closely spaced in latitude, the B x time series are band-pass filtered between 100 s and 1,000 s. Figure 5 shows the filtered B x data from (a) ABK and (b) TRO. 2. Using a 40 min sliding window, incremented by 10 min, the power spectrum and cross-phase spectrum between the two B x time series were computed, following Waters et al. (1991). The spectra for the sample considered here are shown in Figures 5c and 5d. A window width of 40 min was deemed suitable as it will contain at least two wave periods for the lowest frequency perturbations being considered and is similar to widths used in previous studies (e.g., Berube et al. (2003)). 3. A defining feature of a FLR is a relatively stable oscillation with time, so in order to neglect variations on small time scales, hourly averages of the power and cross-phase spectrum are calculated for each frequency bin. This is illustrated in Figures 5e and 5f. Furthermore, frequency bins with an average power below the corresponding median hourly power are excluded and plotted as black. Cross-phase values with a small power are essentially random and add noise to the spectrum. 4. The next step follows from the automated method of Berube et al. (2003). The t statistic, defined as the mean cross phase divided by the cross-phase standard deviation, is calculated for each hourly frequency bin. The resulting spectrum is shown in Figure 5g for the sample considered. The t statistic is low for bins consisting of large cross-phase fluctuations, and consequently, it is useful for identifying stable, nonrandom, peaks in cross phase. All bins where t ≤ 1 are plotted as black for the t statistic spectrum shown in Figure 5g. Values where t ≤ 1 have a standard deviation greater than the mean cross phase, and the significance of the associated cross phase is low (Berube et al., 2003;Press, 1992). This condition will be detailed further in the next step.
5. The final step of the method identifies the frequency bin associated with the peak hourly cross-phase value (i.e., the frequency of the peak cross-phase values in Figure 5f ), taking each hour interval individually. The frequencies for this sample are plotted in Figure 5h for each hour interval, shown by the horizontal black and gray lines. In order to identify whether the peak in cross phase at the frequency is associated with an FLR, the value of the t statistic corresponding to the frequency bin for the corresponding hour interval is consulted (Figure 5g). Provided the t statistic meets the criteria t > 1, the frequency is validated as corresponding to a FLR and plotted as black in Figure 5h. If the criterion is not met, it is interpreted that a FLR is not present. The frequency is considered not valid and plotted as gray.
Therefore, this method enables the automated identification of FLRs, and the field line eigenfrequencies, from simultaneous north-south (B x ) observations at two latitudinally separated ground magnetometer stations. The method was tested against a number of representative samples of magnetometer time series; an example of which is shown in Figure 5, and no biases in the method could be identified. For the sample shown in Figure 5, the cross-phase spectrum shows a clear enhancement in cross phase at a frequency of approximately 4 mHz (Figure 5c), which can be visually identified as evidence of an FLR. The automated method identifies valid frequencies across this time interval (Figure 5h), in agreement with a visual identification. In addition, Figure 5h shows calculated frequencies using the time-of-flight approximation, assuming the T96 magnetic field model and the empirical mass density model of Sandhu et al. (2016a). The solar wind dynamic pressure, IMF B y , IMF B z , and Dst index values at the observation times are inputted into the T96 magnetic field model. A comparison of the frequencies identified using the automated method is in encouraging agreement with the time-of-flight estimates.
In order to further verify the validity of the automated identification method, a comparison to the results of Kale et al. (2009) was performed. Kale et al. (2009) included an analysis of IMAGE ground magnetometer data from the Nurmijärvi (NUR) and Oulujärvi (OUJ) stations, obtained between 28 October and 7 November 2003, corresponding to the Hallowe'en 2003 geomagnetic storm. Kale et al. (2009) identified FLR signatures through a visual inspection of the cross-phase spectra, determining the frequency of the oscillations. The automated method described here was applied to the same data set, with the resulting frequencies observed to be in close agreement with the results of Kale et al. (2009). This provides strong evidence in support of the validity of the automated FLR detection method.

IMAGE Analysis Results
The automated FLR identification method is used in this study to examine the average variations in FLR eigenfrequencies. Multiple station pairings were manually selected, with suitable separations and covering a range of latitudes. The magnetic latitude locations of the station pairs are shown in Figure 6a, with the midpoint magnetic latitude indicated by the red points. The choice of station pairs was influenced by the station separation. Closer station spacing results in smaller cross-phase peak values at the resonant frequency, whereas greater spacing results in increased noise in the cross-phase spectrum due to lower coherence Waters et al., 1995). Previous studies suggest an optimum station separation of approximately 110 km in latitude (Chi et al., 2013;Menk et al., 2004). Although ideal station spacing is restricted by the station locations of the IMAGE array, the station pairs used in this study are reasonable. Furthermore, latitudinal coverage is limited at the higher latitudes of the range considered, due to the IMAGE station locations, although it remains sufficient. For each station pair, the automated FLR identification method was applied to the simultaneous B x observations, for all data obtained within the time period of 2001-2012. This provides a database of FLR frequency observations for each station pair, which are presented in the following section.
Using the IMAGE data set, covering 2001-2012, the automated FLR identification method is applied to observations for multiple station pairs, as described in section 3.1. Figure 6 shows the results from the analysis. Figure 6c shows the frequencies, binned for MLT and averaged, and Figure 6d shows the corresponding number of detected FLRs, identified via the automated method. Each row of data shown in Figures 6c and 6d represents observations from the station pair in the same row of Figure 6a. Clear dependences in the average frequency distribution with MLT and latitude are apparent from Figure 6c. For comparison, the time-of-flight calculations using the empirical mass density model of Sandhu et al. (2016aSandhu et al. ( , 2016b are included in Figure 6b, where each row shows estimates for field lines at the corresponding midpoint AACGM (Altitude Adjusted Corrected Geomagnetic) latitude shown in Figure 6a.

A Comparison of Time-of-Flight Calculations and IMAGE Observations
The analysis of the IMAGE data set has provided observations of fundamental eigenfrequencies for field lines spanning latitudes within 50 to 75 ∘ magnetic latitude and over all MLT. In the context of previous work, the average frequencies observed in this study are in close agreement with the results of Poulter et al. (1984). For example, both this study and Poulter et al. (1984) observe frequencies of approximately 4 mHz at footprint latitudes of 65 ∘ at 06 MLT. The average FLR frequencies observed here are generally lower than the frequencies observed by Samson and Rostoker (1972). The difference is most significant for field lines in the dayside sector. For example, field lines at latitudes of 65 ∘ at 12 MLT have an average frequency of approximately 5 mHz, whereas Samson and Rostoker (1972) report frequencies of approximately 7-8 mHz. A possible cause of the discrepancy may be due to the relatively limited temporal extent of the data set used by Samson and Rostoker (1972), covering only 1 year compared to the 11 years covered by the IMAGE data set used in this study. A discussion of the other biases is considered at the end of this section, which will also contribute to deviations in the results compared to previous studies.
The validity of the refined time-of-flight estimates can now be assessed through a comparison with the statistical analysis of FLRs observed using IMAGE magnetic field data. The distribution of observed FLR frequency with footprint latitude and MLT is shown in Figure 6c, where results are based on IMAGE data obtained between 2001 and 2012 and station pairs have been ordered from top to bottom in decreasing midpoint latitude. Figure 6d shows the corresponding number of detected FLRs for each MLT bin. It can be seen from Figure 6c that the field line eigenfrequency tends to decrease with increasing footprint latitude, in agreement with the time-of-flight calculations (Figure 2). In order to present a more quantitative comparison of the time-of-flight calculations and IMAGE analysis, the latitudinal dependences are explored further in Figure 7. are discussed later, it is difficult to compare the effects of using the empirical mass density model (Sandhu et al., 2016a) compared to the radial power law mass density model in the time-of-flight calculations. However, a slight difference is apparent for 1800 MLT (Figure 7e), where the green profile appears to provide estimates closer to the IMAGE observations compared to the gray profile. This suggests that the time-of-flight calculations are improved through the use of a more realistic mass density model in the dusk sector. In this region, the mass density is strongly influenced by the plasmaspheric bulge and enhanced heavy ion population. These features are encapsulated in the empirical mass density model, unlike the simple radial power law mass density model, such that the resulting time-of-flight calculations demonstrate better agreement with the IMAGE observations in this region. Figures 7a-7d also show the number of detections for each station pair binned for frequency, using a bin size of 1 mHz, as indicated by the color of the bins. Figures7a-7d show that FLRs are detected at all frequencies considered, displaying significant variability in the oscillation frequency. The causes of the variability are discussed later.
In addition to latitudinal features, Figure 6b also shows MLT asymmetries in the frequency of standing Alfvén waves observed in the IMAGE data set. Frequencies peak in the morning MLT sector, in agreement with the time-of-flight calculated frequencies (Figure 2). Figure 8 compares the MLT variations in field line eigenfrequency, for the time-of-flight calculations (green profiles and gray profiles), and the average IMAGE observations (red profile), where the color coding used is as in Figure 7. A representative sample of station pairs is shown. Each panel corresponds to a different station pair from those shown in Figure 6a, and the midpoint latitude is labeled on the right of each panel for reference. The results shown in Figure 8 provide further evidence that the time-of-flight calculations are improved through using the empirical mass density model, as opposed to the radial power law mass density model. The estimated frequencies corresponding to the empirical mass density model (green profile) tend to have values closer to the average IMAGE frequencies (red profile), than estimated frequencies using the radial power law mass density model (gray profile). Similarly to Figure 7, Figure 8 shows the number of detections binned for MLT, using a bin size of 1 h, and frequency, using a bin size of 1 mHz, as indicated by the color of the bins. Large variability in the detected frequencies is apparent over all MLT, which is further discussed later in this section. It is also noted here that the occurrence of detections varies with MLT (clearly shown in Figure 6d), with decreased detections on the nightside compared to the dayside. This is an expected feature and is due to low ionospheric conductances decreasing the magnitude of the phase shift of the FLR. Figures 7 and 8 show that there is a large amount of variation in the frequencies observed in the IMAGE data set. In addition, there are deviations between the observed frequencies, the calculated frequencies using the time-of-flight technique, and frequencies reported by previous studies. The spread in values and differences in values will have contributions from biases in both the models used in the time-of-flight analysis, as well as the automated FLR identification method, which will now be considered.
The time-of-flight calculated frequencies are biased by the choice of magnetic field model and plasma mass density model. The T96 magnetic field model is known to overestimate the ring current contribution to the magnetic field (Zhang et al., 2010), which would correspondingly underestimate the frequencies in the time-of-flight analysis. The plasma mass density model used in this study Sandhu et al. (2016bSandhu et al. ( , 2016a, while an improvement over previous models, also has limitations that introduce biases to the time-of-flight calculations. In particular, Sandhu et al. (2016b) shows that, although efforts were made to remove radiation belt contamination that overestimates the average ion mass, up to 4% of CODIF observations within a bin are expected to suffer from penetrating energetic radiation belt particles. An overestimated average ion mass would result in an overestimated mass density, which correspondingly underestimates the field line eigenfrequencies using the time-of-flight approach. However, the fraction of contaminated data is sufficiently small that this bias is not considered to be significant. Although the mass density model used in the analysis has limitations and scope for further improvement, the results present notable improvement compared to previous models. The model includes the contribution of ion composition, accounts for a localized enhancement in electron density, and incorporates spatial variations with L and MLT.
The results from the analysis of IMAGE observations involve the use of an automated FLR identification method. It is assumed that this method selects only the fundamental mode; however, Obana et al. (2008Obana et al. ( , 2015 present observations of quarter mode waves occurrent at dawn and dusk MLT sectors when there is a strong hemispheric asymmetry in ionospheric conductivity. They can be identified from ground magnetometer observations by a cross mode waves, the automated FLR identification method used in the analysis of the IMAGE data could include detections of quarter mode perturbations. The quarter mode waves have a frequency equal to half the fundamental eigenfrequency of the field line, which would act to decrease the average frequencies observed near dusk and dawn. Bulusu et al. (2015) report that quarter mode waves may be observed in up to 55% of events in the dusk sector at geosynchronous orbit. However, due to the high damping of the quarter mode waves (Obana et al., 2015), this factor is not expected to significantly bias the average frequencies reported in Figures 6-8, and no obvious decrease at dawn and dusk from the background MLT variation is apparent.
Another significant factor that contributes to the large spread in detected eigenfrequencies from the IMAGE observations is the high variability of field line properties throughout the time period considered. The magnetic field configuration, as well as the plasma mass density distribution along the field line, changes in response to solar wind conditions and geomagnetic activity. In particular, the field lines map from the ionospheric footprint to different regions of the magnetosphere for different conditions. Therefore, the eigenfrequency of a field line at a given footprint latitude and MLT can take a wide range of values, as reflected in Figures 7 and 8. It is noted here that the field line configuration (field line length and variation of field strength along the field line) is expected to be the dominant factor relative to the mass density variation along the field line, due to the dependence of the Alfvén velocity (proportional to the magnetic field strength but inversely proportional to the square root of the mass density). The variability in field line configuration is largely defined by the solar wind dynamic pressure, the convective electric field, and the ring current, which act to compress and inflate the magnetosphere as the factors change. The contracting and stretching of the field lines, moving up to several R E in the equatorial plane, act to blur latitudinal variations in the eigenfrequencies measured at the field line footprints. For example, the plasmapause radial position is controlled by the strength of the convective electric field and can vary from 2 to 8 R E (Chappell, 1972), with a typical location of 5 to 7 R E (Darrouzet et al., 2013). For a field line that moves to within the plasmasphere, the relatively high plasma mass density would result in a decrease in the measured eigenfrequency. The blurring of the latitudinal variations, as a consequence of plasmapause motion and geomagnetic conditions, is apparent from Figure 7, which shows relatively weak variations with footprint latitude compared to the time-of-flight calculations. In contrast, the local time variations of the IMAGE observations appear to be well preserved ( Figure 8). This is expected to be due to relatively limited azimuthal perturbations compared to the radial perturbations imposed on the magnetosphere. We also note here that the presence of the plasmapause may contribute to the lower mean eigenfrequencies observed in the IMAGE data set compared to some previous studies (e.g., Samson & Rostoker, 1972).
It is also important to note that the mean frequencies reported from the IMAGE data set may not be comparable to the estimated time-of-flight frequencies if the detection efficiency of FLRs using the automated technique possesses a bias. For example, if the automated FLR detection technique identifies FLRs more easily during storm times, compared to quiet times, then the mean frequency observed from the IMAGE data set would be biased toward storm times and not representative of the same conditions as described by the time-of-flight calculations. Although an examination of a representative sample of magnetometer data did not suggest any detection bias with geomagnetic activity, this will be assessed in detail in a future study.
At higher footprint latitudes, some deviations between the estimated and observed frequencies are apparent. The time-of-flight estimates using the empirical mass density model (green profile) appear to be consistently lower compared to the average IMAGE observations, evident from both Figures 7 and 8. Field lines at high footprint latitudes are highly variable in terms of the magnetic field configuration as well as the plasma mass density along the field line, due to changes in solar wind conditions and geomagnetic activity. Therefore, it can be expected that the mean frequency of detections may not closely correspond to the time-of-flight calculations that are based on typical quiet time conditions at spring equinox (see section 2).
In addition, it is noted that the station pair coverage of IMAGE ground magnetometer stations is relatively sparse for the higher magnetic latitudes considered here, in comparison to station pairs located at magnetic latitudes of approximately 65 ∘ . Furthermore, due to the restricted station coverage, station spacing deviates slightly from ideal separation. As a consequence, the reliability of the average frequencies observed in the IMAGE data set is reduced for the higher latitude station pairs shown in Figures 6-8.
Although there are various important factors that introduce biases for both the time-of-flight calculations and the analysis of the IMAGE data set, the results presented in Figures 7 and 8 clearly demonstrate the improvement in estimates through the use of a more realistic empirical mass density model as opposed to a simple radial power law model. Furthermore, the analysis of IMAGE observations has provided a statistical understanding of how field line eigenfrequencies vary spatially and the variability of the eigenfrequencies, extending upon previous literature.

Conclusions
In this paper, we improve upon current time-of-flight estimates of the fundamental FLR eigenfrequencies for closed geomagnetic field lines. By implementing a more realistic mass density model (Sandhu et al., 2016a) combined with the T96 magnetic field model, the eigenfrequencies were estimated for 5.9 ≤ L < 9.5 across all MLT sectors. This represents an improvement on the estimates of Wild et al. (2005), where the T96 magnetic field model was used in conjunction with a relatively simplistic power law mass density model. The results showed significant deviations in calculated frequencies through the use of a more realistic mass density model. Furthermore, the effect of including the equatorial enhancement in number density was examined.
In order to test the validity of the refined estimates, the calculated frequencies were compared to the average frequencies of resonant field line oscillations observed by ground magnetometers in the IMAGE array. The determination of the average frequency, and its variations with footprint magnetic latitude and MLT, involved the implementation of an automated FLR identification routine, which is based on the cross-phase method (Waters et al., 1991). A comparison of the time-of-flight calculations and average frequencies observed by IMAGE indicates encouraging agreement in the values, as well as agreement in latitudinal and local time dependences. Furthermore, by using the more realistic empirical mass density model, the time-of-flight estimates appear to provide values closer to the IMAGE observations, than calculations using a simple radial power law mass density model. Therefore, it can be concluded that the time-of-flight estimates here provide valid frequency estimates, relative to previous estimates using a simple mass density model.