*Article* **Measurement of Aspect Angle of Field-Aligned Plasma Irregularities in Mid-Latitude E Region Using VHF Atmospheric Radar Imaging and Interferometry Techniques**

**Jenn-Shyong Chen 1,\*, Chien-Ya Wang <sup>2</sup> and Yen-Hsyang Chu <sup>3</sup>**


**Abstract:** Multireceiver and multifrequency radar imaging were carried out with the 46.5 MHz MU radar in Japan (34.85◦N and 136.10◦E) to examine the aspect sensitivity of field-aligned plasma irregularities (FAIs) in the mid-latitude ionosphere E region. A radar beam was directed to geographic north and at 51◦ zenith angle, which was normal to the geomagnetic field line around 100 km height. Nineteen receivers and five carrier frequencies were used for radar imaging to retrieve the power distribution in the radar volume, and then the aspect angle along the geomagnetic field line was calculated according to the angular power distribution. Retrieval algorithms such as Fourier, Capon, and norm-constrained Capon (NC-Capon) were employed, in which the NC-Capon was applied to FAIs for the first time and found to be more suitable for the present study. The aspect angles estimated by the NC-Capon ranged between 0.1◦ and 0.4◦ mostly, and averaged around 0.2◦, which were the same order to the previous measurements with radar interferometry (RI), made for equatorial electrojet irregularities and the lower mid-latitude sporadic E region. For comparison, RI-estimated aspect angles were also investigated and found to be close to that of NC-Capon, but distributed over a wider extent of angles.

**Keywords:** radar imaging; aspect angle; field-aligned plasma irregularities; mid-latitude E region; norm-constrained Capon; VHF radar

### **1. Introduction**

Aspect sensitivity of refractive irregularities in the atmosphere describes a dependence of radar echo intensity on angular direction. In horizontally structured atmosphere, the radar echoes are usually strongest near the zenith and drop off with off-zenith angle, which is evident in observations throughout the low and middle atmosphere with the radar at MF, HF, and VHF bands. For the field-aligned plasma irregularities (FAIs) in the ionosphere, however, the radar echoes are strongest at the beam direction normal to the geomagnetic field line, and decay very quickly with off-perpendicular angle along the geomagnetic field line. Aspect angle is a measurement of aspect sensitivity, which can be defined by the half width of half power, or the standard deviation of Gaussian fitting in the angular power distribution; the former is about <sup>√</sup>*ln*<sup>4</sup> (ؒ1177) times the latter if the angular power distribution is characterized by a Gaussian function. The order of aspect angle in the lower atmosphere is within several degrees [1,2]. By contrast, it can be as small as several times 0.1◦ or less in FAI echoes, depending on latitude and height [3–6]. It has been addressed in [3–6] that the FAIs' aspect angle in the direction of geomagnetic field line relates partly with the values of electron-neutral collision frequency (νe), electron gyrofrequency (Ωe), ion-neutral collision frequency (νi), ion gyrofrequency (Ωi), the ratio νe/Ωe, among others,

**Citation:** Chen, J.-S.; Wang, C.-Y.; Chu, Y.-H. Measurement of Aspect Angle of Field-Aligned Plasma Irregularities in Mid-Latitude E Region Using VHF Atmospheric Radar Imaging and Interferometry Techniques. *Remote Sens.* **2022**, *14*, 611. https://doi.org/10.3390/ rs14030611

Academic Editor: Fabio Giannattasio

Received: 26 December 2021 Accepted: 17 January 2022 Published: 27 January 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and moreover, is also subjected to nonlinear coupling process of unstable waves in the plasma irregularities.

Because of highly localized FAI echoes in the radar volume, and much, much smaller aspect angle in the ionosphere than in the lower atmosphere, many of the approaches/techniques proposed for measurement of aspect angle in the lower and middle atmosphere, such as those used in [1,2,7–13] (i.e., beam swing, spaced antenna, spectral width, etc.), are difficult to apply to FAIs. A proper approach to measurement of FAIs' aspect angle is the radar interferometry (RI) technique that was first applied to equatorial ionosphere by Kudeki and Farley [3] and observed the aspect angles of electrojet irregularities ranging between 0.1◦ and 0.4◦ generally. A later study for equatorial spread F irregularities found an even smaller value of aspect angle, i.e., a few times 0.01◦ [4]. Additionally, RI was applied to the lower mid-latitude sporadic E region (~25.6◦N and ~120.8◦E geography; ~16.33◦N and ~167.34◦W in geomagnetic coordinate) but with a different process and definition [14], in which the aspect angles of layer-type and clump-type plasma irregularities were estimated, and were found to be the same order of magnitude as that of equatorial electrojet irregularities.

In this study, we attempted to apply the approach of coherent radar imaging (CRI) to estimate the aspect angle of mid-latitude E region FAIs. CRI has been employed for study of aspect sensitivity in the lower and middle atmosphere [2,15,16]. The technical requirement of CRI is the use of separate antennas as independent receiving channels to collect the radar echoes. The multiple-channel echoes received can be employed to retrieve the in-beam angular power distribution through a beamforming process. The algorithm of beamforming can be linear-based, e.g., the Fourier method, or adaptive, such as the Capon [15] and norm-constrained Capon methods [17–20], and others. CRI can be regarded as an advanced version of RI.

To meet the requirement of CRI, the Middle and Upper (MU) Atmosphere radar, a pulsed radar operated at 46.5 MHz and maintained by Kyoto University, Japan, was employed in this study. The MU radar is located at 34.85◦N and 136.10◦E in geographic coordinates, and 26.36◦N and 153.78◦W in geomagnetic coordinates. A total of 25 receiving channels can be assigned, and moreover, frequency hopping of five carrier frequencies between radar pulses is available. Such flexible capabilities of operation allow radar image processing not only in 1D range imaging and 2D angular imaging, but also 3D spatial imaging. To this end, the retrieval methods used consistently, such as the Fourier, Capon, and norm-constrained Capon methods, were examined to find their suitability for the present study. As far as we know, the norm-constrained Capon method has not yet been applied to full 3D radar imaging. In addition to radar imaging, several sets of subarrays aligned approximately with the direction of the geomagnetic field line are also available to test the RI technique in measurement of aspect angle. Therefore, comparison of measurements is possible to validate the usability of radar imaging and RI, or evaluate which kind of method is more suitable for estimating the aspect angle of FAIs with the MU radar configuration. To complete the study, the effect of signal-to-noise ratio (SNR) on the estimate of aspect angle was also examined.

The rest of this paper is organized as follows. In Section 2, measurement techniques are described, including experimental setup and data processing of radar imaging and RI. Section 3 provides the results and discussion, and Section 4 states the conclusions.

#### **2. Measurement Techniques**

Radar experimental setup, radar imaging, and radar interferometry are described in Sections 2.1–2.3, respectively.

#### *2.1. Radar Experimental Setup*

Figure 1 shows the antenna array configuration of the MU radar, in which each plus symbol denotes a single antenna with an independent transceiver. The diameter of the array is about 103 m and the array can be partitioned into 25 subarrays, as denoted by the letters from A1 to F5. The separation of adjacent subarrays is 19.55 m, and central carrier frequency is 46.5 MHz. A detailed description of the radar characteristics can be found through the website http://www.rish.kyoto-u.ac.jp/English/MU/index.html (accessed on 25 December 2021). Some operational parameters related to the present study are described below and listed in Table 1.

**Figure 1.** Configuration of MU radar array. Each plus symbol denotes a single antenna with an independent transceiver. The diameter of the array is about 103 m, and the subarrays denoted by the letters from A1 to F5 are individual receiving arrays. The 19 subarrays in the interior and the seven subarrays encircled are employed in radar imaging. The geomagnetic south pole (Smag) is about 7◦ west of geographic north.

**Table 1.** Radar parameters.


Full array was used for transmission of radio waves, and the radar beam was steered to geographic north and at the zenith angle of 51◦. As a result, the radar beam bearing was perpendicular to the geomagnetic field line around the height of 100 km. Moreover, five carrier frequencies were employed in sequential radar pulses that were operated repeatedly. In reception, the whole array and the 19 subarrays in the interior of the array were assigned to 20 receiving channels to collect the radar echoes from five different carrier frequencies. This multireceiver and multifrequency operation allows us to execute 1D range imaging, 2D angular imaging, and 3D spatial imaging. The imaging here means the retrieval of power distribution within the radar volume.

The sampling range interval was between 120.0 and 196.8 km. The corresponding sampling height interval was between about 75.518 and 123.850 km, which was estimated by projecting the sampling range at vertical direction and ignoring the curvature of the Earth surface. Sampling time was 0.0075 s for each carrier frequency returns (=5 × IPP = 5 × 0.0015 s). In calculation of covariance function for radar imaging, 256 raw data points were taken for ensemble average, resulting in an estimate of 1.92 s (=256 × sampling time = 256 × 0.0075 s).

#### *2.2. Radar Imaging*

The retrieval methods used for radar imaging in this study were Fourier, Capon, and norm-constrained Capon. We compared their effectiveness for deriving aspect angle in the 2D and 3D imaging processes. For 2D angular imaging, only multireceiver data are needed. For 3D spatial imaging, both multireceiver and multifrequency data are required. These methods of radar imaging were examined and employed in several studies in the research field of atmosphere [15–25]. A brief mathematical description of the radar imaging employed in this study is provided in Appendix A.

Two-dimensional angular imaging yields the power density (termed brightness distribution hereafter) as a function of angle for each sampling gate, and there could be more than one local maximum in the brightness distribution. Nevertheless, only the one with the maximum brightness level was used in the study. That is, the mean location of the brightness center encircled by the highest contour level was determined first, and then Gaussian-function fitting was applied to the brightness values around the estimated mean location of the brightness center to obtain a brightness width, i.e., the Gaussian-fitted standard deviation. The fitting was executed only in the direction along the geomagnetic field line, the aspect angle of which we surveyed. It was estimated that the geomagnetic south pole was about 7◦ west of geographic north in the field of view of the MU radar in the observation. Therefore, it was expected that the brightness distribution of FAI echoes is field-aligned and usually has a slant of about 7◦ in the 2D imaging surface viewed from the radar beam direction [25]. Notice that the imaging surface is the transverse plane of the radar beam, and the origin is at the radar beam bearing. Abscissa and ordinate of the imaging surface are zonal and elevation directions relative to the radar beam bearing, respectively (Figure 2b,c; explained later).

Since the aspect angle of FAIs could be as small as only 0.1~0.5◦, the imaging was processed with an angular step of 0.05◦ in both zonal and elevation directions. To confirm the Gaussian-fitting result, we also used the peak-find process built in the computing software (i.e., MATLAB) to obtain the half-power peak width. The half-power peak width or the Gaussian-fitted standard deviation is an indication of the aspect angle of the FAIs, in which the former is about <sup>√</sup>*ln*<sup>4</sup> (ؒ1177) times the latter if the angular brightness distribution around the mean location of the brightness center, in the geomagnetic field-aligned direction, is characterized by a Gaussian function.

It should be noticed that the aspect angle estimated from 2D angular imaging is a mean value within the range interval of the sampling gate. If one expects to observe the spatial structure and estimate the aspect angles at different ranges in the radar volume, 3D spatial imaging is one of the solutions. In the 3D radar imaging executed in this study, the range step of imaging was 37.5 m, giving 17 range locations in the sampling gate of 600 m length, i.e., at the range locations of −300, −262.5, −225, ... , 300 m relative to the range center of the sampling gate. A smaller range step is certainly available, but it consumes more computation time and may not be helpful for improving the range resolution due to the finite number of carrier frequencies used in observation. For each range location, the Gaussian fitting and peak-find process can be performed for the 2D-angular brightness distribution to obtain two brightness widths in the geomagnetic field-aligned direction, which denote two possible values of aspect angles at that range location. The results from 19 subarrays and 7 subarrays (indicated by the dashed circles in Figure 1) are compared and shown later.

**Figure 2.** One-, two-, and three-dimensional radar imaging. (**a**) Left: radar echoes; middle: 1D range imaging; right: one example of range-Doppler spectrum; (**b**,**c**) are the brightness distribution obtained from 2D and 3D radar imaging with Fourier, Capon, and norm-constrained Capon methods, respectively. Two examples are shown. Brightness values in (**b**,**c**) are self-normalized and displayed in dB scale. SNR is given in the Fourier panel title, present in linear scale, and following the SNR value is the data time in the format hh:mm:ss LT, on 26 July 2019. The numbers in the Capon column of (**c**) are the range positions of the slices of brightness distribution in the sampling gate, which are relative to the central range of the sampling gate.

#### *2.3. Radar Interferometry*

Radar interferometry (RI) has been widely used in atmospheric studies such as mean angle of arrival of radar echoes [26], lightning detection [27], meteor radar wind [28], among others. Use of RI in measurement of the aspect angle of FAIs was first reported by Kudeki and Farley in 1989 [3] for the equatorial electrojet region. In their approach, normalized cross spectra were calculated for all receiver pairs and then a suitable fitting of a parabola to the absolute values (termed coherence) of the cross spectra was made to retrieve aspect angles at various spectral lines (i.e., Doppler shifts). Lu et al. [6] extended the study of [3] and provided the following normalized cross spectrum:

$$S\_{mn}(kD\_{mn}, w) = \varepsilon^{jkD\_{mn}(\theta)\_w} \exp\left[-\frac{1}{2}k^2 D\_{mn}^2 \theta\_{rms}^2(w)\right],\tag{1}$$

where *k* is wave number, *w* is Doppler shift, *Dmn* is the length of the baseline between receiver antennas *m* and *n*, 〈*θ*〉*<sup>w</sup>* is the mean angular location of the irregularity pattern at the Doppler shift *w*, and *θrms*(*w*) denotes the angular spread value relative to 〈*θ*〉*w*, which can be regarded as the aspect angle at *w*. Gaussian function is then employed to fit the coherence of the cross spectra to derive the parameter *θrms*.

In this study, we did not examine the aspect angle at different Doppler shifts. Therefore, the cross-correlation function *Smn*(*kDmn*) was estimated in time domain. Unfortunately, there was no set of subarrays with the baseline parallel to the direction of the geomagnetic field line exactly. As shown in Figure 1, the three subarrays (A2, F5, D2) and the two sets of five subarrays (A2, F4, F5, D4, D2; A2, A4, F5, C4, D2) had the alignments of about 7.60◦ east of geographic north, and thus had about 14.85◦ difference from the geomagnetic southern pole. Another set of subarrays, F3, F4, F5, C4, and C3, was aligned in the direction of about 28.07◦ west of geographic north, and thus had a difference of about 20.82◦ from the geomagnetic south pole. These sets of subarrays were used for RI in this study.

#### **3. Results**

The aspect angles retrieved from radar imaging and RI are shown and discussed in this section. Moreover, the effect of SNR on the estimate of aspect angle is also examined.

#### *3.1. Radar Imaging in 1D, 2D, and 3D*

First, the benefits and shortages of the Fourier, Capon, and NC-Capon methods in 2D and 3D radar imaging were investigated, as shown in Figure 2. Figure 2a shows the radar echoes collected by the full array of the radar, where the original range-time intensity of FAIs is present in the left panel and 1D range-imaging result is displayed in the middle panel. In the range imaging, the radar echoes collected from the full array and five carrier frequencies were processed with the standard Capon method, and the range step was 5 m. It is clearly seen in the range-imaged brightness distribution that these FAI echoes were positive range rate, indicating a northern component of motion. This direction of motion can be revealed from the typical range–Doppler spectrum of the echoes provided in the right panel of Figure 2a. As shown, the Doppler velocities spread out over a wide range, but were negative mostly and centered on −100 m/s, indicating a motion of away from the radar site according to our definition of Doppler velocity direction in calculation. Since the tilted radar beam was directed to geographic north, apparently a northern component of motion existed in the FAIs structure. Certainly, the zonal motion might also contribute a part to the Doppler spectrum, but could not be as large as 100 m/s because the radar beam was directed to geographic north and had a very narrow beamwidth (~3.6◦). Notice that some visible spectral power at the Doppler velocities of around 200 m/s and in the range interval between 164 and 166 km were aliases of some negative Doppler velocities that were smaller than the negative Nyquist velocity (about −215 m/s).

Figure 2b shows 2D angular brightness distributions of two examples retrieved from different methods. The value given in the title of the Fourier panel is SNR in linear scale, and following the SNR value is the data time on 26 July 2019 (hh:mm:ss LT). Obviously, the first example shows that the Fourier and Capon methods were not good for retrieving the brightness distribution (present in dB scale). The Fourier result had a coarser angular resolution and suffered from serious contamination of sidelobes arising from the weighting vector used in retrieval (refer to the Appendix A). The Capon result, on the other hand, provided an unreliable distribution, which is known to be related to the sensitivity of the method to small errors in the received signals. Nevertheless, the Capon method is still workable for many of the echoes, as shown in the example displayed in the second row of Figure 2b. By contrast, the NC-Capon method exported more stable and reliable outcomes, and moreover, a larger δ value yielded a higher angular resolution. A comparison of the results from different δ values is given later. Readers can also find the meaning of the δ value in Appendix A.

Figure 2c shows two examples of the angular brightness distributions in 3D imaging retrieved from different methods. The upper two rows provide two slices of brightness distribution at the range position of 150 m and 0 m in the sampling gate, and the bottom row presents the brightness distribution at the central position (0 m) of the sampling gate in another example. In general, the features resulting from different methods were similar to that in 2D imaging.

Figure 3 shows the scatter plots of brightness width (degree) vs. SNR (dB). The outcomes of Gaussian fitting and peak-find methods, 2D and 3D imaging, and Capon and NC-Capon methods are present, where the δ value in the NC-Capon method was 60. In this study, we define the half width of half power as representative of aspect angle, and therefore the brightness width obtained from Gaussian fitting was multiplied by <sup>√</sup>*ln*4. This correction was also made for the results of the RI technique shown later. The brightness centers with angular locations within 2◦ of the radar beam direction were adopted, and moreover, the scatter plots of 3D imaging included only the brightness widths from the slices of angular distributions at the range positions of −75, −37.5, 0, 37.5, and 75 m in the sampling gate. The reason for this selection of outcomes is shown later in Figure 4. Several findings from Figure 3 are listed below:


Based on the consequences of what is shown in Figures 2 and 3, the NC-Capon method is definitely more suitable than the standard Capon method for retrieval of brightness width and thereby aspect angle. In the rest of this paper, therefore, only the results of the NC-Capon and peak-find methods are present. Readers can find additional maps of brightness distribution via the web link provided in [29].

As mentioned, in the scatter plot of 3D imaging, only the outcomes from the five slices of angular distributions at the range positions of −75, −37.5, 0, 37.5, and 75 m in the sampling gate were included. This was based on the features shown in Figure 4, where the profiles of mean brightness widths within five sampling gates are shown for three δ values: 20, 40, and 60. The left and right panels demonstrate the results from the outcomes with SNR > 1 and >20, respectively. The horizontal bar denotes the width of two times the standard deviation of brightness width. It can be seen clearly in the left panel that in each sampling gate, the mean value of brightness width varied with the range and had a smaller value around the central range of the sampling gate (indicated by g1 to g5). If a higher threshold of SNR was used, i.e., SNR > 20 as given in the right panel, the dependence of brightness width on the range in the sampling gate existed still although the dependence was mitigated significantly. One cause of this characteristic could be the volume weighting effect. In this study, range-weighting and beam-weighting effects were not corrected in either 2D or 3D imaging. This is under the consideration that the echoes of the FAIs are highly localized, and therefore correction of the weighting effect could produce overcorrected brightness values around the edge of the radar volume. To avoid the influence of the weighting effect on estimate of brightness width, we used the outcomes with the brightness centers located around the center of the sampling gate, that is, the range interval within −75 and 75 m, and the angular interval between −2◦ and 2◦.

**Figure 3.** Scatter plots of brightness width (aspect angle) vs. SNR, retrieved from Capon and NC-Capon methods: (**a**,**b**) show 2D radar imaging results obtained from Capon and NC-Capon methods, respectively. The parameter of δ is 60 in NC-Capon method. (**c**,**d**) are similar to (**a**,**b**) but show 3D radar imaging results. Both outcomes of Gaussian-fit and peak-find methods are present.

Figure 4 shows that mean brightness width varied with the δ value used in the NC-Capon method. In view of this, a dependence of the outcome on δ value should be examined further to find a suitable value of δ for the following data analysis. Figure 5a shows the results from the 2D NC-Capon method, where the upper two panels are the scatter plots of brightness width vs. SNR resulting from the δ values of 10 and 30. Obviously, the brightness width at larger SNR decreased with increased δ value. To reduce the influence of SNR, therefore, the outcomes with SNR > 10 were used to estimate mean brightness width, and the results are shown in the lower panel of Figure 5a. As seen, the change of mean brightness width was tiny when the δ value was larger than 40, although the standard deviation of brightness widths increased slightly with δ value. Similar behavior can also be found in the 3D NC-Capon results, as shown in Figure 5b.

More issues of the NC-Capon method that deserve inspection are computation time and brightness value at different δ values. Figure 6 provides an example of 3D NC-Capon showing the variations of brightness sum and computation time with δ value, where the values on ordinate are self-normalized. Note that the brightness sum in panel (a) was the summation of brightness values in the radar volume just from a single estimate; the other estimates had similar scenarios (not shown). On the other hand, the computation time shown in panel (b) was recorded from fifty estimates. As seen, the self-normalized brightness sum decreased

with increased δ value, and had tiny changes at the δ values larger than ~40. Meanwhile, the computation time also decreased quickly with increased δ value.

**Figure 4.** Profiles of brightness widths obtained from the NC-Capon method with three different δ values (20, 40, and 60). The brightness widths for statistical estimates are collected with the threshold of (**a**) SNR > 1 and (**b**) SNR > 20; g1 to g5 indicate the range center locations of five sampling gates.

Based on the investigation results shown in Figures 5 and 6, we can use the value larger than 40 for δ to save much computation time, but should not be too large to avoid magnifying the divergence of brightness widths. In view of this, the consequences of δ = 60 are presented in the remaining discussion of the 2D and 3D NC-Capon methods.

**Figure 5.** (**a**) Upper: brightness width (aspect angle) vs. SNR for 2D radar imaging, retrieved from the NC-Capon method with δ = 10 and 30, respectively. Bottom: brightness width vs. δ value; dots are mean values and short vertical segments denote two times the standard deviation of brightness widths. The brightness widths are selected from SNR > 10 and angular center location within 2◦. The value of standard deviation is indicated above the segment. (**b**) The same as (**a**), but for 3D radar imaging, in which the restriction of range interval between −75 and 75 m is added to select the outcomes for statistical estimates.

**Figure 6.** (**a**) Sum of brightness values vs. δ value, resulting from a single estimate. (**b**) Computation time of 50 estimates vs. δ value. Ordinate values are self-normalized. The computation is performed with NC-Capon method for 3D radar imaging.

#### *3.2. SNR Effect*

As seen in the scatter plots of brightness width vs. SNR shown in Figure 5, a lower SNR led to a larger brightness width. In many previous studies, finite correction of SNR effect for the coherence of the covariance function had been suggested to improve the consequences. Usually, this can be made through multiplying the coherence by (1 + 1/SNR1) 1/2(1 + 1/SNR2) 1/2, in which SNR1 and SNR2 are the respective SNRs of the two data sets used in the analysis. Such a correction will increase the coherence. However, a smaller SNR magnifies the coherence greatly, and could lead to a coherence larger than 1 sometimes, which is irrational and may cause the coherence matrix to be more singular, or lead to

improper Gaussian fitting in the RI approach. In view of this, only the coherence with SNR > 1 was corrected, and moreover, the corrected coherence was assumed to be 0.99 when it was larger than 1.

Figure 7a,b provides the comparisons of original and SNR-corrected results of 2D and 3D imaging, respectively. The original scatter plots of 2D and 3D imaging show a significant decrease of brightness width with SNR for the echoes of SNR < ~15 dB. After SNR-correction, the dependence of brightness width on SNR was mitigated appreciably in the 2D imaging result, as shown in Figure 7a. By contrast, SNR-corrected coherence yielded less alleviation of SNR-dependent features in the 3D imaging result, as shown in Figure 7b.

**Figure 7.** Scatter plots of brightness width vs. SNR. (**a**) Results of 2D radar imaging with normconstrained Capon method, δ = 60. The left and right subpanels show the outcomes from the original and SNR-corrected coherence functions, respectively. (**b**) The same as (**a**) but for 3D radar imaging.

Based on the consequences shown in Figure 7, it is suggested that the outcomes with SNR > 10 dB or even >15 dB can be adopted to calculate the brightness width of the echoes. On the condition of SNR > 10 dB, most of the aspect angles, represented by brightness widths, spread from ~0.15◦ to ~0.4◦. This order of magnitude is consistent with the RI measurements of aspect angles at the equator [3,6] and lower mid-latitude E region [5,14].

#### *3.3. Radar Interferometry*

As mentioned in Section 2, four sets of subarrays in the MU radar array, of which baselines are aligned approximately with the direction of the geomagnetic field line, are available for radar interferometry (RI). To compare with the results of radar imaging, however, expression (1) without the factor of Doppler shift was employed. Figure 8 exhibits the results, where the scatter plots of brightness width vs. SNR are shown. Several features can be found in Figure 8:


**Figure 8.** Scatter plots of brightness width vs. SNR resulting from the radar-interferometry approach proposed by [3]. Four sets of receiving subarrays are used, as indicated in the panels (**a**–**d**). Left and right subpanels of each panel show the outcomes from the original and SNR-corrected coherence functions, respectively.

For feature (3), it can be understood that (b) and (c) have more subarrays than (a), and the baselines of (b) and (c) are parallel to the direction of the geomagnetic field line more than that of (d). Accordingly, the subarrays in (b) and (c) are more suitable for retrieval of aspect angle of FAIs. As for the features (1) and (2), the aspect angles have more outliers and spread in a broader interval of angle as compared with 2D and 3D imaging results; one of the causes could be the certainty that the baselines of subarrays did not align with the geomagnetic field lines exactly, and another cause was the smaller length of baselines between subarrays. The maximum length of baseline was about 70 m in the subarrays used in Figure 8a–c, and about 80 m in Figure 8d, which were much smaller than that used in [3] and [6] (~310 m).

#### *3.4. Statistical Comparison of Aspect Angles*

Since different methods may result in different brightness widths (i.e., aspect angles), we compared their brightness widths, as shown in Figure 9. In panel (a), the results from 2D and 3D imaging with 19 subarrays and 7 subarrays (legend " ... rx19" and " ... rx7", respectively), and RI (legend "RI: ... ") are shown. The 7 subarrays used here are indicated by the dashed circles in Figure 1. The two sets of subarrays for RI are those used in Figure 8c,d. Correction of the SNR effect has been made for coherence and the threshold of 10 dB was given for SNR to select the outcomes for statistical calculation. The δ value was 60 in the NC-Capon method. As a result, the profile of mean brightness width denoted by RI:A2A4-D2 in panel (a) had smaller values than the others below the range of 162 km, while that denoted by RI:F3F4-C3 had the largest values above the range of 162 km. In addition, both mean profiles of RI were more range-variant than other profiles. On the other hand, the standard deviations of RI aspect angles were mostly larger than the others above 162 km, as observed in the right panel of Figure 9a. In the study, we also checked the RI results from the set of subarrays E3, E4, F5, B4, and B3, of which baseline alignment turns aside from the direction of the geomagnetic field line more than that of the subarrays present here, and found that the resultant aspect angles spread over a much broader interval of angle (not shown). It is thus demonstrated that the aspect angles resulting from RI are highly dependent on the baseline alignment of the subarrays. In contrast to the RI results, the mean profiles of 2D imaging, denoted by "2Dd60rx19" and

"2Dd60rx7", had the smallest variation through the range, and moreover, the standard deviations of brightness widths of 2D imaging were generally the smallest among the results. There was not much difference between the mean profiles of "2Dd60rx19" and "2Dd60rx7", but the standard deviation marked by "2Dd60rx7" had values slightly larger than that of "2Dd60rx19", suggesting that a smaller number of subarrays may deteriorate the outcomes of 2D imaging. On the other hand, 3D imaging provided consequences approximately between 2D imaging and RI results. In general, the outcomes of 3D imaging resulting from 19 subarrays were closer to that of 2D imaging.

**Figure 9.** Profiles of mean and standard deviation of brightness widths estimated from various approaches. (**a**) Results from 2D and 3D radar imaging with the norm-constrained Capon method (δ = 60), and from the radar interferometry (RI) used by Kudeki and Farley [3]. In radar imaging, the locations of echo centers within 2◦ and between −75 and 75 m in the sampling gate are adopted for statistical calculation. In RI, the outcomes shown in Figure 8c,d are used. Only the data with SNR larger than 10 are taken. (**b**) Results from the RI approach proposed by Wang et al. [14].

Another outcome of aspect angle estimated from the definition of Wang et al. [14] is provided in Figure 9b. The method of [14] was also based on RI but three non-coplanar receiving antennas/subarrays were used to derive a mean echo center in the radar volume. These locations of echo centers obtained were generally geomagnetic field-aligned. The aspect angle was then defined to be half the angle subtended by the field-aligned length of echo center locations. There are many associations of three receiving subarrays available in the MU radar array. Among them, the results from the three sets of subarrays—A2E2C2, A2E4B4, and F5D4C4—are shown, which are typical and represent larger, medium, and smaller triangle configurations of receiving subarrays, respectively. Note that the abscissa scale in panel (b) is twice that in panel (a). Obviously, the estimated aspect angles depend on the baseline length between receiving antennas, in which a larger triangle configuration yielded smaller and more concentrated values of aspect angles (solid curves). The smallest triangle configuration (F5D4C4) provided the largest mean and most divergent outcomes (circles). In quantitative analysis, the mean aspect angles of the medium triangle configuration (A2E4B4) were around 0.2◦, which were closer to those of 2D radar imaging. The largest triangle configuration (A2E2C2), however, produced mean aspect angles less than 0.2◦ and close to 0.1◦ at some heights; these were the smallest values among all results. Nevertheless, the deviations of aspect angles yielded by the largest and medium triangle configurations mostly ranged between 0.05◦ and 0.2◦, which were comparable to those in panel (a) but still lager than the results of 2D radar imaging (<0.05◦).

Based on the comparison made in Figure 9, 2D radar imaging can yield the brightness width/aspect angle of FAIs varying with range more smoothly and converging on smaller interval. It should be known, however, that the outcome from 2D imaging is an average throughout the range of the sampling gate (600 m in the present study); 3D imaging, on the other hand, is likely to be capable of inspecting the variation of aspect angle or irregularity structure within the sampling gate. A problem for 3D imaging is that its computation time is much longer than that of 2D imaging, especially when using the NC-Capon method. A test shows that the computation time of 3D imaging can be about 150 times that of 2D imaging. The iteration number in the computation of the NC-Capon method is unpredictable because it depends on the characteristics of the individual data set. In addition, the more subarrays employed and more slice numbers given in the processing, the more time consumed in computation, although using more subarrays and slice numbers for radar imaging may export a higher resolution of brightness distribution.

#### **4. Conclusions**

In this study, radar imaging and radar interferometry (RI) techniques were applied to the measurement of the aspect angle of field-aligned plasma irregularities (FAIs) in the mid-latitude ionospheric E region. This was achieved by means of multireceiver and multifrequency capabilities built in the Middle and Upper (MU) Atmosphere radar (34.85◦N and 136.10◦E) of Japan. The brightness distribution retrieved from 2D and 3D radar imaging with three retrieval methods, i.e., Fourier, Capon, and norm-constrained Capon, were investigated. It was shown that the norm-constrained Capon method could produce more reliable brightness distribution for most circumstances and thereby yield a more trustworthy value of aspect angle.

By taking the radar data of SNR > 10 dB for 2D and 3D radar imaging, the estimated mean aspect angles ranged from 0.2◦ to 0.3◦, and the standard deviations of aspect angles were less than 0.1◦ and could be as small as 0.01◦ at some heights for 2D radar imaging. These values are consistent with those measured in the equatorial electrojet (EEJ) region and in the lower mid-latitude E region with RI technique. With the MU radar, the RI results are also acceptable for certain receiving configurations, but the derived aspect angles spread over an angular range wider than that of radar imaging. This is attributed to a much shorter baseline used for RI in the MU radar array, as compared with the radar used for the EEJ region, and/or different approaches/definitions of aspect angle made for the lower mid-latitude E region.

In conclusion, 2D radar imaging with the norm-constrained Capon method is preferable to other methods, examined in this study, for measurement of aspect angle with the MU radar configuration. Three-dimensional radar imaging is also recommended for further studies of variation in aspect angle and irregularity structure in the radar volume. Nevertheless, 3D radar imaging with the norm-constrained Capon method consumes a much longer time in computation, which is not advantageous to a great amount of data analysis for the ionosphere. In the future, variations of aspect angle with time and altitude can be investigated with 2D radar imaging for the mid-latitude ionosphere. The results can assist in understanding the spatial and temporal characteristics of plasma irregularities in the ionosphere within the radar volume and at the temporal resolution of instrument integration time.

**Author Contributions:** Conceptualization, J.-S.C. and Y.-H.C.; methodology, J.-S.C., C.-Y.W., and Y.-H.C.; software, J.-S.C. and C.-Y.W.; validation, J.-S.C., C.-Y.W., and Y.-H.C.; formal analysis, J.-S.C. and C.-Y.W.; investigation, J.-S.C. and C.-Y.W.; data curation, J.-S.C. and C.-Y.W.; writing—original

draft preparation, J.-S.C.; writing—review and editing, J.-S.C., C.-Y.W., and Y.-H.C.; visualization, J.- S.C. and C.-Y.W.; supervision, J.-S.C. and Y.-H.C.; project administration, J.-S.C. and Y.-H.C.; funding acquisition, J.-S.C. and Y.-H.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Ministry of Science and Technology, and China Medical University, Taiwan (ROC), under grants MOST110-2111-M-039-001 and CMU110-MF-69, respectively.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data sharing is not applicable to this article.

**Acknowledgments:** Experiments were conducted with the support of the MUR International Collaborative Research Program, Research Institute for Sustainable Humanosphere, Kyoto University, Japan.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Radar Imaging**

The equations of radar imaging using multiple receivers and multiple frequencies are provided briefly below. Given the signals from N receivers and M carrier frequencies,

$$\begin{array}{lcl} \mathbf{S}(t) = & \begin{bmatrix} S\_{11}(t)S\_{12}(t) \dots S\_{1N}(t) \\ S\_{21}(t)S\_{22}(t) \dots S\_{2N}(t) \end{bmatrix} \\ \vdots \\ \mathbf{S}\_{M1}(t)S\_{M2}(t) \dots S\_{MN}(t) \end{array} \tag{A1}$$

where *t* is time, *Sij*(*t*) is the signal received by receiver *j* at carrier frequency *i*, and the superscript *T* denotes the transpose of matrix.

First, covariance function of **R**(*t*) is estimated:

$$\mathbf{R}(t) = \mathbf{S}(t)\,\mathbf{S}^{H}(t) = \begin{bmatrix} \mathbf{V}\_{11}\mathbf{V}\_{12}\cdots\mathbf{V}\_{1M} \\ \mathbf{V}\_{21}\mathbf{V}\_{22}\cdots\mathbf{V}\_{2M} \\ \vdots\vdots\ddots\vdots \\ \mathbf{V}\_{M1}\mathbf{V}\_{N2}\cdots\mathbf{V}\_{MM} \end{bmatrix} \,\tag{A2}$$

where the superscript *H* denotes the conjugate and transpose (Hermitian) operator, and **V***pq* is an *N* × *N* matrix composed of covariance functions of all receiver pairs at the carrier frequencies *p* and *q*. As a result, **R**(*t*) is an *MN-by-MN* matrix.

Second, the so-called brightness (power density), *B*, can be estimated with the following expression:

$$B(\mathfrak{a},r) = \frac{1}{\mathcal{W}^2(\mathfrak{a},r)} \mathbf{w}^H \mathbf{R} \mathbf{w}\_{\prime} \tag{A3}$$

where *r* is range, and *a* = [sin*θ*sin*ϕ*, sin*θ*cos*ϕ*, cos*θ*], *θ* and *ϕ* are zenithal and azimuthal angles, respectively. *W*(*a*, *r*) is the spatial weighting function formed by radar beam and range weighting functions. The weighting vector **w** is given as

$$\mathbf{w} = \begin{bmatrix} e^{-j(2k\_1r - k\_1\mathbf{a}\bullet\mathbf{D}\_1)} \ e^{-j(2k\_1r - k\_1\mathbf{a}\bullet\mathbf{D}\_2)} & \dots \ e^{-j(2k\_1r - k\_1\mathbf{a}\bullet\mathbf{D}\_N)} \\ e^{-j(2k\_2r - k\_2\mathbf{a}\bullet\mathbf{D}\_1)} \ e^{-j(2k\_2r - k\_2\mathbf{a}\bullet\mathbf{D}\_2)} & \dots \ e^{-j(2k\_2r - k\_2\mathbf{a}\bullet\mathbf{D}\_N)} \end{bmatrix} \tag{A4}$$
 
$$\begin{bmatrix} \vdots & \vdots & \ddots & \vdots & \vdots \\ e^{-j(2k\_Mr - k\_M\mathbf{a}\bullet\mathbf{D}\_1)} \ e^{-j(2k\_Mr - k\_M\mathbf{a}\bullet\mathbf{D}\_2)} & \dots \ e^{-j(2k\_Mr - k\_M\mathbf{a}\bullet\mathbf{D}\_N)} \end{bmatrix}^T$$

where *D*<sup>j</sup> is the location coordinates of receiver *j*, and *ki* is the wave number of carrier frequency *i*. Expression (A3) is a linear beamforming, called the Fourier method in this study. The weighting pattern of (A4) contains appreciable sidelobes that will contaminate the imaging result of (A3); cases of 2D imaging are shown in Figure 2b,c. On the other hand, the estimator of the standard Capon method is:

$$B(\mathbf{a},r) = \frac{1}{\mathcal{W}^2(\mathbf{a},r)} \frac{1}{\mathbf{e}^H \mathbf{R}^{-1} \mathbf{e}'} \tag{A5}$$

where the superscript −1 is the inverse of matrix, and vector **e** is identical to (A4). In obtaining (A5), a constrained condition, min(B = **w***H***Rw**) subject to **e***H***w** = 1, is given to derive an optimal weighting vector:

$$\mathbf{w}\_{\mathbb{C}} = \frac{\mathbf{R}^{-1}\mathbf{e}}{\mathbf{e}^{H}\mathbf{R}^{-1}\mathbf{e}}.\tag{A6}$$

Substituting (A6) into (A3) yields (A5). The estimator (A5) suppresses the contamination of the sidelobes in the weighting pattern of (A4); however, it is sensitive to small errors in the received signals, which may cause the covariance matrix R to be singular or badly scaled, leading to an inaccurate estimate of R<sup>−</sup>1.

To mitigate the sensitivity of (A5) to small errors in the received signals, the normconstrained Capon method uses two constraints:

$$\min(\mathbf{B} = \mathbf{w}^H \mathbf{R} \mathbf{w}) \text{ subject to } \mathbf{e}^H \mathbf{w} = \mathbf{N} \times \mathbf{M}, \tag{A7}$$

$$\delta \mid \mathbf{w}^H \mathbf{w} \mid \le \delta \text{ (N} \times M\text{)},\tag{A8}$$

where *δ* is a user-defined value. Expression (A7) is basically the same constraint used in derivation of (A6) except for the multiplier *N* × *M*. To meet the demand of norm constraint (A8), a diagonal-loading value *σ* is applied to (A6):

$$\mathbf{w}\_{\rm NC} = \frac{\left(\mathbf{R} + \sigma \mathbf{I}\right)^{-1} \mathbf{e}}{\mathbf{e}^{H} \left(\mathbf{R} + \sigma \mathbf{I}\right)^{-1} \mathbf{e}} \text{ ( $N \times M$ )},\tag{A9}$$

where **I** is the identity matrix. A suitable **w**NC can be obtained by giving a proper value of *δ* and then increasing the value of *σ* from 0 until (A9) satisfies the constraint (A8). The brightness is thus estimated with the following expression:

$$B(\mathbf{a}, r) = \mathbf{w}\_{\rm NC}^H (\mathbf{R} + \sigma \mathbf{I}) \mathbf{w}\_{\rm NC}. \tag{A10}$$

Note that selection of a proper *δ* value is another issue of the NC-Capon method; it could depend on the data and radar system characteristics. The effect of *δ* value on the outcomes of brightness width (aspect angle) were examined and shown in Figures 5 and 6.

In this study, the spatial weighting function *W*(*a*, *r*) was ignored because of highly localized FAI echoes. Moreover, in use of these equations for 1D range imaging, only one receiver is needed and so the receiver-indexed terms and parameters are either omitted or reduced to unit value. For 2D angular imaging, only one carrier frequency is required and thus the frequency-indexed terms and parameters are either left out or replaced by one. Notice that the zenith orientation is toward the bearing of the radar beam because of a tilt beam used, and then the range goes along the radar beam direction and the angular surface is the transverse of the radar beam direction.

One more note is that in multiple frequency operation, the *M* carrier frequencies are used with sequential radar pulses and repeat in observation, therefore the time difference between the echoes from the first and *M*-th pulses is *M* times that of the interpulse period (IPP). If the target moves very fast, such as meteor heads, a significant change of location occurs during the time period of *M* × IPP so that range imaging of the target could deteriorate or fail.

#### **References**


## *Article* **Magnetic Field and Electron Density Scaling Properties in the Equatorial Plasma Bubbles**

**Paola De Michelis 1,\*, Giuseppe Consolini 2, Tommaso Alberti 2, Roberta Tozzi 1, Fabio Giannattasio 1, Igino Coco 1, Michael Pezzopane <sup>1</sup> and Alessio Pignalberi <sup>1</sup>**


**Abstract:** The ionospheric plasma density irregularities are known to play a role in the propagation of electromagnetic signals and to be one of the most important sources of disturbance for the Global Navigation Satellite System, being responsible for degradation and, sometimes, interruptions of the signals received by the system. In the equatorial ionospheric F region, these plasma density irregularities, known as plasma bubbles, find the suitable conditions for their development during post-sunset hours. In recent years, important features of plasma bubbles such as their dependence on latitude, longitude, and solar and geomagnetic activities have been inferred indirectly using their magnetic signatures. Here, we study the scaling properties of both the electron density and the magnetic field inside the plasma bubbles using measurements on board the Swarm A satellite from 1 April 2014 to 31 January 2016. We show that the spectral features of plasma irregularities cannot be directly inferred from their magnetic signatures. A relation more complex than the linear one is necessary to properly describe the role played by the evolution of plasma bubbles with local time and by the development of turbulent phenomena.

**Keywords:** plasma turbulence; ionospheric irregularities

#### **1. Introduction**

One of the most interesting features of the equatorial ionospheric F region is the existence of plasma density irregularities, which find the suitable conditions for their development during post-sunset hours [1]. In the last few years, different terminologies have been used to term the equatorial ionospheric irregularities. Often, these are identified as equatorial spread-F, in virtue of the spread of the experimental trace either in frequency or in amplitude affecting the ionosonde observations, which is the marking of the existence of electron density irregularities in the reflecting layer (e.g., Booker and Wells [2]). With the advent of in situ satellite observations, they started using terms such as "depletions", "bite-outs," or "plasma holes", because the irregularities were recorded mostly at locations where the electron density is depleted with respect to the background ionosphere (e.g., McClure et al. [3], Dyson and Benson [4]). Then, terms such as "plumes" or "wedges" appeared to characterize the morphology of turbulent regions and the development of such irregularities, as observed by radars (e.g., Scannapieco and Ossakow [5], Woodman and La Hoz [6]). In the last few decades, to identify the generation process of these irregularities, the term "bubble" has taken hold more and more (e.g., Burke [7]). It is worth noting that recent studies (e.g., Kil et al. [8]) introduced also the expression "plasma depletion shell" to highlight the three-dimensional structure of bubbles. The formation of these irregularities is recognized to be driven by the Rayleigh–Taylor instability mechanism [6,9,10] that generates in the ionosphere a situation similar to that occurring when a heavy fluid flows over

**Citation:** De Michelis, P.; Consolini, G.; Alberti, T.; Tozzi, R.; Giannattasio, F.; Coco, I.; Pezzopane, M.; Pignalberi, A. Magnetic Field and Electron Density Scaling Properties in the Equatorial Plasma Bubbles. *Remote Sens.* **2022**, *14*, 918. https://doi.org/ 10.3390/rs14040918

Academic Editor: Michael E. Gorbunov

Received: 8 January 2022 Accepted: 10 February 2022 Published: 14 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

a lighter one. The absence of sunlight, which leads to a different rate of recombination in the ionospheric layers, causes a density gradient between the upper and lower ionospheric F region with a plasma density in the upper F region higher than in the lower one. Such steep vertical density gradients realize the conditions for the growth of the Rayleigh–Taylor instability. These plasma density irregularities expand vertically and successively move along the magnetic field lines on each side of the geomagnetic equator. Typically, they also drifted eastward from vertical polarization electric fields associated with zonal neutral winds in the F region.

The result is the formation of magnetic field-aligned plasma-depleted regions, which are usually known as equatorial plasma bubbles. They are unstable to plasma density perturbations and polarization electric fields so that secondary irregularities can be generated, which are mainly the result of cascading processes. Nowadays, many important characteristics of these irregularities are known. They have different scale sizes, from hundreds of kilometers down to a few decameters [11,12], are mainly observed within a narrow band of ±20◦ magnetic latitude [9], are observable in a wide range of altitudes [6,9] up to around 1000/1500 km [13–15], and their spatial and temporal distributions are strongly controlled by various geophysical parameters such as solar [16] and geomagnetic activity [17], season [16,18], local time, and longitude [8,19–21]. What drives the scientific community to investigate the plasma density irregularities is the knowledge of their role in the propagation of electromagnetic signals. They can be responsible for the disturbance of radio waves and are one of the most important sources of disturbance for the Global Navigation Satellite System (GNSS) signal. Indeed, in the worst case, plasma density irregularities can cause a loss of lock event, a condition for which a GNSS receiver can no longer track the signal sent by the satellite, with a consequent degradation of the positioning accuracy. In the last few years, the relevance of the GNSS in our society has substantially increased, as many critical infrastructures and economies are dependent on the positioning, navigation, and timing services of GNSS, so that our society is nowadays vulnerable to damages due to the malfunction of these systems. For this reason, one of the research priorities in the space weather community is to improve the knowledge of these plasma density irregularities, trying to better understand their features and generation mechanisms.

An interesting aspect of these irregularities, which requires more study, is their possible turbulent nature. Indeed, the equatorial ionosphere has proper physical conditions for the establishment of a fully developed turbulence [22,23]. Through numerical simulations, it has been possible to reproduce the ionospheric irregularities and study the corresponding dynamics under collisional and inertial regime flows [24,25]. It has been found that the plasma irregularities may be characterized by energy spectra with slopes ranging between −5/3 and −3 associated with the occurrence of energy and/or enstrophy cascades, respectively [26,27]. Over the years, these numerical models have been improved trying to describe the ionospheric two-dimensional plasma turbulence in a more realistic way. It has been found that after sunset, at apex heights between 600 and 1000 km, large-scale irregularities present a dynamics typical of inertial regime flow [28] and that although the foundations of the Kraichnan's theory could be undetermined in the topside equatorial ionospheric F region, the numerical simulation results seem to confirm the probable existence of large plasma density irregularities. These are characterized by velocity fluctuation spectra following power laws with spectral indices similar to those predicted for energy and enstrophy cascades in two-dimensional turbulent flows [23]. In the last twenty years, thanks essentially to the rapid increase in the performance of computers, it has been possible to improve the models describing the evolution of plasma depletions. Thus, new models capable of reproducing the three-dimensional turbulent structures of equatorial plasma bubbles with high spatial and temporal resolutions have been developed (see for example Yokoyama et al. [29], Yokoyama [30]). The results of such simulations confirm observations of plasma density data recorded, for example, on board DE-2 satellite. In this case, it has been found that plasma density is characterized by an energy spectrum with a slope around −5/3 for scale sizes larger than 1 km in all regions characterized

by strong plasma density irregularities [28]. Nowadays, the numerical simulations have become useful tools both to study the plasma bubbles and their nonlinear evolution, and to understand some of their features, which cannot be fully understood from theoretical predictions. Therefore, it is extremely important to combine the numerical simulation results with real observations.

In order to detect and characterize plasma density irregularities, different methods, approaches, instruments, and data can be used such as, for example, ionosondes, airglow imager, coherent/incoherent scatter radars [31], as well as Global Positioning System (GPS) Total Electron Content (TEC) [32,33] in the topside ionosphere. The ground-based observations are mostly fixed in a single location and cannot give the spatial variations of plasma density irregularities in a large area. Conversely, observations from Low Earth Orbit (LEO) satellites such as for example DMSP (Defense Meteorological Satellite Program), ROCSAT-1 (Republic of China Satellite), C/NOFS (Communications/Navigation Outage Forecasting System), and COSMIC (Constellation Observing System for Meteorology, Ionosphere, and Climate) permitted climatology studies and provided the latitude-local time distributions of plasma density irregularities. The recently launched Global-scale Observation of Limb and Disk (GOLD) mission, which is located at a geo-stationary orbit, seems to be a promising mission to study bubbles, offering the opportunity to observe the Earth's complete disk continuously from the geostationary orbit [34,35]. A different method for the detection of ionospheric irregularities based on the diamagnetic effect was suggested a few years ago. Plasma density irregularities are indeed associated with magnetic field perturbations [36,37] that can be used as a proxy for plasma irregularities. There are different ways to describe the so-called diamagnetic effect but, perhaps, the easiest one is to consider that the sum of magnetic and plasma pressures has to be constant in a quasi-stationary state. This means that when the plasma pressure is reduced, the magnetic pressure and consequently the magnetic field strength must increase. The equatorial plasma bubbles identify regions of reduced plasma density; consequently, they have to be characterized by a magnetic field strength higher than the ambient one. That was confirmed using simultaneous observations of magnetic field and electron density recorded by the CHAllenging Minisatellite Payload (CHAMP) satellite. These measurements revealed the increase of a few nT in the main field intensity inside the plasma density depletions [38] Park et al. [39] quantitatively investigated the actual balance between plasma and magnetic pressures within the plasma bubbles, showing that a dominant part of the magnetic pressure change could be explained in term of plasma density change under specific assumptions. Based on the high correlation between the variations of plasma density and those associated with magnetic field strength, an extensive survey of magnetic signatures related to plasma density irregularities was used for statistical studies [15,40,41] devoted to the analysis of the global distribution of equatorial plasma irregularities in the topside ionosphere. These studies revealed features of the magnetic signatures that closely reflected those of the plasma bubbles previously obtained using different methods as GPS scintillation [19], in situ plasma density measurements [42,43], and radio wave propagation [44]. Furthermore, magnetic field measurements recorded on board the CHAMP satellite with different time resolution (50 Hz and 1 Hz) permitted addressing some properties of the plasma density irregularities related to their different spatial sizes. Evidently, these results do not reflect directly the features of plasma depletions but their associated magnetic signatures. However, the features obtained using the magnetic signatures of the plasma depletions are consistent with those obtained by directly using the electron density measurements. For this reason, the diamagnetic effect, as an indirect way of sampling plasma density structures, can be extremely useful in studying plasma irregularities.

While it is correct to assume that some plasma bubbles features can be inferred using the corresponding magnetic signatures, the assumption that plasma bubbles scaling properties can be inferred using the associated magnetic field scaling properties may be not necessarily correct as well. The purpose of this study is to understand if the hypothesis of a linear relation between the scaling properties of the electron density and of the magnetic field is valid. So, needless to say, the clarification of this point is of crucial importance.

#### **2. Data Description**

To investigate the scaling features of the electron density irregularities and of their magnetic signatures at low latitude, we use measurements on board the Swarm A satellite during its near-polar orbit around the Earth at an altitude of approximately 460 km [45]. Swarm A is one of the satellites constellation launched by ESA for Earth observation. Specifically, the ESA mission consists of three identical satellites named Alpha, Bravo, and Charlie (A, B, and C), which were launched on 22 November 2013 into a near-polar orbit. Swarm A and C fly side-by-side at an altitude of 462 km (initial altitude) and at an inclination angle of 87.35◦, whereas Swarm B flies at a higher orbit of 511 km (initial altitude) and at an inclination angle of 87.75◦.

We consider measurements collected from 1 April 2014 to 31 January 2016 at low and mid latitudes (|*Lat*| < 40◦) under conditions of high solar activity. It is known indeed that the equatorial plasma bubbles occur mainly during periods characterized by high extreme ultraviolet solar flux [41]. In the selected time interval, the value of the solar radio flux index (F10.7) is equal to (100 ± 30) sfu. This is not an extremely high value for the solar activity, but the selected period, which covers a part of the descending phase of Solar Cycle 24, is the only period of moderate/high solar activity during which Swarm is flying. Electron density and magnetic field vector data at a rate of 1 Hz are selected from the ESA dissemination server (http://earth.esa.int/swarm accessed on 12 February 2022). To highlight the magnetic signatures associated with plasma bubbles, we remove from the magnetic field measurements the contribution due to the magnetic field of internal origin, as modeled by CHAOS-6 [46]. CHAOS-6 is the latest generation of the CHAOS series of global high-resolution geomagnetic field models, which spans a period between 1997 and 2016. It is derived primarily from magnetic satellite data (Swarm, CHAMP, Ørsted, and SAC-C), although ground-based activity indices and observatory monthly means are also used. It is able to model the dominating core field, the crustal field, and the magnetic fields from large-scale magnetospheric currents, but there is no explicit representation of the ionospheric field or fields due to magnetosphere—ionosphere coupling currents. One of the applications of this model is the removal from the observed magnetic field of the large contributions from the prime sources when the weaker effects such as ionospheric currents, plasma irregularities, or even ocean tidal motions want to be investigated. We use the CHAOS-6 package to compute the geomagnetic field of internal origin, which is successively subtracted to the magnetic field observed by Swarm. Lastly, we focus on the 18:00–24:00 local time (LT) sector and select periods characterized by low/moderate geomagnetic activity. To select these periods, we used the *Kp* index, which is designed to measure the solar particle radiation through its magnetic effects, and today, it is considered a proxy for the energy input from the solar wind to Earth. It is a global geomagnetic activity index based on 3-hour measurements from ground-based magnetometers around the world. The *Kp* index ranges from 0 to 9, where a value of 0 means that the geomagnetic activity is absent, and a value of 9 means that an extreme geomagnetic storming is occurring. Data are distributed by GFZ (Helmholtz Centre Potsdam) [47] and they are redistributed by various data centers and databases. To select periods characterized by low/moderate geomagnetic activity, we consider all those time intervals characterized by *Kp* ≤ 3. According to previous results [40], the probability of observing plasma bubbles is maximum under these conditions.

To select the plasma bubbles events, we use the Swarm Level-2 Ionospheric Bubble Index (IBI) product, which was evaluated using plasma density and magnetic field data with a cadence of 1 Hz. A full description of the data file and explanations regarding the IBI data can be found in the L2-IBI product description file (https://earth.esa.int/eogateway/ documents/20142/37627/Swarm-Level-2-IBI-product-description.pdf/3e9f6c3a-1ffc-ea5 3-0161-a18b63f90c6f, accessed on 13 February 2022). This product is composed of three

parameters: Bubble Index, Bubble Probability, and Bubble Flag. The Bubble Index can assume three different values: 1, when the analyzed data are affected by a plasma bubble, 0 when there is not a plasma bubble, and −1 when conditions do not allow unequivocally identifying a plasma bubble in the data [48]. Furthermore, when IBI = 1, there is an additional flag, the Bubble Flag (BF), indicating the quality of plasma bubbles detection that is related to exceeding a certain probability threshold. When BF = 1, we are in the presence of a high correlation between plasma density depletion and magnetic field signature. Here, we chose the condition IBI = 1 and BF = 1 to identify the plasma bubbles. However, we note that this automatic detection algorithm cannot distinguish equatorial plasma bubbles from plasma blobs, which are localized plasma density enhancements. This means that the events identified by the L2-IBI processor can contain also blob events. However, the occurrence probability of equatorial plasma bubbles is significantly higher than that of blobs mainly during solar maximum years.

#### **3. Method and Data Analysis**

The scaling features of the electron density and magnetic field at different temporal scales are evaluated using the so-called generalized *qth*-order structure function, namely *Sq*. For a signal *X*(*t*), *Sq* is given by the following equation:

$$S\_q(\tau) = \langle | \left| X(t+\tau) - X(t) \right| \left| \prime \right\rangle,\tag{1}$$

where *t* is the time, *τ* is the time delay, and ... stands for a statistical average. When *Sq*(*τ*) has a power-law behavior as a function of *τ*:

$$S\_q(\tau) = \tau^{\gamma(q)},\tag{2}$$

the analyzed signal is characterized by a scale invariance, and *γ*(*q*) is the so-called *qth* order scaling exponent, which defines the scaling nature of the increments of the signal *X*(*t*). A signal possesses a global scale-invariance if *γ*(*q*) = *H q*; that is, it is a linear function of the order *q*. Conversely, a signal is defined as multifractal, that is, it has a local scale-invariance, if *γ*(*q*) is a nonlinear convex function of *q*. Here, we investigate the second-order scaling exponent, *γ*(2), which is related to the power spectral density (PSD) exponent, *<sup>β</sup>* (PSD *<sup>f</sup>* <sup>−</sup>*β*), via the following relation:

$$
\beta = \gamma(2) + 1.\tag{3}
$$

The use of the structure function analysis, rather than the standard Fourier analysis, as a method to detect scaling features is motivated by the fact that the analysis of scaling/spectral features in the real time space, using signal increments, is less affected by the choice of a specific decomposition basis, as it is instead the case for the Fourier analysis. Furthermore, being geophysical signals, such as magnetic field and electron density fluctuations in the ionosphere, generally characterized by power-law spectra, i.e., *PSD*(*f*) *<sup>f</sup>* <sup>−</sup>*β*, with spectral exponents *β* ranging between 1 and 3, one deals with nonstationary signals with stationary increments [49]. In this case, the Fourier analysis, which requires stationarity, might produce unreliable results especially when evaluating spectral properties over short time intervals. Conversely, the structure function analysis works with signal increments that are expected to be nearly stationary, and thus, the corresponding results over short time windows/intervals are more reliable. Furthermore, the structure function analysis is a companion method that is widely used in the framework of turbulence-related studies.

The measurements collected by instruments on board a satellite require some precautions to be correctly used. The satellite crosses regions characterized by different physical processes, which change in space and time. Consequently, the collected measurements, such as in our case electron density and magnetic field strength, present properties with a local character. The scaling properties of these real-time series also acquire a local character, and it is necessary to apply methods of analysis capable of evaluating their local scaling

features describing the different regions explored by the satellite. The structure function approach allows doing this. Among the different methods to evaluate the structure functions, here, we apply the local detrended structure function analysis (DSFA) method discussed in De Michelis et al. [50,51] and already successfully applied in a series of previous works showing results that are well in agreement with the present literature on the scaling and/or spectral features of turbulent fluctuations in the ionosphere. In detail, we apply the analysis on a moving window of *N* = 301 points using a time series obtained from the original one after removing its linear trend, if present. This gives us the opportunity to remove from the initial signal all the possible large-scale variations that can affect a correct estimation of the scaling features. The choice of a moving window of 301 points is the result of an accurate series of test over a wide set of simulated multifractional Brownian motions [52–54], which are characterized by a local dependence (nonstationarity) of the scaling/spectral features. Indeed, this number of points is the right compromise to get an estimation of local scaling exponents sufficiently precise (with a typical error ≤ 7%), being 1 order of magnitude larger than the longest investigated timescale (*τ* = 30 s). In each moving window, we evaluate the second-order structure functions associated with the electron density and magnetic field strength time series, respectively. These are evaluated for time delay values (*τ*) ranging from 1 to 40 s, and the estimated second-order scaling exponents are associated with the position of the satellite at the center of the fixed window. We assume that the results obtained in the time domain are valid in the space domain too (see e.g., [55,56] and references therein). Thus, considering that the Swarm A orbital velocity is around 8 km/s, the results obtained in the time domain between 1 and 30 s are valid in the spatial domain between ∼8 and ∼ 250 km, i.e., *δr Vsτ* where *Vs* is the satellite velocity. We remark that the previous conversion between the timescale and the corresponding spatial scale is based on the assumption of the Taylor's hypothesis according to which instantaneous spatial averages coincide with temporal averages calculated from the recorded signal. Given this hypothesis as valid, which can be applied for frozen turbulence, we assume that the temporal evolution of turbulent structures is longer than that required for the satellite to cross them [22,55–57]. Hence, by evaluating the scaling properties of the magnetic field strength and electronic density associated with the plasma bubbles, we can investigate the properties of these plasma density irregularities from both a magnetic and plasma point of view in a range of scales between some kilometers and a few hundred kilometers.

The structure function analysis is initially done on the entire time series of around <sup>5</sup> · <sup>10</sup><sup>7</sup> points (22 months with a time resolution of 1 Hz) and only later, the scaling exponent values of the electron density (*Ne*) and magnetic field strength (|**B**|) associated with the occurrence of plasma bubbles are selected.

#### **4. Results**

Figure 1 reports a map of the values of the second-order structure function exponent, *γ*(2), which was obtained analyzing *Ne* and |**B**| inside the plasma bubbles. These two different parameters are characterized by similar values of *γ*(2), supporting the previous findings by Lühr et al. [15]. However, to look for any difference, we start by plotting in Figure 2 the corresponding histograms. The second-order structure function exponents range from 0.4 to 2 with the highest probability occurring at (1.0 ± 0.2) in both cases.

Taking into account the relation between the second-order scaling exponent and the Fourier power spectral density exponent (*β* = *γ*(2) + 1), we notice that the spectral features obtained in the case of electron density (*β* 2) are in agreement with previous results. For example, Livingston et al. [58] analyzing the spectra of electron density data recorded from the Atmospheric Explorer-E satellite (AE-E) in the interval 20:00–02:00 LT found a distribution of spectral slopes in the range between 1 and 3 centered around 1.9, which is a mean value that agrees with our results and is close to those of Dyson et al. [59], who used data recorded by the OGO satellite, Basu et al. [60], who analyzed AE-E measurements, and Le et al. [61], who analyzed ROCSAT-1 burst mode measurements at 1024 Hz sampling. Concerning the spectral scaling exponents obtained for the magnetic field strength, there

are not many previous analyses with which to compare them. Lühr et al. [15] used magnetic field data recorded by CHAMP to investigate the spectral properties of the signal and found values of the spectral indices ranging between 1.4 and 2.6 with a peak around 1.9. Again, these results agree with ours. Although the mean values of the two distributions reported in Figure 2 are equal, their shapes are different, in which the *γ*(2) distribution for *Ne* is skewed. Indeed, the skewness *<sup>λ</sup>*<sup>3</sup> values of *<sup>γ</sup>*(2) distributions of *Ne* and <sup>|</sup>**B**<sup>|</sup> are *<sup>λ</sup>Ne* <sup>3</sup> 0.45 and *λ*|**B**<sup>|</sup> <sup>3</sup> 0.08, respectively.

**Figure 1.** Values of the second-order scaling exponent, *γ*(2), associated with electron density (**top**) and magnetic field strength (**bottom**) inside the plasma bubbles in the geographic latitude–longitude plane. Data are selected according to the following conditions: 18:00 ≤ LT≤ 24:00, ±40◦ geographic latitude, *Kp* ≤ 3.

Figure 3 reports the values of the second-order structure function exponent obtained analyzing *Ne* and |**B**| inside the plasma bubbles according to QD magnetic latitude. For *Ne*, the highest occurrence of *γ*(2) is found for *γ*(2) < 1 at a latitudinal distance from the equator of around 10◦ in magnetic latitude in both hemispheres, while in the case of |**B**|, the highest occurrence of *γ*(2) is found for a value around 1 occurring always in the same regions. Thus, the peak of the highest occurrence of plasma bubbles is associated with a range of *γ*(2) values that is smaller for *Ne* than |**B**|.

**Figure 2.** Histograms of the second-order scaling exponent values associated with electron density (**left**) and magnetic field strength (**right**) inside the plasma bubbles, respectively. Data are selected according to the following conditions: 18:00 ≤ LT ≤ 24:00, ±40◦ geographic latitude, *Kp* ≤ 3.

To go into the observed differences between the two distributions, we repeat our analysis considering two-hour time intervals from 18:00 LT to 24:00 LT. Figure 4 reports the distribution of *γ*(2) values relative to *Ne* and |**B**| in the three different selected time intervals. A double peak distribution of *γ*(2) clearly appears analyzing the scaling properties of |**B**| in the range 20:00 ≤ LT < 22:00. The *γ*(2) distribution has one peak in the next local time interval (22:00 ≤ LT ≤ 24:00) and a rather flat shape in the previous one (18:00 ≤ LT < 20:00). In the case of *γ*(2) values associated with *Ne*, the data distribution has only one peak in the second and third local time intervals while a distribution with a rather flat shape appears in the first local time interval. The emergence of a bimodal distribution for the magnetic field *γ*(2) exponent could be an indication of a non-unique relation between the scaling features of electron density and magnetic field fluctuations in plasma bubbles. Furthermore, as clearly shown in several works (see, e.g., [62]) (and references therein), in the second LT interval here considered, plasma bubbles are still evolving at Swarm altitude, so that the bimodality could be a reflection of the nonstationarity of the plasma bubble structures.

**Figure 3.** Dependence on QD magnetic latitude of the second−order scaling exponent values associated with *Ne* (**left**) and |**B**| (**right**) inside the plasma bubbles.

**Figure 4.** Histograms of the second-order scaling exponent values associated with electron density (**left**) and magnetic field strength (**right**) inside the plasma bubbles, respectively. Data are selected according to the following conditions: 18:00 ≤ LT < 20:00 (light blue), 20:00 ≤ LT < 22:00 LT (blue) and 22:00 ≤ LT ≤ 24:00 (purple). In all cases, the geographic latitude is ±40◦ and *Kp* ≤ 3.

To get more insights on the bimodal character of the *γ*(2) distribution that appears, analyzing the scaling properties of |**B**| inside the plasma bubbles between 20:00 LT and 22:00 LT, we check how the *γ*(2) values of the two peaks are distributed in latitude and longitude. Taking into account that the two different populations are characterized by values of *γ*(2) < 1 and *γ*(2) > 1, respectively, Figure 5 shows their distributions in the latitude–longitude plane. The two geographic distributions seem similar, although the values *γ*(2) < 1 seem to prefer the areas closest to the equator. This feature is well highlighted on the bottom of Figure 5, where the *γ*(2) values are sorted according to magnetic latitude. Here, the red line clearly shows how the plasma bubbles around the magnetic equator (±4◦) are often characterized by a magnetic field strength with a secondorder scaling exponent *γ*(2) < 1. The values of the scaling exponents of the magnetic field strength seem to depend on magnetic latitude at least during the hours following their formation. To better investigate this point, we report in the left panels of Figure 6 the joint probability distributions of *<sup>γ</sup>*(2) values relative to *Ne* (*γNe* (2)) and |**B**| (*γ*|**B**|(2)) along with the mean values of *<sup>γ</sup>*|**B**|(2) exponent (white circles) for fixed *<sup>γ</sup>Ne* (2) at the three different selected two-hour time intervals. There is a large spreading of the observed distributions with most of the values concentrated around [*γNe* (2), *<sup>γ</sup>*|**B**|(2)] = [1, 1], which indicates that the relation between electron density and magnetic field scaling exponents can be assumed to be linear only in a zero-order approximation. However, looking at the trend of the mean values of *<sup>γ</sup>*|**B**|(2) fixed *<sup>γ</sup>Ne* (2), we can realize how the relation between the two scaling indices significantly departs from a linear trend for *γNe* (2) < 0.9 at LT ≥ 20 : 00. This result suggests that a more complex relation might exist between plasma density and magnetic field variations. Furthermore, the departure from a linear relation increases with the midnight approaching. We can conjecture that the observed departure from linearity occurs when the fluctuations of the electron density show scaling features similar to those expected for convective turbulence, i.e., *β* ∼ 5/3 [63]. In the right panels of Figure 6, we report the conditioned probably distribution density functions of *<sup>γ</sup>*|**B**|(2) values for fixed *<sup>γ</sup>Ne* (2) = 1. The observed distributions are multimodal, thus supporting the hypothesis that a linear relation between magnetic field intensity and plasma density spectral features has to be considered valid only on average being a very crude approximation. The multimodal structure of the observed distribution could also be the counterpart for different physical processes at the origin of the observed fluctuations, or it could be due to the role played by other physical quantities, such as for example the plasma temperature fluctuations, that could modify the relation between electron density and magnetic field second-order scaling

exponents. Moreover, a very interesting result stands in the clear bimodal character of the distribution of the magnetic field *γ*(2) exponent in the first LT sector. Indeed, there is not a one-to-one correspondence between the magnetic field *γ*(2) exponent and that of the electron density fluctuations, i.e., *<sup>γ</sup>*|**B**|(2) < *<sup>γ</sup>Ne* (2) ≡ 1. A possible explanation could be that, as shown in Gruzinov et al. [63], we are still in the early stage of a rising bubble so that the spectral features are still not stationary.

**Figure 5.** Values of the second-order scaling exponent associated with magnetic field strength inside the plasma bubbles in the geographic latitude–longitude plane during the time interval 20:00 ≤ LT < 22:00 LT for *γ*(2) < 1 and *γ*(2) > 1, respectively. On the bottom, the *γ*(2) < 1 and *γ*(2) > 1 values are sorted according to magnetic latitude.

**Figure 6.** On the left, the joint probability density functions of second-order scaling exponents for *Ne* and for |**B**| in the three different selected two-hour time intervals. White circles refer to the mean values of the magnetic field *<sup>γ</sup>*|**B**|(2) exponent at fixed electron density *<sup>γ</sup>Ne* (2) exponent values. The black line is the bisector line for which *<sup>γ</sup>*|**B**|(2) = *<sup>γ</sup>Ne* (2). On the right panels, the conditioned probability distribution density functions of *<sup>γ</sup>*|**B**|(2) values for fixed *<sup>γ</sup>Ne* (2) = 1.

We conclude our analysis reporting in Figure 7 the probability density functions of the location of the plasma bubbles as a function of both the magnetic latitude and LT. These distributions suggest the tendency of plasma bubbles to move away from the equator as the time goes by. This could indicate that the plasma bubbles observed near the midnight sector have not been generated in that LT sector but could have been generated near sunset, have moved at high altitudes, and successively have started to descend in altitude, moving along the magnetic field lines. Indeed, the flux tube intersecting the 460 km Swarm A orbit at ∼10◦ magnetic latitude has an apex height around 600 km at the equator where the plasma bubbles are generated. This is a crucial point in the interpretation of the LT evolution of the scaling features of magnetic field and density fluctuations inside the plasma bubbles. It makes you understand how, moving away from the sunset, Swarm A crosses older plasma bubbles at locations where the turbulence had time to fully develop.

**Figure 7.** Dependence on magnetic latitude and local time of location of the plasma bubbles.

#### **5. Discussion**

Our findings seem to suggest that a non-trivial linear relation exists between the spectral features of the plasma density and magnetic field in the plasma irregularities which depends on the local time and latitude. Actually, the fact that the magnetic field and electron density scaling properties are not simply linearly related is somehow expected. Usually, when we study the plasma bubbles dynamics, we start with a few assumptions.

First, following the derivation of the diamagnetic current in Alken et al. [64], i.e., considering the standard magnetohydrodynamic equations for a plasma in a steady-state configuration, that is, *∂<sup>t</sup>* → 0, and neglecting gravity, we can write that the plasma pressure force (∇*p*) is counter-balanced by the Lorentz force (**j** × **B**), i.e.,

$$
\nabla p = \mathbf{j} \times \mathbf{B},
\tag{4}
$$

where **B** is the ambient magnetic field. Combining Equation (4) with the Maxwell's equation describing the Ampère's circuital law in absence of displacement currents

$$
\nabla \times \mathbf{B} = \mu\_0 \mathbf{j},\tag{5}
$$

where *μ*<sup>0</sup> is the permeability of free space, we can write

$$
\nabla p = \frac{1}{\mu\_0} (\nabla \times \mathbf{B}) \times \mathbf{B}.\tag{6}
$$

By using the vector calculus identity

$$(\nabla \times \mathbf{B}) \times \mathbf{B} = (\mathbf{B} \cdot \nabla)\mathbf{B} - \nabla \left(\mathcal{B}^2 / 2\right),\tag{7}$$

we can rewrite Equation (4) as

$$
\nabla(p + \frac{B^2}{2\mu\_0}) = \frac{1}{\mu\_0}(\mathbf{B} \cdot \nabla)\mathbf{B}.\tag{8}
$$

The term on the left of Equation (8) is the sum of plasma pressure force and magnetic pressure one, while on the right, the term represents magnetic tension due to the curvature of the field lines. Taking into account that the plasma structures are much smaller than the geomagnetic field curvature radius, it is possible to assume a linear field line geometry. In this case, we can neglect the magnetic curvature term, (**B** · ∇)**B** ∼ 0, in the momentum balance so that a change in plasma pressure is, in first approximation, compensated by a change in the magnetic pressure (second assumption). Thus, according to Lühr et al. [37] and taking into account that *p* = *kBNe*(*Te* + *Ti*), we can write that across the irregularity walls, the following relation holds:

$$
\Delta \left( \frac{B^2}{2\mu\_0} \right) + k\_B \Delta [N\_t(T\_t + T\_i)] = 0 \tag{9}
$$

where Δ denotes a finite scale difference (i.e., the gradient in Equation (4) computed at the bubble intensity scales), *kB* is the Boltzmann constant, *Ne* is the electron density, and *Te* and *Ti* are electron and ion temperatures, respectively.

Third, considering that the strength of the ambient magnetic field is around four orders of magnitude greater than the magnetic field generated by the equatorial plasma bubbles, Equation (9) can be linearized as:

$$
\Delta B \approx -\frac{\mu\_0 k\_B}{B} \Delta \left[ \mathcal{N}\_\varepsilon (T\_\varepsilon + T\_i) \right]. \tag{10}
$$

Fourth, assuming a constant plasma temperature, i.e., assuming a local thermodynamic equilibrium *Te* + *Ti* = *Teq*, according to the latter expression, a high correlation must exist between Δ*B* and Δ*Ne* across plasma bubble walls, and plasma density irregularities can be detected reliably through magnetic field data, being valid the following relation

$$
\Delta B \approx \mathcal{C}\_0 \,\Delta N\_\varepsilon \,\tag{11}
$$

where *C*<sup>0</sup> is a constant. Thus, the scaling properties of the electron density and magnetic field data should be equal, leading to a universal character of the statistical properties of fluctuations as

$$
\Delta B \approx \mathbb{C}\_0 \, \Delta \mathcal{N}\_{\mathfrak{k}} \sim \pi^{\gamma(1)},\tag{12}
$$

that can be generalized as

$$
\langle |\Delta B|^q \rangle \approx \langle |\Delta \mathcal{N}\_{\mathfrak{e}}|^q \rangle \sim \pi^{\gamma(q)}.\tag{13}
$$

Conversely, our analysis find a discrepancy between the theoretical prediction by Lühr et al. [37] and the experimental results.

As observed by Alken et al. [64], studying the global flow of the gravity-driven and diamagnetic currents, the diamagnetic effect formula proposed by Lühr et al. [37] is strictly valid only when the plasma has reached equilibrium and the ambient field has no curvature. The Earth's magnetic field has a large dipole moment and the assumption to ignore the magnetic field line curvature could be one of the possible reasons responsible for the observed discrepancies. Furthermore, it is however true that the plasma bubbles are unstable systems that evolve rapidly in time. This means that the assumption of a plasma in equilibrium that, in first approximation, may be only valid in the early hours after sunset, is no longer true in the following hours when the turbulence inside plasma bubbles is fully developed [51]. This point is supported by the discrepancies from a linear relation between magnetic field and plasma density scaling exponents for increasing LT, as shown

in Figure 6, and by the geographic location of plasma bubbles in LT, as shown in Figure 7. We remark that Equation (10), being a linear approximation of Equation (4), is expected to be valid in the first phases of the plasma bubbles rising, i.e., for LT near the dusk–night sector, as highlighted in Figure 6. There is also an additional issue related to the assumption of a constant plasma (electron and ion) temperature that has to be considered. This might be true in the first instants of plasma bubbles formation, i.e., when plasma bubbles rise just after the sunset and the instability growth is still quasi-linear. Conversely, when these structures start evolving, the emerging of turbulent fluctuations can strongly affect both plasma density and temperature fluctuations which could passively respond to turbulent fluctuations. This means that the assumption of negligible temperature fluctuations could be no longer valid for hours approaching midnight. This is a crucial point that will be investigated in a successive work.

Last but not least, we note that Equation (10) neglects the role of other terms that might affect the evolution of plasma bubbles, such as for instance the emergence of **E** × **B** convective turbulence and/or Hall MHD compressive turbulence, which might be at the basis of some discrepancies between density and magnetic field intensity fluctuations. In particular, we observe that some of the features of the observed fluctuations of the magnetic field, such as the spectral slope near *k*−<sup>2</sup> are in good agreement with the emergence of 2D **E** × **B** turbulence of electron density, which is indeed isomorphic to the viscous convection of an ordinary fluid driven by temperature gradients [63]. These extra elements require the investigation of the dynamics of the magnetic field components with a further detail than that discussed here in the case of the magnetic field intensity.

#### **6. Summary and Conclusions**

In order to understand whether it is possible to study the dynamical features of plasma density irregularities by using independently either the magnetic field or the electron density measurements, this study focused on determining the relationship between the spectral features of electron density and magnetic field strength inside plasma bubbles. In the past, many important features of the equatorial plasma bubbles have been obtained indirectly from the analysis of their magnetic signatures, using the so-called diamagnetic effect. For example, it has been possible to investigate their dependence on latitude, longitude [8,19–21], and solar and geomagnetic activity [16,17] by using magnetic field data recorded by LEO satellites. However, the study of the dynamical features of these plasma density irregularities by using only the magnetic field might not be correct, since it is necessary to assume that the scaling properties of the electron density and magnetic field fluctuations are equal.

In order to address this point, we studied the scaling properties of both the electron density and the magnetic field inside plasma bubbles using measurements on board the Swarm A satellite from 1 April 2014 to 31 January 2016. These were obtained using the so-called generalized *qth*-order structure function analysis evaluated by applying the local detrended structure function analysis (DSFA) method [50,51]. The considered method, along with the time resolution of the electron density and magnetic field data, permitted investigating spatial scales ranging from some kilometers to a few hundred kilometers. At this range of spatial scales, our findings suggest that it is not correct to hypothesize that the spectral features of the plasma irregularities can be inferred from their magnetic properties. For example, this type of assumption was done in the past by Lühr et al. [15], who using 50 Hz magnetic field measurements from CHAMP obtained the spectral properties of the magnetic data and considered them also valid to investigate the spectral features of the electron density irregularities.

Our findings suggest that a more complex relation may exist between the spectral features of the electron density and magnetic field in the plasma irregularities, which depends on the local time and latitude. Indeed, the results of this research support the idea that there are some relevant discrepancies between the spectral features of electron density fluctuations and magnetic field intensity ones. The discrepancies depend on the evolution of plasma bubbles with local time and on its turbulent nature. In the early stage of plasma bubble formation, i.e., when plasma bubbles rise just after the sunset and the instability growth is still quasi-linear, we find a quasi-linear relationship on average between the scaling properties of electron density and those of magnetic field. Conversely, when these structures develop, the evolution of turbulent fluctuations can strongly affect also other physical quantities such as the plasma temperature fluctuations, which could modify the spectral features. This means that the assumption of negligible temperature fluctuations could no longer be valid for hours approaching midnight. In other words, the electron temperature fluctuations could play an important role when turbulence is fully developed.

A better comprehension of the plasma bubbles dynamics and of the turbulence processes that characterize their time evolution may benefit from the use of very high-resolution vector magnetic field and plasma density measurements such as those that should be available from the future NanoMagSat mission [65]. This future satellite, which should fly at an altitude of 500 km and an inclination angle of 60◦, should allow a quick full local time coverage and for the first time joint measurements of electron density and magnetic signals at very high frequencies. The mission should give the opportunity to investigate plasma variability at very short spatial scales ranging from a few meters to some kilometers.

**Author Contributions:** Conceptualization, P.D.M.; methodology, P.D.M. and G.C.; investigation, P.D.M., G.C. and T.A.; data curation, R.T., F.G., M.P., A.P. and I.C.; writing—original draft preparation, P.D.M.; writing—review and editing, all. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received financial support from Progetto INGV Pianeta Dinamico funded by MIUR (Fondo finalizzato al rilancio degli investimenti delle amministrazioni centrali dello Stato e allo sviluppo del Paese, legge 145/2018)-Task A1-2021 under Grant codice CUP D53J19000170001.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Publicly available data were analyzed in this study. ESA-Swarm data are available from the ESA dissemination server (http://earth.esa.int/swarm, accessed on 13 February 2022). OMNI data are available at NASA CDAWeb site (https://cdaweb.gsfc.nasa.gov/ index.html/, accessed on 13 February 2022). Elaborated data are available on request from the corresponding author.

**Acknowledgments:** The results presented rely on data collected by ESA-Swarm mission. We thank the European Space Agency that supports the Swarm mission. Swarm data can be accessed at http://earth.esa.int/swarm (accessed on 12 February 2022). The authors kindly acknowledge V. Papitashvili and J. King at the National Space Science Data Center of the Goddard Space Flight Center for permission to use the OMNI data and the NASA CDAWeb team for making these data available. This work was supported in part by Progetto INGV Pianeta Dinamico funded by MIUR (Fondo finalizzato al rilancio degli investimenti delle amministrazioni centrali dello Stato e allo sviluppo del Paese, legge 145/2018)-Task A1-2021 under Grant codice CUP D53J19000170001

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:



#### **References**


## *Article* **Pressure-Gradient Current at High Latitude from Swarm Measurements**

**Giulia Lovati 1, Paola De Michelis 2,\*, Giuseppe Consolini <sup>3</sup> and Francesco Berrilli <sup>4</sup>**


**Abstract:** The pressure-gradient current is among the weaker ionospheric current systems arising from plasma pressure variations. It is also called diamagnetic current because it produces a magnetic field which is oriented oppositely to the ambient magnetic field, causing its reduction. The magnetic reduction can be revealed in measurements made by low-Earth orbiting satellites flying close to ionospheric plasma regions where rapid changes in density occur. Using geomagnetic field, plasma density and electron temperature measurements recorded on board ESA Swarm A satellite from April 2014 to March 2018, we reconstruct the flow patterns of the pressure-gradient current at highlatitude ionosphere in both hemispheres, and investigate their dependence on magnetic local time, geomagnetic activity, season and solar forcing drivers. Although being small in amplitude these currents appear to be a ubiquitous phenomenon at ionospheric high latitudes characterized by well defined flow patterns, which can cause artifacts in the main field models. Our findings can be used to correct magnetic field measurements for diamagnetic current effect, to improve modern magnetic field models, as well as to understand the impact of ionospheric irregularities on ionospheric dynamics at small-scale sizes of a few tens of kilometers.

**Keywords:** ionosphere F region; high latitude; pressure-gradient current; diamagnetic current; swarm measurements

#### **1. Introduction**

The Earth's magnetosphere is characterized by different plasma regions, each one with particular density and temperature distributions [1]. Different currents flow in the magnetosphere characterized by patterns and intensities changing in time and space. Some of these currents are created by plasma space inhomogeneity and driven by plasma pressure gradients. The neutral sheet current, the magnetopause current and the ring current are examples of pressure-gradient currents, which flow in the Earth's magnetosphere [2].

This phenomenon also occurs in the ionosphere where the pressure-gradient current is among the weaker ionospheric current systems arising from plasma pressure variations. Indeed, due to the coupling between geomagnetic field and plasma pressure gradient, electrons and ions drift in opposite directions generating an electric current whose intensity is of the order of a few nA/m2. This current is also called diamagnetic because it produces a magnetic field, which is oriented oppositely to the ambient magnetic field, causing its reduction. The magnetic reduction associated with the pressure-gradient current can be revealed in measurements made by low-Earth orbiting (LEO) satellites flying close to ionospheric plasma regions where rapid changes in density occur. Anyway, identifying diamagnetic current by using its magnetic signature is not easy due to the weak intensity of generated magnetic perturbation, that is about 10,000 time smaller than the ambient geomagnetic field. This is why the studies investigating this current are quite recent, since high-accuracy satellite magnetic field measurements are available. Due to its origin, the

**Citation:** Lovati, G.; De Michelis, P.; Consolini, G.; Berrilli, F. Pressure-Gradient Current at High Latitude from Swarm Measurements. *Remote Sens.* **2022**, *14*, 1428. https://doi.org/ 10.3390/rs14061428

Academic Editor: Michael E. Gorbunov

Received: 27 January 2022 Accepted: 9 March 2022 Published: 15 March 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

diamagnetic current can be revealed at both low and high latitudes, and more generally in all those regions where the plasma pressure gradients are greatest.

Several studies focus on low latitude, in the equatorial belt, where irregular plasma density structures of various scale sizes, from centimeters up to hundreds of kilometers, can be observed after sunset and where exists a region with increased plasma density, known as Appleton or equatorial ionization anomaly. About twenty years ago, Lühr and collaborators [3] were the first to reveal the presence of these currents in the ionosphere. By using the high-resolution magnetic field data in combination with plasma density observations recorded by CHAMP (CHAllenging Minisatellite Payload) satellite [4], they revealed the magnetic signatures of plasma density variations at low latitudes and reconstructed the strength and geographic location of diamagnetic field caused by the Appleton anomaly. They estimated the strength of this field up to 5 nT at a variety of local times, showed how the diamagnetic effect could severely affect the correct estimation of the ionospheric current distributions from space, and how it could be capable of producing artifacts in the new generation main field models [5].

In the following years some other papers dealt with this subject. Interesting is the study carried out by Alken et al. [6] where the global gravity and diamagnetic current systems in the ionospheric topside F2 region were modeled by using the Thermosphere Ionosphere Electrodynamics General Circulation Model (TIEGCM) [7] with empirical density, wind and temperature inputs. The results obtained from this model permitted to reconstruct the diamagnetic current structure at low latitude, which was found to be prominent near the equatorial ionization anomaly, and to calculate the associated magnetic field. Although the obtained results showed some discrepancies with what found using the diamagnetic effect formula proposed by Lühr et al. [3], the study confirmed the need to pay attention to all those ionospheric current systems that, although smaller than others well-known ionospheric current systems, could play a role in the ionospheric dynamics and help in the development of more accurate magnetic models.

Successively, by considering 10 years of CHAMP measurements (2000–2010) and two years of Swarm data (late 2013–2015), Alken [8] reconstructed the gravity and diamagnetic currents flowing at low latitude in the topside F2 region. He found that the strength of these currents depended on season, local time and solar activity. These currents were strongest during spring and fall and their magnetic signatures ranged from 1 to 7 nT, depending on solar activity level. Another important feature of these currents was their presence until midnight. For this reason, the assumption made by many magnetic models according to which the ionospheric currents were negligible during the nighttime and geomagnetically quiet periods, was not correct.

Recently, awareness of the existence of a large amount of plasma density irregularities at high latitude, especially at polar latitudes, has led to the hypothesis that diamagnetic currents had to exist also in these regions, where they could play a role in the ionospheric dynamics. Here, the pressure-gradients are mainly related to the convection pattern and to electron precipitation, that characterizes both the polar cap and the auroral oval [9]. In particular, soft electron precipitation caused by the direct entry of magnetosheath electrons in the cusp region and from electrons of ionospheric origin that are energized at intermediate altitudes by wave-particle interactions with dispersive Alfvén waves, are supposed to generate density anomalies with elevated electron temperature in the F region ionosphere [10]. At the present time, few studies have been published on the diamagnetic currents at high latitudes: we know only two papers where these currents are investigated analysing their magnetic signatures. The first paper is by Park et al. [9] and the second one by Laundal et al. [11]. Using magnetic data from CHAMP satellite and a combination of plasma density measurements from both the Digital Ion Drift Meter and the Planar Langmuir Probe on board CHAMP, Park et al. [9] presented a climatological study of the diamagnetic signatures at high latitude due to ionospheric plasma irregularities, finding that ionospheric irregularities did not have a uniform zonal spatial distribution, but exhibited a maximum at pre-midnight ionosphere and dayside cusp. The geographic location

of the irregularities also depended on geomagnetic activity and season. In this paper, the authors do not reconstruct the diamagnetic currents but use their magnetic signatures in order to investigate the occurrence rate of ionospheric irregularities. Successively, analyzing two years of magnetic and plasma density measurements from Swarm constellation, Laundal et al. [11] showed that the magnetic field variations recorded at satellite altitude were often associated with plasma pressure variations in the polar cap. They found that superposed on the magnetic perturbations due to the auroral electrojet currents was an irregular pattern, anticorrelated with plasma density measurements, that was the signature of the diamagnetic current. They suggested that this current could dominate the disturbances in the magnetic field strength at small-scale size of a few tens of kilometers.

Here, using geomagnetic field, plasma density and electron temperature measurements on board Swarm constellation, we reconstruct the flow patterns of the diamagnetic current at high-latitude ionosphere in both hemispheres, and investigate their dependence on geomagnetic activity, season and solar forcing drivers.

#### **2. Data and Method**

We study the pressure-gradient current at ionospheric high latitude using measurements on board Swarm constellation [12]. This constellation consists of three identical satellites launched into a near polar orbit at the end of 2013. Two of these satellites (Swarm A and Swarm C) fly side-by-side (1.4◦ separation in longitude) at an altitude around 460 km (initial altitude), whereas the third satellite (Swarm B) flies at higher orbit, around 510 km (initial altitude). These satellites are equally equipped and simultaneously acquire magnetic field, plasma density and electron temperature measurements, thanks to the presence in each satellite of a Vector Field Magnetometer (VFM) and a Electric Field Instrument (EFI) [12,13].

In this study, we use magnetic field components, electron density *Ne* and electron temperature *Te* measurements acquired by Swarm A [14], at a rate of 1 Hz, during a period of four years from 1 April 2014 to 28 February 2018. Data are selected from the ESA dissemination server (http://earth.esa.int/swarm accessed on 14 March 2022). Starting from April 2014, the selected time interval excludes the period where Swarm constellation was stabilizing its orbit and makes possible to assume that Swarm A is approximately at the same altitude for all the period under analysis. The choice of February 2018 is due to the lack of the Auroral Electrojet index data beyond this period.

Since the ionospheric currents are mainly organized by Earth's magnetic field, the study is done using a magnetic coordinate system, the non-orthogonal Quasi-Dipole (QD) system of coordinates [15–17]. We select data with a quasi-dipole magnetic latitude *λQD* ≥ |50◦| in order to focus on mid- and high-latitude ionospheric topside F2 region, and we use the Magnetic Local Time (MLT) [18] instead of universal time (UT) in order to consider the position respect to the Sun. To report our data on the *λQD*-MLT plane, we realize maps of our parameters with a 1◦ × 1◦ binning (1◦ longitude corresponds to 4 min in MLT). Time series belonging to each bin are initially filtered removing values beyond 3*σ* from the mean value in order to take off possible spikes. The value representative of each bin corresponds to the average value of the filtered time series.

We use the Auroral Electrojet AE index to study the dependence of the pressuregradient currents on the geomagnetic activity. This index is a good proxy of the global geomagnetic activity at high latitude, well describing the auroral zone magnetic activity due to the increase of the ionospheric currents flowing below and within the auroral oval. The AE time series has been downloaded from the Space Physics Data Facility of the NASA Goddard Space Flight Center (OMNI data set http://cdaweb.gsfs.nasa.gov/) where it is available only until 28 February 2018.

To analyse the dependence on the solar activity of the flow patterns of the diamagnetic current at high latitude we consider two different proxies: the Mg II core-to-wing ratio and the F10.7 index. The Mg II index refers to a broad absorption feature with narrow emission peaks in the core. This emission doublet, near 280 nm, originates in the Sun's chromosphere, while the line wings part originates in the photosphere, showing much less variability. Therefore, the ratio of line core intensity to wing intensity provides a good estimate of solar variability, avoiding degradation effects [19]. For the selected period Mg II is provided by GOME-2 (Global Ozone Monitoring Experiment) mission instruments [20]. F10.7 index, instead, is the radio emission at 10.7 cm, that originates from regions of intense magnetic field, characterized by structures like plages, networks, and sunspots [21]. It is recorded by the Dominion Radio Astrophysical Observatory (Canada) with three flux determinations for each day [22]. The time series containing these two proxies have been downloaded from the Space Physics Data Facility of the NASA Goddard Space Flight Center (OMNI data set http://cdaweb.gsfs.nasa.gov/ accessed on 14 March 2022) in the case of F10.7, and from UVSAT Bremen University (GOME-2A, http://www.iup.uni.bremen.de/UVSAT/Datasets/mgii accessed on 14 March 2022) in the case of MgII.

To reconstruct the flow patterns of the pressure-gradient current at high latitude we can observe that for time scales longer than about 1 min the ionospheric electrodynamics is approximately in a steady-state with divergence-free current density (-*J*), and time-invariant electric field (- *E*). From the generalized Ohm's law, the current density -*J* can be expressed as the sum of different contributes due to the neutral wind (-*Jw*), the global electrostatic electric field (*σ*- *E*, with *σ* equal to the conductivity tensor), gravitational force (-*Jg*), and pressuregradient force (-*JP*). At large spatial scale the electric force and the neutral wind dominate the pressure-gradient and the gravitational forces. Thus, the first two forces are responsible for the main ionospheric currents while the other two generate currents that can usually be neglected. However, these currents can become important on smaller spatial scale, and can have non-negligible effects on low Earth orbit satellite because of their magnetic perturbations. By assuming that it is possible to ignore the effects of neutral collisions, the current density due to pressure-gradient force can be expressed as in Equation (1) [15,23].

$$f\_P^\* = -\frac{\nabla P \times \vec{B}}{B^2} \tag{1}$$

where *P* = *NekB*(*Ti* + *Te*) is the plasma pressure, *kB* is the Boltzmann constant, *Te* is the electron temperature, *Ti* is the ion temperature, and - *B* is the ambient magnetic field, which is produced by all sources within and outside the solid Earth up to the so-called magnetopause. The plasma pressure-gradient current, which is mainly in the zonal direction due to the strong vertical gradient of the plasma density, is naturally responsible of variations in the magnetic field strength. In first approximation, the magnetic field variations associated with this current can be evaluated by considering the balance between the plasma and magnetic pressures. According to Lühr et al. [3], in a quasi-stationary plasma where the structures are much smaller than the geomagnetic field curvature radius, the signature in the magnetic variation strength due to plasma pressure gradient is equal to:

$$
\Delta B = -\frac{\mu\_0 k\_B}{B} \Delta \left( N\_\ell \left( T\_\ell + T\_{\bar{l}} \right) \right) \tag{2}
$$

where *μ*<sup>0</sup> is the susceptibility of free space, and *kB* is the Boltzmann constant. Here, the main assumption is that the measured magnetic fluctuations are spatial on the time scale for the satellite to traverse them. This interpretation is valid when the wave frequency in the plasma rest frame is much less than the Doppler frequency due to spacecraft motion across the horizontal spatial structure. Assuming a local thermodynamic equilibrium, the magnetic field perturbations are anticorrelated with variations in plasma density: a plasma density depletion is associated with a magnetic field enhancement while a reduction in the magnetic field is related to an increase in plasma density. The current density driven by pressure gradient flows around the density enhancement in a direction that produces a magnetic field capable of reducing the ambient magnetic field.

To reconstruct the flow patterns of the plasma pressure-gradient current at Swarm altitude, it is necessary to estimate the plasma pressure. This requires to know the temperature and electron density distributions within the region of interest, being *P* = *NekB*(*Te* + *Ti*), according to the equation of state. Swarm satellite measures only the electron density and electron temperature along its orbit but, at its altitude, the temperature of oxygen ions (*O*+), which can be assumed as the main ionic species present, is of the same order of magnitude as the electron temperature and we make the assumption that *Te Ti* [11,24]. Thus, in first approximation, we assume that plasma pressure is equal to *P* = 2*NekBTe*. From *P* values it is possible to compute the pressure gradient (∇*P*) on the horizontal plane. In order to approximate the geometry of the system, we consider the gradient in spherical coordinates:

$$\nabla P = \frac{\partial P}{\partial r}\vec{r} + \frac{1}{r}\frac{\partial P}{\partial \theta}\vec{\theta} + \frac{1}{r\sin\theta}\frac{\partial P}{\partial \phi}\vec{\phi} \tag{3}$$

where *θ* is the magnetic co-latitude, and *φ* is the longitude, which corresponds to MLT. Since Swarm A flies approximately at a constant altitude, we have *r const* and in addition we are not able to evaluate the gradient along this direction. Therefore, we estimate the pressure gradient in the plane perpendicular to the radial direction:

$$\nabla P\_{\perp} = \frac{1}{r} \frac{\partial P}{\partial \theta} \vec{\theta} + \frac{1}{r \sin \theta} \frac{\partial P}{\partial \phi} \vec{\phi} \tag{4}$$

However, considering Equation (1), we remark that the neglected *<sup>∂</sup><sup>P</sup> <sup>∂</sup><sup>r</sup>* term provides a small contribution to the plasma pressure-gradient current being the magnetic field nearly parallel to it. Indeed, at high latitudes the magnetic field is approximately aligned with the radial direction and, compared to this term, the other two components result negligible. Therefore, in this approximation, the lack of information on *<sup>∂</sup><sup>P</sup> <sup>∂</sup><sup>r</sup>* does not affect significantly our analysis.

In order to estimate *<sup>∂</sup><sup>P</sup> ∂θ* and *<sup>∂</sup><sup>P</sup> ∂φ* , we consider the matrix of pressure binned values, having a number of rows equal to QD-latitude bins and a number of column equal to MLT bins. Having selected latitudes from |50◦| upwards and a binning of 1◦ × 1◦, it results in a 40 × 360 matrix. The gradient is evaluated through second order accurate central differences [25]. Once obtained *<sup>∂</sup><sup>P</sup> ∂θ* and *<sup>∂</sup><sup>P</sup> ∂φ* , we complete the gradient estimation in spherical coordinate using Equation (4) considering *r* as fixed at the altitude of Swarm A, and *θ* equal to the mean QD-latitude of each bin under consideration. Finally, the current density due to pressure-gradient force is estimated according to Equation (1).

#### **3. Results**

#### *3.1. Variation of Pressure-Gradient Current with Geomagnetic Activity Level*

The analysis of the current density spatial distribution due to pressure-gradient force at high latitude in the Northern and Southern hemispheres starts investigating the dependence on the geomagnetic activity level. Using the AE index we select two different levels of geomagnetic activity. Quiet periods correspond to values of *AE* ≤ 80 nT while disturbed ones to *AE* ≥ 120 nT [26,27]. The two selected activity levels are chosen in order to clearly separate active and quiet periods according to the distribution function of AE index values, which is bimodal with a crossover around 100 nT. For this type of analysis we consider the whole dataset, that is, the electron density, electron temperature, and magnetic field recorded by Swarm A from 1 April 2014 to 28 February 2018.

Figures 1–3 report a polar view of the large-scale spatial distribution of electron density (*Ne*), electron temperature (*Te*), and plasma pressure (*P*) in the Northern and Southern hemispheres during quiet and disturbed periods, respectively. Data are mapped into grid binned at 1◦ × 1◦ in QD latitude and MLT coordinates. The two selected geomagnetic activity levels correspond to different spatial distributions of the plasma parameters in the two hemispheres. Focusing on plasma pressure, we notice that its spatial distribution shows a clear dependence on magnetic latitude and magnetic local time (see Figure 3). It is characterized by a marked dayside/nigthside asymmetry, which mainly reflects the dayside/nigthside asymmetry of the electron density and temperature spatial distributions. Indeed, as reported in Figure 1, the electron density is generally almost twice in the dayside compared to the nightside due to the dependence of the electron density on the upper atmosphere ultraviolet ionization.

**Figure 1.** Polar view of the average spatial distribution of electron density (*Ne*) for quiet (**left** column) and disturbed (**right** column) geomagnetic activity conditions, for the Northern and Southern hemispheres, in QD-latitude (*λQD* > 50◦ and *λQD* < −50◦ for the Northern and Southern hemisphere, respectively) and MLT reference system. Maps are obtained using data recorded by Swarm A from 1 April 2014 to 28 February 2018. The concentric dashed circles represent contours of magnetic latitude separated by 10◦.

Moreover, a marked dayside/nigthside asymmetry also characterizes the spatial distribution of the electron temperature (see Figure 2) which is, even in this case, due to the increase of ultraviolet photoionization. However, the spatial distribution of the electron temperature also depends on the particle precipitation caused by magnetosphereionosphere coupling. Indeed, soft particle precipitation heats the F region ionosphere and is a likely source of the observed correlation of the anomalous density with elevated particle temperature. Thus, the magnetosphere-ionosphere coupling is responsible of the temperature enhancement both in the cusp region around noon [28], where particles of solar origin are directly injected into the ionosphere, and at auroral and sub-auroral latitudes

where particles arrive directly from the geomagnetic tail regions and plasmasphere during disturbed periods (see Figure 2).

Not all the features of the electron density and electron temperature are recognizable in the plasma pressure spatial distribution. For example, the so-called main ionospheric trough (MIT), which corresponds to a prominent plasma depletion at subauroral regions during the night hours in a latitudinally limited band between 60◦ and 70◦ (see Figure 1), is not clearly recognizable in the pressure data. It is due to the fact that there is an electron temperature enhancement owing to the joint action of particle precipitation and decreased collisional cooling in correspondence with the main ionospheric trough [29,30]. In the nightside, at auroral and sub-auroral latitudes this process is dominant and becomes more and more important with the increase of the geomagnetic activity level (see Figure 2). Due to the anticorrelation between the electron density and electron temperature, the plasma pressure does not show a particular behaviour in this region. Conversely, Figure 3 reveals a plasma pressure depletion at low latitudes, between 50◦ and 65◦, in the nightside, between 21:00 and 06:00 MLT, whose position moves progressively to a lower latitude as the level of geomagnetic activity increases.

**Figure 2.** Polar view of the average spatial distribution of electron temperature (*Te*) for quiet (**left** column) and disturbed (**right** column) geomagnetic activity conditions, for the Northern and Southern hemispheres, in QD-latitude (*λQD* > 50◦ and *λQD* < −50◦ for the Northern and Southern hemisphere, respectively) and MLT reference system. Maps are obtained using data recorded by Swarm A from 1 April 2014 to 28 February 2018. The concentric dashed circles represent contours of magnetic latitude separated by 10◦.

**Figure 3.** Polar view of the average spatial distribution of plasma pressure for quiet (**left** column) and disturbed (**right** column) geomagnetic activity conditions, for the Northern and Southern hemispheres, in QD-latitude (*λQD* > 50◦ and *λQD* < −50◦ for the Northern and Southern hemisphere, respectively) and MLT reference system. Maps are obtained using data recorded by Swarm A from 1 April 2014 to 28 February 2018. The concentric dashed circles represent contours of magnetic latitude separated by 10◦.

Another interesting feature of the plasma pressure is its marked increase in the cusp region around noon, which becomes more marked in the disturbed periods in both hemispheres. This increase is associated with the observed electron temperature and density enhancements.

Lastly, Figure 3 shows two other remarkable features, which mainly reflect *Ne* spatial distribution features. There is an enhancement of plasma pressure from noon across the polar cap to the nightside during disturbed conditions in both hemispheres. This plasma pressure enhancement corresponds to the characteristic tongue of ionization (TOI) of *Ne* [31,32], which is a large-scale feature of the F-region polar ionosphere. Furthermore, there is a plasma pressure increases around dusk in disturbed conditions reflecting the wellknown storm-enhanced density (SED) plumes, which are prominent ionospheric electron density increases at the dayside mid and high latitudes [33].

Figure 4 reports the flow patterns of the plasma pressure-gradient current (black arrows) at high latitudes in the Northern and Southern hemispheres during the two different geomagnetic activity levels. Only for graphical reasons, the black arrows are mapped into

grids binned at 2◦ × 4◦ in QD latitude and MLT coordinates, where 4◦ magnetic longitude corresponds to 16 min. Below the flow patterns of the reconstructed current there is the corresponding plasma pressure map, which permits us to better visualize the regions where the current develops. As expected, the current tends to flow around the plasma pressure enhancements and is stronger in region where the plasma pressure changes rapidly. Some well defined patterns can be recognized in the plasma pressure-gradient current regardless of geomagnetic activity level.

**Figure 4.** Flow patterns of pressure-gradient current (black arrows) superimposed on the plasma pressure spatial distribution for quiet (**left** column) and disturbed (**right** column) geomagnetic activity conditions, for the Northern and Southern hemispheres, in QD-latitude (*λQD* > |50◦|) and MLT reference system. Maps are obtained using data recorded by Swarm A from 1 April 2014 to 28 February 2018. For graphical reasons, the current's vector field is mapped into grids binned at 2◦ × 4◦, where 4◦ magnetic longitude corresponds to 16 min. The concentric white circles are plotted in 10◦ intervals, corresponding to QD-latitudes of |80◦|, |70◦| and |60◦| starting from the centre, respectively.

A vortex exists around the cleft region, whose position shifts equatorward as magnetic activity increases, following the natural shift of center of the precipitation region. Around the cleft region the current flows counter-clockwise in the Northern hemisphere and clockwise in the Southern one. This current flows in a direction which reduces the ambient geomagnetic field. At high latitude the magnetic field is nearly vertical and the horizontal pressure-gradient current produces a magnetic field which is, in first approximation, parallel to the main field but with opposite direction. Thus, the current produces a

reduction of the magnetic field inside the plasma region around which the current flows. The existence of this current structure is indirectly confirmed by Park et al. [9] analysing the diamagnetic signatures of high-latitude ionospheric irregularities. In both hemispheres they found depletions in magnetic field strength of about 1 nT concentrated in the cusp region that they associated with the presence of local plasma irregularities with scale sizes below some hundred kilometers. Anyway, this region can be also characterized by large anomalies in the electron density correlated with intense small-scale magnetic fluctuations, which are associated with incident Alvén waves producing a local heating. This means that in this region an opposite scenario to that implicitly assumed may be possible: the observed pressure gradients may be also produced by magnetic disturbances [10].

At high latitudes (*λQD* > |75◦|) the current flows counter-clockwise in the Northern hemisphere and in opposite direction in the Southern one. The spatial distribution of the current seems to identify, especially during disturbed periods, another region characterized by an enhancement of plasma density, which is due to the propagation of plasma from the cusp to the polar cap. This region, that is visible in both hemispheres, is known to be characterized by the presence of plasma instability and the formation of ionospheric irregularities known as polar cap patches, where the density is at least twice that of the background [34]. However, the reconstructed flow patterns of the pressure-gradient currents capture plasma pressure variations which are not necessary due to the presence of plasma density irregularities as for example polar cap patches and auroral blobs. At lower latitudes in the Northern and Southern hemispheres the flow pattern of the pressuregradient current identifies another region where the ambient magnetic field can be reduced by this current. In first approximation, this region corresponds to the auroral zone from 15:00 MLT to 09:00 MLT, passing by 24:00 MLT. Thus, the pressure-gradient currents characterize also the auroral zone and equatorward of the auroral zone on the nigthside.

In general, the flow patterns of the pressure-gradient current are similar in both hemispheres and what one finds is their shift towards lower latitudes with the increase of the geomagnetic activity.

Figure 5 reports the distributions of the pressure-gradient current values for disturbed and quiet geomagnetic activity conditions for the Northern and Southern hemispheres. The comparison among the four distributions shows that the pressure-gradient current intensity is slightly higher in the Southern hemisphere than in the Northern one, while there is no great difference between quiet to disturbed periods. At high latitudes the mean values of these currents are quite low, around 1 order of magnitude less than the same diamagnetic currents observed at low latitudes [35]. Our findings are in agreement with previous studies. Analyzing two years of Swarm measurements Laundal et al. [11] estimated the small variations in the magnetic field strength due to pressure-gradient currents at polar latitudes by evaluating the magnetic field intensity variations anticorrelated with plasma density variations. As result of their statistical analysis, they showed that the polar caps, the auroral zone, and equatorward of the auroral zone on the nightside were the regions where the magnetic field variations were well explained by plasma pressure gradients. Their findings are slightly different from the results by Park et al. [9], where the highest occurrence rates between the magnetic field and the electron density variations are found in the auroral oval. In both papers the flow pattern of the pressure-gradient current is not calculated, nor its intensity. The existence of this current is invoked to justify the variations of magnetic field correlated with those of plasma density and climatological analyses of their correspondence are proposed. This means that, to make a comparison with our results, it is necessary to consider that the pressure-gradient currents flow around the plasma pressure enhancement regions, which are often characterized by the occurrence of ionospheric irregularities and are the regions identified in the previous papers. It is also important to notice that both Laundal et al. [11] and Park et al. [9] assume that the variations in plasma pressure are dominated by electron density, under the hypothesis of a local thermodynamic equilibrium, but looking at Figure 2 we notice that also the temperature plays a role in the plasma pressure, mainly in the auroral and sub-auroral latitudes.

**Figure 5.** Distributions of the pressure-gradient current values for disturbed and quiet geomagnetic activity conditions for the Northern and Southern hemispheres. Solid line represents the continuous probability density curve and it is plotted above the more transparent relative histogram.

To evaluate the pressure-gradient current we have used the magnetic field measurements recorded by Swarm A. It contains all the different contributes coming from its sources both internal (core and crust) and external (magnetospheric and ionospheric currents) to the Earth. In order to investigate the possible effect on the pressure-gradient current flow patterns of remote currents, in particular of those current systems that flow about 300 km below Swarm A and are expected to produce a smooth magnetic signature on the measured magnetic field [36], we evaluate the pressure-gradient current considering in this last case the magnetic field strength associated with internal sources. The internal part of the magnetic field is estimated using the CHAOS model [37]. Thus, the pressure-gradient current is evaluated and compared with that obtained using the ambient magnetic field. Figure 6 reports the percentage difference between the current evaluated in the two ways. This procedure should leave only the contribute of external sources to the pressure-gradient current. The effect on the pressure-gradient currents at Swarm altitude of the remote current systems is very low: it does not exceed 0.3%. However, it is not uniform, depending on the geographic location, magnetic local time, and geomagnetic activity level. Our findings show that the features of the reconstructed pressure-gradient currents at Swarm altitude depend only partially on geomagnetic field variations due to external disturbances, mainly they depend on the main field structures and on the features of pressure gradients.

**Figure 6.** Difference in percentage between the pressure-gradient current evaluated using the magnetic field measured by Swarm A and that evaluated using the magnetic field of internal origin, in the case of quiet (**left** column) and disturbed (**right** column) geomagnetic activity conditions, for the Northern and Southern hemispheres, in QD-latitude (*λQD* > 50◦ and *λQD* < −50◦ for the Northern and Southern hemisphere, respectively) and MLT reference system.

#### *3.2. Variation of Pressure-Gradient Current with Season*

In order to investigate the seasonal dependence of the pressure-gradient current, we consider the period between 5 May and 5 August of each year to identify the Northern local summer and the Southern local winter, and the period between 5 November and 5 February of each year to identify the Northern local winter and Southern local summer. Thus, by portioning the entire dataset according to summer and winter solstices, we can study the dependence of the pressure-gradient current on solar illumination.

Figure 7 reports the obtained plasma pressure distribution according to the selected seasons for both hemispheres, together with the resulting flow patterns of the pressuregradient current. The current is described using black arrows whose lengths are proportional to the intensity. Differently from previous figures, we use different scales for the plasma pressure in the two seasons. Indeed, the plasma pressure is almost three times higher in summer than in winter. Looking at Figure 7 we notice that from winter to summer, there is an increase of plasma pressure in the dayside ionosphere at all MLTs, especially between the noon and the dusk at middle latitudes. Plasma pressure increases even in the post-nightside sector at low latitudes in the Southern hemisphere, and at high latitudes in

the polar cap. Lastly, an increase is observable in the polar cusp during summer in both hemispheres.

**Figure 7.** Flow patterns of pressure-gradient current (black arrows) superimposed on the plasma pressure distribution during summer (**left** column) and winter (**right** column) periods, in the Northern and Southern hemispheres. Maps are in QD-latitude (*λQD* > |50◦|) and MLT reference system. The binning window is 1◦*λQD* × 4 MLT for the pressure map, while for the current's vector field the window is 2◦*λQD* × 16 MLT for graphical reasons. The concentric white circles are plotted in 10◦ interval, corresponding to QD-latitudes of |80◦|, |70◦| and |60◦| starting from the centre.

The plasma pressure distributions reported in Figure 7 show that summer-winter asymmetry is more pronounced in the Southern hemisphere than in the Northern one. This is due to the different summer-winter asymmetry of *Ne* spatial distribution between the two hemispheres, which is a combination of two different effects, the first produced by the solar radiation and the second by the F region annual anomaly [38,39]. Indeed, *Ne* values are higher during summer than during winter, as a result of a larger solar radiation but, at the same time, they are higher during December solstice than during June one. Since in the Southern hemisphere the increase in *Ne* values happens during summer season, it amplifies the difference with respect to local winter. Conversely, the F region annual anomaly has the effect of reducing the asymmetry between local seasons in the Northern hemisphere.

The dependence of the plasma pressure on the seasons also involves a dependence of the pressure-gradient current on the solar illumination. The flow patterns of the currents and their strengths change from winter to summer. The pressure-gradient currents are visible especially in summer where the structures, which have previously identified by studying the dependence on the geomagnetic activity levels, can be easily recognizable. Even in this case, we can observe three regions characterized by the presence of these currents, one around cleft region, another at high latitudes around the magnetic poles, and the third confined at auroral latitudes. These regions should identify the zones characterized by the occurrence of ionospheric irregularities.

As regards the current intensity in the two hemispheres during the two selected seasons, we report the distributions of the pressure-gradient current values in Figure 8. We notice that the current intensity is higher during summer than winter, when more plasma is produced by sunlight and plasma pressure gradients are greater. The current intensity is greater in the Southern hemisphere than Northern one.

#### *3.3. Variation of Pressure-Gradient Current with Solar Activity*

To evaluate the level of solar activity during the selected period of Swarm observations we use, as previously described, the daily solar flux at 10.7 cm (F10.7) and the Mg II coreto-wing ratio as proxies. The values of F10.7 and Mg II from 1 April 2014 to 31 December 2017 are reported in the lower panel of Figure 9. Here, the vertical dotted lines indicate the transition from one year to the next, the horizontal lines identify the yearly mean values, and the shaded parts cover the values between the yearly mean values ± one standard deviation year by year for each index. The analysed period corresponds to the decreasing phase of the solar cycle, and is characterized by a decrease in the F10.7 and Mg II yearly mean values.

For what concerns the pressure-gradient current intensity, the most remarkable features emerging from Figure 9 are a decrease in maximum and minimum values of the current intensity with the decrease of solar activity as well as a decrease of the intensity around the cusp. Moreover, regardless of the solar activity level, the pressure-gradient current is substantially more intense in the Southern hemisphere than in the Northern one at almost all magnetic latitudes and MLTs. In the Southern hemisphere, besides the three main regions already identified, another region of enhanced current is visible. It is at a latitude between 70◦ and 80◦ in the nightside. This structure is well visible in 2014 and 2015 when the solar activity is higher.

**Figure 9.** From top to bottom: Polar view of the average spatial distribution of the pressure-gradient current intensities, year by year, for the Northern and Southern hemisphere, respectively; F10.7 (blue) and the Mg II (red) indices during the 4 years of Swarm observations. The vertical dotted lines indicate the transition from one year to the next, the horizontal lines identify the yearly mean values, while the shaded parts cover the values between the yearly mean values ± one standard deviation year by year for each index.

#### **4. Discussion and Conclusions**

The aim of the present research was to examine the flow patterns of the pressuregradient currents and their strengths at high-latitude ionospheric topside F2 region in both hemispheres, and investigate their dependence on magnetic local time, geomagnetic activity, season and solar forcing drivers. All this has been achieved thanks to the geomagnetic field, electron density and electron temperature variations observed by Swarm A satellite from 2014 to 2018. The flow patterns of the pressure-gradient currents have been estimated under the hypothesis that it is possible to ignore the effects of neutral collisions and consider the electron and ion temperature equal. Our study is the first to show polar flow patterns of the pressure-gradient currents and their strengths at Swarm altitude. From the observations the following conclusions can be drawn:

• During geomagnetically disturbed periods (*AE* ≥ 120 nT) the plasma pressure gradients are particularly large around cleft region, where the electron density is changing rapidly. The pressure-gradient current flows around this plasma pressure enhancement region in both hemispheres. The existence of this flow pattern agrees with an increased probability of finding magnetic field variations correlated with plasma density ones around cleft region [9]. Anyway, we remark that additional contributions

to the observed magnetic field variations and plasma pressure gradients can be due to incident Alfvén wave [10];


The findings of this study confirm that the strength of the pressure-gradient currents is weak in the high-latitude ionospheric topside F2 region, at Swarm A altitude (460 km), but these currents appear preferentially at the same geographic locations regardless of geomagnetic activity, season and solar activity. This means that the magnetic field variations associated with these flow patterns can induce disturbances on the ionospheric magnetic field measurements and cause artifacts in main field models.

**Author Contributions:** Conceptualization, G.L., P.D.M. and G.C.; methodology, software and formal analysis, G.L.; investigation, G.L., P.D.M. and G.C.; writing original draft preparation, G.L. and P.D.M.; writing review and editing, all Authors; supervision, P.D.M. and F.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** Publicly available dataset were analyzed in this study. Elaborated data are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **On Turbulent Features of E** *×* **B Plasma Motion in the Auroral Topside Ionosphere: Some Results from CSES-01 Satellite**

**Giuseppe Consolini 1,\*, Virgilio Quattrociocchi 1,2, Simone Benella 1, Paola De Michelis 3, Tommaso Alberti 1, Mirko Piersanti <sup>2</sup> and Maria Federica Marcucci <sup>1</sup>**


**Abstract:** The recent Chinese Seismo-Electromagnetic Satellite (CSES-01) provides a good opportunity to investigate some features of plasma properties and its motion in the topside ionosphere. Using simultaneous measurements from the electric field detector and the magnetometers onboard CSES-01, we investigate some properties of the plasma **E** × **B** drift velocity for a case study during a crossing of the Southern auroral region in the topside ionosphere. In detail, we analyze the spectral and scaling features of the plasma drift velocity and provide evidence of the turbulent character of the **E** × **B** drift. Our results provide an evidence of the occurrence of 2D **E** × **B** intermittent convective turbulence for the plasma motion in the topside ionospheric F2 auroral region at scales from tens of meters to tens of kilometers. The intermittent character of the observed turbulence suggests that the macro-scale intermittent structure is isomorphic with a quasi-1D fractal structure, as happens, for example, in the case of a filamentary or thin-tube-like structure. Furthermore, in the analyzed range of scales we found that both magnetohydrodynamic and kinetic processes may affect the plasma dynamics at spatial scales below 2 km. The results are discussed and compared with previous results reported in the literature.

**Keywords:** plasma turbulence; auroral ionosphere; **E** × **B** plasma motion

## **1. Introduction**

In the geospace environment, such as the Earth's ionosphere and magnetosphere, plasma dynamics can display multiscale features and chaos. As it occurs in several space and astrophysical plasmas, the chaotic nature of the plasma dynamics in the geospace environment is due to the occurrence of turbulent processes. This happens, for instance, at high-latitude ionospheric regions, where plasma dynamics was shown to be characterized by turbulence as revealed, for example, by the irregular and chaotic character of the electric and magnetic field fluctuations in the ULF (Ultra Low Frequency) and ELF (Extra Low Frequency) spectral ranges [1–3]. In addition, previous studies of the ionospheric electric, magnetic and electron density fluctuations have shown that the fluctuations are characterized by power-law spectral densities, scaling features and non-Gaussian statistics of spatial and temporal increments at small/short scales (see, e.g., [4–9] and references therein). In detail, in the auroral regions, it has been found by Kintner et al. [10] that the origin of electric field broadband spectra and the features of its fluctuations are probably caused by the occurrence of intermittent turbulence [10,11], due to the sporadic fast interactions between localized coherent plasma structures [4,12]. In any case, most of the observed fluctuations of fields and plasma parameters in the high-latitude ionospheric

**Citation:** Consolini, G.;

Quattrociocchi, V.; Benella, S.; De Michelis, P.; Alberti, T.; Piersanti, M.; Marcucci, M.F. On Turbulent Features of E × B Plasma Motion in the Auroral Topside Ionosphere: Some Results from CSES-01 Satellite. *Remote Sens.* **2022**, *14*, 1936. https:// doi.org/10.3390/rs14081936

Academic Editor: Michael E. Gorbunov

Received: 12 February 2022 Accepted: 9 April 2022 Published: 17 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

regions are characterized by very complex multiscale features that might affect the plasma dynamics [4,13,14].

The high variability of the electric and magnetic field fluctuations can affect the plasma dynamics through the term **E** × **B**, which can be strongly influenced by the multiscale and turbulent character of the field fluctuations. The plasma **E** × **B** drift has been widely investigated in the ionospheric equatorial region, where the variations of the vertical plasma velocity drift play a relevant role in generating ionospheric irregularities [15–18]. Conversely, it is believed less relevant in the auroral region, especially at the typical altitudes where auroral emission occurs [19]. This is because at this altitude the presence of neutrals can interfere with the **E** × **B** drift disrupting the coherent plasma motion perpendicular to the magnetic field direction due to the ion-neutral collisions. This mechanism is at the origin of the well-known Hall current [19]. However, at higher altitudes where the neutrals' density is lower, the **E** × **B** drift may play a relevant role in the dynamics of the plasma also in the auroral regions, and can be relevant to the formation of plasma irregularities.

Recently, the Chinese Seismo-Electromagnetic Satellite (CSES-01) [20,21], launched on February 2018 and equipped with a large set of instruments, including an electric field detector (EFD) [22,23], a flux gate (HPM) and a search-coil (SCM) [24] magnetometer, has had the opportunity to explore the Southern auroral oval region during a time interval characterized by the occurrence of a magnetospheric substorm. During the crossing of the southern auroral region, the satellite has measured simultaneously the magnetic and electric field fluctuations [25]. CSES-01, which has been originally designed to study possible correlations between seismic events and iono/magnetospheric perturbations, flies in the topside F2 ionosphere at ∼500 km of altitude on a Sun-synchronous orbit. Thus, it can also provide the opportunity to explore the plasma dynamics at an altitude higher than that at which the typical auroral phenomena occur. Generally CSES-01 is operative at geographic latitude |*Lat*| < 65◦, and thus the observation of high-latitude phenomena is sporadic and limited to a few periods in which the Earth's dipole tilt angle is higher than |15◦|, corresponding to solstices. On 11 August 2018 from 21:41 UT to 21:45 UT, CSES-01 was in operational mode during a crossing of the Southern polar F2 topside ionosphere. Its measurements allowed us [25] to investigate some interesting features of the electric field fluctuations in the auroral region and to establish the very complex structure of electric field fluctuations probably due to the occurrence of intermittent turbulence [4,12].

This study sits out to investigate the spectral and scaling features of the plasma **E** × **B** drift velocity using magnetic and electric field measurements from the HPM, SCM and EFD experiments onboard CSES-01 during the crossing of the Southern auroral region occurred on 11 August 2018. In detail, this case study seeks to examine the multiscale character of the plasma drift velocity **E** × **B**.

#### **2. Data Description**

We consider magnetic and electric field data recorded by HPM, SCM and EFD instruments [22,24] onboard CSES-01 during a crossing of the Southern polar ionosphere on 11 August 2018 from 21:42:30 UT to 21:45:00 UT. During this time interval CSES-01 partially crosses the Southern auroral oval region in its ascending orbit near 02:00 MLT when auroral precipitation is going on. This time interval is, indeed, characterized by a moderate/high auroral substorm activity, as indicated by the value of the auroral electrojet index, AE > 500 nT [25]. This period has already been investigated in previous papers [25,26], to which the reader may refer for more information.

We use magnetic field data with a resolution of 1 sample/s from HPM and with a resolution of 10 k samples/s from SCM and electric field measurements with a resolution of 5000 samples/s. Magnetic field data from HPM and SCM are joined and successively reduced at the same resolution of EFD, i.e., 5000 samples/s, in order to compute the plasma **E** × **B** drift velocity, i.e.,

$$\mathbf{v}\_D = \frac{\mathbf{E} \times \mathbf{B}}{B^2},\tag{1}$$

where **E** and **B** are the electric and magnetic field, respectively. Data are in the geographic coordinate system, [*x*, *y*, *z*], where the *X*-component is pointing **to** the geographic North, the *Y*-component **to** East and the *Z*-component to Earth's center. Furthermore, the electric field components are corrected subtracting the term **E***<sup>s</sup>* = **v***<sup>s</sup>* × **B**, where **v***<sup>s</sup>* is the satellite velocity, which is due to the satellite motion in the Earth's magnetic field. Lastly, the velocity field is downsampled to the resolution of 250 samples/s by averaging the original 5000 samples/s data on consecutive boxes of 20 points. Thus, considering the satellite velocity, which is of the order of *vs* 8 km/s, and the chosen resolution of 250 samples/s, we are capable of exploring velocity fluctuations up to scales of the order of ∼30 m. This length scale is well below the ion inertial length scale, which is of the order of 2 <sup>÷</sup> 3 km for *<sup>O</sup>*<sup>+</sup> and 10 ÷ 20 m for *e*−. This means that we are capable of exploring also the sub-ionic scales.

Figure 1 shows the electric and the magnetic field data collected by CSES-01 and used in our analysis. The values of the electric field along the three different directions are similar, ranging between −0.3 V/m and 0.15 V/m. Conversely, the magnetic field is stronger along the *z*-axis as expected at high latitude being the magnetic field nearly vertical.

**Figure 1.** Electric (solid lines on the left vertical axes) and magnetic (dashed lines on the right vertical axes) field measurements collected by CSES-01 EFD and HPM instruments for the time interval under consideration. The three components are in the geographic coordinate system.

Figure 2 reports the plasma drift velocity, **v***D*, computed using Equation (1). Due to the quasi-dipole magnetic field configuration, the velocity field of the plasma in the crossed polar region is mainly horizontal in the [*x*, *y*]-plane. Although the mean flow direction is along the *Y* component, the plasma drift velocity shows large fluctuations along the *x* and *y* axes, as confirmed by the values of the root mean square along the three axes **v***RMS <sup>D</sup>* ≡ [0.83, 1.34, 0.35] km/s.

Anyway, the geographic coordinate system is not the optimal reference system for our analysis. Usually, the analysis of the fluctuation field is done in the parallel and perpendicular directions to the magnetic field. For this reason, we change coordinate system and compute the parallel and perpendicular components to the local magnetic field direction of the plasma **E** × **B** drift velocity. The component of the drift velocity, which is parallel (*<sup>v</sup>* ) to the magnetic field direction, is practically zero being <sup>|</sup> *<sup>v</sup>* /*v*<sup>⊥</sup> | <sup>10</sup>−<sup>6</sup> (we remark that it is expected to be zero by construction). The two perpendicular components, which describe the horizontal plasma motion, are reported in Figure 3. They are chosen so that *v* (1) <sup>⊥</sup> and *<sup>v</sup>* (2) <sup>⊥</sup> are almost along the direction of North (*x*) and East (−*y*) axis, respectively. In detail, for each point we consider the local magnetic field direction (i.e., the geomagnetic field direction that is essentially oriented in the *z*-direction) and define two unitary vectors

in the plane perpendicular to geomagnetic field direction, with components mainly along *x*- and *y*-directions. This defines the local reference frame. Successively, we compute the components of the **E** × **B** along these two perpendicular unitary vectors.

**Figure 2.** Drift velocity (**v***<sup>D</sup>* <sup>=</sup> (**<sup>E</sup>** <sup>×</sup> **<sup>B</sup>**)/*B*2) in the geographic coordinate system.

**Figure 3.** The two components of drift velocity (**v***<sup>D</sup>* <sup>=</sup> (**<sup>E</sup>** <sup>×</sup> **<sup>B</sup>**)/*B*2) perpendicular to the local magnetic field direction.

The effective dimension *d*∗ of the drift velocity fluctuation field can be measured by applying the covariance analysis, which allows us to estimate the eigenvalues of the covariance matrix of the drift velocity. Indeed, knowing the eigenvalues of covariance matrix, we can compute the effective dimension, *d*∗, of the fluctuating field as

$$d^\* = \frac{1}{\lambda\_{\text{max}}} \sum\_{i} \lambda\_i \tag{2}$$

where *λ<sup>i</sup>* are the eigenvalues and *λmax* is the maximum eigenvalue. In our case, we get an effective dimension *d*<sup>∗</sup> 1.35. This quantity provides an indication of the dimension of the phase-space spanned by the field fluctuations. In our case, being 1 < *d*∗ < 2, we expect that the drift velocity fluctuation field can be confined in a 2D space. Consequently, we could assume to deal with a 2D plasma motion. Anyway, we remark that the effective dimension *d*∗ is not to be confused with the fractal dimension associated with the plasma dynamics under the drift velocity.

Before proceeding in the analysis, we check the consistency of the obtained average velocity field with the macroscopic overall convection pattern, that can be reconstructed from measurements of the Super Dual Auroral Radar Network (SuperDARN) radars operating in the Southern polar hemisphere. SuperDARN is an international scientific radar network consisting of high frequency (HF) coherent scatter radars that is operated under international cooperation. The SuperDARN data are mainly used to monitor the dynamics of the ionosphere and upper atmosphere in the high- and mid-latitude regions through the production of global plasma convection maps. However, these data can also be used to study many others geospace phenomena including field aligned currents, geomagnetic storms, magnetospheric substorms, magnetic reconnection, and interhemispheric plasma convection asymmetries (see, e.g., [27–29] and references therein). Figure 4 reports a comparison between the ionospheric convection map as reconstructed from SuperDARN observations relative to the time interval 21:44-21:46 UT, corresponding to the central part CSES-01 track and the convection drift velocity obtained by **E** × **B** from CSES-01 measurements.

Global maps of ionospheric convection can be derived from the SuperDARN data using the technique developed by Ruohoniemi and Baker [30] and Shepherd and Ruohoniemi [31]. In this technique, the SuperDARN line-of-sight velocities are fitted to an expansion of the high-latitude electrostatic potential in spherical harmonic functions. In order to constraint the solution in regions where no SuperDARN data are available, additional data from statistical models are used. The convection map shown in Figure 4 has been computed using the statistical model by Thomas and Shepherd [32] for adding data where SuperDARN data are missing. This statistical model (named TS18) has been derived for the Southern and Northern hemispheres using measurements from all the mid-latitude, high-latitude, and polar radars available in the years from 2010 to 2016. The TS18 is capable of producing climatological patterns of ionospheric convection organized by solar wind, interplanetary magnetic field (IMF), and dipole tilt angle conditions. The TS18 model inputs used for the computation of the map in Figure 4 correspond the average values of IMF clock angle (223◦), solar wind electric field (*Esw* = 2.2 mV/m), and dipole tilt angle (−17.7◦) for the time interval of interest (11 August 2018 from 21:42:30 UT to 21:45:00 UT). In Figure 4, the colored vectors represent the fitted SuperDARN velocity vectors. The convection map has been generated by using the Radar Software Toolkit 4.2 [33]. Lastly, the green arrows are designed to facilitate the reader in locating the direction of the plasma convection motion. On the bottom of the same figure, we report a zoom of the plasma convection pattern, which permits us to identify the CSES-01 trajectory during the analyzed period. The trajectory is the blue line, while the color lines refer to the reconstructed plasma **E** × **B** drift velocity (**v***D*). The length of the lines is proportional to the intensity of the drift velocity while the color (from red to light violet) indicates the time. An overall agreement between the observed convection pattern (green arrow) and the reconstructed **E** × **B** drift velocity (colored lines) is found. This confirms the correct reconstruction of the convective plasma **E** × **B** drift velocity by using CSES-01 measurements. The only small discrepancy is observed in correspondence with the higher latitude part of CSES-01 trajectory, where there is a plasma flow reversal. This can be due to a localized vortex structure near the lower latitude boundary of the convection cells not resolved in the SuperDARN convection map.

In Figure 5, we report the motion of the velocity vector tip in the plane perpendicular to the magnetic field direction. The motion of the vector tip shows a certain degree of chaoticity, which may be an indication of the occurrence of turbulence in plasma motion.

**Figure 4.** On the **top**: the instantaneous convection cells as reconstructed from SuperDARN observations in Antarctica. Green arrows show the overall plasma convection. On the **bottom**: a zoom of the region of the CSES-01 trajectory. Colored lines refer to the velocity vector field in the *XY*-plane. Colors (from red to light violet) indicate the time.

**Figure 5.** Motion of the velocity vector tip in the plane perpendicular to the magnetic field. The color is associated with the universal time UT (see the color bar).

#### **3. Methods**

The analysis of the spectral and scaling features of the plasma **E** × **B** drift is done using the standard Fourier power spectral analysis and the structure function analysis applied in fluid and MHD turbulence studies. While the Fourier analysis is a standard method to detect the presence or not of specific modes in the fluctuation field, the structure function analysis, as clearly discussed in Frisch [34], is one of the most powerful methods to investigate scaling features for fully developed turbulence. This approach is, indeed, a powerful method to detect the occurrence anomalous scaling features, i.e., what is generally called *intermittency* in fluid and MHD turbulence.

In the case of fully developed turbulence, a *q*th-order structure function, *S*(*q*)(*δr*), is defined as the *q*th-order moment of the longitudinal velocity increments, *δvl*, at the spatial scale *δr*, i.e.,

$$S^{\langle q \rangle}(\delta r) = \langle (\upsilon\_l(r + \delta r) - \upsilon\_l(r))^q \rangle \equiv \langle (\delta \upsilon\_l)^q \rangle. \tag{3}$$

Looking at the definition of Equation (3) for a fixed spatial scale *δr*, the *q*th-order structure function represents the moment of order *q* of the distribution of the velocity increments at that scale. In the case of a turbulent flow, this quantity is expected to scale in the inertial range according to the following expression,

$$S^{(q)}\left(\delta r\right) \sim \delta r^{\mathbb{Z}(q)},\tag{4}$$

where the scaling exponents *ζ*(*q*), as stated by the usual Kolmogorov K41 theory of turbulence, are expected to be *ζ*(*q*) = *q*/3 [34]. In particular, following the K41 theory and according to the Navier–Stokes equation for an inviscid fluid (i.e., in the limit of an infinite Reynolds number) the third order structure function (*S*(3)(*δr*)) is expected to scale with *ζ*(3) = 1. From the theory, it is possible to write

$$S^{(3)}(\\\delta r) = -\frac{4}{5}\epsilon\delta r.\tag{5}$$

where is the energy transfer rate along the cascade in the inertial range. Conversely, the *ζ*(*q*) = *q*/3 dependence of the structure function scaling exponents is essentially a conjecture that implies that a global scale-invariance should exist. However, the observed scaling exponents can show departures from the linear dependence due to intermittency

effects. The observed dependence of scaling exponents on the moment order is well described by a convex function of *q*, so that the global scale invariance is missing. By the way, in the case of intermittency effects, Equation (5) is also expected to hold.

Sometimes, the previous definition of the structure function is generalized considering the absolute value of the velocity differences, *δvl*, i.e.,

$$S^{(q)}(\delta r) = \langle |\delta v\_l(\delta r)|^q \rangle,\tag{6}$$

which allows both a better analysis of the dependence of scaling exponents on moment order (see, e.g., [35,36]) and the investigation of non-integer *q*. We remark that, although structure function analysis is generally defined for increments along the flow direction, it has also been generalized to the case of transverse directions, as well as, for time increments, assuming Taylor's hypothesis to be valid [37,38].

A simple way to investigate the emergence of intermittency or the validity of the existence of a global scale invariance is to compute the relative scaling of the *q*th-order structure function on the *p*th-order one, which is expected to scale as

$$S^{(q)}(
\delta r) = \left[ S^{(p)}(
\delta r) \right]^{\gamma\_p(q)},\tag{7}$$

where *γp*(*q*) = *q*/*p* in the case of a global scale invariance as the one predicted by K41 theory. We can note how *γ*(*q*) ≡ *ζ*(*q*) if Equation (5) holds and *p* = 3 [35,39]. This method has been introduced by Benzi et al. [35] and it is named as *extended self-similarity* (ESS), which is also valid in the case of low Reynolds number turbulence.

Another interesting quantity is the so-called *generalized kurtosis* Γ*q*(*δr*) [37], which is defined as

$$\Gamma\_q(\delta r) \equiv \frac{\mathcal{S}^{(q)}(\delta r)}{\left[\mathcal{S}^{(2)}(\delta r)\right]^{q/2}}.\tag{8}$$

It is expected to be constant in the inertial range when global scale invariance holds. Conversely, departures from a constant value are generally referred as the evidence for the occurrence of *intermittency*. We remark that the emergence of intermittency manifests also in the lack of a scale-invariant shape of the probability density functions (PDFs) of the velocity increments *δv*(*δr*) with the scale *δr*.

We apply the described analysis to the drift velocity field.

#### **4. Results**

We start our study of the properties of the plasma **E** × **B** drift velocity by investigating its spectral features. Figure 6 reports the trace, *S*(*f*), of the power spectral density (PSD) of the components of the drift velocity field, which are perpendicular to the local magnetic field direction. The trace of PSD is defined as the sum of the PSDs of the two perpendicular components. The observed shape of PSD is quite consistent with the spectra expected in the case of two dimensional (2D) **E** × **B** convective turbulence in a quasi-steady state [40] or two-dimensional magnetohydrodynamic (MHD) turbulence in absence of Alfvén effect [41]. Indeed, for example, numerical simulations by Gruzinov et al. [40] revealed that when ionospheric plasma density irregularities break up into fingers leading to the formation of stable shear flows, the omnidirectional energy spectrum scaling exponent approaches to *α* −2 after reaching a quasi-steady state rather than *α* −5/3 as expected from the K41 theory prediction in the case of isotropic turbulence. In both cases [40,41], the spectral scaling exponents are steeper than the K41 theory prediction, as also occurs in our case being *α* −1.7. Consistent with the literature, the observed PSD of the components of the plasma drift velocity field seems to support the idea of the occurrence of a strong 2D turbulence in a low−*β* plasma medium [42], where *β*, which is the ratio of plasma to magnetic pressure, is a parameter classifying plasma conditions being usually *β* 1 for the Earth's ionospheric plasma.

Let us now move to the analysis of the scaling features by the generalized structure functions, *S*(*q*)(*τ*),

$$S^{(q)}(\tau) = \langle |v\_i(t+\tau) - v\_i(t)|^q \rangle,\tag{9}$$

where *τ* is the timescale, i.e., the time delay used to compute the velocity increments. The analysis is done in the plane perpendicular to the geomagnetic field local direction by considering separately the two components and then averaging the obtained structure functions. The averaging procedure is justified by the fact that the observed scaling features in the two separate perpendicular directions are isotropic, i.e., we do not observe significant differences. The timescale can be associated with a spatial scale *δr* by considering the satellite velocity *vs*, i.e., *δr* = *τvs*. Here, we investigate the timescale in the range from 4 ms to 3 s, approximately corresponding to a range of spatial scale *δr* ∈ [0.03, 23.50] km. Clearly, this assumption, which is analogous to the Taylor's hypothesis, is strictly valid if the transit time is faster than the evolution time of the velocity field at the investigated time scale, an assumption that is nearly valid in the low-frequency range, that is, where we can assume that the temporal fluctuations are principally due to Doppler-shifted and stationary spatial variations (see, e.g., [11,25,43] and references therein). We will return on this point in Section 5.

**Figure 6.** The trace, *S*(*f*), of the PSD of the perpendicular components of the **E** × **B** drift velocity. The two power laws are that expected for K41 theory (∼*<sup>f</sup>* <sup>−</sup>5/3) and that observed for 2D **<sup>E</sup>** <sup>×</sup> **<sup>B</sup>** convective turbulence (∼*<sup>f</sup>* <sup>−</sup>2) in a quasi-steady state [40].

In the case of scale-invariant signals, the generalized *q*-order structure functions, *S*(*q*)(*τ*), are expected to depend on the scale *τ* according to a power law, i.e.,

$$S^{(q)}(\tau) \simeq \tau^{\gamma(q)},\tag{10}$$

where *γ*(*q*) are the scaling exponents, which for a signal characterized by a global scale invariance are expected to depend linearly on the moment order *q*, i.e., *γ*(*q*) = *αq*. In particular, the K41 theory of homogeneous turbulence predicts *γ*(*q*) = *q*/3.

Figure 7 shows the generalized structure functions, *S*(*q*)(*τ*), of the velocity field components perpendicular to the magnetic field direction. Good scaling regimes are observed. In particular, the scaling of the structure functions is observed over near two orders of magnitude from *τ* ∼ 0.04 s to *τ* ∼ 4 s. Furthermore, the third-order structure function, *S*(3)(*τ*), displays a range of time scales where a reasonable agreement with the expected linear scaling predicted by the K41 fluid turbulence theory can be found. This is very clear by plotting the compensated structure functions, *τ*−1*S*(*q*)(*τ*) for *q* = 2, 3 and 4 as shown

in Figure 8. This result suggests that to compute the scaling exponents of the structure functions, we can apply the ESS method introduced by Benzi et al. [35], which is based on the investigation of the relative scaling of the structure functions as reported in Equation (7).

**Figure 7.** The generalized structure functions, *S*(*q*)(*τ*), of the velocity field components perpendicular to the magnetic field direction. The black line refers to a linear scaling, i.e., *<sup>S</sup>*(3)(*τ*) *<sup>τ</sup>*.

**Figure 8.** The compensated generalized structure functions, *τ*−1*S*(*q*)(*τ*), of the velocity field components perpendicular to the magnetic field direction for *q* = 2, 3 and 4. The black lines refer to power law fits.

Figure 9 shows relative scaling of the *q*th-order structure functions, *S*(*q*)(*τ*), on the 3rd-order one for a selected number of moment orders *q*. For all the moment orders, the *q*th-order structure functions show a power-law dependence on the corresponding 3rdorder one in agreement with Equation (7). In other terms, ESS is a property of the observed velocity fluctuations perpendicular to the magnetic field direction.

**Figure 9.** The relative scaling of *q*th-order structure functions, *S*(*q*)(*τ*), of the velocity field components perpendicular to the magnetic field direction versus the corresponding 3rd-order one. The black dashed lines are power-law fits.

The relative scaling exponents *γ*(*q*) are reported as a function of the moment order *q* in Figure 10. A clear departure on the linear scaling is observed, being the dependence of the scaling exponents on moment order of a convex function. The departure from a linear scaling does not support the occurrence of a global scale-invariance, suggesting instead the emergence of anomalous scaling features, i.e., *intermittency*.

**Figure 10.** The relative scaling exponents *γ*(*q*) as a function of the moment order *q*. The dashed line is the expected trend of *γ*(*q*) for the K41 theory of turbulence (*γ*(*q*) = *q*/3 being *ζ*(3) = 1). The blue solid line is a nonlinear best fit done using the Meneveau and Sreenivasan *P-model* [44]. The solid green line is the *γ*(*q*) trend for the She-Leveque model of Equation (12) [45].

To characterize and quantify the degree of intermittency, that is how the scaling exponents *γ*(*q*) deviate from the Kolmogorov prediction, some intermittency models were proposed, which attempt to explain in particular the anomalous scaling exponents. A model capable of explaining the anomalous scaling exponents is the multifractal model by Meneveau and Sreenivasan [44], known as the *P-model*. It is a simple model able to describe the energy-cascading process in the inertial range whose scheme is based on the generalized two-scale Cantor set with equal scales, but unequal weights. Thus, we compare the observed trend of the *γ*(*q*) exponents with the one predicted by *P-model*. In detail, we fit the trend of *γ*(*q*) as a function of moment order *q* using the following expression,

$$\gamma(q) = 1 - \log\_2\left(p^{\frac{q}{5}} + (1 - p)^{\frac{q}{5}}\right),\tag{11}$$

where *p* is the weight of the measure repartition in the multiplicative process of the twoscale Cantor set. The *P-model* is an excellent approximation of the observed trend of *γ*(*q*) and we get *p* = [0.25 ± 0.01]. Being the value of *p* significantly different from one-half (*p* = 0.5) we may surely assert that the drift velocity fluctuation field has a multifractal structure and that the observed energy repartition at the different scales is an anomalous multiplicative process as the one described by the *P-model*.

In the framework of fluid turbulence, other models have been proposed to explain the convex trend of *γ*(*q*). Among these, the She-Leveque (SH) model [45] relates the anomalous scaling with the dimension of the dissipative structure in turbulence. In detail, in the case of fluid turbulence, the SH model predicts for *γ*(*q*) the following behavior:

$$\gamma(q) = \frac{q}{9} + 2\left[1 - \left(\frac{2}{3}\right)^{\frac{q}{9}}\right].\tag{12}$$

The previous expression, which is specialized for the case of 3D fluid turbulence, was generalized later by Politano and Pouquet [46] (see also ref. Biskamp and Müller [47]). The generalized expression of SH model assumes the following form:

$$\gamma(q) = (1-x)\frac{q}{3} + \mathcal{C}\_0 \left[ 1 - \left( 1 - \frac{x}{\mathcal{C}\_0} \right)^{\frac{q}{3}} \right],\tag{13}$$

where *C*<sup>0</sup> is the co-dimension of the intermittent structures and *x* is the scaling exponent of the dynamic timescale associated with the most intermittent structure (see [46,47]). Thus, in the case of 3D fluid turbulence, one gets *C*<sup>0</sup> = 2 and *x* = 2/3.

The SH-model (see Equation (13)) can be also written as,

$$\gamma(q) = \left(1 - \frac{\mathbb{C}\_0}{3}\right)\frac{q}{3} + \mathbb{C}\_0\left(1 - \beta^{\frac{q}{3}}\right) \tag{14}$$

where *β* = 1 − *x*/*C*<sup>0</sup> is a parameter related to the degree of intermittency. Indeed, for *β* → 1, there is no intermittency, while if *β* ≤ 1, intermittency is present. In our case, we get *β* 2/3.

Table 1 shows the values for the first 4 scaling exponents in comparison with Ruiz-Chavarria et al. [48] results from fluid turbulence, with Benzi et al. [49] for 3D convective turbulence, Biskamp and Schwarz [50] simulations on 2D MHD turbulence (in this case data refer to Elsässer variables), and SH 3D fluid model. The observed values are very well in agreement in spite of the different physical scenarios, i.e., **E** × **B** convective 2D plasma motion, 3D fluid turbulence [48], 2D MHD simulations [50] and theoretical models [45]. A very good agreement is found with exponents from fluid and convective turbulence. In particular, the good agreement with the SH model [45] suggests that the co-dimension of the intermittent structures is *C*<sup>0</sup> 2.


**Table 1.** Observed scaling exponents *γ*(*q*) and comparison with results from models and literature.

Data from Ref. [50] refer to scaling exponents of Elsässer variables.

Figure 11 reports the generalized kurtosis Γ4(*τ*) as a function of the timescale *τ*. Γ<sup>4</sup> is greater than 3 for all the considered timescales and exhibits an increasing trend for decreasing values of *τ*. These features provide a strong indication of the non-Gaussian character of the velocity increments (being Γ<sup>4</sup> > 3) that is related to the occurrence of intermittency. Another interesting feature of the trend of Γ<sup>4</sup> is the presence of a local maximum at *<sup>τ</sup>* 0.2 s. This timescale is near the one corresponding to the O<sup>+</sup> ion inertial length, *<sup>η</sup>*, assuming a density *<sup>n</sup>*O<sup>+</sup> ∈ [1, 2] · <sup>10</sup><sup>5</sup> cm−<sup>3</sup> (this range of values for the oxygen ion density *n*O<sup>+</sup> is estimated from electron density *ne*<sup>−</sup> [51] assuming quasi-neutrality and that oxygen ions are the more relevant ions at CSES-01 altitude). Indeed, if we consider the Taylor's hypothesis valid, we have *τO*<sup>+</sup> *<sup>η</sup>* = *η*/*vs* 0.3 s being *vs* ∼ 7.8 km/s the satellite speed.

We compute the probability density functions (PDFs) of the velocity increments in the time scale interval from 4 ms to ∼0.5 s. Figure 12 shows the evolution of the standard deviation normalized PDF with the timescale. PDF collapsing does not occur since the shape of the PDFs evolves across different timescale *τ*, supporting the emergence of intermittency in the velocity field.

**Figure 11.** The generalized kurtosis Γ4(*τ*). The solid curve is a guide for eye. The vertical dashed line indicates the expected timescale *τO*<sup>+</sup> *<sup>η</sup>* 0.3 s corresponding to the O<sup>+</sup> inertial length, *<sup>η</sup>*, assuming a density in the range [1, 2] <sup>×</sup> <sup>10</sup><sup>5</sup> cm<sup>−</sup>3.

**Figure 12.** The PDFs of the velocity increments in the time scale interval from 4 ms to ∼0.5 s. Data are rescaled by the corresponding standard deviation.

We remark that the PDFs are not Gaussian being, indeed, characterized by a leptokurtic shape, which can be well fitted by the following function:

$$p(\mathbf{x}) = \frac{N\_0}{\left[1 + \frac{1}{\kappa} \left(\frac{\mathbf{x}}{\mathbf{x}\_0}\right)^2\right]^\kappa} \exp\left(-\left|\frac{\mathbf{x}}{\mathbf{x}\_c}\right|\right),\tag{15}$$

where *x* = *δv*(*τ*)/*στ*, *x*<sup>0</sup> and *xc* are characteristic scales, *N*<sup>0</sup> is a normalization factor and *κ* is the exponent governing the tail behavior. In the case of the smallest-scale PDF at *τ* = 4 ms we get *x*<sup>0</sup> = [0.130 ± 0.003], *xc* = [2.06 ± 0.03] and *κ* = [0.6 ± 0.1].

To better characterize the evolution of the PDFs as a function of the timescale, we estimate the Kullback–Leibler (*KL*) distance between the PDFs. In detail, we compute the following quantity:

$$KL(\mathbf{r} \mid \mathbf{r}\_0) = \int\_{-5}^{+5} p(\mathbf{x}, \mathbf{r}) \log\_2 \frac{p(\mathbf{x}, \mathbf{r})}{p(\mathbf{x}, \mathbf{r}\_0)} d\mathbf{x} \tag{16}$$

where *x* = *δv*/*σ* and *τ* and *τ*<sup>0</sup> are the timescales of the two considered PDFs. Here, we set *τ*<sup>0</sup> = 4 ms, i.e., the minimum available timescale.

To quantify the significance of the measured *KL*-distance we compute a critical threshold via Monte Carlo simulation. In detail, we generate a series of 1000 random samples *xi* containing the same number of values of the actual time series and following the PDF *p*(*x*, *τ*0) of Equation (15). Successively, we compute the reciprocal *KL*-distance between the distributions of this set of 1000 random samples and evaluate the cumulative distribution of the values of *KL*-distance. The value corresponding to the 95% level of the cumulative distribution of the *KL*-distance is set as the critical threshold to discriminate the significance of the measured *KL*-distance. This threshold value is *KL*<sup>∗</sup> = <sup>1</sup> × <sup>10</sup><sup>−</sup>3.

Figure 13 shows a clear increasing of *KL*-distance with the time scale well above the 95% significance threshold, which is the evidence for the absence of PDF-collapsing. In other words, the PDFs evolve with the timescale of increments so that there is not an invariant shape of the PDFs. This is another counterpart of the absence of a global scale-invariance of increments at different timescales, i.e., of the emergence of intermittency. Furthermore, we observe a short intermediate range of timescales from 0.1 s to 0.6 s where

there is a plateau in *KL* distance. This plateau is well in agreement with the range of timescales *τO*<sup>+</sup> *<sup>η</sup>* associated with the *O*<sup>+</sup> inertial length *η*.

**Figure 13.** The Kullback–Leibler (*KL*) distance between the PDFs. The vertical dashed line is in correspondence of the timescale *τO*<sup>+</sup> *<sup>η</sup>* <sup>∼</sup> 0.35 s associated with the *<sup>O</sup>*<sup>+</sup> inertial length *<sup>η</sup>* for a density of <sup>∼</sup>105 cm<sup>−</sup>3. The solid horizontal line indicate the 95% critical threshold, *KL*<sup>∗</sup> <sup>=</sup> <sup>1</sup>×10<sup>−</sup>3, below which the *KL*-distance between PDFs is not significant (grey region).

#### **5. Discussion**

In this study we provided a first analysis of the scaling features of plasma **E** × **B** drift velocity in the auroral regions using data collected by the CSES-01 for a case study during a crossing of the Southern polar ionosphere. In detail, our study investigates the drift velocity fluctuations in a range of timescales *τ* ∈ [0.004, 3] s, approximately corresponding to a range of spatial scale *δr* ∈ [0.03, 23.50] km, assuming that Taylor's hypothesis is valid.

The reconstructed 3D convective velocity field is anisotropic, the vertical velocity being very small. This is due to the geomagnetic field configuration in the polar regions, which is mainly vertical. We found a general good consistency of the reconstructed convective velocity field and the overall large scale plasma convection as observed by SuperDARN observations.

Regarding the spectral features, we found evidence of the anisotropic character of the fluctuation field and of the existence of a large interval of frequencies (wave-numbers via Taylor's hypothesis) where the spectral features follow a power law. The observed spectral exponent of the PSDs in the plane perpendicular to the main geomagnetic field is well in agreement with that expected in the case of 2D **E** × **B** convective turbulence (|*α*| - 2) in a quasi-steady state [40]. Furthermore, the observed spectral features are also in agreement with the occurrence of 2D MHD turbulence and in strong low-*β* 2D plasma turbulence. Indeed, low-*β* plasma, due to the strong mean magnetic field the plasma motion, is mainly confined in the perpendicular plane, so that the turbulence is essentially 2D. The presence of a strong mean field prevents the bending of magnetic field lines in the mean field direction [50]. Similar spectral features have been found by Kelley and Kintner [52] and Cerisier et al. [53] in the high-latitude ionosphere for the electric field fluctuations, which are strongly correlated to the **E** × **B** velocity in that region, and in inertial range of strong turbulence in low-*β* plasmas [42].

We remark that the previous considerations on the spectral properties are based on the assumption of the validity of the Taylor's hypothesis, which is assumed to be true in the investigated frequency range. Indeed, according to Tchen et al. [42], the conversion of Eulerian observations into Lagrangian ones using the Taylor's hypothesis is strictly valid only if the flow velocity is very high, or in the case of frozen turbulence, i.e., in the lack of a nonlinear dispersion relation. In other words, the assumption of a linear relationship between *f* and *k* is restricted to those cases in which the plasma drift velocity is very high or, if this is not the case, for large wave numbers, as it occurs in the inertial range.

From a general point of view, a difference between 3D and 2D MHD turbulence is the occurrence of an inverse magnetic helicity cascade, which is capable of generating self-organization, i.e., large-scale coherent structures [41]. This scenario seems to be very well in agreement with the occurrence of large scale plasma motion, which manifests in the formation of the two macroscopic convection cells observed in the polar ionosphere (see Figure 4 top panel). The convective motion is clearly forced by the magnetic reconnection phenomenon at the magnetopause, which enhances the plasma convection in the magnetosphere–ionosphere system. In other words, it seems that the nature of the observed turbulence could be isomorphic to a forced convective turbulence.

The analysis of the scaling features of the velocity field increments evidences the occurrence of anomalous scaling properties, supporting the intermittent character of the velocity fluctuation field. The emergence of intermittency is also supported by generalized kurtosis Γ4, which shows an increasing trend for decreasing timescale *τ*, and by the absence of collapsing of the PDFs of the velocity increments at different timescales as shown in Figure 12 and by the evolution of *KL*-distance with the timescales, as shown in Figure 13.

In spite of the different character of the observed turbulence with respect to the ordinary 3D fluid turbulence, the observed anomalous scaling properties are consistent with what is found in the case of the 3D fluid turbulence [48] and theoretical derivations [44,45]. This is a very interesting result that requires a deeper theoretical investigation, but that in any case seems to support the similarity of 2D and 3D plasma turbulence [41], although we have to remark that in the case of 2D MHD turbulence, a more intermittent character is expected with respect to 3D MHD turbulence [50]. Furthermore, the comparison with SH model suggests that the co-dimension of the intermittent structures is *C*<sup>0</sup> 2 and *β* 2/3.

Being *C*<sup>0</sup> 2, in a first approximation considering the case of a 3D system, the fractal dimension of the intermittent structures is *dF* = *D* − *C*<sup>0</sup> = 1 (where *D* = 3 is the system dimension), which suggests that the *macro-scale* structure is isomorphic with a quasi-1D fractal structure, i.e., a flux tube or a line. This is quite well in agreement with what can be observed in Figure 4 when the satellite enters in the convection region, where the structure of the **E** × **B** drift velocity resembles a quasi-cylindrical configuration, i.e., a structure of co-dimension *C*<sup>0</sup> = 2 and a dimension 1.

Another interesting result from our analysis stands in the behavior of the generalized kurtosis Γ<sup>4</sup> and the *KL*-distance as a function of timescale *τ*. As expected from the occurrence of intermittency and anomalous scaling properties we observe an increase of the generalized kurtosis Γ<sup>4</sup> for decreasing timescales and a departure of the shape of the PDFs from the one at the smallest investigated timescale *τ*<sup>0</sup> = 4 ms. However, around the typical timescale associated with the oxygen inertial length we observe a maximum in the generalized kurtosis Γ<sup>4</sup> and a plateau in the increasing trend of the *KL*-distance. This result suggests that in the investigated range of scales, the turbulence associated with **E** × **B** convective motion could be more complex than the standard single fluid 2D MHD turbulence. In other words, the multi-ion character of ionospheric plasma could manifest in contiguous multiple regimes of different physical processes, both MHD and kinetic, in the range of timescales from *τ* = 5 ms to *τ* = 1s [13,14].

#### **6. Summary and Conclusions**

In this work, we have investigated the scaling features of the **E** × **B** drift motion in the topside F2 polar ionosphere for a case study relative to a crossing of the Southern auroral region by the Chinese Seismo-Electromagnetic satellite (CSES-01) during a geomagnetically disturbed period.

We have found a clear evidence for the occurrence of a turbulent character of the **E** × **B** drift velocity in the analyzed region. The observed turbulence seems to be isomorphic to the 2D **E** × **B** convective turbulence [40]. Furthermore, the drift velocity fluctuations display an intermittent character as evidenced by the scaling analysis and the behavior of the generalized kurtosis Γ4(*τ*). The observed intermittency features of the **E** × **B** drift velocity suggests the hypothesis that intermittency can be related to a macro-scale structure, which is consistent with a quasi-1D fractal structure (thin-flux-like or filamentary structures).

Another interesting result is the emergence of the occurrence of both MHD and kinetic processes in the range of investigated timescales, which can affect the dynamics of the plasma motion at the smallest spatial scales, i.e., below 2 km. This is clearly shown by the trend of the generalized kurtosis Γ<sup>4</sup> and the Kulback–Leibler distance, which show some deviations from the general trend at a scale consistent with the O<sup>+</sup> inertial length. Our findings support the view of very complex dynamics and the idea that it is necessary to consider the multi-species character of the ionospheric plasma medium in order to properly understand and correctly model the small scale ionospheric dynamics, not being sufficient the use of a single fluid theory as is the MHD theory [13,14].

Clearly, further research is required to disentangle and characterize the features of the MHD and kinetic processes involved at the smallest spatial scales. A better comprehension of the observed turbulence processes may benefit from the use of very high-resolution measurements such as those available from the future NanoMagSat mission [54].

**Author Contributions:** Conceptualization, G.C. and V.Q.; methodology, G.C. and T.A.; investigation, V.Q., G.C., S.B., T.A. and P.D.M.; data curation, M.P. and M.F.M.; writing—original draft preparation, V.Q., G.C. and P.D.M.; writing—review and editing, all. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is supported by the contract ASI LIMADOU Scienza+ n. 2020-31-HH.0 and Italian PNRA under contract PNRA18 00289-A Space weather in Polar Ionosphere: the Role of Turbulence

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** This work made use of the data from the CSES mission (http://www. leos.ac.cn/, accessed on 1 October 2020), a project funded by China National Space Administration and China Earthquake Administration in collaboration with Italian Space Agency and Istituto Nazionale di Fisica Nucleare.

**Acknowledgments:** This work is in the framework of the CSES-LIMADOU Collaboration (http: //cses.roma2.infn.it). We acknowledge the Italian Space Agency (ASI) and the Italian National Project for Anctartic Research (PNRA) for supporting this work in the framework of contract ASI "LIMADOU Scienza+" n◦ 2020-31-HH.0 and PNRA18 00289-A "Space weather in Polar Ionosphere: the Role of Turbulence". V. Quattrociocchi and M. Piersanti thank the Italian Space Agency for the financial support under the contract "LIMADOU-2 fase B2/C/D/E1". The authors acknowledge the use of SuperDARN data. SuperDARN is a collection of radars funded by national scientific funding agencies of Australia, Canada, China, France, Italy, Japan, Norway, South Africa, United Kingdom and the United States of America.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **References**

