Use of the International Reference Ionosphere 2012 model to calculate emission frequency scale of the ionospheric Alfvén resonator

We apply the IRI-2012 version of the International Reference Ionosphere model to calculate the difference between frequencies of adjacent harmonics (frequency scale) of the ionospheric Alfvén resonator (IAR) emission. The calculated values are compared with the frequency scale data obtained from search-coil magnetometer measurements. To obtain satisfactory results, it appeared necessary to modify the IRI-2012 model by replacing the vertical profile of ionospheric parameters adopted in the standard model with the profile elongated along the magnetic field lines. Subsequent improvement was obtained through the model correction by using local f0F2 measurements. Finally, our results showed strong correlation between estimated and measured values of frequency scale with the root-mean-square error of 0.14 Hz, or 15%.


Introduction
One of the achievements in the last century magnetospheric physics was the theoretical prediction (Polyakov & Rapoport 1981), and then the experimental detection (Belyaev et al. 1987) of resonance structures in the 0.2-10 Hz band. Dynamic spectra of the oscillations look like a series of fan-shaped bands (Fig. 1). Interpretation of emissions proposed by the above authors and adopted in many subsequent studies is based on the concept of resonant Alfvén oscillations in the ionosphere. Hence is the abbreviation IAR (ionospheric Alfvén resonator) often used in the literature. The IAR formation region is located in the cavity formed along the magnetic flux tube between the ionosphere bottom (conductive E region) and the Alfvén speed maximum in the transition region from the ionosphere to the magnetosphere at about 2000-4000 km. In this area, Alfvén oscillations are trapped, and form a standing wave structure (Polyakov & Rapoport 1981).
The further study of the IAR advanced in two directions. On the one hand, a lot of theoretical papers appeared, in which the authors tried to take into account the main features of the ionospheric plasma environment, creating increasingly sophisticated theoretical models. Thus, Trakhtengerts & Feldstein (1984) and Lysak (1991) attracted feedback instability developing in the auroral ionosphere to explain the IAR phenomenon. Initially, the theoretical models of the resonator were built in a simplified geometry with magnetic field lines entering the ionosphere vertically (Lysak 1999), then the models used the geometry with sloping field lines (Sydorenko et al. 2008), and then the dipole field geometry . To formulate the boundary conditions at the lower boundary of the cavity, at first only Pedersen conductivity was involved, and then the Hall conductivity (Lysak & Yoshikawa 2006). Streltsov & Lotko (2008) included coupling to the slow mode in the basic model; Lysak et al. (2013) and Waters et al. (2013) included the coupling to the fast mode. Pokhotelov et al. (2000) found frequency dispersion of the IAR eigenmode due to the Hall conductivity. Alpatov et al. (2005) calculated the characteristics of the Alfvén waves propagating between the Earth's surface and the upper ionosphere using a numerical model. They analyzed two types of external wave sources: pulse and monochromatic oscillations. In their paper, they also gave analytical approximations of numerical results, including a formula for the bounce period.
On the other hand, the morphological studies of the IAR emissions, based both on the ground-based and on satellite observations, have been continued. Daily, seasonal, and cyclic variations in the properties of the IAR (Belyaev et al. 1997;Yahnin et al. 2003) were examined. Their relation to various geomagnetic and geophysical processes (Guglielmi et al. 2006;Potapov et al. 2008;Dovbnya et al. 2010Dovbnya et al. , 2012 and to the state of the ionosphere (Yahnin et al. 2003;Hebden et al. 2005;Parent et al. 2010;Potapov et al. 2014) was studied. A number of authors made suggestions on using the observations of the electromagnetic emission spectral structure within the 0.5-10 Hz range for sounding the ionospheric parameters and specifying ionospheric models (Demekhov et al. 2000;Yahnin et al. 2003;Parent et al. 2010). However, those were only proposals, not supported by numerical estimates.
The purpose of this paper is to quantitatively estimate the closeness of the relationship between the on-ground measured spectral characteristics of IAR emissions and ionospheric parameters within some models of the IAR and of the ionosphere. To do this, as a first stage, we chose the simplest IAR model and the newest IRI-2012 model (Bilitza et al. 2014) of the ionosphere. In some ways, this article is a continuation of Potapov et al. (2014), where we examined in detail the properties of the IAR emissions and for the first time tried to link the changes in ionospheric parameters predicted by IRI-2012 model with the frequency structure of emission patterns. In the present paper, we attempt to explain and remove the discrepancy between observations and model calculations reported in Potapov et al. (2014).

Observations
As observational data of IAR emissions, we used magnetic measurements performed during the 15-month period, March 2010 through May 2011, at the Mondy observatory located on the border between Russia and Mongolia. The observatory geographic coordinates are 51.6°N, 100.9°E; McIlwain parameter is L = 2.1. The station is equipped with a LEMI-30 induction magnetometer. The magnetometer frequency range is 0.001-30 Hz. Within the 1-30 Hz range the device has a flat frequency response with a conversion factor of 20 mV/nT. The noise level does not exceed 0.2 pT/Hz 1/2 at 1 Hz and 0.2 pT/Hz 1/2 at 10 Hz frequency. ADC has 24 digits. Sampling rate is 60 Hz. The timing is provided by the GPS receiver. As the observatory does not possess any ionosounding equipment, we used the data from the Digisonde Ò DPS-4 installed in Irkutsk, about 260 km from the observatory. The Digisonde (Reinisch et al. 1997) is an ionospheric radar that uses high-frequency (HF) radio waves for the remote sensing of the ionosphere. The latitude difference between the two points is less than 0.6°; Irkutsk is to the north. As the main ionospheric data we used 15-min values of critical frequency f 0 F 2 and altitude h m F 2 of the F 2 layer maximum obtained by the standard method.
A special computer code was developed to process digital magnetic records within the 1-10 Hz range. The program calculates the dynamic spectrum of the signal in two components X and Y and builds daily ''frequency-time'' spectrograms. The signal intensity at a given frequency at a given time is displayed in color. The 200 s width Gabor window was used with a 4096-point overlap. The obtained spectrograms were analyzed visually. Most often, 2-4 spectral bands could be digitized. The frequency scale Df values were measured on spectrograms by using a manual on-screen cursor clicking technique for each pair of the adjacent bands opposite every hour mark given that the spectral bands (harmonics) of IAR emission were evident for that hour. (A similar technique was used in Hebden et al. 2005.) Then, the Df values measured between different harmonic pairs were averaged providing one Df value for each hour. Figure 2 shows a typical example of the measured frequencies of spectral bands versus harmonic number for 11 December 2010. Linear regression lines and confidence limit lines (99%) are shown.
For further analysis, we selected 13 daily spectrograms so that they more or less equally represented the change in the IAR emission properties throughout the year. This gave us 99 one-hour experimental values of the frequency scale. The average K p index for the selected days was about 2_ and the three-hour K p did not exceed 4_.

Calculations based on the standard IRI-2012 model
If we are interested in calculating only Df, not estimating harmonic frequencies f j (here j is the number of harmonics), this allows us to avoid analysis of Alfvén wave reflection from the cavity walls and to consider a simple model of standing Alfvén waves reflected ideally from the walls, with plasma and magnetic field in the cavity being inhomogeneous. It is clear that our time-of-flight model is the WKB approximation that is considered usually to be applicable for the wavelengths much shorter than the system scale. But sometimes the WKB approach gives satisfactory results for wavelengths comparable with the system size. In case of the IAR emission, we have long wavelengths for low harmonics, but the space between the harmonics is almost constant for all harmonics, as one can see from Figure 2 (see also Hebden et al. 2005). This is a sign indicating that in this particular case the WKB can be applied, at least as a zero approximation.
In such a case we come to a straightforward expression for Df: here T 0 is the Alfvén travel time between the resonator walls, or the bounce period in terms of Alpatov et al. (2005); l bottom and l top are the positions of the resonator lower and upper walls, respectively; A(l) is the Alfvén speed where the magnetic field B, the electron number density N e , and the effective ion mass m eff vary along the axis of the resonator. The effective mass of the ions at each height level is determined by adding the masses of the ions with a weight corresponding to the relative concentration of ions at a given altitude.
In contrast to Lysak (1991); Demekhov et al. (2000) and Plyasov et al. (2012), we do not develop a theory of the IAR to compare theoretical and experimental spectral characteristics of the IAR emission. Instead, we compare the observational data with the results of Df calculations based on formula (1). For the Alfvén speed computation, we use on-line available models DGRF/IGRF (http://omniweb.gsfc.nasa.gov/ vitmo/igrf_vitmo.html) for the magnetic field, and the International Reference Ionosphere IRI-2012 (http://omniweb.gsfc.nasa. gov/vitmo/iri2012_vitmo.html) for the ionospheric plasma (Bilitza et al. 2014). Because IRI-2012 is limited in altitude by l * = 2000 km, we supplemented this model by extrapolation formulae for number density profiles of different sorts of ions. These formulae are similar to those used in Lysak (2004). Thus, for heights l ! l * , we have the following expressions for ion l dependence and electron number densities: i denotes number densities in cm À3 of oxygen (i = 1), helium (i = 2), and nitrogen (i = 3) singly charged ions at a height of l * = 2000 km according to IRI-2012 model; N 0 H þ is the number density in cm À3 of protons at the same height. For the effective mass density m eff , the first three terms (i = 1, 2, 3) are multiplied by the oxygen, helium, and nitrogen atomic masses, respectively, and the last term by the hydrogen atomic mass. The values h 1 , h 2 , and h 3 are chosen so that Alfvén speed profile A l ð Þ ¼ 21:8B l ð Þ= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N e l ð Þm eff l ð Þ p below l * matches the profile above l * without break; p = 1; R E = 6371 km is the Earth radius; B and m eff are measured in nT and a.m.u. units, respectively; Alfvén speed A(l) in km/s. In total, for 13 selected days of observations, we have performed 99 runs of the frequency scale Df calc tied to the corresponding hours of observations, and compared the results with the measured values Df meas . We took l bottom = 100 km, l top was taken equal to the Alfvén speed maximum altitude, and the step in altitude was Dh = Dl = 50 km. (The step size, as well as the effect of h 1 and h 2 choice, is discussed in Sect. 6.) The comparison result is presented in Figure 3; it shows that the correlation between the calculated and measured values is rather high (the Pearson's correlation coefficient r = 0.90), but the regression line does not run through zero, although it has the~45°slope. The root-mean-square error is about 0.4 Hz. The main cause for such a discrepancy can be as follows: the IRI-2012 model gives variations in the ionosphere variables along a vertical profile, but actually IAR captures Alfvén waves  propagating along the inclined field lines of the geomagnetic field. This leads to the Df overestimation, first, because the resonator length turns to be shorter than its actual length, second, due to the fact that the wave trajectory along a vertical profile runs higher than a field line conjugated with the observatory, i.e. it runs in the ionosphere region, where the Alfvén speed is higher, leading to a further increase in the calculated Df values.

Calculations based on the IRI-2012 model with
a profile elongated along the magnetic field line conjugated with the observatory Our next step was to calculate the frequency scale from the model based on the inclined profile. To do this we used multiple computations of the IRI-2012 results for a chain of ground sites along the Mondy meridian. The procedure is in tracing along the IGRF field line. First, we calculate the required parameters of the ionospheric plasma for 100 km altitude above the Mondy observatory (1st point of profile). Then, we move along the field line that is magnetically conjugated with the 1st point. The 2nd point position is at the crossing of this field line with the 150 km altitude level. Using coordinates of on-ground site under the 2nd point found from IGRF, we find plasma parameters in the 2nd point from the IRI-2012 model. The profile 3rd point is at the crossing of the field line with the 200 km altitude level, and so on, until we reach the 2000 km level. To obtain the plasma parameters at higher levels, we used the extrapolation (2) with the difference that l is now a distance along the magnetic field line, not along the vertical profile. The same is true for l in relation (1). Therefore, step Dl used for numerical integration in (1) is now Dl = Dh/sin I h , where I h is the magnetic field inclination at the altitude h.
In such a way, we transferred from the vertical profile of plasma parameters to the profile aligned along the local geomagnetic field line conjugated with the measurement site. Thus, geometrically, the path is one-dimensional, but it becomes curved. Since the IRI model assumes a smooth change in parameters, these parameters along the new profile also keep continuously changing. Figure 4 shows the results. This time, as we have expected, the regression line running through the origin of coordinates almost coincides with the line of the perfect match Df calc = Df meas . The correlation between the measured and the calculated frequency scales this time is stronger. The rootmean-square error of Df estimate is now about 0.21 Hz; it is approximately one-half of its previous value obtained from IRI-2012 with the vertical profile. However, the scatter of the points on the plot is quite large. Unlike Figure 3, the error is random, not systematic. There are a few possible sources of the error, both in calculating Df calc and in measuring Df meas . One of them is a measurement error due to the blurring of spectral band lines on spectrograms, as well as owing to the band distortions caused by the ionospheric profile inhomogeneity (Ermakova et al. 2011), or by an interaction between the IAR emission and field-line resonances (Bösinger et al. 2004). The calculation error may be caused by the IRI-2012 model non-perfection: for certain geophysical conditions, it may give inaccurate values of the ionospheric parameters above Mondy.

Correction of results by using local DPS-4
measurements of the f 0 F2 critical frequency The latter possibility can be tested. To do this, we compared the Df calc deviation from Df meas with the difference between the F 2 layer critical frequency f 0 F 2 estimated by the IRI-2012 model and that measured by the DPS-4 digital ionosonde for the whole series of 99 observational hours. Figure 5a shows the result. We see that the relationship is linear and has a relatively high correlation, 0.72. The regression equation is Df calc À Df meas = 0.215 ( f 0 F 2meas À f 0 F 2calc ). This regression means that, for example, when IRI-2012 underestimates the critical frequency by 1 MHz, our model with the inclined profile (based on IRI-2012) overestimates the frequency scale by 0.215 Hz on average. Hence, to be corrected, the computed frequency scale Df calc shown in Figure 4 should be replaced with Df calc * = Df calc À 0.215 ( f 0 F 2meas À f 0 F 2calc ). By plotting Df calc * against Df meas , we obtain our final result shown in Figure 5b. The root-mean-square error of Df estimate is now 0.14 Hz, reduced by half as compared to the non-corrected value. Taking into account the Df meas % 0.95 Hz mean value, we obtain the relative root-mean-square error of 15%. Again, the remaining scatter of points around the line of perfect match may have various reasons. If we exclude the measurement error, the main reasons are the two mentioned as follows. First, this is a mismatch of sites of magnetic measurements and of the ionosphere radio sounding. Two hundred sixty km between them may be a substantial distance, especially if there are horizontal ionospheric irregularities. Second, this is an additional non-perfection of the model used. The point is that, by using the local f 0 F 2 measurements we corrected only the profile's lower part. But its part above the F 2 layer maximum may be represented by both the IRI-2012 model and the extrapolation we have adopted only with a significant error.

Summary and discussion
So, we have shown that the modified IRI-2012 version corrected by local f 0 F 2 measurements permits to estimate the IAR frequency scale with a very high correlation between the estimated and measured values and with the~15% relative mean-square error. We used a method that can be termed a sequential elimination technique. We started with using the standard IRI-2012 model. The result revealed this model's deficiency when using it for our purposes due to the vertical profile effect implemented in the model. Therefore, we excluded the effect by replacing the vertical profile with an inclined one elongated along the magnetic field lines. This required the model modification. The new result was much more satisfactory. The residual standard error looked purely random until we compared it with the error of the critical frequency estimation by the IRI-2012 model. We have adjusted our model, taking into account the results of local measurements of the critical frequency by on-ground instrument. This has improved the result significantly. In principle, we could continue to adjust our model, excluding other extraneous factors. For example, some seasonal variation can be traced in the residual error of the model's last version. We could exclude it, thus reducing the mean-square error. However, this will likely not lead to the model improvement, because the seasonal variation in Df calc * is weak, and there is no certainty that, with the next sample of the IAR emission measurements, the seasonal variation in the residual error will not have a very different kind. Now, we discuss some details of our model. First, it should be noted that the step size, Dh = 50 km, was selected as a compromise between the calculation accuracy and the speed of computation. Figure 6a shows two altitudinal profiles of the Alfvén speed calculated with a different step size: Dh = 50 km (thick gray line) and Dh = 10 km (dash line). We see that the significant differences between the profiles take place only at the lowest levels of the ionosphere, from 100 to 150 km. At the same time, it should be noted that the main part of time-of-flight is accumulated in the middle of the wave path, where the Alfvén speed is close to its minimum. That is why the Dh decrease down to 10 km gives less than 1% change in Df calc . For the 7 March 2011, 2100 UT event, we have Df calc = 1.30561 Hz with Dh = 10 km, and Df calc = 1.29357 Hz with Dh = 50 km.
The graphs in Figures 6b and 6c show other possible sources of variations in the Df calc values, apart from those given by the model IRI-2012. These sources are the changes in scales h 1 and h 2 that determine slope of oxygen (h 1 ) and helium (h 2 ) ion number density above the l * = 2000 km level (see Sect. 3). (Note that h 3 = 500 km in all our calculations.) Remember that we chose the h 1 and h 2 parameters so that the curve of the Alfvén speed profile had no break at 2000 km. We see in Figure 6b that the variations in h 1 have very slight effect on the speed profile. Naturally, they almost do not affect Df calc as well. Therefore, even if we choose a not so entirely correct value of h 1 , it does not introduce large errors in Df calc . On the contrary, variations in h 2 change the profile noticeably as one can see from Figure 6c. In this case, we see clearly that the ''correct'' choice is h 2 = 400 km. Nevertheless, the difference between Df calc computed with three versions of h 2 is not large. This may be because the two tendencies, the change in the maximum speed and variation in the maximum value height, can compensate each other.
On the one hand our results confirm again the existing theory of IAR; on the other they open some ways to search for using magnetic observations of IAR emission to obtain information on the ionosphere, including its upper regions. These data are hidden in deviation of the model estimations from the measured values. However, to extract the information an inverse problem must be solved.
The high correlation between the IRI forecast and the measurement of IAR frequency structure obtained with the simplest model, evidences, in our view, that, in this case, more complex models will not give significant improvement. Of course, using a simple IAR model can be justified only to assess the frequency scale Df of the emission; this method cannot be applied to estimate the eigenmode frequency or the amplitude of spectral bands. To do this, more sophisticated methods are needed. Indeed, the eigenmode frequency depends Fig. 5. Correction of results by taking into account the difference between the predicted and measured f 0 F 2 : (a) Df calc deviation from Df meas versus the difference between F 2 layer critical frequency f 0 F 2 estimated by the IRI-2012 model and that measured by DPS-4 Digisonde Ò ; dashed line is the regression line running through the origin of coordinates; thin lines are 99% confidence limit lines; (b) the same as in Figure 3, but with correction taking into account the Df calc À Df meas difference dependence on the IRI-2012 model's inaccuracy in f 0 F 2 estimates. not only on the Alfvén speed inside the resonance cavity and on the cavity size, but also on the processes of wave reflection from the lower and upper boundaries of the resonator as well. These processes include the influence of the Hall and Pedersen conductivities and coupling of the Alfvén mode to the fast and slow modes, see Lysak & Yoshikawa (2006) for a review.
The above example of applying the IRI model to IAR emission suggests that this model has the potential to a much wider application than purely ionospheric problems. We have shown that, after a certain model modification, IRI can be successfully used when analyzing the wave fields, primarily of the Alfvénic type. These fields include, besides IAR emissions, pearls, other pulsations of the Pc1 frequency band, and the fluctuations associated with the field-line resonances. We encourage IRI developers to develop a mode version with the profile inclined along the magnetic field lines, so that researchers can apply it to study ULF waves in the magnetosphere-ionosphere system.