Identifying ULF Wave Eigenfrequencies in SuperDARN Backscatter Using a Lomb‐Scargle Cross‐Phase Analysis

The eigenfrequencies of standing Alfvén waves on closed magnetospheric field lines can be estimated using the cross‐phase technique. These eigenfrequencies can be used to monitor the plasma mass density distribution along the field line. So far, this has only been applied to ground‐based magnetometer data. The Super Dual Auroral Radar Network (SuperDARN) radars offer some benefits over magnetometers. They provide greater spatial resolution and coverage, as well as direct sensing above the E region ionosphere, which screens ultralow frequency (ULF) waves from the ground. However, there are significant data quality issues. These include the uncertain origin of radar backscatter, uneven sampling of data due to data gaps, and inaccurate fitting to the autocorrelation functions. Artificial backscatter from an ionospheric heater has been used to remove the uncertainty in backscatter location. We have developed a Lomb‐Scargle cross‐phase analysis for application to discontinuous radar data. The First Principles Fitting Methodology has been used to improve the fitted data products derived from the autocorrelation functions. Using these techniques, we have shown that it is possible to measure eigenfrequencies with SuperDARN data, and we have verified an example using ground‐based magnetometer data. Finally, we have demonstrated that the eigenfrequency signature in this example was caused by a broadband source of energy.


Introduction
Closed magnetospheric field lines can host standing Alfvén waves and therefore have a natural set of frequencies, or eigenfrequencies. These eigenfrequencies are determined by the plasma mass density distribution, the topology of the magnetic field along the flux tube, and the polarization of the standing wave. Therefore, knowledge of the magnetospheric eigenfrequencies can be useful in determining the properties of near-Earth space. Standing waves are usually assumed to be toroidal or poloidal in polarization. These polarizations refer to plasma displacement in the azimuthal (east-west in ionosphere) or radial (north-south in ionosphere) directions. However, in a realistic nonaxisymmetric magnetosphere, the standing waves can have an intermediate polarization (e.g., Elsden & Wright, 2018). In this study, we are most interested in the toroidal mode because they tend to have larger-scale sizes, be more stable, and have a consistent poleward phase. Mager and Klimushkin (2006) have shown that over time the poloidal oscillations will change their polarization to being toroidal.
In the classical picture of a field line resonance (FLR) as described by Southwood (1974) and Chen and Hasegawa (1974), a field line is driven by an evanescent mode originating from a surface wave on the magnetopause, which couples to a field line with a matching Alfvén wave eigenfrequency. Such a discrete, driven FLR is characterized by a peak in wave power and a 180 • phase shift around the resonance location (e.g., Plaschke et al., 2008;Walker et al., 1979). More generally, the ensemble of magnetospheric field lines provide an Alfvén eigenfrequency continuum, which gives rise to a continuum of standing Alfvén waves if the source has a wide-frequency band. The resonant characteristics at each frequency in this continuum are equivalent to those of a FLR of the same frequency. These characteristics can be used to measure the local eigenfrequencies.
Eigenfrequencies can be measured from the ground using the cross-phase technique of Waters et al. (1991). This utilizes two latitudinally spaced magnetometers and calculates the cross-phase spectrum between 10.1029/2018JA025859 them. A minimum or maximum (depending on the order of subtraction of the spectra) in the cross-phase spectrum occurs at the eigenfrequencies. The technique requires the presence of a broadband source as described above to excite the field lines at each station at their local eigenfrequency (Menk et al., 2004;Waters et al., 1991). Claudepierre et al. (2010), Ellington et al. (2016), and Elsden and Wright (2018) have all shown that such a broadband source can be used to excite field lines using simulations. Many studies have utilized the cross-phase technique to study magnetospheric ultralow frequency (ULF) waves (e.g., Clausen et al., 2008;Kawano et al., 2002;Sandhu et al., 2018;Waters et al., 1995). If multiple harmonics of a field line are known, then the plasma mass density distribution along that field line can be estimated (e.g., Denton et al., 2001;Price et al., 1999;Takahashi & Denton, 2007). Recently, Wharton et al. (2018) have shown it is possible to estimate the higher harmonics of the field line automatically using the cross-phase technique. However, determining the plasma mass density distribution is not the focus of this paper. However, the technique has not been widely applied to high-frequency (HF) coherent radars, such as the Super Dual Auroral Radar Network (SuperDARN; Chisham et al., 2007;Greenwald et al., 1995). The electric field component of ULF waves causes an oscillating ExB drift of the ionospheric plasma. This can be observed in SuperDARN line-of-sight Doppler velocity measurements. The line-of-sight velocity can therefore be used as a proxy for the magnetic field oscillations. Fenrich and Waters (2008) used a similar method to diagnose the phase coherence between a FLR observed by the Kodiak SuperDARN radar and Advanced Composition Explorer solar wind data, but as far as the authors are aware, the cross-phase technique has not been applied to HF coherent radar data for the purpose of measuring eigenfrequencies.
The SuperDARN radars offer some advantages over using ground-based instruments. The SuperDARN network provides significantly better spatial coverage and resolution than ground-based magnetometers and can directly observe the F region ionosphere. Therefore, the radars detect these velocity perturbations before they have been either attenuated or rotated by the E region ionosphere (Hughes, 1983;Hughes & Southwood, 1976). Many studies have used SuperDARN to study ULF waves (e.g., Bland & McDonald, 2016;Wright et al., 1998;Yeoman et al., 2000).
However, there are disadvantages to using SuperDARN data. Backscatter can be divided between that which scatters off the ground and irregularities in the ionosphere. Milan et al. (1997) showed that a significant proportion of ground backscatter returns from behind the radars. The uncertainty in the associated ground location of the scatter also increases with range . If the returning signal is actually coming from behind the radar, then the velocity measurements will correspond to a different magnetic latitude to that expected from the front. Therefore, interpretation of the eigenfrequencies from the cross-phase analysis will be incorrect. SuperDARN data can also be unevenly sampled, unlike ground magnetometer data. By uneven, we mean both the sampling rate is inconsistent and the time series contains significant data gaps, which is more likely for SuperDARN. This results in errors in the spectral analysis when using the fast Fourier transform (FFT), which in turn affects the cross-phase. A further issue is how the SuperDARN software derives physical parameters such as line-of-sight velocity. When SuperDARN transmits, it uses a nonredundant sequence of pulses to build up an autocorrelation function (ACF). The standard SuperDARN processing software, FitACF, then uses a fitting routine on the ACF to return useful parameters including backscatter power and line-of-sight velocity. In the FitACF software, the weighting used to fit the ACF lags is based on ad hoc estimates of the ACF lag variance. This can result in incorrect uncertainties in the fitted parameters (Reimer et al., 2018). The lag-filtering criterion can also result in useful data being discarded by the algorithm. Also, because the magnitude and phase of the ACF are fitted directly, aliasing errors can arise in the fitting process. This generates incorrectly fitted data.
In this study, we have attempted to overcome these problems in order to apply the cross-phase technique. We use intervals of data from the Þykkvibaer radar located on Iceland when the European Incoherent Scatter (EISCAT) ionospheric modification (heating) facility (Rietveld et al., 1993) was operating at Tromsø. It was chosen over the Hankasalmi radar due to its eastward viewing direction, required to observe oscillations from the toroidal mode. The heater can artificially generate field-aligned irregularities in the F region ionosphere. These are easily detectable with coherent radars such as SuperDARN. This removes the uncertainty of where the backscattered signal is returning from. The irregularities generated by the heater also have a very low spectral width, which makes them ideal for observing ULF waves (Yeoman et al., 2006).
The uneven sampling of the data can be dealt with using a Lomb-Scargle (LS) spectrogram (Lomb, 1976;Scargle, 1982). LS is a modification of the FFT that can deal with uneven sampling and can offer a greater spectral resolution. Bland et al. (2014), Bland and McDonald (2016), and Shi et al. (2018) have shown LS can be used to analyze ULF waves in SuperDARN data. LS can also yield a phase equivalent to the phase of an FFT spectrogram (Hocke, 1998), which has not been applied to SuperDARN before. We have used these methods to create a LS cross-phase technique to determine eigenfrequencies in SuperDARN data.
To improve the estimation of physical parameters from the ACFs, we have adopted the First Principles Fitting Methodology (FPFM) of Reimer et al. (2018). This uses analytical expressions for the ACF variance to fit the ACFs correctly. Reimer et al. (2018) have shown this method produces more fitted data with more accurate velocity errors.
Finally, we compare our SuperDARN eigenfrequency examples to ground magnetometer estimations of eigenfrequencies and explore the wave populations responsible for the eigenfrequency signatures.

Instrumentation
For this study, we have applied the cross-phase technique to the Þykkvibaer radar from the SuperDARN (Chisham et al., 2007;Greenwald et al., 1995). This is a global network of HF coherent radars designed to monitor the motion of ionospheric irregularities in the E and F region ionosphere. The radars are frequency agile, operating between 8 and 20 MHz and can sample at multiple ranges. In the common mode of operation, the Þykkvibaer radar sweeps over 16 beams, or look directions, and has a range resolution of 45 km. The time resolution in common mode is 60 s. However, the radars can operate with a higher spatial and temporal resolution by reducing the number of beams employed, using smaller range gates and increasing the distance to the first range gate from 180 km to the desired value. An example of such high-resolution viewing is later shown in Figure 1.
To deal with the uncertainty in the return direction of the scattered signal, this study makes use of the EIS-CAT ionospheric modification facility (the heater; Rietveld et al., 1993) located at Ramfjordmoen, Tromsø. The locations of the Þykkvibaer radar and EISCAT heater are illustrated in Figure 1. The colored patch shows the viewing area of Þykkvibaer during a typical heating campaign when the radar was run in a high-resolution mode. The blue dots show the locations of magnetometers used in section 5. The velocity oscillations observed in this mode will be east-west so are associated with the toroidal mode. In this mode, only beams 13-15 were operating. The heater can transmit radio waves with an effective radiated power of 300 MW vertically into the ionosphere (Stubbe, 1996). Here the O-mode polarization can couple to the upper hybrid resonance and feed energy into the ionospheric plasma through a wave-wave interaction (Robinson, 1989). This creates field-aligned irregularities in the F region, which are observable by the SuperDARN radars . We refer to these artificially induced irregularities as "artificial backscatter." This technique involves unmodulated F region heating during quiet magnetospheric conditions. The action of the heater is thus restricted to increasing the backscatter cross section, which the F region ionosphere offers to HF coherent radars (Yeoman et al., 2006). Other experiments have also made use of the generation of artificial backscatter (e.g., Baddeley et al., 2002;Wright & Yeoman, 1999a; for investigating ULF waves. Characteristic properties of artificial backscatter are the strong backscattered power and extremely low spectral widths (e.g., Wright & Yeoman, 1999a;Yeoman et al., 2006), which increases the precision of the velocity measurements.
Finally, ground-based magnetometers from the International Monitor for Auroral and Geomagnetic Effects (IMAGE) array (Luhr, 1994) were used in section 5 for validating results in section 4.

Cross-Phase on Unevenly Spaced Data
The cross-phase technique utilizes two latitudinally spaced ground magnetometers to determine the eigenfrequency of the field line at the midpoint of the station locations (Waters et al., 1991). The eigenfrequencies occur at the frequencies where the cross-phase spectrum is a minimum. It makes two primary assumptions: The length of the field line where the poleward station is located is longer than the field line where the equatorward station is located, and the field lines are driven by broadband power. Such power might be supplied by a spectrum of fast mode waves originating from the magnetopause, for example, which would create a continuum of resonating field lines across the two magnetometers. The phase difference can become positive if the wave velocity gradient changes direction.
The LS periodogram was developed for spectral analysis of unevenly sampled data (Lomb, 1976;Scargle, 1982). Uneven sampling also refers to situations where the data are consistently sampled in time but contain numerous, large data gaps, as is typical of SuperDARN time series. Bland et al. (2014) have shown that the LS periodogram can be used on this kind of unevenly sampled data, and it has been successfully used to study ULF waves in SuperDARN data (Bland et al., 2014;Bland & McDonald, 2016;Shi et al., 2018). However, cross-phase requires complex spectra for its calculation (equation (1)). F Pol ( ) and F Eq ( ) are the complex spectra of the time series data derived from the poleward and equatorward stations, respectively. N 0 is the number of data points in the sample, is the angular frequency under consideration, and the asterisk represents the complex conjugate. Therefore, we need a complex version of LS, in an equivalent form to FFT.
We start with the LS periodogram, defined by Scargle (1982) with equation (2). From this, Scargle (1982) developed the LS periodogram, P( ), given by equation (3). The parameter , the time translation parameter, is defined in equation (4). t j is the set of time values of the data points in the sample, and x(t j ) is its amplitudes. For the SuperDARN time series, x(t j ) would be the line-of-sight velocities.
In the case of even sampling, it can be shown that the LSP reduces to the classical periodogram (equation (2)) that results from the FFT.
Equation (3) produces only real values. We need the LS periodogram in an equivalent form to the FFT with phase and amplitude information. The phase and amplitude of the LS periodogram were derived by Hocke (1998) by considering the LS as a method for least squares fitting harmonic functions to data as was done by Lomb (1976). (Note there is a typographical error in the amplitude in Hocke, 1998). Equation (3) can be written as equation (5), where a and b are defined by equations (6) and (7).
Substituting equation (2) into equation (5) gives the amplitude (equation (8)) in terms of a and b.
Finally, Hocke (1998) showed that the phase can be determined with equation (9). Using equations (8) and (9), we can create a LS equivalent of the FFT, which is given by equation (10). Only a and b need calculating to work out the phase and amplitude.
Equation (10) is the direct LS equivalent of FFT and is identical under even sampling conditions. Spectra calculated with this formula can be used in conjunction with equation (1) to create cross-phase spectra on unevenly sampled data. It should be noted that the time stamps of the data from each beam used to create the two sets of spectra do not have to be identical to each other if a common frequency grid is used to calculate the spectra. Another advantage of using a LS method is we can oversample the frequency axis to improve the frequency resolution. However, this comes at the expense of increased computation time and gives diminishing returns as the oversampling increases. For this study, we have used an oversampling factor of 4, as in Bland et al. (2014).

The First Principles Fitting Methodology
The signal received by the radar can be decomposed into three types: signal from the desired range, noise, and self-clutter (cross-range interference). All of these signals contribute to the variance of the ACF. In the FitACF software, the self-clutter is not accounted for. An ad hoc expression of the variance was used in FitACF when fitting to the magnitude and phase of the ACF, but this expression does not account for the variance introduced by noise and self-clutter. Therefore, empirical filtering conditions were used to remove lags that might contain self-clutter, so-called "bad lags." More detail on this can be found in Reimer et al. (2018). Reimer and Hussey (2015) showed that the self-clutter can be analytically estimated. In a following paper, Reimer et al. (2016) showed that the variance of the ACF can be calculated analytically. This removes the need for the ad hoc estimates of the variance used in FitACF. Reimer et al. (2018) then used these techniques to create a FPFM to create a new fitting routine, LMFIT2 (Levenburg-Marquardt Fitter version 2). LMFIT2 calculates the variance analytically and fits directly to the real and imaginary components of the ACF instead of the amplitude and phase. This removes aliasing errors that arise when fitting to the phase. Reimer et al. (2018) showed that their method produced more fitted data. This was because FitACF was removing the bad lags and so reduced the quantity of fitted data. Using LMFIT2, the velocity errors now represent the fitting error, or the uncertainty, in the velocity. This provides the means to filter data using the velocity error.
We have chosen to use the LMFIT2 fitting routine to improve the data product we apply the cross-phase method to. If the phase of the waves measured in the line-of-sight velocity data have large errors, this will affect the cross-phase results. This could easily arise from incorrect velocity measurements. LMFIT2 provides more reliable velocity values with an associated uncertainty. Figure 2 shows a comparison between data processed with FitACF and LMFIT2. The top row shows backscattered power, or signal-to-noise ratio, of the radar backscatter; the second row line-of-sight velocity; and the third row the velocity error. The left and right columns show the data processed with FitACF and LMFIT2, respectively. The increase in data quantity is clear for the LMFIT2 processing, and the lower velocity errors are clear in the bottom right panel.
Because Reimer et al. (2016) have shown that the velocity errors represent a statistical uncertainty, we have chosen to use the velocity error as a filter to remove poorly fitted data. For all the examples in section 4, the median velocity error for all the data in both beams has been used as a threshold: That is, all data with a velocity error greater than this median have been removed. If the median velocity error is calculated to be less than 10 m/s, then the threshold is set at 10 m/s instead.

Artificial Backscatter Examples
In order to view the artificial backscatter above the heating facility, radio waves propagating from the Þykkvibaer radar must travel 1.5 hops to reach it. Therefore, the radio wave propagation is highly sensitive to the ionospheric conditions, so even when the heater is operating, it is not always possible to observe the artificial backscatter from Þykkvibaer. We selected time periods of heater operation that were long enough to perform the LS cross-phase analysis on. For example, Robinson et al. (1997) ran the heater for 4-min intervals during their experiments, a time period too short to observe a full cycle of Pc5 waves below 4.17 mHz. However, other experiments were run for longer, and it is these examples we have used.
The first example comes from an experiment run on 15 October 1998. This interval of heating has already been analyzed by Yeoman (1999b, 1999c), , and Wright et al. (2004). During this interval, the Þykkvibaer radar was operating beams 11-15, giving a time resolution of ∼11 s in the situation where there are no data gaps. During all of the heater experiments, the Þykkvibaer radar was operating in high time-resolution modes with a reduced number of beams operating.
Range-time-power plots for this interval from beams 14 and 15 of the Þykkvibaer radar are shown in Figures 3a and 3b, respectively. A band of backscatter with high signal-to-noise ratio is present in the range gates centered on range gate 37. This is 1.5-hop scatter. There is also a weaker band at ∼60 range gates. This is 2.5-hop scatter to the same location . This represents observations of the same scatter but with radio waves propagating with a higher elevation angle. To observe 1.5-hop scatter, 1-hop ground scatter must also be present. However, due to the mode of operation, the distance to the first range gate is 1,470 km, beyond the 1-hop ground scatter location. Hence, it is not observed. Figure 3c shows filtered velocity data from range gate 37 on beams 14 and 15. The data have been processed with an 800-s boxcar filter to remove frequencies below the Pc5 range, and spurious velocities outside of −300 < v < 300 are removed. Finally, the mean is removed. Figure 3c shows the presence of waves of many frequencies, ideal for cross-phasing.
Next, the dynamic LS cross-phase spectrum is calculated between range gates 37 on each beam using equation (10). This calculation is only performed if at least 50% of the possible data are available (this is the data quality flag at the top of Figure 3). A 10-min sliding window was incremented every minute to provide good temporal resolution to the spectrogram. The small size of the window and incrementation does make neighboring values highly dependent on one another, but this decision was taken due to the small sample size of 40 min where the scatter could be seen in this example. This calculation was repeated for range gate pairs 35, 36, 38, and 39, but it was found that range gate pair 37 gave the clearest signal in the cross-phase. The spectrogram was then smoothed in both time and frequency with a box 16 min by 3 mHz wide to remove random fluctuations and retain more persistent features. This is shown in Figure 3e. Clear bands of cross-phase minima can be seen in Figure 3e, which may be indicative of eigenfrequencies.
Finally, the B x component of the magnetic field as measured by the Tromsø magnetometer (colocated with the heater) is shown in Figure 3f. This shows very little wave activity on the ground. The average Sym H index was 5.90 nT, so geomagnetic conditions were quiet. The period displayed in Figure 4 was also analyzed by . The data from the two range gates plotted in Figure 4c are in phase and show a wave with a narrow band dominant frequency, which can also be seen in Figure 4f. Unfortunately, due to more than 50% of the data being missing in some time windows, it was unreasonable to calculate the cross-phase spectrum (hence the white spaces in Figures 4d and 4e). In such cases, we should also be wary of edge effects that could arise from the smoothing in Figure 4e. This is an inevitable problem with using SuperDARN data. However, we can see an interesting effect in Figure 4c. The velocities jump repeatedly between approximately 0 m/s −1 and large negative values. This is most likely the received signal being a mixture of artificially induced ionospheric scatter and ground scatter. The propagation path to the artificial backscatter is inconsistent so the radar keeps picking up ground scatter. Such cases make automated determination of cross-phase signatures more complicated.
The period displayed in Figure 5 was analyzed by Wright et al. (2004). It shows a much stronger backscattered signal and a dominant frequency of ∼3.3 mHz, or ∼5 min, from 08:45 to 09:35 UT in Figures 5c and 5f. This is coincident with much weaker cross-phase minima at these times in Figure 5e than the example in Figure 3e than the example in Figure 3. These signals are dominated by an identical narrow band signal at both range gates, and this may reduce the sensitivity of the cross-phase technique. The lack of phase difference between the range gates will not give a minimum in the cross-phase spectrum in Figure 5e. The cross-phase minima become stronger as time progresses, and the signals in Figure 5c become noisier. However, the increase in cross-phase minima may be because of the edge effects previously described.

A Natural Scatter Example
Ideally, the cross-phase technique should be applicable to all natural ionospheric scatter. However, high-resolution beams with a time resolution of less than 10 s are optimal for this, but in the common operational mode, the time resolution is 60 s. This gives a Nyquist limit of 8.3 mHz. In common mode, which is the primary mode of operation of SuperDARN, the higher harmonics will be aliased to lower frequencies as they are generally above this limit (e.g., Takahashi & McPherron, 1982;Wharton et al., 2018). However, the fundamental frequency should still be observable, so the common mode data can be used for this purpose. Without the use of the EISCAT heater, exact determination of the location of the scatter is difficult. Milan et al. (1997) state that ground scatter signals can return from behind the radars, but that ionospheric scatter is always thought to originate from the front. This is because of the angle the magnetic field makes with the ground behind the radar means the propagating rays are unlikely to meet the orthogonality condition to return to the radar. Therefore, for this experiment, we have constrained ourselves to look at ionospheric scatter only.
To find a suitable example, we devised a simple algorithm to count the number of ionospheric data points per hour, assuming the standard ground scatter condition. This algorithm was applied to data from 2008 on the Þykkvibaer radar, and the eastward beams were searched for instances of high numbers of ionospheric data points. It was found best to search closer to the radar as half hop scatter was more common. Due to the convergence of the beams closer to the radar, the chosen beams for cross-phasing were not adjacent to one another. Menk et al. (2004) have shown that the ideal separation is 110 km, and we have tried to maintain this resolution. This however was performed with ground magnetometers. With the greater spatial resolution of SuperDARN, it may be possible to test this determination more thoroughly in a future experiment.
A suitable example of a cross-phase signature was found on 5 March 2008 between range gates 10 and 9 on beams 12 and 15, respectively. The separation of these range gates was 106 km when using a virtual height in the projection model of 400 km. The longitudinal separation was 0.35 • , sufficiently small that we can ignore the presence of phase differences with longitude. This example is shown in Figure 6, in the same format as Figures 3 -5, but without the magnetometer data, as there was no suitable colocated ground station. Figures 6a and 6b illustrate a consistent band of ionospheric scatter in both beams. A variety of wave activity can be seen in Figure 6c. Figures 6d and 6e show this is the case, with a clear minimum in cross-phase occurring at ∼3.5 mHz from 20:00 to 22:00 UT, consistent with the expected fundamental frequency at 66.64 • , the Corrected Geomagnetic latitude of Tromsø (Fenrich et al., 1995;Sandhu et al., 2018). This example shows that it is possible to apply the cross-phase technique to natural ionospheric scatter.

Comparison to Ground Magnetometer Data
Next, we wanted to verify the signatures seen in the examples in section 4 by comparing them to ground magnetometer data. The clearest example was considered to be that in Figure 3. It was compared to the expected eigenfrequencies seen in ground magnetometer data as deduced by the cross-phase algorithm of Wharton et al. (2018). This uses an automated smoothing and filtering algorithm to determine periods of significant cross-phase. Details of how this varies compared to other algorithms are presented in Wharton et al. (2018), but its main advantage is that it can detect the higher harmonics. To do this comparison, B x (geographic north-south) data were used in the algorithm to account for the ionospheric rotation of the magnetic vector (e.g., Hughes & Southwood, 1976). This allowed direct comparison to the eastward-looking beams from Þykkvibaer. Sørøya and Masi were the IMAGE magnetometer stations chosen for this, as the L shell value of their midpoint was approximately the same as that of the heater. The algorithm was run for the month of October 1998. A "slice" was taken through Figure 3e at 12:20 UT and is shown in Figure 7a. This shows the smoothed cross-phase with its associated standard error as a function of frequency (the error bars). The magnetometer eigenfrequency data from the algorithm for a 3-hr-wide Magnetic Local Time (MLT) window centered on 14.74 MLT (the MLT of the 12:20 UT slice) is shown in Figure 7b. This data can be used to estimate the eigenfrequency distribution at the time this slice was taken. Figure 7a shows a series of troughs, the clearest being at ∼4 mHz, followed by three troughs between 8 and 16 mHz and a large trough at 35 mHz. These troughs are the dark red bands in Figure 3e. These troughs are clear despite a standard error of approximately ±20 • in the cross-phase. Figure 7b shows a very clear peak in occurrence at ∼4 mHz, which is assumed to represent the fundamental. This coincides very closely with the lowest-frequency trough in Figure 7a. This suggests that we are observing the fundamental frequency in the SuperDARN data in this example. There are also small peaks in occurrence in Figure 7b, of which some coincide with troughs in Figure 7a. Some of these troughs could be harmonics, but due to the low occurrences for those peaks in Figure 7b, they should not be considered significant measurements.
The large trough at 35 mHz in Figure 7a is particularly clear and may represent a higher harmonic. However, Figure 7b does not confirm this. It may be that the higher harmonics are easier to measure in SuperDARN data than magnetometer data, which could be due to a dependence on scale size. However, it is difficult to validate this claim, and a large-scale study would be needed to see if this peak appears consistently.

The Source of the Cross-Phase Signature
We then examined the variation in wave activity with magnetic latitude. This was undertaken to try and understand what kind of wave activity might be driving the cross-phase signature. B x data from a range of longitudinally aligned magnetometers from the IMAGE array are shown in Figure 8 (names given in the caption). The wave activity in the magnetometers is relatively weak in this example, but greatest in amplitude at the BJN magnetometer. There is also a clear phase difference moving north to south across BJN from LYR down to SOR. This implies a FLR occurred across the stations centered on BJN. FLRs display characteristic amplitude and phase variations with latitude (Samson et al., 1971;Walker et al., 1979). Those are a peak in amplitude at the latitude where the local eigenfrequency matches the FLR frequency, and a 180 • phase change across this latitude, where that phase changes most rapidly at the same latitude. Hence, a latitudinal array of magnetometers can be used to search for these properties. This is what we see in Figure 8. Figure 9 gives us further insight into the wave activity during the example in Figure 3. Figure 9a shows the LS spectrum of BJN (Bear Island) on the dates shown with four times oversampling (e.g., Frescura et al., 2008). The blue vertical line shows the frequency of the FLR as measured at Bear Island. The red line shows the limit of the 800-s second boxcar filter. This filter should have suppressed any frequencies below this limit, although we can see some power still remains. Peaks below 800 s, or 1.25 mHz, should be treated with caution. The green line represents the fundamental eigenfrequency derived from the SuperDARN measurements and the magnetometer algorithm from Figure 7. Figure 9b shows the equivalent spectrum at Tromsø, and Figure 9c shows the spectrum of range gate 37 on beam 15 from the Þykkvibaer radar during the event. Colored vertical lines have the same meaning as in Figure 9a. From these three plots, we can see that the estimated eigenfrequency at Tromsø from the cross-phase method does not contain the most power, neither does it contain the Bear Island FLR. It is worth noting that we expect the amplitude peak to be weaker at Tromsø than at Bear Island due to the lower latitude but, nevertheless, this smaller peak is not obvious.
This is illustrated further in Figures 9d and 9e, which show the amplitude and phase variation with Altitude Adjusted Corrected Geomagnetic latitude at the frequency of the Bear Island FLR. Tromsø and the heater (black horizontal line) lie outside of the peak in power of the Bear Island FLR and its associated phase change. Both the frequency of the Bear Island FLR of ∼2 mHz and the estimated eigenfrequency at Tromsø of ∼4 mHz are consistent with the expected fundamental frequencies at these latitudes according to previous studies (e.g., Mathie et al., 1999;Menk et al., 2004;Sandhu et al., 2018;Waters et al., 1995).
From this we can deduce that the Bear Island FLR was not responsible for the cross-phase signature seen in the magnetometer and radar data because it was at the wrong frequency. Instead, a much weaker source of power was responsible. Waters et al. (1991), Menk et al. (2004), and Clausen et al. (2008) all state that the field lines need only be driven by a weak source with a wide enough band of power in frequency space to drive the field lines at the two magnetometers (or range gates) at their natural frequencies. Such a source of power would create a continuum of resonating field lines across the magnetometers, which is the basis of the broadband assumption mentioned in sections 1 and 3.1. We can see that Figure 9 is consistent with this requirement.

Conclusions
In this study, we have shown it is possible to apply the LS cross-phase technique to SuperDARN ionospheric radar data. This opens additional possibilities of using SuperDARN to study resonant frequencies and magnetospheric plasma densities. Measured eigenfrequencies have been used in conjunction with suitable magnetic field models to estimate plasma mass densities with ground magnetometers in previous work (Berube et al., 2006;Sandhu et al., 2018;Waters et al., 1996), with Wharton et al. (2018) showing that the higher harmonics can be routinely detected and used to determine the shape of the plasma mass density distribution.
In this work, we have focused on resolving some common issues with SuperDARN data. These were the uncertain location of radar backscatter, uneven sampling of the data, and inaccurate fitting of the ACFs. We have resolved these problems with the following solutions: using artificial backscatter from an RF ionospheric heater, developing a LS cross-phase analysis that is suitable for use on unevenly sampled data, and the FPFM of Reimer et al. (2018).
We presented a few examples of eigenfrequencies derived from SuperDARN data, and validated one example by comparing it to ground magnetometer data. We also showed that this signature was caused by a broadband source of power that triggered a continuum of resonating field lines across the station pair/range gates.
However, the expansion of this method is not straightforward. Detecting toroidal eigenfrequencies requires east-west looking beams, however, most SuperDARN radars look northward. It should be possible to expand the cross-phase method to other directions and look at other polarisations, but they may have different driving mechanisms. The method also requires a high time resolution across adjacent beams in order to detect the higher harmonics, but most radars do not currently run these modes. The common mode limits investigation to the fundamental frequency only. Nevertheless, all of these issues can be overcome with subsequent investigation.
In future, we intend to expand this work and create an automated method for determining the eigenfrequencies in a large data set. A significant statistical study will overcome the sporadic nature of scatter seen in the SuperDARN data.