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

2.1 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₀ is the number of data points composing the spectrum, which should 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.

(2)

2.2 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.

2.3 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.

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.

3.1 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., 2008, 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 measurements were binned by MLT. The MLT of the midpoint of the magnetometers when each spectrum was taken was used. The bin width was 3 hr, giving eight bins. These occurrence histograms are shown in Figure 4, which only includes results for the MUO-PEL pair of magnetometers. The histograms show the distributions of the estimated eigenfrequencies. Black lines are a five-bin moving average. It should be noted that one measurement does not correspond to one standing Alfvén wave. The presence of one standing wave can yield many measurements, so longer-lasting features or features with a greater spread in frequency are represented by more measurements in the histograms.

Figure 4 shows the technique detects more standing waves in daylight hours. There is a higher occurrence of measurements in Figures 4d, 4e, and 4f than in Figures 4a, 4b, and 4h. This was also observed by Chi et al. (2013). Observations are still apparent at night but are fewer in number and less structured in frequency. Data in the daylight hours show clear discrete peaks in frequency, those being most prominent in Figures 4d–4f. 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.

3.2 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 Figures 5c–5e (the station pairs at the lowest latitudes), there is a clear grouping of measurements into bands of frequency during the daylight hours. At higher latitudes, the eigenfrequency estimates decrease and the distinct groups of frequencies become harder to resolve. This can be seen in Figures 5a and 5b. There are also measurements for the LYR-NAL station pair (Figure 5a), despite their proximity to the open-closed field line boundary. Station pairs on open field lines would not be expected to identify standing waves. Johnsen and Lorentzen (2012) found from a 15-season statistical study of meridian scanning photometer data that the mean location of open-closed field line boundary was 75.4° magnetic latitude. This means that the LYR station will be on open field lines approximately half of the time. Note that data for the first 9 days of the month were missing for the Hankasalmi station, which will have reduced the occurrence values for the OUJ-HAN pair by approximately a third.

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.

4.1 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 the Dst model of O'Brien and Moldwin (2003). The dipolar L shells of the four lowest-latitude station pairs are shown as dotted lines. 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.

4.2 Solving the MHD Wave Equation With a Power Law Plasma Model

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₀ is the ambient magnetic field, ρ is the plasma mass density, μ₀ 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.

(3)

The magnetic field was modeled using the T96 model (Tsyganenko, 1995, 1996) to calculate the strength of the magnetic field, B₀, and the azimuthal separation of the field lines, h_α. This is calculated from the normal separation of the field lines in the azimuthal direction, where the equatorial footprints of two closely spaced field lines have a small separation. The position of the midpoint of the magnetometer is used for the field line tracing. For the model arguments, 15 September 2010 was used for the date and the time was calculated from the MLT. Average values for September 2010 of the following quantities were used, those being: P_dyn=1.479 nPa, Dst =− 6.635 nT, B_y=−1.385 nT, and B_z=0.145 nT. These values were derived from the OMNI database, where the SymH index was used in place of the Dst index for ease of calculation. 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 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.

(4)

(5)

(6)

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.

4.3 Comparing the Model Values to the Noon Data at L = 5.44

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.

Figure 7e shows how the modeled eigenfrequencies, f_R, vary as a function of p (black lines). The fundamental frequency is the same for all values of p to provide a consistent reference level. The colored lines show four possible and unknown harmonics derived from Figure 4e: f_a–f_d. Their frequencies are written on the left for convenience. The surrounding colored bars represent their error bars, where the error represents the full width half maximum of those frequencies. The dashed vertical black lines and horizontal double-headed black arrows will be explained in section 5.2.

5.1 Variation of Eigenfrequencies With MLT and Latitude

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.

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 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, 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.

5.2 The Absence of Even Harmonics

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).

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.

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⁻⁴ (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⁻³ inside the plasmasphere and R⁻⁴ 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⁻⁴ 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. (1995, 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₁, f_b=f₂, …) we calculate (f₂−f₁)/f₃=0.426 and (f₃−f₂)/f₃=0.406; f₄ 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₂−f₁)/f₃ and (f₄−f₃)/f₃ with the assumption that only the odd modes are present, so f_a=f₁, f_b=f₃, and f_c=f₅. This gives (f₂−f₁)/f₃=0.358, which (using table 2 of Takahashi & Denton, 2007) gives a value for 3 < p < 4. (f₄−f₃)/f₃=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 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.

5.3 Modeling an Equatorial Bulge

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₇−f₅<f₅−f₃. There is also an uncertainty in the value of p.

Some authors have observed a peak in plasma mass density at the equator based on satellite measurements (e.g., Denton et al., 2004, 2006, 2009; Takahashi & Denton, 2007; Takahashi et al., 2004). Sandhu et al. (2016a, 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.

(7)

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³, p = 4.45, and ρ_b=−0.0003 amu/cm³. 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³. Just using a power law model, ρ_eq=216 amu/cm³ 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.