Looking for a proxy of the ionospheric turbulence with Swarm data

The present work focuses on the analysis of the scaling features of electron density fluctuations in the mid- and high-latitude topside ionosphere under different conditions of geomagnetic activity. The aim is to understand whether it is possible to identify a proxy that may provide information on the properties of electron density fluctuations and on the possible physical mechanisms at their origin, as for instance, turbulence phenomena. So, we selected about 4 years (April 2014–February 2018) of 1 Hz electron density measurements recorded on-board ESA Swarm A satellite. Using the Auroral Electrojet (AE) index, we identified two different geomagnetic conditions: quiet (AE < 50 nT) and active (AE > 300 nT). For both datasets, we evaluated the first- and second-order scaling exponents and an intermittency coefficient associated with the electron density fluctuations. Then, the joint probability distribution between each of these quantities and the rate of change of electron density index was also evaluated. We identified two families of plasma density fluctuations characterized by different mean values of both the scaling exponents and the considered ionospheric index, suggesting that different mechanisms (instabilities/turbulent processes) can be responsible for the observed scaling features. Furthermore, a clear different localization of the two families in the magnetic latitude—magnetic local time plane is found and its dependence on geomagnetic activity levels is analyzed. These results may well have a bearing about the capability of recognizing the turbulent character of irregularities using a typical ionospheric plasma irregularity index as a proxy.

www.nature.com/scientificreports/ index values obtained over the years, the existence of a power law suggests scale invariance, meaning that the same physical process is responsible for the formation of irregularities at different scales 7 . Of course, the observed values of the spectral index may depend on different factors like, e.g., latitude or altitude. Furthermore, since different spectral indices may be proxies of several underlying physical processes, it is important to discuss them in connection with the selected ionospheric region.
In this work, we investigate the scaling features of the electron density fluctuations at high latitudes, in both hemispheres, recorded on-board one of the three ESA-Swarm satellites 12 under two different geomagnetic activity conditions. This work is a generalization of a single case study recently published 13 where, for the first time, it was suggested a possible link between the high values of Rate Of change of the electron Density Index (RODI) and some specific features of the electron density fluctuations for a period of 7 days around the 2015's St. Patrick's Day storm. In this study, we consider an extended dataset covering about 4 years (April 2014-February 2018), where we identify quiet and active geomagnetic conditions. The findings suggest the existence of two families of electron density irregularities characterized by different scaling properties and different values of RODI. The analysis focuses on the middle-and high-latitude regions and for the first time climatological maps of the probability densities of the two families in the magnetic latitude-magnetic local time plane are presented as a function of geomagnetic activity for both the hemispheres. Thus, our findings confirm what has been achieved in the case of a single storm and, at the same time, demonstrate that the existence of the two families is independent on the geomagnetic activity level. What changes is the weight of one family compared to the other as a function of the geomagnetic activity level and their probability density distributions in the magnetic latitude-magnetic local time plane. The results also confirm the possibility that RODI, which quantifies the relevance of ionospheric irregularities, may provide information on the features of the electron density fluctuations and on the possible physical mechanisms at their origin, as for instance, turbulence phenomena. This is of great importance in the framework of space weather. Indeed, ionospheric plasma density irregularities can strongly influence the quality of electromagnetic signals that propagate in space and consequently can influence the performance of those systems which are based on their propagation as for example the global navigation satellite systems (GNSS) and the global positioning system (GPS) on which many of our infrastructures depend either directly or indirectly 14 . The identification of regions characterized by turbulence processes could be the keystone for predicting where the degradation and loss of GPS signals have a higher probability to occur.

Data
The present work focuses on the analysis of the scaling features of the electron density fluctuations in the midand high-latitude F-region ionosphere under two different conditions of geomagnetic activity. For this reason, we consider 1 Hz electron density measurements collected by the Langmuir Probes (LPs) on-board ESA Swarm A satellite, during a period of about 4 years (from April 2014 to February 2018). Using the Auroral Electrojet (AE) index, which is a good proxy of the geomagnetic activity level at mid/high latitudes, we select two distinct datasets corresponding respectively to geomagnetically quiet (AE < 50 nT) and active (AE > 300 nT) periods. The electron density measurements, retrieved from the ESA ftp repository (ftp://swarm -diss.eo.esa.int) are successively used to evaluate RODI, which is defined as the standard deviation of the Rate of change Of electron Density (ROD) in a sliding window of fixed size, in our case of 10 s 13,15,16 . We focus our attention on mid-and high-latitude ( > |50| • ) regions in the Northern and Southern Hemisphere analyzing data in the Quasi-Dipole (QD) latitude 17 and magnetic local time (MLT) coordinates.

Method of analysis
One of the most striking features of turbulent signals is the scale-invariant nature of the fluctuations shown at different spatial and temporal scales. Indeed, turbulence is an inherently multiscale phenomenon, showing fluctuations either in space or in time (or both) and resulting in a noisy and chaotic dynamics. These multiscale fluctuations are due to the fact that in turbulent media the disturbances/excitations occur at macroscopic scales that are far from the microscopic scales at which dissipation takes place. According to the general picture of the Richardson's cascade, turbulent media develop an inertial range where the mix of excited spatial and temporal scales is characterized by specific symmetry features (a phenomenon known as scale invariance). In that range, the energy injected at the largest scale is transferred to the smallest ones (where dissipation occurs) without a significant dissipation. The manifestation of scale invariance in a fluctuation field is the absence of any characteristic length or timescale for the emerging fluctuation structures, i.e., the occurrence of self-similarity in a statistical sense. In the case of a turbulent signal, a simple way to investigate the scaling features of the moments of the distribution of signal increments at different spatial/temporal scales is represented by the so-called q th -order structure function at scale l, namely S q (l) . Here, we compute the generalized q th -order structure function of the electron density N e as a function of the time delay τ , i.e., where t is the time and ... stands for a statistical average.
In presence of scale invariance S q (τ ) has a power-law behavior as a function of τ: where γ (q) is the so-called scaling exponent, which allows us to characterize the scaling nature of the increments of the signal under investigation (in our case N e ).
Here, we investigate the first-and second-order scaling exponents. The first-order scaling exponent γ (1) , also known as Hurst exponent (H) 18 , is a proxy of the persistent character of the analyzed time series increments, i.e., www.nature.com/scientificreports/ of the long-term memory of the time series 19 . It can assume any value in the range between 0 and 1. When the time series are characterized by values of H greater than 0.5, this means that the increments of the time series tend to cluster in a direction, and positive increments are more probably followed by positive increments while negative increments are more probably followed by negative ones. In this case the time series are characterized by a long-term positive autocorrelation. Conversely, when the scaling features of series are characterized by values of H lower than 0.5, the increments of the series tend to switch between positive and negative values, namely a positive increment is more likely to be followed by a negative one and vice-versa. In this case, the time series have the tendency to return to a long-term mean. Lastly, a value of the Hurst exponent equal to 0.5 indicates an absence of correlation in the increments (which are representative of fluctuations at the investigated time scale) of a time series that happens, for example, in the case of a Brownian time series. The second-order scaling exponent, γ (2) , provides us information on the spectral features of the signal under investigation. In particular, for signals whose Fourier power spectral density (PSD) functions follow a power law behavior, the power spectral density exponent, β , is related to the second-order structure function scaling exponent as follows This is a consequence of the Wiener-Khinchin theorem relating the power spectral density of a signal to its autocorrelation function 20 . The study of the second-order scaling exponent is important for understanding the different physical processes (turbulence processes, stochastic processes, avalanche processes and so on) at the origin of a measured signal.
From a practical point of view, it is worth highlighting that the electron density time series recorded by Swarm satellite, like many real time series, are generally non-stationary, reflecting the diverse physical processes occurring in the different regions crossed by the satellite along its orbits. However, stationarity may occur locally, i.e., the scaling features may acquire a local character as in the case of a multifractional Brownian motion [21][22][23] . To evaluate the local scaling features of the electron density time series, we apply a local structure function method based on the detrended structure function analysis proposed by De Michelis et al. 24 . According to this method we use overlapping moving windows to remove possible large-scale variations in the time series, which can affect the correct estimation of the scaling features. In detail, we fix an overlapping moving window of N points and consider the portion of the time series falling in this window. We detrend this portion of the time series using a linear function and successively apply the classical structure function analysis. Here, we use a moving window of size equal to N = 301 points that means a moving window of 301 s. Lastly, the same analysis is repeated on the new window obtained by shifting the previous window of 1 point, i.e., of 1 s.
Taking into account the size of the moving window, we can only evaluate the first two structure functions and consequently only the first two scaling exponents ( γ (1) ≡ H and γ (2) ). Indeed, to obtain a good estimation of a scaling exponent the corresponding structure function must be computed on a window characterized by a size at least of the order of 10 q . The structure functions are evaluated for different values of the time delay τ chosen in the range between 1 s and 40 s. After evaluating the first-and second-order structure functions in each moving window, we can estimate the first-and second-order scaling exponents whose values are associated with the position of the satellite in the instant of time corresponding to the center of the considered time window. The error associated with the first-and the second-order scaling exponents is approximately 7% and < 10%, respectively. These errors have been estimated via both a simulation of multifractional Brownian motions as described in Consolini et al. 23 and a fitting procedure.
Concerning the stationarity of the analyzed time series, we can observe that the reduction of the time interval over which measurements are acquired corresponds to decrease the possibility of observing fast variations. Furthermore, since in several cases the spectral features associated with the actual measurements are characterized by power-law spectral densities with a spectral exponent β in the range from − 3 to − 1, the corresponding time series are generally non-stationary with stationary increments in a range of time scales bounded below and above 25 . Thus, in our case the structure function analysis can be considered to be reliable. Moreover, although the analysis is done in the time domain, we can assume that it is also valid in the spatial one. Indeed, the analysis is applied to a range of time scales corresponding to those ranges where low-frequency temporal fluctuations are mainly due to Doppler-shifted and stationary spatial variations (see e.g., Refs. 26,27 and references therein). This means that the scaling properties, that we find and that are valid in the range 1-40 s, can be considered valid also in the spatial domain which ranges between ∼ 8 and ∼ 300 km, considering the orbital velocity of Swarm satellite that is ∼ 7.6 km/s.
Using the two first scaling exponents, H and γ (2) , it is possible to introduce an intermittency parameter, I, computed according to the following expression, This parameter provides a measure of the inhomogeneity with which the energy is distributed in the cascade process from the largest to the smallest scales. More precisely, intermittency is a measure of the relevance of anomalous scaling features related to a multifractal character of the signal increments. It is a measure of the departure from a simple scaling, i.e. www.nature.com/scientificreports/ being d 2 γ (q)/dq 2 ≤ 0 ∀ moment orders p and q . In other words, the scaling exponents do not follow a linear dependence on q, but show a convex trend. High values of intermittency in a time series are generally due to a non uniform distribution of the energy at all the time scales and manifest in temporal bursts. On the other hand, if the time series refers to a spatial crossing, the occurrence of intermittency can be related to a multifractal structure of the energy distribution in space as it occurs in the case of intermittency in turbulence.

Results and discussion
To limit a considerable loss of points due to the moving window approach necessary to evaluate the first-and second-order scaling exponents, we evaluate them using a single electron density time series covering the whole selected period, 47 months with a time resolution of 1 s. In this way, the time series relative to the first-and second-order scaling exponents have only 300 points less than the original one. Indeed, due to the moving window applied to evaluate the scaling exponents, the scaling exponent time series begin and stop 150 points after and before the original electron density time series, respectively. Clearly, data refer to the entire globe and to different values of AE index. Thus, we select two different datasets, consisting of values of the scaling exponents associated with AE < 50 nT and with AE > 300 nT, respectively. Successively, we divide the QD-Latitude-MLT plane into a regular grid of 1 • × 4 min and compute polar view maps applying a weighted Gaussian interpolation scheme.  www.nature.com/scientificreports/ scaling exponent lower than 0.9 (black line). Considering the relationship between the second-order scaling exponent and the spectral exponent of the power spectral density (Eq. 4), this means that the N e fluctuations characterized by an antipersistent behavior are often associated with spectral exponents lower than 2. These spectral features of the electron density fluctuations agree with results by De Michelis et al. 13 . However, we can notice that the second-order scaling exponent is characterized by values which tend to be slightly different in the high-latitude regions. For example, in the polar cap region, the second-order scaling exponent is characterized by a value around 0.66, meaning electron density fluctuations whose power spectral density has a spectral exponent β ≃ 5/3 . In the auroral regions, the second-order scaling exponents are, on the other hand, characterized by a spectral exponent β ∼ 1.9 − 2 . Conversely, the regions outside the auroral oval are characterized by N e fluctuations with a persistent character that are often associated with spectral exponents greater than 2. These features, which depend on the magnetic latitude and magnetic local time, are found during quiet and active periods and therefore seem to be a well-defined characteristic of N e fluctuations of the high latitude ionosphere.
Moving to the intermittency features, polar maps show that the region characterized by the highest values of intermittency tends to contract towards the pole or to expand towards the equator according to the level geomagnetic activity, i.e., moving from quiet to active periods. This suggests that the low-latitude edge of the auroral oval can be a region where turbulence processes may play a key role in the plasma dynamics and where the energy cascade is unevenly transferred from large to small scales. This assertion can be made on the basis of the relation existing between the spectral features of passive scalar quantities, as it is expected to be the electron density, and the ones corresponding to the velocity field in the case of turbulent media 28,29 . Thus, under the assumption of a definite link between the electron density fluctuations and the plasma velocity ones, we can observe that in these regions the energy dissipation associated with the turbulent fluctuations is strongly inhomogeneous being indeed localized in space.
Our previous analysis on the St. Patrick's Day storm 13 suggested a possible relationship between the high values of RODI and some specific features of N e fluctuations. To further investigate this aspect we repeat the analysis done by De Michelis et al. 13 for each of the two different levels of geomagnetic activity in both the hemispheres. Specifically, we evaluate the joint probability distributions between RODI and the first-order scaling exponent, the second-order scaling exponent and intermittency, respectively. Corresponding results are reported in Figs. 2 and 3 in the case of quiet and active periods in the Northern and Southern Hemisphere, respectively. The study confirms from a statistical point of view that the joint probability distributions identify two families which are characterized by different scaling properties of N e fluctuations. High values of RODI are generally associated with N e fluctuations characterized by an antipersistent ( H < 0.5 ) behavior and a second-order scaling exponent γ (2) < 1 . Conversely, low values of RODI are mainly associated with N e fluctuations characterized by a persistent ( H > 0.5 ) behavior and a second-order scaling exponent γ (2) > 1 . However, it should be noted that some of the low values of RODI can be also associated with antipersistent N e fluctuations. These two families can be identified in both hemispheres and their existence is independent on the geomagnetic activity level. To clarify the meaning of persistency in our context we can say that a persistent signal tends to generate large scale fluctuations (structures), while antipersistency has to be associated with the formation of small scale structures. Thus, the generation of irregularities at small scales should be associated with the antipersistent character of the analyzed signal ( H < 0.5 ). The value of RODI that reasonably seems to indicate the separation between the two families is the one for which Log(RODI) = (3.25 ± 0.05) . This threshold value is obtained by determining the average position of the minima between the power density functions (PDFs) of the two families along the y axis ( Log(RODI) ) in the range between 3 < Log(RODI) < 3.5 and considering all the analyzed quantities and cases. It is highlighted by the horizontal black dashed lines in each panel of Figs. 2 and 3. This value is slightly different from what was previously achieved by analyzing the St. Patrick's Day storm 13 where we found a value around Log(RODI) = 3.5 . The difference is due to a more rigorous methodology in determining the locus of minimal superposition between the two families that it has been possible thank to a wider statistics.
In order to establish whether the existence of the two families has any latitudinal dependence, we consider progressively smaller portions of the polar region. As an example, in the first column of Fig. 4 we report the joint probability distributions between RODI and the first-order scaling exponent (first row), the second-order scaling exponent (second row) and intermittency (third row), respectively obtained considering QD-latitudes higher than 50 • in the Northern Hemisphere during quiet period. The next 2 columns display the plots relative to the same joint probability distributions evaluated considering different values of the minimum latitude. What Fig. 4 points out is the existence of a dependence on latitude of the relative weight between the two different families. The family of electron density fluctuations with an antipersistent behavior and a second-order scaling exponent γ (2) < 1 associated with high RODI values is more definite at high latitudes. The family of electron density fluctuations with a persistent behavior and a second-order scaling exponent γ (2) > 1 associated with low RODI values disappears for QD-latitudes higher than 65 • −70 • .
Thus, these figures suggest that high values of RODI are mainly associated with N e fluctuations characterized by scaling properties which recall those resulting from fluid-type turbulence processes. Conversely, other processes (a turbulence of different origin or other instability mechanisms) could be at the origin of the N e fluctuations associated with the other family. Moreover, both families seem to show a similar degree of intermittency which suggests that the associated fluctuation fields are non homogeneous. Similar results are obtained analyzing the disturbed period for the Northern Hemisphere and the quiet and disturbed periods for Southern Hemisphere.
To visualize the distributions of these two different families in the QDLat-MLT plane we try to separate them. Figure 5 reports a sketch of the identification of the two families, indicated in red and purple, done using the joint probability distribution between RODI and the first-order scaling exponent (H) evaluated for magnetic latitudes higher than 50 • in the Northern Hemisphere for the quiet period. The selection between the two families is done using as a threshold for the RODI a value Log(RODI) = 3.25 . The obtained spatial distributions of the probability densities of the two families in the QDLat-MLT plane associated with high (top panels) and low (bottom panels)  Fig. 6. The same quantities relative to the Southern Hemisphere are reported in Fig. 7.
To compute the spatial distribution of the probability density function (PDF) reported in Figures 6 and 7 we use a kernel method based 30 on a 2D Gaussian kernel, k(u x , u y ) , defined as being ǫ the kernel width (here of the order of 1 degree), while u 2 x + u 2 y is the squared distance between the point where the probability is computed and the dataset points (i.e., u 2 x + u 2 y = (x i − x) 2 + (y i − y) 2 where {(x i , y i )} refers to the dataset). Thus, the PDF is given by www.nature.com/scientificreports/ where N is the number of points of the dataset. Furthermore, because the spatial distribution of the satellite measurements is not uniform in the QDLat-MLT plane, the obtained spatial distributions have been renormalized taking into account the corresponding ones of the satellite observations. The location of these two different families in the QDLat-MLT plane clearly shows that the first population, which is characterized by antipersistency, β < 2 and high values of RODI, is mainly located at high latitude inside the auroral oval and the polar cap, while the second population, which is characterized by persistency, β > 2 and low values of RODI, is located at lower latitudes, mainly outside the auroral oval. With the increase of the geomagnetic activity, the distribution of the first population tends to expand towards lower latitudes, in perfect analogy with what happens to the auroral oval regions.

Summary and conclusions
Plasma irregularities and turbulence processes characterize both the low-and high-latitude ionosphere. We used observations from Swarm constellation to study the scaling features of the electron density fluctuations in the mid-and high-latitude F-region ionosphere during geomagnetically quiet and active periods. The aim was to investigate the possibility that RODI, which quantifies the relevance of ionospheric irregularities, may provide information on the properties of the electron density fluctuations and on the possible physical mechanisms at their origin, such as, for instance, turbulence phenomena. It is known that due to the solar wind-magnetosphere-ionosphere coupling, the high-latitude ionosphere is characterized by irregularities (plasma density fluctuations and structures) ranging from hundreds of kilometers down to meters. Processes such as for example plasma convection, particle precipitation, currents both parallel and perpendicular to the geomagnetic field, plasma instabilities and neutral fluid dynamics, such as gravity waves, have been proposed to account for the generation of these irregularities at large and intermediate scales, which mainly characterize the cusp region, the polar cap and the nightside sector auroral ionosphere. Under certain conditions these plasma irregularities can become unstable, producing scale-invariant structures via turbulent phenomena. Indeed, phenomena such as the � E × � B gradient drift instability (GDI) or the current convective instability (CCI), can both directly and indirectly produce a cascade process through which smaller  www.nature.com/scientificreports/ size irregularities are generated 7,31 . These scale-invariant features of plasma density fluctuations associated with irregularities manifest in different spectral features i.e., in different power spectral exponents, and persistency features.
One of the main findings of our statistical analysis is that we can clearly identify two different families of plasma density fluctuations characterized by different mean values of scaling exponents and RODI. This finding suggests that, in general, two main different classes of physical phenomena can be at the origin of the different scaling features. Furthermore, a clear different localization of the two families in the QDLat-MLT plane is found. The first population, which is characterized by antipersistency, β = γ (2) + 1 < 2 and high values of RODI, is mainly located inside the auroral oval where particle precipitation dominates. The second population, which is characterized by persistency, β > 2 and low values of RODI, is located at low latitude, mainly outside the auroral oval. These findings suggest that, in general, the first population could be due to turbulent mechanisms generated by GDI. Indeed, the observed average value of the spectral density exponent β ∼ 1.7 is in agreement with previous findings 2,9,32 and with theoretical predictions for this turbulent mechanism 7 . Furthermore, the polar cap and, in particular, its trailing side are regions where due to the large-scale convection GDI may be an active phenomenon. At the same time, we cannot exclude that other mechanisms could be at the origin of the observed spectral features, as for instance the occurrence of coherent structures and nonlinear wave interactions 32 . However, we remark that coherent structures are also a product of an intermittent turbulent cascade that is a common feature of strong turbulence, so that the observed spectral features could be due to the occurrence of an intermittent strong turbulence generated by GDI.
For what concerns the second population, this is characterized by spectral features ( β ∼ 2 ) that remind random turbulent fluctuations in regions where a density enhancement is observed in correspondence with plasma coming from the inner magnetospheric regions (see, e.g. 32,33 ). The intermittency degree associated with this population seems to be slightly lower than that of the first one, suggesting that the observed density fluctuations could be due to mechanisms different from strong/fully developed turbulence. To support this hypothesis www.nature.com/scientificreports/ we remind that electron density could be considered as a passive scalar in turbulent medium. Indeed, it has been reported in the literature 34 that passive scalar quantities can show a certain degree of intermittency also in the case of purely stochastic velocity fields, so we can speculate that the low degree of intermittency of this population with respect to the first one could be the counterpart of a more stochastic nature of the driving field, i.e., the electrostatic field 7 . However, we cannot exclude that other mechanisms such as CCI and Kelvin-Helmoltz instability could be responsible for the observed features. Thus, a careful discussion of this second population needs further theoretical/numerical investigations, including also the possible effects of neutrals 35 . Anyway, a detailed discussion of the different instability mechanisms at the origin of the observed scaling features requires the investigation of other physical quantities, such as magnetic and electric field fluctuations. Some preliminary analyses of the magnetic field scaling features as a function of different interplanetary magnetic field orientations 36 and geomagnetic activity levels 24,37 have shown that there is not a simple relationship between the scaling properties of these two physical quantities, suggesting that a more detailed analysis should be performed in order to identify and characterize the instability mechanisms at work. This is postponed to future works, being the aim of the present study limited to the analysis of the link between RODI and the electron density scaling features. An implication of our analysis is the possibility to consider RODI as a proxy to identify different families of irregularities and the different nature of physical processes at the base of their origin. Specifically, our study showed that RODI could be used to detect plasma density irregularities caused by fully developed turbulent phenomena.
This is a positive outcome for space weather purposes. Indeed, ionospheric plasma density irregularities can strongly influence the quality of propagating electromagnetic signals in terms of delay, distortion and signal loss, adversely affecting systems as GNSS and GPS on which many of our infrastructures depend either directly or indirectly. Recently, it has been found 38 that at high latitudes plasma density irregularities, and in particular those associated with large plasma gradients, play an important role in the GPS signal loss. In fact, the GPS signal loss events seem to occur with high probability when the electromagnetic signals propagate inside regions characterized by large electron density gradients. The climatological maps of GPS signal loss based on Swarm measurements reveal that they occur in the Northern Hemisphere mainly around the cusp region and along nightside auroral latitudes 38 , that are precisely the regions identified by our analysis in terms of families of electron density fluctuations characterized by antipersistency, a second-order scaling exponent γ (2) < 1 and high values of RODI. This suggests that GPS signal loss events could be linked to the occurrence of turbulence processes characterizing the plasma density fluctuations and structures. The identification of regions characterized by turbulence processes could be the keystone for predicting the ones where the degradation and loss of GPS signals have a higher probability to happen.