Observing the cold plasma in the Earth ' s magnetosphere with the EMMA network

We illustrate a semi-automated procedure to detect the field line resonance (FLR) frequencies and the derived equatorial plasma mass densities in the inner magnetosphere from ULF measurements recorded at the European quasi-Meridional Magnetometer Array (EMMA). FLR frequencies are detected using the standard technique based on cross-phase and amplitude ratio spectra from pairs of stations latitudinally separated. Equatorial plasma mass densities are then inferred by solving the toroidal MHD wave equation using the TS05 Tsyganenko magnetic field model and assuming a 1/r dependence of the mass density along the field line. We also present a statistical analysis of the results obtained from 165 non-consecutive days of observations at 8 station pairs covering the range of magnetic L-shells 2.4-5.5 and encompassing a wide range of geomagnetic conditions. The rate of FLR detection maximizes around local noon at each pair of stations, reaching the highest values (~95%) around L = 3. A clear diurnal modulation of the FLR frequency is observed at all L values. At the lowest latitudes, the variation is characterized by a rapid decrease in the early morning hours, a stagnation in the middle of the day, and an increase in the evening hours. At higher latitudes, a continuous and more pronounced decrease of the FLR frequency is observed during all daytime hours reflecting a permanent state of recovery of flux tubes depleted by events of enhanced magnetospheric convection. Consistently, the radial profiles of the inferred equatorial mass density show a density increase from morning to afternoon which gets more pronounced with increasing distance and with the level of the preceding geomagnetic activity. The results also confirm the formation of the plasmapause at geocentric distances that decrease as the disturbance level increases. Mean mass density distributions in the equatorial plane are also shown in 2-D maps for different geomagnetic conditions, as well as for a representative stormy day.


INTRODUCTION
Field line resonances (FLRs) naturally occur through− out the magnetosphere and are a useful tool to sound the dynamics of the cold plasma in plasmasphere and plas− matrough [Waters et al., 2006;Menk et al., 2014]. In fact, the resonance frequencies of a geomagnetic field line depend on the field line length, and on the magnetic field intensity and plasma mass density ρ along the field line. Therefore, assuming a model for the magnetos− pheric field and a functional form for ρ along the field line, the mass density at a given point of the field line (usually at the equatorial point) can be inferred from the observed FLR frequency.
In the past, several techniques were developed to de− termine FLR frequencies from ground−based ULF meas− urements. They mostly used latitudinal variation of ULF waves properties such as spectral power and phase or polarization, which show characteristic features at the resonance frequency [Sugiura and Wilson, 1964;Sam− son, 1973;Miletits et al., 1990;Ziesolleck et al., 1993;Fenrich et al., 1995]. Estimates of the FLR frequency were also obtained by comparing the spectral power of the northward (H) component of magnetometer data with the eastward (D) one . Menk and Waters [2013] made an historical review of these techniques.
In this work we refer to the well established gradi− ent method introduced by Baransky et al. [1985] and further developed by Waters et al. [1991]. The gradient method works using data collected by a pair of magne− tometers slightly separated in latitude and approxi− mately aligned along the same magnetic meridian. The method takes advantage of the fact that the field lines connected to the two stations have similar frequency response, in both amplitude and phase, but slightly shifted in frequency. In particular, the resonance fre− quency generally decreases with increasing latitude, ex− cept near the plasmapause [Menk et al., 2004] or at very low latitudes ( λ < ~ 35°) [Menk et al., 2000;Kawano et al., 2002]. The latitudinal separation between magne− tometers must be enough to distinguish the resonance frequency of the two field lines but at the same time ensure sufficient coherency between the two signals. Typically, a separation of 1−3 degrees is required, de− pending on the latitude and the quality of the signals. Baransky et al., [1985] proposed the idea to evalu− ate the resonance frequency of the field line whose footprint is halfway between the stations, by analyzing the H−component amplitude spectra. They showed that, at this frequency (f r ), the amplitude ratio A E /A P between the equatorward (E) and poleward (P) station should cross unity with positive slope or, equivalently, the am− plitude difference should cross zero. Waters et al., [1991] pointed out that a more reliable determination of f r is provided by the so−called "cross−phase technique" which identifies f r as the frequency where the phase dif− ference φ E − φ P maximizes.
The FLR frequencies f r detected as described above are assumed to be eigenfrequencies of the toroidal os− cillations of the given geomagnetic field line. Theoret− ically, such frequencies can be computed for an arbitrary magnetic field geometry by solving the fol− lowing wave equation developed by Singer et al., [1981]: (1) where s is the curvilinear coordinate along the field line, ξ α is the field line displacement along a normal direc− tion α (azimuthal direction for the toroidal mode), h α is a scale factor proportional to the distance to an adja− cent field line in the α direction, B is the magnetic field, V A is the Alfvén velocity = B/(µ o ρ) 1/2 , and ω = 2πf r is the angular eigenfrequency. ξ α = 0 at the ionospheric footprints of the field line are used as boundary condi− tion. In order to solve equation (1), both magnetic field B(s) and mass density ρ(s) models are required.
It is common use to assume the following power law for the mass density distribution: (2) where r is the geocentric distance, r eq is the geocentric distance of the equatorial point of the field line, ρ eq is the mass density at the equator, and m is the power law index.
The proper value of m has been investigated by many authors [e.g., Menk et al.,1999;Takahashi et al., 2004;Vellante and Förster, 2006;Reinisch et al., 2009]. In par− ticular, Takahashi et al., [2004], by performing a multi− harmonic analysis on CRRES satellite observations, found an optimal value of m ∼ 0.5 in the range of mag− netic shells 4 < L < 6. Vellante and Förster, [2006], using a plasmaspheric physical−numerical model, concluded that a good choice for midlatitude field lines (2.3 < L < 3.4) is m = 1 for a large variety of solar and geomagnetic conditions. They also pointed out that significantly higher m values would be necessary for L < 2, due to the increased contribution of ionospheric mass loading. For the present study we assumed m = 1 which, accord− ing to these previous results, is a reasonable choice for the investigated latitudinal range (2 < L < 6). In any case, in this L−range, different choices of m in the range 0−6 would produce variations in the inferred equatorial den− sity smaller than the typical indetermination due to the FLR frequency uncertainty Berube et al., 2006;Vellante and Förster, 2006].
As regards the geomagnetic field model, a dipole ap− proximation has been usually adopted for deriving the equatorial density from FLRs detected at low and mid latitudes (L < 4) [e.g., Menk et al., 1999;Berube et al., 2005;Villante et al., 2006;Chi et al., 2013;Wang et al., 2013]. Vellante et al., [2014a], however, found that the use of the dipole model at these latitudes may result in an error in the inferred density appreciably larger than what is usually assumed, up to 40% at L = 2 when com− pared with the results obtained using the International Geomagnetic Reference Field (IGRF) model.
At higher latitudes, the currents associated with the interaction of the geomagnetic field with the solar wind may give a significant contribution and even the IGRF model may become inappropriate. Warner and Orr [1979] were the first to point out that a more realistic model is needed to describe the geomagnetic field at high latitudes. Nowadays, the Tsyganenko's models [see Tsyganenko, 2013 for a review] are widely used in plasma mass density evaluation [Waters et al., 1996;Berube et al., 2006;Takahashi et al., 2006;Kale et al., 2009;Maeda et al., 2009]. Huang et al., [2008] evaluated the performance of some versions of the Tsyganenko models during magnetic storm times. They compared the output of the T96 [Tsyganenko, 1996], the T02 [Tsyga− nenko, 2002a,b] and the TS05 [Tsyganenko and Sitnov, 2005] models with the magnetic field observations of the GOES satellite at geosynchronous orbit. They con− cluded that TS05 predicts the geostationary field best for all conditions, while the other models have larger errors for the main and the recovery phases of a storm. On the basis of these results, we adopted the TS05 model for the present study. By applying the gradient method to pairs of stations of a latitudinally extended network, it is possible to mon− itor the radial dependence of the equatorial plasma mass density along the longitudinal sector identified by the array. Such a possibility is offered by EMMA, the Euro− pean quasi−Meridional Magnetometer Array, which was established in 2012 in the framework of the PLASMON FP7 European project [Lichtenberger et al., 2013]. Section 2 introduces the EMMA network and describes in detail the automated or partially automated algorithms imple− mented to infer the equatorial plasma mass density. Sec− tion 3 shows the results of a statistical analysis conducted applying the procedure described in section 2 to selected periods encompassing a wide range of geomagnetic con− ditions. Section 4 presents conclusions.

THE EMMA ARRAY AND THE MASS DENSITY INFERENCE IMPLEMENTATION
EMMA currently consists of 27 stations and was set up by unification and extension of previously existing European magnetic arrays: SEGMA (South European GeoMagnetic Array), MM100 (Magnetic Meridian 100) and Finnish part of IMAGE (International Monitor of Auroral Geomagnetic Effects). Figure 1 shows the location of the stations in geo− graphic coordinates. As can be seen, the array is suit− able for determining FLR frequencies by means of the gradient technique for an extended range of L values (1.6−6.2). Therefore, EMMA is in general able to mon− itor the plasma mass density both in the plasmasphere and in the plasmatrough.
This monitoring system is operative on a near real time base. 1 sec data are indeed collected at two main servers (one in Hungary and the other in Italy) which continu− ously monitor the operational status of the whole net− work, and cyclically (every 15 min) run a procedure (FLRID) to detect the fundamental FLR frequencies from several station pairs (http://geofizika.canet.hu/plasmon/). A separate automated process (FLRINV) converts the de− tected FLR frequencies into estimates of equatorial plasma mass densities. This is done by numerically solving the wave equation (1) using the T02 Tsyganenko magnetic field model. Solar wind parameters which are necessary for the T02 model are taken online from the NOAA Space Weather Prediction Center which provides the latest 2 hours of magnetic and plasma data of the ACE satellite located at the L1 libration point (http://services.swpc.noaa.gov/text), while real−time Dst data are taken from the Dcx server of the University of Oulu, Finland (http://dcx.oulu.fi/DstDcxDxtData/Real− Time/Dxt/DxtRT.txt).
After having carefully analyzed the first results de− rived from the original code, we realized that, under particular conditions, the automated procedure fails in detecting the correct FLR frequencies. At high latitudes, sometimes, the selected frequency is not the fundamen− tal one but a higher harmonic, which leads to underes− timate the inferred plasma density. The fully automated algorithm also fails sometimes when applied to station pairs near the plasmapause.
A modified and supervised version of FLRID has been developed and is presented in section 2.1. This version is a transitional code meant to provide checked and trusted FLR frequencies and at the same time to gather infor− mation useful to improve the original code. The final goal is to reach a reliable fully automated version.

FLRID -IDENTIFICATION OF FLR FREQUENCY
FLR frequencies are detected by performing a cross−spectral analysis of ULF magnetic signals, and the capability to observe FLR signatures is directly linked to the spectral method chosen to derive the spectra and to the parameters used in applying that method. The gradient technique is generally imple− mented through a Fast Fourier Transform (FFT) [e.g., Berube et al., 2003] and so does FLRID, although some authors tried alternative methods [e.g., Boudouridis and Zesta, 2007].
We developed an interactive version of FLRID with a graphical user interface (GUI) to support the visual in− spection of the cross−spectra and facilitate the change of the default parameters. It consists of a control panel called FLR FINDER, and two interactive windows to se− lect/deselect frequencies directly from graphs. Figure 2 shows a screenshot of FLR FINDER while FLRID is pro− cessing data from TAR and BRZ stations on 6 June 2013.
At each run magnetic data from two stations at 1 Hz resolution, organized in daily files, are analyzed. Then, the FLR frequency detection proceeds as follows.
1) A synchronization of the data from the two sta− tions is automatically performed, if a time delay is found from cross−covariance analysis of the D−component signals that are expected to be in− phase at such small station separation. A dedi− cated control in FLR FINDER allows one to man− ually change the time delay.
2) The data are differentiated, and the cross−spec− trum and the power spectra are evaluated in a moving 2−hour time window which is shifted by half an hour in successive steps. This length of the time window was found to be the minimum length able to confidently detect the low reso− nance frequencies at high latitudes. Shorter time windows could be used at lower latitudes (plas− masphere region), but for the present study we decided to choose a fixed time window for all pairs in order to have the same temporal resolu− tion. In any case, both the time window length and time step can be modified in the specific section of FLR FINDER called "Spectral analy− sis". 3) For each 2−hour time window the spectra are evaluated for a variable number of overlapping sub−intervals that are successively averaged. The number of sub−intervals depends on the latitude of the stations considered and for the present study we varied it from 1 (highest latitudes) to 7 (lowest latitudes). The final spectrum is smoothed over a number of samples that de− creases with latitude from 11 to 9. All these pa− rameters are set in the section "Spectral analysis" of FLR FINDER. Figure 3 shows the screenshot of a 2−hour window analysis that is shown at each half hour step. Spectral density of the two stations, coherence, amplitude ratio and cross− phase are shown simultaneously. 4) The cross−phase and the amplitude ratio are fit− ted with a smoothed spline to avoid rapid oscil−  lations (magenta curves in Figure 3). FLR FINDER has also a section called "FIT" to set the fit method and related parameters. Then FLRID searches for maxima in the smoothed cross− phase curve. To identify a maximum as an FLR frequency, we required the peak to exceed the baseline by at least 10° and the coherence at the selected frequency to be higher than 0.1. The uncertainty associated to the frequency is the half width of the peak at 80% of its height. This uncertainty is represented as a horizontal bar in the cross−phase panel of Figure 3. 5) The expected behavior of the amplitude ratio q is used as an additional criterion for the selection. The cross−phase identified frequency f r must lie between a minimum q − and a maximum q + of q, and q + must be greater than q − by a given threshold which we set to 5%. The analysis of q allows one to obtain also an alternative estimate of the resonance frequency (f r ′), as the frequency where q = M = (q − q + ) 1/2 [Green et al., 1993]. The deviation between the two estimates | f r ′− f r | provides an additional measure of the uncer− tainty in the f r estimate. The finally selected fre− quencies are indicated by red dashed lines. 6) Sometimes a selected frequency appears clearly unreasonable (with respect to the typical values expected for that station pair, or with respect to the previous/next values). This may be due, for example, to particular noise conditions affect− ing one or both stations which generate spectral features matching, just by chance, the selection conditions. These frequencies can be manually discarded in an interactive way. 7) After the whole day is analyzed, the program generates a second window containing the dy− namic amplitude ratio and the dynamic cross− phase spectra. Figure 4 shows such an output for the run of Figure 2. The scale of the amplitude ratio colorbar is log 10 (q/M * ), where M * is the me− dian of M. Selected frequencies are overplotted onto the dynamic color spectra. In this particular DEL CORPO ET AL.
6 case three harmonics are clearly visible for most of the day. From the bottom panel it is possible to select in an interactive way the data points to be used for the inversion (or for other scopes) by drawing a curve around them. In the typical sit− uation we select points which confidently corre− spond to the first harmonic. A text file containing the selected frequencies (and related information) is then automatically generated. Typically, an FLR frequency is identified by a posi− tive maximum in the cross−phase, but across the plasmapause, the steep density variation may cause a reversal in the slope of the Alfvén velocity (or, equiv− alently, of the FLR frequency) radial profile, leading to a cross−phase reversal [Menk et al., 2004;Kale et al., 2007]. In this case the FLR frequency can be identified by a cross−phase minimum, as well as by a reversal in the amplitude ratio behavior. These cases are consid− ered by re−running the same algorithm but by revers− ing the order of the two stations in the "General settings" window ( Figure 2).
The output shown in Figure 4 may also be gener− ated directly (skipping the examination of the single spectral results showed in Figure 3), by deselecting the "deep analysis" button in the control panel in Figure 2, resulting in significant time saving. This is actually the most common way we use FLR FINDER. Only critical cases require a deeper analysis. For example, when the field line of the station pair maps near the plasmapause, or during highly dynamic geomagnetic conditions when the diurnal variation of the FLR frequency may show unusual patterns.
We point out that all the parameters we chose for the selection derive from a series of trials at each sta− tion pair and are by no means mandatory values, we just found them to give a good compromise between quantity and quality of the selected frequencies.

FLRINV -FLR FREQUENCY INVERSION
For each selected FLR frequency f r , an automated process (FLRINV) infers the corresponding equatorial mass density ρ eq . It is basically the same code devel− oped during the PLASMON project with minor updates.
We use the following reformulation of the Singer equation (1) in terms of normalized quantities [Vellante et al., 2014a]: . The power law model for the mass den− sity along the field line (equation 2) with m = 1 has been used in the expression of the Alfvén velocity. The boundary conditions become ξ ′ α (0) = ξ ′ α (1) = 0. The magnetic field strength B and the coordinates along the field line s (and hence l and x) are retrieved from the TS05 magnetic field model and by tracing the magnetic field line whose footprint is halfway between the con− sidered stations. The scale factor h α is obtained by com− puting the distance along s from an adjacent field line separated at the apex by an angle Δφ = 0.01°. The eigenvalues λ, solutions of equation (3), are found using the Dormand-Prince method [Dormand and Prince, 1980] implemented as a built−in function on the MAT− LAB environment, and a shooting code. If the detected frequency f r is the fundamental frequency it is suffi− cient to determine only the first eigenvalue. The corre− sponding equatorial mass density is then derived by the following expression [Vellante et al., 2014a]: (4)

FITTING PROCEDURE OF THE RADIAL DENSITY PROFILE
The overall picture of the distribution of the cold plasma in the inner magnetosphere derived by EMMA measurements can be better appreciated by investigat− ing the radial profiles of ρ eq . The equatorial densities derived from a given sta− tion pair correspond to a time−changing radial distance. This effect is more important for high latitude station pairs, and in general during disturbed magnetospheric conditions. A fitting procedure makes it possible to de− termine ρ eq at fixed radial distances.
The criteria adopted in fitting the data can be sum− marized as follows: 1) data points with r eq > r lim are not considered, where r lim = 15 Earth radii (R ⊕ ) is the limit of va− lidity of the TS05 model; 2) only profiles with at least 5 data points are fitted; 3) at each time, plasma density data points are in− terpolated by a smoothing spline curve; 4) radial distance intervals larger than 3 R ⊕ without experimental data points are excluded in the fit− ted profiles. Point 1 ensures that only reliable data points are considered in the fitting procedure. Points 2 and 4 avoid questionable fitted values in case of poorly populated profiles.
The smoothing spline is a tradeoff between a pure 7 PLASMASPHERE MONITORING BY EMMA NETWORK ρ eq = λB eq spline and a linear least square fitting and was found to be the best way to describe the typical radial plasma den− sity profile. The computation is conducted using the Curve Fitting Toolbox built−in functions based on the MATLAB environment. An example is shown in Figure 5. Experi− mental data points, together with their uncertainty, are represented with open blue circles, while the fitted pro− file is shown with the red solid curve. The dashed lines represent the uncertainty in the fitted profile and were obtained by applying the procedure described above to the upper and lower values of the uncertainty interval. As can be seen, the procedure provides a way to de− scribe the radial distribution of the plasma mass density with a continuous smooth curve and is able at the same time to well characterize transition regions like the plasmapause.
The fitted profiles obtained at different local times can be merged in a single image in polar coordinates to obtain a two−dimensional density map in the magnetic equatorial plane. Figure 6 shows an example of these maps taken dur− ing the main and early recovery phase of the storm event of 1 June 2013. Note that this image is not a snapshot of the plasmasphere at a given time. EMMA takes 24 hours to make a full scan of the equatorial mass density, thus the temporal and spatial variations are inevitably merged. Nevertheless, this kind of maps can help us to understand the plasmasphere dynamics during the various phases of a geomagnetic storm. For example, the map represented in Figure 6 clearly shows an outward extension of the plasmasphere on the af− ternoon side which is consistent with the expected for− mation in that local time sector of a drainage plume in the aftermath of the enhanced magnetospheric convec− tion [Grebowski, 1970;Goldstein et al., 2004].

A STATISTICAL SURVEY
The procedure described in Section 2 has been ap− plied to data collected by the EMMA stations in seven selected periods, for a total of 165 days. The periods are summarized in Table 1.
In order to consider the widest possible range of ge− omagnetic disturbance level, the periods have been se− lected taking into account intervals of prolonged quiet geomagnetic conditions as well as highly disturbed events. In particular, the data set contains 13 geomag− netic storm events, 9 of them characterized by a mini− mum Dst less than −100 nT. For each period a variable number of station pairs was used, depending on the data availability and the signal quality. A total number of 41 different station pairs was used for the entire analysis, and 8 stations pairs, which were common to all periods and adequate to examine the radial varia− DEL CORPO ET AL.  tion, were further selected to perform a statistical sur− vey of the EMMA observations. These station pairs are listed in Table 2 and are also highlighted in red in Fig−  ure 1. The L−parameter of the midpoint between the stations is shown in the third column of Table 2. The analysis involved the identification (and inversion) only of the fundamental resonance frequency which here− after is simply indicated as the resonance frequency. Figure 7a shows the rate of FLR occurrence as a function of local time for each station pair. The occur− rence maximizes around 12:00 LT, reflecting the ex− pected behavior of the Q−factor of the magnetospheric resonator [Newton et al., 1978;. Chi et al., [2013] found similar behaviors using one year of Mid−continent Magnetoseismic Chain (McMAC) [Chi et al., 2005] observations in the range 1.6 < L < 3.3.

FLR OCCURRENCE
The detection rate depends also on the L−shell con− sidered. As an example, the L−dependence of the FLR occurrence at 12:00 LT is shown in Figure 7b. The FLR detection rate reaches a maximum around L ∼ 3. Chi et 9 PLASMASPHERE MONITORING BY EMMA NETWORK  al., [2013] found a similar behavior with a maximum for L ∼ 2.4. Several factors contribute to the L−depen− dence of the FLR occurrence. The amplitude of the driv− ing compressional wave generally decreases with decreasing L due to the increasing distance from the source which is typically located in the remote parts of the magnetosphere [Yumoto, 1986]. Also, at low lati− tudes (L < ∼2), where a significant portion of the field lines lies in the ionosphere, the Q−factor of the magne− tospheric resonator rapidly decreases [Yumoto et al., 1995;Menk et al., 2000]. But other factors can coun− teract the tendency of decreasing FLR detectability with decreasing latitude. For example, the presence of the plasmapause, typically located at 3 < L < 5, can alter or completely remove any FLR signature [Milling et al., 2001;Menk et al., 2004]. Also, the low resonance fre− quency at the highest latitudes may be comparable to the spectral frequency resolution, making harder its de− tection.
The maximum detection rate for each pair of sta− tions is in general greater than those found by Chi et al., [2013] and are in agreement with the results of Berube et al., [2003] and Waters et al., [1994]. As suggested by Chi et al., [2013], the discrepancy between their results and those of others can be attributed to more stringent criteria adopted by the FLR detection algorithm and to the period of low solar activity (2006−2007) in which their observations were made.
Although the FLR occurrence rapidly falls off during nighttime, it does not completely vanish. This feature observed also by Chi et al., [2013], has not been inves− tigated in detail, and represents an interesting subject for future research.
As previously stated, a cross−phase reversal may occur across the plasmapause and a negative peak is identified in the cross−phase analysis. Figure 8 shows the percentage of such occurrences with respect to the total number of FLRs detected for each station pair. As can be seen, the occurrence is generally low. This is mainly due to the limited efficiency of the technique near the plasmapause. Indeed, in order to a cross−phase reversal to occur, the radial mass density profile in the plasmapause region must be quite steep (steeper than r − 8 , [Kale et al., 2007]), and the field lines connected to the two stations must both lie in the plasmapause [Milling et al., 2001;Menk et al., 2004]. The first re− quirement is hampered by the effect of increased heavy ion contribution in the region immediately external to the plasmapause which gives rise to plasmapause mass density profiles less steep than electron density pro− files [Fraser et al., 2005]. Also, spatial integration ef− fects, which are intrinsic to ground magnetometer observations [Poulter and Allan, 1985], may contribute to make hard the plasmapause detection by the FLR technique. Nonetheless, Figure 8 clearly shows an in− crease of the occurrence by a factor of 5 between L = 3.2 and L = 3.6. This result suggests that the inner edge of steep plasmapauses is more frequently located inside this range.

DIURNAL VARIATION
A diurnal variation in the FLR frequency has been reported in several papers. At low latitudes (L < 2) the FLR frequency has been observed to exhibit an early local morning decrease [Green et al., 1993;Waters et al., 1994;Vellante et al., 2002] as well as an evening increase [Green et al., 1993] which are attributed to corresponding increase/decrease in the O + concentra− tion at dawn/dusk [Poulter et al., 1988]. Conversely at high latitudes (L > 5), an inverted U−shaped ("arch") temporal variation with maximum around (or slightly before) noon has been reported and attributed to an opposite variation in the field line length [Waters et al., 1995;Mathie et al., 1999]. Figure 9 shows the diurnal variation of f r obtained from our dataset. In order to reduce the level of fluc− tuations we generated, for each station pair, hourly values of f r from the original 30−min values and me− DEL CORPO ET AL. dian values from each hourly bin were finally com− puted. Median values obtained from a number of points (days) less than 20 are excluded in Figure 9. Error bars are not shown in order to better follow the different profiles.
A decrease of f r with increasing L is generally ob− served, except for a few cases in nighttime hours. These cases could be indicative of a plasmapause crossing, but could also be an effect of the lower sample size in these hourly bins.
The diurnal variation at all L values is characterized by a general decrease in the daytime hours which gets steeper and steeper with increasing L. At L = 2.4 and 2.6, the decrease is steeper in the early morning hours (03−07 LT), similarly to what found at low latitudes in previous investigations. A clear increase in the evening hours (a local time sector rarely covered by FLR ob− servations) is also observed at all station pairs. The FLR frequency at TARBRZ (L = 2.9) also shows an anom− alous behavior (when compared to that of the adjacent pairs) between 04−07 LT, with a minimum value of ~13 mHz at 04 LT. We attribute this behavior to the pres− ence in our dataset of quarter−wave mode events which can occur when one end of the field line is sun− lit and the other one is in darkness [Obana et al., 2008[Obana et al., , 2015. In such cases, indeed, the detected resonance frequency is approximately half of the fundamental eigenfrequency of a field line with fixed ends. We ac− tually noticed several of these events during the selec− tion phase (see for example, for the same station pair, the jump in frequency between 02:30 UT and 03:30 UT in Figure 4) and decided to keep them to highlight pos− sible systematic features in the diurnal variation and for a future study. In the inversion process we consid− ered however only the frequencies selected during day− time hours; the interval was actually extended to include one hour before sunrise and one hour after sunset on the ground, which usually appear not af− fected by the frequency change typical of the quarter wave phenomenon.
The diurnal variation at the highest latitudes (L = 4.8 and 5.5) is quite different from the "arch" shape observed at similar latitudes by Waters et al., [1995] and Mathie et al., [1999]. The present results seem to be more consistent with a daytime refilling from the ion− osphere rather than with a variation of the field line length. To test this hypothesis we evaluated the ex− pected LT variation of f r at MUOPEL (L = 5.5) using the TS05 model to describe the magnetic field geometry and two different models of the plasma density distri− bution. Average values of solar wind/Dst parameters 11 PLASMASPHERE MONITORING BY EMMA NETWORK  were considered for the TS05 model (see caption of Figure 10). The resulting temporal variation of the length of the field line traced from the MUOPEL mid− point is shown in the top panel of Figure 10. A mini− mum value of ~13.4 R ⊕ is reached at ~ 10 LT and a maximum value of ~14.7 R ⊕ occurs at ~ 22 LT. In a first run, the plasma density was assumed not to vary in LT but to depend only on the radial distance r. In particular we used an r −4 dependence for the equatorial plasma density (taken from the empirical plasmatrough electron density model of Sheeley et al., [2001]), and an r −1 dependence along the field line. The fundamental toroidal mode frequency was then com− puted using the Singer equation. The results are shown in the middle panel of Figure 10. The typical "arch" structure appears, reflecting the variation of the field line length shown in the upper panel.
In a second run we used the full plasmatrough elec− tron density model of Sheeley et al., [2001] which in− cludes a LT dependence modeled by a cosine function. The resulting LT variation of f r is shown on the bottom panel of Figure 10. The similarity with the correspon− ding experimental observation (bottommost curve in Figure 9) is remarkable. Both the model and the ex− perimental variation show a maximum at ~ 04 LT and a minimum at ~ 18 LT. The calculated frequencies are higher than those observed by a factor 1.5−2. This can be explained if one considers that the Sheeley et al. model refers to the electron density while our observa− tions depend on the mass density. The observed differ− ence between the observed and calculated frequencies would be consistent with an average ion mass of 2−4, which agrees with typical values of this parameter in the plasmathrough [Takahashi et al., 2006]. Lastly, the observed daytime frequency decrease appears more pronounced than the calculated one. It might indicate that in our data set, there is a greater weight of inter− vals characterized by strong ionospheric refilling of flux tubes depleted by previous disturbance activity (recovery phase following plasmasphere erosion) [Car− penter and Lemaire, 1997].

GEOMAGNETIC ACTIVITY DEPENDENCE
It is well known that the level of plasma density at a given radial distance may significantly depend on the previous history of geomagnetic activity. In order to examine this effect on our observations, we used two indicators of the preceding disturbance activity. The first indicator, Dst (min) , is the minimum of the Dst index in the 24 hours preceding the observation. The second indicator, Kp (max) , is the maximum of the Kp index in the same time interval. Many authors in the past used these or similar quantities in the attempt to quantify the effect of the magnetic disturbances on the magnetospheric plasma [e.g., Carpenter and Anderson, 1992;Sheeley et al., 2001;Moldwin et al., 2002;Taka− hashi et al., 2006].
We divided the original data set into three subsets for each indicator of the disturbance level and com− pared the corresponding distributions of f r and ρ eq . In terms of Dst (min) the data set was divided into the fol− lowing subsets: 1) Dst (min) > −20 nT, quiet magnetosphere; 2) −50 nT < Dst (min) ≤ −20 nT, moderately dis− turbed magnetosphere; 3) Dst (min) ≤ −50 nT, highly disturbed magnetos− phere. In terms of Kp (max) the subsets were arranged as follows: 1) Kp (max) ≤ 2 + , quiet magnetosphere; 2) 3 ≤ Kp (max) ≤ 4 + , moderately disturbed magne− tosphere; 3) Kp (max) ≥ 5 − , highly disturbed magnetosphere. Figure 11 shows the dependence of f r and ρ eq on the geocentric distance r eq for the three disturbance lev− els defined above. The figure shows the results for the Kp (max) −derived subsets, but very similar results (not shown) can be obtained by performing the same analysis on the Dst (min) −derived subsets.
For each level of activity there are three panels. The top panel shows the number of observations for each station pair, the middle (bottom) panel shows the de− pendence of f r ( ρ eq ) on the geocentric distance r eq . The vertical arrows on the bottom indicate the nominal L values of each station pair (i.e., the equatorial field line distance if only sources internal to the Earth, IGRF model, were considered). Open circles are individual observations, while open squares are the median val− ues for each station pair. Error bars indicate first and third quartile values. The results relative to two dis− tinct local times, 07:00 LT (blue) and 17:00 LT (red), are presented simultaneously to highlight morning− afternoon differences.
Individual f r and ρ eq observations are widely scat− tered at each station pair and the dispersion generally increases with increasing L and disturbance level. The figure also shows the variability of r eq for any given station pair (as predicted by the TS05 model) for the different disturbance levels, i.e., how variable is the distance of the region monitored by a given ground measurement point. Note also that r eq values are typ− ically larger than the nominal L values, and the devi− ation increases with increasing L and disturbance level. Such an effect reflects the field line stretching due to the ring current contribution modeled by TS05 (see for example, Berube et al., [2006] and Vellante et al., [2014b] for an analogous effect modeled by T02).
As expected from the diurnal variation of f r exam− ined in the previous section, plasma densities in the afternoon are generally greater than those in the morning and the difference tends to increase with in− creasing distance and level of previous geomagnetic activity. As already mentioned, this is related to the daytime refilling of the flux tubes which, especially for the highest L, are typically in a permanent state of recovery.
There is also some evidence in the radial profiles of a plasmapause formation at a geocentric distance which decreases as the disturbance level increases. The plasmapause signature is more evident for disturbed conditions (Figure 11c) when a peak in the frequency radial profile at 07 LT appears at ~3.8 R ⊕ in corre− spondence to a knee in the mass density radial profile, which is consistent with the negative cross−phase oc− currence shown in Figure 8. Altogether the results suggest that, on average, the plasmapause forms at r eq > 5 R ⊕ during quiet times, shifts inward to 3.8 R ⊕ < r eq < 5 R ⊕ during moderate disturbed conditions and reaches 3.4 R ⊕ < r eq < 3.8 R ⊕ for highly disturbed con− ditions. Obviously, these have to be considered only as average values for the considered disturbance lev− els; for example, we found evidence of a plasmapause formation at r eq values as low as 2.5 R ⊕ for severe storms. The plasmapause signature is less evident at 17 LT with respect to 07 LT probably because the day− time refilling of the depleted flux tubes tends to reduce the difference between plasmaspheric and plasma− trough density levels.
The same fitting procedure described in Section 2.3 can be applied to the median data shown in Figure 11. To exclude statistically unreliable values from the fit− ting procedure, only median values obtained by at least 10 data points were considered. The polar plots obtained by merging the fitted profiles provide infor− mation about the average configuration of the plas− masphere and the plasmatrough for different geomagnetic activity conditions. Figure 12 shows the resulting averaged equatorial mass density as a func− 13 PLASMASPHERE MONITORING BY EMMA NETWORK FIGURE11. H (a) Field line resonance frequency (fr) and plasma mass density (ρ eq ) as a function of geocentric distance (r eq ) for quiet conditions (Kp (max) ≤ 2 + ) and two distinct local times: 07:00 LT (blue points) and 17:00 LT (red points). Open squares con− nected by straight lines are the median values for each station pair. Error bars span the range between the first and third quartile. Top panel shows the occurrence for each station pair. (b) The same but for moderately disturbed conditions (3 − ≤ Kp (max) ≤ 4 + ). (c) The same but for highly disturbed conditions (Kp (max) ≥ 5 − ). tion of the radial distance and magnetic local time for quiet (panel a) and highly disturbed (panel b) periods, respectively.
Compared to Figure 11, the maps in Figure 12 pro− vide a global overview of the plasma distribution in the inner magnetosphere more directly. Also, the MLT dependence of the average plasma mass density can be better visualized. During quiet time the plasmas− phere is mostly symmetric, although a slight increase of the density in the afternoon sector is appreciable if we look, for example, at the level curve corresponding to 1000 amu cm −3 . During disturbed conditions a plasmapause signature appears near 3−4 R ⊕ in the morning/pre−noon sector, as evidenced by the rapid change of the colors in the radial direction. Moreover, an evident dawn−dusk asymmetry arises, as clearly indicated by the increasing separation with the mag− netic local time between the level curves at 100 and 1000 amu cm −3 . For r eq > 4 R ⊕ the density may drop by an order of magnitude from quiet to disturbed magnetospheric conditions.

SUMMARY AND OUTLOOK
Based on consolidated techniques, a careful and thorough procedure to derive the equatorial cold plasma mass density in the inner magnetosphere has been de− veloped. The method is meant to work with ground magnetometer meridional arrays like EMMA, a mag− netic chain recently established in Europe.
The plasma mass density is inferred from FLR fre− quencies of the geomagnetic field driven by ULF waves. FLR frequencies are detected by performing a cross− spectral analysis of signals detected by several magne− tometer station pairs. Automated algorithms for FLR detection often fail in the identification of the funda− mental harmonic at high latitudes and near the plasma− pause. In order to improve the reliability of the FLR detection, a check processing is then needed.
The main innovation in our procedure is the intro− duction of an interactive software integrating a graph− ical user interface that supports the visual inspection of the cross−spectra, speeding up the checking process.
Another important innovation concerns the intro− duction of a fitting procedure that enables to evaluate the equatorial density at specific geocentric distances, allowing to separate radial and temporal effects which are merged when analyzing the results from a single station pair. The fitting procedure also facilitates the creation of 2D−maps of the equatorial density that are very useful to figure out the general dynamics of the cold plasma in the inner magnetosphere during per− turbed periods.
The method has been thoroughly presented and ap− plied to a set of data encompassing a wide range of ge− omagnetic conditions. A database of FLR frequencies, and corresponding equatorial plasma densities, obtained from 165 non consecutive days of observations at 8 sta− tion pairs covering the range of L values 2.4−5.5 was created and a statistical analysis was performed.
A clear diurnal modulation of the FLR frequency is observed at all L values. At the lowest latitudes, the vari− ation is characterized by a rapid decrease in the early morning hours, a stagnation in the middle of the day, and an increase in the evening hours, which agrees with past observations at low latitudes [Green et al., 1993;Waters et al., 1994;Vellante et al., 2002] and with theo− retical modeling of the "diurnal breathing" of the plas− masphere [Poulter et al., 1988;Krall and Huba, 2016]. At higher latitudes, a continuous and more pronounced de− crease of the FLR frequency is observed during all day− time hours which reflects a situation of permanent recovery of flux tubes (refilling process from the iono− sphere) depleted by events of enhanced magnetospheric convection [Carpenter and Lemaire, 1997]. The radial profiles of the inferred equatorial mass density show in− deed a density increase from morning to afternoon which gets more pronounced with increasing distance and with the level of the preceding geomagnetic activity. The re− sults also confirm the formation of the plasmapause at geocentric distances that decrease as the disturbance level increases.
The present data set is already suitable to construct an empirical model of the equatorial mass density both in the plasmasphere and in the plasmatrough. This would represent an important reference model which is presently missing, except for the plasmaspheric mass density model by Berube et al., [2005] which is limited to the L−range 1.7−3.2. A dedicated paper is in preparation.
Other important topics which we are investigating using our (possibly extended) database are: a) compar− ison with in situ measurements of the electron number density (e.g., from Van Allen Probes satellites) [Kurth et al., 2015;Zhelavskaya et al., 2016] to get information on the plasma composition for different geomagnetic activity conditions; b) quantification of the plasmas− phere refilling from the underlying ionosphere after ge− omagnetic disturbances.
The final goal is to develop a fully automated algo− rithm which improves the current real−time monitoring system and can be used for space weather purposes. The dataset used in this work is large enough to derive a general description of the variability of the FLRs and plasma mass densities, especially for the dayside region. However, the availability of a reliable automated algo− rithm would increase considerably the dataset enabling a more reliable statistical description. In order to achieve this objective, there are a number of aspects that still need to be developed and implemented: 1) Identification, when possible, of the FLR harmonic structure for a correct interpretation of the har− monic number of the selected frequency (a prob− lem which we often find at high latitudes). The information provided by the knowledge of a set of harmonics can also be used to guess, case by case, a more proper distribution of the mass den− sity along the field line (rather than using a fixed power law form for all cases and L values) and then obtaining a more correct estimation of the equatorial density [Denton et al., 2014]; 2) Simultaneous analysis of more station pairs to help deciding whether to accept or reject a se− lected frequency; 3) Use of more general boundary conditions in the wave equation which take into account realistic conditions of the ionosphere, which may be im− portant in case of asymmetric conductance be− tween the conjugate ionospheres (quarter wave phenomenon) [Obana et al., 2008], and for night− time hours [Ozeke and Mann, 2005]; 4) Development of specific algorithms for the cases when the station pair is near the plasmapause and complex cross−phase structures appear with both positive and negative peaks observed at different frequencies at the same time; 5) Critical revision of the possible sources of error associated with the entire inference procedure and quantitative evaluation of the contribution of each source to the mass density uncertainty. It would be also important to extend the procedure to other magnetic arrays longitudinally separated from EMMA (like, e.g., CARISMA/McMAC) [Mann et al., 2008;Chi et al., 2013] to discriminate between tempo− ral and spatial variations of the plasma density.
We thank the Finnish Meteorological Institute (FMI), the University of Oulu (Finland), the Institute of Geophysics of the Polish Academy of Sciences (IGF-PAS), the Mining and Geological Survey of Hungary (MBFSZ) and the University of L' Aquila for contributing to EMMA. EMMA data used in this study are available from Zenodo (https://www.doi.org/10.5281/zenodo.3387216). We also acknowledge use of NASA/GSFC's Space Physics Data Facility's OMNIWeb service, and OMNI data. A.D.C, M.V, and B.H. thank the International Space Science Institute (ISSI, Bern) for hosting the activities of the team Investigating the Magnetosphere through Magnetoseismology (led by Peter Chi, USA) during which we had fruitful discussions with the other team members.