© A.S. Potapov et al., Published by EDP Sciences 2015

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

Fig. 1. Example of the IAR emission discrete spectral bands. Shown are daily spectra of the magnetic field X- and Y-components from the measurements with LEMI-30 induction magnetometer.

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 (Lysak et al. 2013). 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. 2010, 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).

2. 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₀F₂ and altitude h_mF₂ of the F₂ 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 Δf 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 Δf values measured between different harmonic pairs were averaged providing one Δf 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.

Fig. 2. Scatter plot of harmonic frequency (asterisks) plotted against harmonic number for each hour 08:00 through 12:00 UT, with the least squares regression lines through these points (sky-blue lines) and the confidence limit (99%) shown by thin lines.

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

3. Calculations based on the standard IRI-2012 model

If we are interested in calculating only Δf, 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 Δf:(1)here T₀ 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 speedwhere 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 Δf 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: ; N_e(l) = N_i(l), where denotes number densities in cm⁻³ 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; is the number density in cm⁻³ 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₁, h₂, and h₃ are chosen so that Alfvén speed profile 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 Δf_calc tied to the corresponding hours of observations, and compared the results with the measured values ∆f_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 Δh = Δl = 50 km. (The step size, as well as the effect of h₁ and h₂ 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 ∆f 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 Δf values.

Fig. 3. Calculated ∆fcalc versus measured ∆fmeas values of the IAR frequency scale. Calculations were performed by using the standard IRI-2012 model with a vertical profile. Bold brown line is the linear regression; thin lines are the 99% confidence limits; thick gray line is the line of perfect match Δfcalc = Δfmeas.

4. 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 Δl used for numerical integration in (1) is now Δl = Δh/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 Δf_calc = Δf_meas. The correlation between the measured and the calculated frequency scales this time is stronger. The root-mean-square error of ∆f 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 Δf_calc and in measuring Δf_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.

Fig. 4. Calculated ∆f calc versus measured ∆f meas values of the IAR frequency scale as a result of calculations based on the modified IRI-2012 model with the inclined profile elongated along the local geomagnetic field line. Bold brown line is the regression line running through the origin of coordinates; thin lines are the 99% confidence limits; gray line is the line of perfect match Δf calc = Δf meas.

5. Correction of results by using local DPS-4 measurements of the f₀F2 critical frequency

The latter possibility can be tested. To do this, we compared the ∆f_calc deviation from ∆f_meas with the difference between the F₂ layer critical frequency f₀F₂ 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 ∆f_calc − ∆f_meas = 0.215 (f₀F_2meas − f₀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 ∆f_calc shown in Figure 4 should be replaced with ∆f_calc^* = ∆f_calc − 0.215 (f₀F_2meas − f₀F_2calc). By plotting ∆f_calc^* against ∆f_meas, we obtain our final result shown in Figure 5b. The root-mean-square error of ∆f estimate is now 0.14 Hz, reduced by half as compared to the non-corrected value. Taking into account the ∆f_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₀F₂ measurements we corrected only the profile’s lower part. But its part above the F₂ layer maximum may be represented by both the IRI-2012 model and the extrapolation we have adopted only with a significant error.

Fig. 5. Correction of results by taking into account the difference between the predicted and measured f0F2: (a) ∆fcalc deviation from ∆fmeas versus the difference between F2 layer critical frequency f0F2 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 Δfcalc − Δfmeas difference dependence on the IRI-2012 model’s inaccuracy in f0F2 estimates.

6. Summary and discussion

So, we have shown that the modified IRI-2012 version corrected by local f₀F₂ 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 ∆f_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, ∆h = 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: ∆h = 50 km (thick gray line) and ∆h = 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 ∆h decrease down to 10 km gives less than 1% change in ∆f_calc. For the 7 March 2011, 2100 UT event, we have ∆f_calc = 1.30561 Hz with ∆h = 10 km, and ∆f_calc = 1.29357 Hz with ∆h = 50 km.

Fig. 6. Variations in vertical profiles of the Alfvén speed based on the results of calculations depending on (a) the step ∆h in height, (b) the characteristic scale h1 of the density slope of oxygen ions, and (c) the scale h2 of the density recession of helium ions.

The graphs in Figures 6b and 6c show other possible sources of variations in the ∆f_calc values, apart from those given by the model IRI-2012. These sources are the changes in scales h₁ and h₂ that determine slope of oxygen (h₁) and helium (h₂) ion number density above the l_* = 2000 km level (see Sect. 3). (Note that h₃ = 500 km in all our calculations.) Remember that we chose the h₁ and h₂ 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₁ have very slight effect on the speed profile. Naturally, they almost do not affect ∆f_calc as well. Therefore, even if we choose a not so entirely correct value of h₁, it does not introduce large errors in ∆f_calc. On the contrary, variations in h₂ change the profile noticeably as one can see from Figure 6c. In this case, we see clearly that the “correct” choice is h₂ = 400 km. Nevertheless, the difference between ∆f_calc computed with three versions of h₂ 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 ∆f 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 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.

Acknowledgments

The data supporting Figures 2–5 are available in the Supporting Information Data Table. Authors are grateful to the team of the IRI model. We also thank the team of the ISTP Complex Magneto-Ionospheric Observatory (Irkutsk) for providing magnetic and ionospheric data. The part of work made by A.S.P. (development of the IAR spectrum formation model) was supported by the Russian Science Foundation, Grant 14-37-00027; the part of work made by T.N.P. and B.V.D. (observations and morphological studies) was supported by the Russian Foundation for Basic Research, Grant 13-05-00529. The editor thanks C. Waters and an anonymous referee for their assistance in evaluating this paper.

References

All Figures

	Fig. 1. Example of the IAR emission discrete spectral bands. Shown are daily spectra of the magnetic field X- and Y-components from the measurements with LEMI-30 induction magnetometer.
In the text

	Fig. 2. Scatter plot of harmonic frequency (asterisks) plotted against harmonic number for each hour 08:00 through 12:00 UT, with the least squares regression lines through these points (sky-blue lines) and the confidence limit (99%) shown by thin lines.
In the text

	Fig. 3. Calculated ∆fcalc versus measured ∆fmeas values of the IAR frequency scale. Calculations were performed by using the standard IRI-2012 model with a vertical profile. Bold brown line is the linear regression; thin lines are the 99% confidence limits; thick gray line is the line of perfect match Δfcalc = Δfmeas.
In the text

Fig. 4. Calculated ∆f calc versus measured ∆f meas values of the IAR frequency scale as a result of calculations based on the modified IRI-2012 model with the inclined profile elongated along the local geomagnetic field line. Bold brown line is the regression line running through the origin of coordinates; thin lines are the 99% confidence limits; gray line is the line of perfect match Δf calc = Δf meas.

In the text

Fig. 5. Correction of results by taking into account the difference between the predicted and measured f0F2: (a) ∆fcalc deviation from ∆fmeas versus the difference between F2 layer critical frequency f0F2 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 Δfcalc − Δfmeas difference dependence on the IRI-2012 model’s inaccuracy in f0F2 estimates.

In the text

	Fig. 6. Variations in vertical profiles of the Alfvén speed based on the results of calculations depending on (a) the step ∆h in height, (b) the characteristic scale h1 of the density slope of oxygen ions, and (c) the scale h2 of the density recession of helium ions.
In the text