Cross‐Phase Determination of Ultralow Frequency Wave Harmonic Frequencies and Their Associated Plasma Mass Density Distributions

Latitudinally spaced ground‐based magnetometers can be used to estimate the eigenfrequencies of magnetic field lines using the cross‐phase technique. These eigenfrequencies can be used with a magnetic field model and an assumed plasma mass density distribution to determine the plasma mass density in the magnetosphere. Automating this process can be difficult, and so far, it has not been possible to distinguish between the different harmonics. Misidentification of the harmonic mode will lead to incorrect estimations of the plasma mass density. We have developed an algorithm capable of identifying multiple harmonics in cross‐phase spectrograms, using International Monitor for Auroral Geomagnetic Effects magnetometers. Knowledge of multiple harmonics allows the distribution of plasma mass density to be estimated instead of assumed. A statistical study was performed that showed clear common bands of eigenfrequencies, interpreted as different harmonics. These eigenfrequencies were lowest in the early afternoon and at higher latitudes. There was also a greater occurrence of measurements in the dayside. We then modeled the plasma mass density distribution with a power law characterized by the exponent p and compared the model eigenfrequencies to the data. This suggested that the even modes did not form during the interval of this study. Examination of the harmonic spacing and the high occurrence of the third harmonic supported this suggestion. We attribute the absence of the even modes to the driving mechanisms. Finally, we show that an equatorial bulge in plasma mass density was not present in our study.


Introduction
Magnetic field lines can host standing Alfvén waves if both ends of the field line are bounded by a conducting ionosphere. The eigenfrequencies of a magnetic field line are determined by its length, magnetic field magnitude, polarization, and the distribution of plasma mass density along it. Therefore, these eigenfrequencies can be used with a suitable magnetic field model to investigate magnetospheric plasma populations (e.g., Berube et al., 2006;Menk et al., 1999Menk et al., , 2004Singer et al., 1981;Waters et al., 1995). The polarization of standing Alfvén waves is usually considered to be either toroidal or poloidal, where the displacement of the plasma is either azimuthal or radial, respectively. In reality, the polarization is most likely an intermediate polarization due to the nonaxisymmetry of the magnetosphere (e.g., Claudepierre et al., 2010;Elsden & Wright, 2018;Wright & Elsden, 2016). However, the poloidal polarization has a smaller scale size and will transform into the toroidal mode over time (Mager & Klimushkin, 2006). Therefore, it is less likely to be seen in ground magnetometer data. For this work, we have only looked at toroidally polarized waves.
Ground-based magnetometers can be used to estimate the eigenfrequencies of resonating magnetic field lines. It can be difficult to separate the resonating signal from other signals such as the driving source or surface waves on the magnetopause (e.g., Chi et al., 2013). This was resolved by Baransky et al. (1985), who showed that the eigenfrequency of standing Alfvén waves could be estimated from the cross-amplitude spectra of two latitudinally spaced ground-based magnetometers. Waters et al. (1991) later demonstrated that using the cross-phase spectra was more reliable. Menk et al. (2014) has shown that cross-phase can be used to monitor a wide range of dynamic phenomena including the location of the plasmapause and plumes. Many other studies have utilized cross phase for monitoring eigenfrequencies (e.g., Clausen et al., 2008;Kale et al., 2007;Obana et al., 2008;Sandhu et al., 2018;Waters et al., 1995Waters et al., , 1996. If only the fundamental frequency is known, then the distribution of plasma mass density is unknowable. It must then be assumed in order to calculate the equatorial mass density. The form of this distribution Journal of Geophysical Research: Space Physics 10.1029/2018JA025487 determines the spacing of the different harmonics. The most common form for describing the plasma distribution along the field line is a power law (equation (1)). This simple model has been used by many authors (e.g., Angerami & Carpenter, 1966;Cummings et al., 1969;Denton & Gallagher, 2000;Denton et al., 2015;Menk et al., 2004;Sinha & Rajaram, 1997;Takahashi & McPherron, 1982;Waters et al., 1996); represents the plasma mass density, r represents the radial distance on the field line from the center of the Earth, and p is a free parameter. Subscripts eq represent values in the equatorial plane. If p > 0, there is a minimum in plasma mass density at the equator. If p < 0, this gives a maximum at the equator.
This simple power law may not always be appropriate. Takahashi et al. (2004), Denton et al. (2004Denton et al. ( , 2006, Takahashi and Denton (2007), and Denton et al. (2009) suggested there was a peak in mass density at the equator, based on satellite measurements. A recent statistical study by Sandhu et al. (2016aSandhu et al. ( , 2016b used Cluster data from 2001 to 2012 to deduce the overall mass density distribution between 5.9 < L < 9.5. They fitted equations to these data and deduced that the power law distribution was only applicable away from the magnetic equator. This was due to a bulge in electron number density and an increase in heavy ion concentration at the equator. Sandhu et al. (2016aSandhu et al. ( , 2016b also showed that the plasma distribution was not symmetrical in magnetic local time (MLT), with the bulge being greater in the dusk sector. These studies suggest the simple power law model is insufficient to describe the plasma mass density distribution of the magnetosphere.
In most circumstances, estimations of eigenfrequencies using cross-phase have been undertaken on a case-by-case basis (e.g., Menk et al., 1999), with appropriate features being selected by eye from dynamic cross-spectra. For large scale statistical work, this approach is inefficient. A few authors have attempted to automate this process. Berube et al. (2003) were the first authors to automate searches for eigenfrequencies in dynamic cross-phase spectra. They applied this method to low-latitude magnetometers within 39 ∘ to 50 ∘ Corrected GeoMagnetic latitude. Their method yielded an uncertainty in the eigenfrequencies by comparing results with the amplitude gradient method of Baransky et al. (1985). Sandhu et al. (2018) expanded this work to high-latitude magnetometers from the IMAGE (International Monitor for Auroral Geomagnetic Effects) array (Luhr, 1994). Chi et al. (2013) created a similar method to look at the Mid-continent Magnetoseismic Chain magnetometer array (1.46 < L < 3.42). Their algorithm was designed to select out the fundamental frequency using a method depending on the stability of the estimated frequency. It stated it should not vary by more than 3 mHz in 10 min. This does not rule out higher harmonics, which may also show such stability. Indeed, a group of higher frequencies was found for the Glyndon-Cambridge pair of magnetometers (L = 3.3) at 37 mHz, on top of the fundamental at 7 mHz. Chi et al. (2013) did not use these observations as they gave plasma mass densities unreasonably low. However, they may have been observations of a higher harmonic. Lichtenberger et al. (2013) developed an automated cross-phase technique for the European quasi-Meridional Magnetometer Array, a part of the PLASMON space weather monitoring project. However, no details about the algorithm are given in their paper. The European quasi-Meridional Magnetometer Array website (http://geofizika.canet.hu/plasmon/) states that an automated algorithm, called Field Line Resonance Identification, is under development but implies it is not complete. This algorithm was also used by Jorgensen et al. (2017) for comparing ground-based measurements of plasmaspheric mass densities to the Dynamic Global Core Plasma Model, but again, no details of the algorithm were given.
These algorithms were limited as they could only estimate one harmonic at a time. Takahashi and Denton (2007) used multiple harmonics to determine the plasma mass density distribution from spacecraft data, but it has not been done with ground-based data. Ground-based data are advantageous because of the improved spatial and temporal coverage. In this paper, the methods of Berube et al. (2003) and Sandhu et al. (2018) have been adapted to create an algorithm that can detect multiple harmonics simultaneously. This should be possible as it has been shown that magnetic field lines can support many harmonics at the same time (e.g., Takahashi & McPherron, 1982). Hasegawa et al. (1983) showed that a broadband source could excite all harmonics on a field line that had a frequency within the range of the broadband source. Here we apply this algorithm to five pairs of IMAGE Magnetometers for September 2010 to create average distributions of eigenfrequencies. Variation in the average eigenfrequencies are investigated in latitude and MLT. For one pair, these average eigenfrequencies are then compared with eigenfrequencies derived using the wave equation of Singer et al. (1981) and a power law mass density model. This information is used to deduce which harmonics were present in the magnetosphere at this time. Finally, we consider the impact of adding a bulge in equatorial plasma density on the structure of the whole plasma mass density distribution.

Cross-Phase Estimation of Eigenfrequencies
The cross-phase technique can estimate the eigenfrequency of the magnetic field line at the midpoint of two magnetometer stations with a few assumptions (e.g., Clausen et al., 2008). Firstly, the magnetic field line is constantly driven by a broadband set of quasi-random fast mode waves. Secondly, the natural frequency of the magnetic field at the poleward station is slightly lower than the magnetic field at the equatorward station. This is primarily due to the field line length increasing with geomagnetic latitude. An exception to this occurs at the plasmapause. The sharp increase in plasma mass density moving equatorward causes the eigenfrequencies to decrease across this region. This would appear as a peak in cross-phase of opposite sign (e.g., Kale et al., 2007), where peaks in cross-phase are explained below. Finally, the amplitude spectra are assumed to have a single peak. Figure 1 illustrates the principles of the cross-phase technique. Model spectra are shown in Figure 1a. Subtracting the equatorward amplitude spectrum from the poleward amplitude spectrum yields a negative gradient in the amplitude difference spectrum centered on where the amplitude difference is 0 ( Figure 1b). The frequency at which this occurs is assumed to be the eigenfrequency (Baransky et al., 1985). Waters et al. (1991) found the phase spectra more reliable. Figure 1c shows model phase spectra. Subtracting them in a similar manner to the amplitudes gives the cross phase, shown in Figure 1d. Baransky et al. (1989) showed that the peak in cross-phase occurs at the eigenfrequency. Whether this is positive or negative depends on the order of subtraction. For example, Obana et al. (2015) subtracted spectra so that eigenfrequencies occurred at minima in cross phase, whereas Chi et al. (2000) subtracted spectra so that the eigenfrequencies occurred at the maxima in cross phase. The cross spectrum is then calculated using equation (2). F pol ( ) and F eq ( ) are the complex Fourier spectra of the magnetic time series from the poleward and equatorward stations, respectively. N 0 is the number of data points composing the spectrum, which should Journal of Geophysical Research: Space Physics 10.1029/2018JA025487 be equal for each spectrum. The argument of the cross spectrum gives the cross-phase, and the modulus gives the cross-power. Calculating the cross-phase in this way prevents aliasing errors when calculating the cross-phase. Also, if a station pair straddles the plasmapause, the cross-phase peak will be positive and the algorithm described below will remove those measurements.

Data Selection
This study was performed with magnetometers from the IMAGE magnetometer network. All of these magnetometers had 10-s time resolution or a Nyquist frequency of 50 mHz. The chosen pairs of magnetometer stations are shown in Table 1, along with their Corrected GeoMagnetic latitudes, longitudes, and meridional separation. The pairing criteria were chosen to minimize longitudinal separation to reduce any phase differences that could arise from azimuthally propagating waves. Menk et al. (1999) showed that a greater latitudinal separation will result in a larger cross-phase measurement but will also reduce the coherency between the signals. Menk et al. (2004) deduced the optimum separation was 110 km. Therefore, the pairings were chosen to fit these criteria as closely as possible, although choice is limited by the station locations. Menk et al. (1999) also notes that the phase change with latitude will be faster for narrower resonating structures. In that circumstance, a closer separation would be better. However, it is not possible to know this information beforehand. For the station pairings chosen, the meridional separation generally decreases with increasing latitude. This should not be an issue, as dipolar L shells are crossed more rapidly at higher latitudes so structures are likely to be narrower in latitude.

Automated Eigenfrequency Searching
The algorithm that searches for eigenfrequencies in magnetometer data has been adapted from Berube et al. (2003). It is demonstrated in Figure 2 for the MUO-PEL pair of magnetometers for 8 September 2010. Firstly, the B x (north-south) data are processed to linearly interpolate across data gaps. Secondly, the data are detrended using an 800-s boxcar window to remove very low frequencies that would otherwise dominate the spectrum (Ponomarenko et al., 2003). These processed data are shown in Figure 2a.
Next, the dynamic cross-spectrum is calculated using equation (2). The data are split into 40 min time windows that are incremented every 5 min. The Fast Fourier Transform for each station in the pair are calculated for each sliding window. The two spectra are then used to calculate the cross-spectrum, where the modulus and argument give the cross-power and cross-phase, respectively. The power ratio is defined as the power spectrum (amplitude spectrum multiplied by its complex conjugate) of the poleward station divided by the power spectrum of the equatorward station. Figures 2b, 2c, and 2d show the dynamic cross-power, cross-phase and power ratio, respectively. Figure 2c shows that a peak in the cross-phase is observed at approximately 10-12 mHz and spans from ∼02:00 UT to ∼12:00 UT. This is a target feature of the algorithm due to its consistently high cross-phase.

10.1029/2018JA025487
The cross-phase is then averaged in both time and frequency. This is to remove cross-phase values that are random and not part of a consistent feature. Each point in the dynamic cross phase is replaced with the local average of all the points in a box 60 min by 2 mHz centered on that point. This averaging is shown in Figure 2e. This is a similar approach to Berube et al. (2003) and Chi et al. (2013), though the box sizes are different. Berube et al. (2003) used a 10-mHz by 30-min box and Chi et al. (2013) used a 3.4-mHz by 50-min box. They were using magnetometers at lower latitudes with a much higher Nyquist frequency. Higher, more widely spread frequencies would therefore be expected. This study utilizes magnetometers at much higher latitudes, so the frequencies are expected to be much lower. Therefore, using a box too broad in frequency may result in features being unresolvable. It was found that this size was optimal for these data.
Next, the dominant frequencies are identified. For each 40-min window in Figure 2e, the mean and standard deviation over all frequencies are calculated. Only cross-phase values that are less than the mean minus one standard deviation are kept, and these values are shown in Figure 2f. Berube et al. (2003) used the mean plus two standard deviations with the cross-phase defined the opposite way around, but this threshold was considered too harsh for our data set and it was revised. Henceforth, these surviving data points are called the dominant averages.
A statistical t test is then applied to the dominant averages. This is the mean cross-phase divided by the standard deviation of the mean's constituents and multiplied by minus 1. These values are shown in Figure 2g.
If the t statistic of the dominant average is greater than 1, it is kept. This leaves the final cross-phase values, which are shown in Figure 2h.
The last step in the algorithm is to assign an uncertainty to these frequency values. A method based on the one used by Berube et al. (2003) was applied as follows. For each final measurement, its column was searched for the nearest frequency where the power ratio is equal to 1 and has a negative gradient. Interpolation across values is usually required between frequencies as the power ratio is rarely exactly equal to 1. The difference in frequency between the final cross-phase value and the nearest frequency where the power ratio is 1 defines the uncertainty in frequency. These uncertainties are shown as the black, vertical error bars in Figure 2h. Any final cross-phase values with an uncertainty of more than 1 mHz were removed. Figure 2 shows that the algorithm was successful in picking out the desired feature. It has identified multiple measurements that would be considered the dominant feature if selecting by eye. Very few measurements that might be considered noise have survived. Sixty examples from two different magnetometer pairs in the form of Figure 2 were examined by eye to see whether the algorithm was finding sensible features. The results were satisfactory in all cases. It should also be pointed out that for time periods where no visually identified cross-phase features were present, the algorithm did not produce any measurements, meaning the algorithm will only identify cross-phase values if there is a feature there to measure. A further point is to notice that the time interval of the identified feature is during a quiet period of geomagnetic activity ( Figure 2a). The cross-phase technique relies on the assumption of a broadband driver, that is, there is power at a wide range of frequencies. Large amplitude waves at a specific frequency will break this assumption, so the cross-phase technique will become ineffective. This is why the feature exists where the B x amplitudes are smallest. The intense wave activity after 14 UT is unlikely to be locally resonant and may be driven by substorms.
Further examples of the algorithm are shown in Figure 3. These four examples only use subplots a, c, and h (without error bars) from Figure 2, which correspond respectively to i, ii, and iii in Figure 3. The error bars were removed in the these plots as they were all less than 1 mHz. Figure 3ai shows data from PEL-MUO on 1 September 2010, with a very clear feature that has been identified in Figure 3aiii along part of its length. This example shows that the algorithm does not always select out the entire feature but it is efficient at removing noisy values. This was considered an acceptable compromise. Multiple frequencies can be seen identified at ∼06:00 UT, showing an example where the algorithm does not just detect the fundamental harmonic. Figure 3bi shows data from the same pair (PEL-MUO) on the 11 September 2010. Three or even four harmonics can be seen simultaneously at ∼09:00 UT in Figure 3bii. Figure 3ci shows data from a station pair (MAS-SOR) at a higher L shell. Two clear harmonic structures can be seen in Figure 3cii, lasting for almost 12 hr in the case of the higher frequency harmonic. Finally, Figure 3di shows data from a lower-latitude pair with two very clear harmonic bands in Figure 3dii. The average cross-phase for this example is much higher than the other examples. This is consistent with both Berube et al. (2003) and Chi et al. (2013), who showed examples of cross-phase spectrograms with cross-phase values much higher than 60 ∘ , which is the saturation limit on the color scale in all subplots of Figure 3. Data in both of these studies were from lower L shells, suggesting that resonant features at lower latitudes have a more rapid change in cross-phase with latitude. However, this effect is left to a future study. It can be seen from these examples that the algorithm is consistently selecting out features that might be identified by eye.

Distribution in MLT
This algorithm was applied to all B x data obtained during September 2010 for the magnetometer pairs given in Table 1. September was chosen as it was at equinox, so would remove complications from field lines that could be crossing the terminator and produce quarter wave modes (Allan & Knox, 1979;Obana et al., 2008Obana et al., , 2015. Although Chi et al. (2013) deduced that quarter wave modes were insignificant in their statistical study after comparing equinoctial and nonequinoctial results, the precaution was taken anyway. September 2010 was a geomagnetically quiet month, with no sudden storm commencements.  The frequencies observed at night are also similar to those observed in the day, though less distinct.
Assuming the lowest frequency peak is the fundamental mode (first harmonic), then subsequent peaks might be the second, third harmonics, and so forth. However, at this stage, we cannot be certain of the harmonic numbers of each peak. Therefore, we adopt the terminology of f a for the fundamental (an assumption we retain throughout this paper), then f b , f c , … , where the values of a, b, c, … will be discussed in more detail in section 5. Figures 4d-4f show that the algorithm identified f b more often than the fundamental. This does not mean that f b contains the most power, only that it had a larger phase difference between the two magnetometer stations so the algorithm was more likely to detect it. Binning the measurements by cross power (not shown) shows that the cross power decreases with frequency and that the fundamental contains the most power.
Asymmetry can also be seen in both frequency and occurrence. The frequency of the peaks is lower in the dusk sector (Figure 4f ) than the corresponding peaks in the dawn sector (Figure 4d). There is an asymmetry in occurrence between Figures 4d and 4f and also Figures 4c and 4g. This asymmetry in dawn/dusk occurrence was observed by Takahashi et al. (2016) in spacecraft data. More measurements were obtained in the morning/dawn sectors than the afternoon/dusk sectors.

Distribution With Geomagnetic Latitude
These data were then reprocessed as follows. Instead of using 3-hr-wide MLT bins, the data for each station pair in Table 1 were split into twenty-four 1-hr-wide MLT bins, and histograms and peak values found again.
Observations for all five station pairs in Table 1 are shown in Figure 5 (where each pair is labeled on the right with its dipolar L value). Figure 5 shows the variation in estimated eigenfrequencies with MLT and having used the five pairs, also in latitude. The color shows the number of measurements (occurrence) in each bin, and the black crosses show the peaks in each histogram. The red lines show peak frequencies that exist over multiple MLT windows. These red lines show how the eigenfrequencies vary with MLT.
For For all of the stations, there are consistent frequencies (see red lines in Figure 5) that last for up to and even longer than 12 hr of MLT. These frequencies maximize in the dawn-to-noon MLT sector.

Calculating Model Eigenfrequencies
In this section, we use the data from the PEL-MUO pair derived in section 3 to determine the plasma mass density distribution. This was achieved by solving the wave equation of Singer et al. (1981; see also Waters et al., 1995Waters et al., , 1996 to produce model frequencies.
Comparing these model frequencies to the data in section 3 allowed us to then deduce the harmonic numbers.

Station Pair Selection for Modeling
The station pair chosen for this modeling work was PEL-MUO, specifically the 12 MLT data set. This was for two reasons: (1) The data for this station pair had the clearest frequency peaks of all the MLT bins from all five station pairs. The fundamental frequency, f a , of these data was 3.71 mHz.
(2) This station pair was chosen due to its position relative to the location of the plasmapause. The plasma population differs on either side of this boundary. Within the plasmapause is cold dense plasma originating from the ionosphere whose motion is dominated by the corotation of the inner magnetosphere. Outside the plasmapause is hotter, less dense plasma whose motion is dominated by the convection of the Dungey cycle. Figure 6a shows the variation in the Dst index throughout September 2010. It shows there were no major storms during this period and that this month was geomagnetically quiet. Figure 6b shows the location of the plasmapause calculated using  Figure 6b shows that the plasmapause remained at a steady position of L ≈ 4 during September 2010 and that the PEL-MUO pair and the MAS-SOR pair were consistently located in the plasmatrough. Selecting one of these pairs therefore removed complications that could arise from the plasmapause passing back and forth over the magnetometers. The PEL-MUO pair was chosen over the MAS-SOR pair due to the better clarity of its frequency peaks. Cummings et al. (1969) calculated the first estimation of eigenfrequencies for a dipolar magnetic field. They found that the eigenfrequencies for toroidal and poloidal modes differ, with the poloidal eigenfrequencies being lower than the toroidal for low harmonic numbers. This is due to the variation along the field line of the normalized separation of adjacent field lines being different for the azimuthal and radial directions. The measurements in this study are from the north-south component of the magnetometer data. Assuming a 90 ∘ rotation of the magnetic field vector through the ionosphere (e.g., Hughes & Southwood, 1976), then pertubations in this component correspond to a toroidal mode in the magnetosphere. In general, resonant waves will not be purely toroidal, but this is a reasonable approximation for a symmetric magnetosphere (Wright & Elsden, 2016). Singer et al. (1981) expanded the theory to magnetic fields with arbitrary geometry by developing a suitable second-order differential wave equation (equation (3)). is the perpendicular displacement of the field line, h is a parameter describing the separation of the field lines (which maximizes at one near the magnetic equator), B 0 is the ambient magnetic field, is the plasma mass density, 0 is the magnetic permeability of free space, ds is an element of length along the field line, and is the frequency of oscillation. Equation (3) can be solved with an appropriate magnetic field model and plasma mass density model to give a theoretical set of eigenfrequencies.

Solving the MHD Wave Equation With a Power Law Plasma Model
The magnetic field was modeled using the T96 model (Tsyganenko, 1995(Tsyganenko, , 1996 to calculate the strength of the magnetic field, B 0 , and the azimuthal separation of the field lines, h . This is calculated from the normal  Wanliss and Showalter (2006) showed that this was an appropriate substitution. The plasma mass density distribution was modeled using the power law form (equation (1)). Equation (3) was used with the T96 model to calculate these densities along the field line.
Equation (3) is solved using a Runge-Kutta-Gill (RKG) numerical shooting method for a given value of the exponent p from equation (1). However, the RKG method can only be used to solve first-order ordinary

Journal of Geophysical Research: Space Physics
10.1029/2018JA025487 differential equations, whereas equation (3) is second order. Therefore, equation (3) was rewritten as a matrix of two first-order ordinary differential equations (equation (6)) using the substitutions given by equations (4) and (5). The definition of the Alfvén velocity, v A , has been substituted in. Variables marked with primes denote derivatives with respect to s.
Equation (6) was solved with the initial conditions of zero displacement at one ionosphere (perfect conductor) and an arbitrary gradient. A range of frequencies were evaluated and those that caused the displacement to be 0 at the other end of the field line were considered to be the eigenfrequencies. Initially, a minimization routine was used to change eq until the fundamental frequency was 3.71 mHz. The frequency ratio is defined as the frequency of a harmonic divided by its fundamental, and it is this quantity that will be used to determine the harmonic number. With this value of eq , the equation was solved again to find all the other eigenfrequencies that caused the displacement to be 0 at the other end of the field line for a range of values of p. Note that eq will vary with p. Figure 7 shows how the eigenfrequencies change with the index p for the MUO-PEL station pair. Figures 7a-7d show how the plasma distribution changes shape along the field line as p is varied from −6 to 6. The length of the field line is L, which here is 71,586 km. For a negative index, the majority of the plasma lies in the equatorial plane, decreasing to a minimum at the ionosphere (the bounds of the plot). For an index of 0, the plasma mass density is constant along the field lines. For a positive index, the majority of the plasma is toward the ionosphere, with a minimum in plasma density at the equator. These distributions have been scaled using eq so that the fundamental frequency is the same as that from the data (3.71 mHz) and that value of eq is shown by the dotted lines in Figures 7a-7d. The value of eq decreases as p increases.

Discussion
The results in section 3 show that the observed frequencies exhibit several spatial dependences in terms of occurrence and frequency. Section 4 also suggests some interesting implications from the harmonic spacing. We now discuss these features. Figures 4 and 5 show a greater occurrence of measurements in the dayside magnetosphere, which is unsurprising. Takahashi and Anderson (1992) showed that there is more ultralow frequency wave power on the dayside than the nightside using data from the AMPTE/CCE spacecraft, though this power would not necessarily contribute to standing Alfvén waves. The greater occurrence in the dayside could be due to the greater ionospheric conductivity, which increases the chance of standing waves forming. The increased occurrence may also be due to a greater number or intensity of sources on the dayside, principally the solar wind. The lower occurrence in the nightside may also be due to the changing polarization of the resonances. The likeliest polarization at night might not be toroidal, in which case, taking observations of the B x component would not be optimal.

Variation of Eigenfrequencies With MLT and Latitude
Another observation was that the nightside eigenfrequencies were similar to the dayside eigenfrequencies ( Figure 4). Field lines are longer on the nightside as they are stretched into the magnetotail, which should reduce the eigenfrequencies. However, the field lines map to different plasma regimes in the night and day

10.1029/2018JA025487
which would also influence the eigenfrequencies. The nightside field lines may be more dynamic too, due to larger changes in field configuration. This would also influence the eigenfrequencies.
Figures 4 and 5 showed that there were more observations in the dawn sector than the dusk sector. This agrees with observations by Lessard et al. (1999). This suggests that field lines in the dawn sector are more likely to resonate. This feature could be the result of asymmetry in the Kelvin-Helmholtz instability (KHI) on the magnetopause flanks. Henry et al. (2017) showed, using Time History of Events and Macroscale Interactions during Substorms data from 2007 to 2013, that there is an asymmetry favoring the dawn side in the occurrence of Kelvin-Helmholtz waves on the magnetopause. However, Takahashi et al. (2016) deduced that the source of this asymmetry was internal to the magnetosphere. This was based on a statistical study of Time History of Events and Macroscale Interactions during Substorms data from 2008 to 2014, which determined that the occurrence rate of Pc5s was independent of the interplanetary magnetic field orientation. Therefore, a driver external to the magnetosphere, such as the KHI, would not be responsible for this asymmetry. Instead, Takahashi et al. (2016) propose that the asymmetry is due to the differing radial distribution of plasma with local time. Thus, the source of this asymmetry is not currently fully resolved.
Figures 4 and 5 show that the eigenfrequencies are lower in the afternoon sector. The higher frequencies in the morning could be due to the asymmetric shape of the plasmasphere, which has been noted before (e.g., Chappell et al., 1971;Chi et al., 2013;Mathie et al., 1999;Poulter et al., 1984;Sandhu et al., 2018;Takahashi & McPherron, 1982;Warner & Orr, 1979). The plasma mass density is higher in the afternoon sector due to the refilling of flux tubes from the ionosphere throughout the day, as well as the average ion mass increasing (Sandhu et al., 2016a(Sandhu et al., , 2016b. Figure 5 shows that all the eigenfrequencies decreased with increasing latitude. This is expected as the field line length increases with latitude. Resolving the different harmonics was also more difficult at higher latitudes. This could be because the outer magnetosphere is more variable than the inner magnetosphere, resulting in a greater eigenfrequency variation. It may also be due to the harmonics decreasing to lower frequencies, making them harder to resolve in the cross-spectrum. The dipolar L shell value for LYR-NAL (17.34) suggests that this station pair will often be on open field lines.
In these circumstances, we would not expect standing waves to form, and thus, we would not expect the cross-phase technique to identify the eigenfrequencies. However, Mathie et al. (1999) showed that field line resonances can exist at sufficiently low frequencies on high-latitude field lines if they are closed. The observations taken by this station pair (Figure 5a) suggest that they spent a significant proportion of September 2010 on closed field lines. Figure 7e shows that the separation of the harmonics decreases as p increases. Increasing p effectively removes plasma from the equatorial plane, increasing the local Alfvén velocity. This in turn increases the eigenfrequencies. The first harmonic (fundamental) is affected most by this change because the maximum displacement for this mode occurs in the equatorial plane. This results in the frequency ratios decreasing with p. The frequency ratios are also not integers or half integers of the first harmonic (Cummings et al., 1969). This is due to the varying magnetic field strength and plasma distribution along the field line (e.g., Denton et al., 2004).

The Absence of Even Harmonics
From Figure 7e, we can see that the spacing of the harmonics in the data are much greater than those in the model. If we assume that all harmonics are present in the data, then matching with the model eigenfrequencies gives p < −6. This is physically unrealistic and implies that the plasma distribution tends to a minimum at the top of the ionosphere.
However, if we assume only odd modes are present in the data and compare again, we get values of p ≈ 4. Specifically, the mean value of p calculated from the third, fifth, and seventh harmonics is p = 4.27 ± 0.39. Figure 7e shows horizontal double-headed black arrows, which show the uncertainty in p for each harmonic. This shows that p could exist in a much wider range of approximately 1 < p < 5.5. However, with each harmonic added, the width of this uncertainty decreases due to the increase in the gradient of the model harmonics with p. Having more harmonics therefore allows us to constrain the value of p with greater certainty. This mean value of p provides a more realistic plasma mass density distribution. The absence of even harmonics has been noted before by Lanzerotti and Fukunishi (1974) in the region L ≈ 4, although they noted that there was a scarcity of measurements from other latitudes at the time.

10.1029/2018JA025487
The lack in occurrence of even harmonics is not expected to be a result of observational biases in the cross-phase technique. Arguably, there is a second harmonic present in Figures 5c and 5d at ∼9 UT, but the occurrence and duration are too low for this feature to be considered significant in the context of the surrounding observations. The phase profiles of odd and even resonating structures are expected to appear the same when viewed from the ground, because both even and odd modes have a node in the electric field in the ionosphere. Therefore, the cross-phase technique should not be able to distinguish between odd and even modes.
These observations suggest that the fast mode driving the resonances is most likely to have a field-aligned structure that is symmetric about the equatorial plane. The field-aligned structure of the fast mode is important in determining whether it will excite a field line eigenmode (Southwood & Kivelson, 1986). A symmetric fast mode may be more likely to form because the displacement of the field line is a maximum at the magnetic equator due to the minimum in the field strength there (Lanzerotti & Fukunishi, 1974;Poulter et al., 1984). Instabilities such as the KHI may therefore be more efficient at driving the odd modes (Lanzerotti & Fukunishi, 1974;Takahashi & McPherron, 1982). Claudepierre et al. (2010) have also conducted simulations showing the effect of solar wind buffeting on the excitation of field line resonances. They found that only the odd modes were excited in their simulations and attributed this to the forcing from the solar wind being symmetric about the equatorial plane. Poulter et al. (1984) suggested even modes would not be able to form efficiently. Even modes can be generated by instabilities in the ring current ions (Southwood, 1976). Figure 6 shows us that the Dst index remained low throughout September 2010. Therefore, the ring current was quiet at L = 5.44 during this period, so it was unlikely even modes would be observed.
We now review values of p from the literature to compare to our estimates. Early studies used whistler measurements (e.g., Angerami & Carpenter, 1966;Carpenter & Smith, 1964) to determine p = 3 inside the plasmasphere and p = 4 in the plasmatrough. These values agreed with the diffusive equilibrium model of Angerami and Thomas (1964) and the collisionless model of Eviatar et al. (1964), respectively.
These results agreed well with satellite observations. Cummings et al. (1969) used a plasma density model of R −4 (i.e., p = 4) in their calculations of harmonic frequencies in a dipole field at geostationary orbit and found a good agreement with data from the ATS-1 satellite. Orr and Matthew (1971) extended this theory to cover a wider range of L shells, making the assumption that the plasma distribution varied as R −3 inside the plasmasphere and R −4 outside, although they struggled to determine which harmonic they were observing. They also showed that the wave period depended only very weakly on the power index p. Chappell et al. (1971) also found that the diffusive equilibrium and collisionless models fitted data from the OGO-5 satellite inside and outside the plasmapause, respectively, but also found that the plasma density varies as R −4 within the plasmapause in the afternoon and dusk sectors.
Other authors that used values of p = 4 in the plasmatrough for their calculations include Warner and Orr (1979), Poulter et al. (1984), Takahashi and Anderson (1992), Waters et al. (1995Waters et al. ( , 1996, Chi et al. (2000), Menk et al. (2004), and Menk et al. (2014). Poulter et al. (1984) used measurements from the STARE radar (Greenwald et al., 1978) and Slope Point Radars (Keys & Johnston, 1979) to show that the value of p was higher in the morning than the afternoon, those values being p = 4.74 and p = 3.08, respectively.
More recent observations have shown that the value of p could be more variable. Menk et al. (1999) used p = 3 in their calculations of eq using the Taylor and Walker (1984) perturbation model, employing data from ground-based magnetometers at low latitudes. They showed that even if the value of p varies significantly (e.g., from 1 to 6), the change in estimated equatorial mass density is only a few percent (this is evident in Figures 7a-7d). Therefore, large changes in the value of p could be expected from only small changes in geomagnetic conditions. Menk et al. (2004) also used p = 3 and p = 4 inside and outside the plasmapause, respectively, and found the value of eq insensitive to the value of p. Denton et al. (2006) found that p ≈ 1 in their study between 4 < L < 5, lower than values deduced by other authors. A value of 1 for p was also used by Takahashi et al. (2010), but these measurement were based on Geostationary Operational Environmental Satellites (GOES) measurements at geostationary orbit. This is at a higher L shell, so the plasma distribution should be expected to be the same as at L = 5.44. However, most authors appear to have found that p = 4 is a reasonable value to describe the plasma mass density distribution in the plasmatrough. Also, to our knowledge, nobody has reported values of p < −6. This supports our hypothesis that only the odd modes are present in the data.
We can further test which harmonics are present by examining the spacing of the harmonic modes. We use a method adapted from Schulz (1996) by Takahashi and Denton (2007). In table 2 of Takahashi and Denton (2007), they converted the frequency quanta (difference in frequency between neighboring harmonics), first calculated by Schulz (1996), into a more useable form. These quanta were calculated for a dipole magnetic field with a power law distribution. Takahashi and Denton (2007) normalized the frequency separation of two harmonics by the frequency of the third harmonic for a range of values of p ( in Takahashi & Denton, 2007). The data employed in Figure 7 from the PEL-MUO pair can be compared to this simple table.
If we make the assumption that all modes are present in the data used in Figure 7e, (so f a = f 1 , f b = f 2 , … ) we calculate (f 2 − f 1 )∕f 3 = 0.426 and (f 3 − f 2 )∕f 3 = 0.406; f 4 was not used for this calculation as it was not considered sufficiently reliable due to a small quantity of data. Comparing these two values to table 2 in Takahashi and Denton (2007) suggests that p < −6. This is the same conclusion arrived at when comparing the data to the model frequencies. Schulz (1996) deduces that the harmonics are approximately evenly spaced. Therefore, we can calculate (f 2 − f 1 )∕f 3 and (f 4 − f 3 )∕f 3 with the assumption that only the odd modes are present, so f a = f 1 , f b = f 3 , and f c = f 5 . This gives (f 2 − f 1 )∕f 3 = 0.358, which (using table 2 of Takahashi & Denton, 2007) gives a value for 3 < p < 4.
(f 4 − f 3 )∕f 3 = 0.338, which gives 5 < p < 6, although the former calculation is more reliable due to better data coverage. These values of p are in agreement when comparing the model odd harmonics in Figure 7 to the harmonics from the data, as discussed above. Examination of the harmonic separation further supports the suggestion that only the odd modes are present in this data set.
A further observation is that Takahashi and Denton (2007) and Takahashi et al. (2010) noted in their studies using GOES that the third harmonic was detected most frequently. If the even harmonics were missing, then the dominant peaks (f b ) in Figures 4c-4g would represent the third harmonic and this would agree with the findings of Takahashi and Denton (2007) and Takahashi et al. (2010).
In summary, the value of p = 4.27 ± 0.39 deduced here agrees with values established by many (if not all) previous authors. It also agrees with values of p usually associated with the plasmatrough, which is where the PEL-MUO station pair was located according to the plasmapause model of O'Brien & Moldwin, (2003; see Figure 6). Based on the high occurrence of f b in Figure 4, comparing the data with model frequencies, and Journal of Geophysical Research: Space Physics 10.1029/2018JA025487 examination of the separation of the harmonic modes, we can conclude that it is most likely that the even modes were not present at L = 5.44 during September 2010. This does not mean that the even modes can never be present in ground magnetometer data (Takahashi & Denton, 2007, show clear detections of even harmonics with GOES at geostationary orbit), only that they were absent during September 2010.
With the harmonic numbers now assumed to be known, variations in the eigenfrequencies with latitude are examined further in Figure 8. Peak values (black crosses) from the 11-12 MLT bin for each station pair in Figure 5 have been plotted as a function of their latitude using circular markers. The color represents the occurrence value of that point. Crosses located on red lines in Figure 5 have been joined between the station pairs, showing how the different harmonics change with latitude. Note that the four frequencies between 10 and 30 mHz for the LYR-NAL pair are anomalous values that survived the selection algorithm and have very low occurrence values.
It can be seen that the spacing between the harmonics increases significantly at lower latitudes and that the size of that spacing is not an integer or half integer of the fundamental frequency. The large spacing between the harmonics agrees with our assessment that only the odd harmonics are present. Schulz (1996) deduced that the harmonic spacing for a power law distribution was approximately equal between the adjacent higher harmonics, though not for the fundamental and the second. Figure 7e shows that this is not the case; that is, f 7 − f 5 < f 5 − f 3 . There is also an uncertainty in the value of p.

Modeling an Equatorial Bulge
Some authors have observed a peak in plasma mass density at the equator based on satellite measurements (e.g., Denton et al., 2004Denton et al., , 2006Denton et al., , 2009Takahashi & Denton, 2007;Takahashi et al., 2004). Sandhu et al. (2016aSandhu et al. ( , 2016b) modeled this plasma distribution using electron density and average ion mass measurements from Cluster. Although Denton et al. (2006) suggests that an equatorial bulge in mass density exists at L > 6, we will now investigate the impact such an equatorial bulge would have on our PEL-MUO data set.
Therefore, we have experimented with a model of the form of equation (7). The left-hand term is the usual power law distribution, and the right-hand term is a Gaussian function to describe the equatorial bulge. It is based on the Gaussian function used to describe the bulge in electron density from Sandhu et al. (2016b); b is the equatorial mass density contribution to the bulge, so the total equatorial mass density would be eq + b .
This model has three free parameters: eq , p, and b . However, if we have at least three harmonics, we can find an exact solution. We define the misfit value as the root-mean-square error of the difference in the harmonics from the data in Figure 4e and those calculated using the model. The three free parameters of the model are adjusted using a Nelder-Mead minimization routine to minimize the misfit value. The harmonic numbers are assumed to be known from the previous analysis, assuming the even harmonics are missing. The fundamental frequency is no longer fixed, although the minimization routine will result in it becoming equal to the fundamental in the data regardless. The result of this analysis gives eq = 216 amu/cm 3 , p = 4.45, and b = −0.0003 amu/cm 3 . A plasma mass density distribution of the form of equation (7) with these values will give eigenfrequencies that match those in Figure 4e. This analysis shows that there was no equatorial bulge present along the field line at L = 5.44 during September 2010.
This analysis can also be done for the power law model, by setting b = 0 amu/cm 3 . Just using a power law model, eq = 216 amu/cm 3 and p = 4.45. Note that the value of p = 4.27±0.39 from section 5.2 was the mean value of p, not the optimized value. The two calculated plasma mass density distributions are equal, leading us to conclude that there was no equatorial bulge. This agrees with Denton et al. (2006), who deduced that the bulge only existed at L > 6, whereas our field line is at L = 5.44.

Conclusions
In this study, an algorithm has been developed to automate cross-phase processing of ground magnetometer data for the purpose of investigating the factors influencing the eigenfrequencies of magnetic field lines. This was modified from previous algorithms by Berube et al. (2003) and Sandhu et al. (2018). This algorithm

10.1029/2018JA025487
can be applied to large data sets and is capable of detecting multiple harmonics from the data, a feature previous search algorithms were not designed for. It has been used in a study to investigate data from five magnetometer pairs during September 2010.
This data set was broken down by MLT, and there were clear bands of common eigenfrequencies throughout the daylight hours. These eigenfrequencies were higher in the morning hours than the afternoon. This was attributed to an increasing plasma mass density as the day progressed, and this had the effect of reducing the eigenfrequencies. More detections were made in the daylight hours, which was attributed to a weaker ionosphere at night that would be less capable of maintaining standing waves. Standing waves may also be more likely in the daylight hours due to increased ultralow frequency wave power (Takahashi & Anderson, 1992). There were also more observations in the morning than the afternoon. Variations of eigenfrequencies with latitude were also considered using data from five station pairs, which showed that the eigenfrequencies of all harmonics decreased with latitude. The algorithm was also capable of identifying features for the NAL-LYR pair of magnetometers on Svalbard, suggesting that closed field lines at these latitudes were also common.
Next, data from a 3-hr-wide MLT bin centered on noon for the MUO-PEL pair of magnetometers was used for comparison to modeling work. Equation (3) (Singer et al., 1981) was solved using an RKG routine, using the T96 magnetic field model (Tsyganenko, 1995(Tsyganenko, , 1996 and a range of power law plasma distributions, in order to deduce the theoretical eigenfrequencies of the magnetic field lines at the MUO-PEL station pair. When using the power law plasma mass density model, comparisons suggested that the even modes were missing in the data; otherwise, unphysical plasma distributions arose. We therefore assumed that even modes do not form efficiently at this location in the magnetosphere. This may be because even modes are excited mostly by ring current instabilities (Southwood, 1976), which may not have been present. Odd modes are more likely to be formed by the KHI as they cause a maximum in displacement at the equator. This suggestion was also confirmed by examining the spacing of the harmonic modes and the high occurrence of the third harmonic.
Finally, we considered the impact on the whole plasma distribution when a bulge in mass density was introduced at the equator. We fitted a plasma mass density model with an additional Gaussian term to represent the bulge then used a minimization routine to determine the models parameters. We concluded that there was no equatorial bulge in plasma mass density and that a simple power law was sufficient to describe the plasma mass density distribution.
In this study, only one station pair at equinox has been studied in detail, and it was deduced that the even modes did not form at this time. It was also established that this station pair mapped to the plasmatrough. However, this assumption may not be valid for other latitudes or at other local times. Work is now underway to investigate this assumption using a much wider range of data. However, this methodology indicates that the most probable harmonic number of a resonating field line can be determined in ground magnetometer data for a specific event, if enough data are used to construct the long-term behavior of the field line. It has also shown that the plasma mass density distribution can be calculated and not just assumed, from having knowledge of multiple eigenfrequencies. Future work will look for examples of times that include even modes to see what impact they have on the overall distribution and what might cause them to occur. The impact of these effects, however, is left to future study.