On the contrary, UofA computed the correlation of 2nd NEC magnetic field component (East-West component) of the magnetic residuals obtained by removing the provided core, external and lithospheric fields from the CHAOS model. They performed the analysis in a 4 min sliding window. The maximum correlation is achieved by first calculating the along-track latitude difference between Swarm A and C, and then fine-tuning this time lag corresponding to optimal correlation. Therefore, UofA provided a time series of the correlation coefficient for each of the events. In order to obtain an average correlation coefficient also for the analysis from UofA, we averaged the correlation time series provided by UofA in the auroral time intervals provided by ISS.

The values of the correlation parameters provided by ISS-JUB, based on filtered and unfiltered data, and the values of the averaged correlation parameter from UofA, for all the events of FAC-MICE are shown in Fig. 8.1. In this graph, we also reported the average values of the AE index, downloaded from NASA cdaweb (https://cdaweb. sci.gsfc.nasa.gov/index.html/). Along the X axis is reported the number of the event, as listed in Table 8.2, identified through the auroral automatic identification procedure from ISS-JUB. The number of events is now 36, because some of the intervals in the FAC-MICE list (Table 8.1) contained two auroral crossings, corresponding to the ascending and descending branches of the orbit. In Fig. 8.1, we can note that the values of the three correlation coefficients for each event are quite similar to each other, and no relation with the geomagnetic activity is observed.

**Fig. 8.1** FAC quality parameters based on cross-correlation analysis of magnetic field data measured by Swarm A and C, for the 36 events of the FAC-MICE dataset listed in Table 8.2. In this figure, also reported, is the average value of the AE index for each of the events

#### *8.4.1 Comparison for the 'High Correlation' Events*

As a first event for the comparison of various FAC estimates, let's consider the event # 15, observed on 2014-07-11 in the interval 19:30–19:45 UT, characterized by very high correlation coefficients and low AE index. This is a conjunction event, when all the three Swarm spacecraft are close together during the auroral crossing, and it was observed in the northern hemisphere at 18 MLT.

The magnetic field residuals from the standard magnetic field models (AUX\_CORE, Lithosphere-MF7, magnetosphere-POMME-6) detected by the three Swarm spacecraft are illustrated in Fig. 8.2, in NEC reference frame. We can see that the magnetic signatures observed by the three spacecraft appear very similar, and are nearly simultaneous, apart for the time delay among Swarm A and C, which is approximately 8 s for this event.

Figure 8.3 shows the comparison among some of the FAC products delivered to ESA for FAC-MICE, for the event # 15. In panel (a) are displayed the single and dual S/C products provided by GFZ, which are equivalent to ESA Level 2 FAC products. The filtered single S/C products, the three S/C products and the unfiltered dual S/C products provided by ISS-JUB are displayed in panels (b), (c), (d) respectively. In panel (e) the products from FMI based on SECS method, in panel (f) the multi-point estimates from RAL, and in panel (g) the single S/C products from UnCA. In all these panels, the dual S/C product from GFZ is also shown for comparison.

The single and dual FAC estimates provided by GFZ, which are equivalent to ESA Level 2 FAC products, are displayed in panel a. The single spacecraft FAC product

**Table 8.2** The FAC-MICE auroral events, as identified through the auroral automatic identification procedure from ISS-JUB. The number of events is now 36, because some of the intervals in the FAC-MICE list (Table 8.1) contain two auroral crossings, corresponding to the ascending and descending branches of the orbit

**Fig. 8.2** The magnetic field residuals detected by the three Swarm spacecraft, in NEC reference frame, for the auroral event on 2014-07-11, 19:30–19:44 UT. Note that this is a conjunction event, when the three spacecraft are close together during the auroral crossing

obtained from Swarm C has been delayed by 7.9 s to take into account the spacecraft separation along the orbit.

The single spacecraft FAC products show several high frequency fluctuations, at time scales of few seconds, which sometimes are similar among the three spacecraft.

The dual spacecraft FAC product appears much smoother with respect to the single spacecraft FAC products. This difference is related to the low-pass band filter adopted for the dual spacecraft FAC product, at 20 s time scale (Ritter et al. 2013). This filter is necessary in the dual spacecraft FAC product to meet the time stationarity assumption, suppressing the FACs carried by kinetic Alfven waves and also to filter out FAC structures smaller than the separation among the 4 point measurements used for FAC calculation (quad size), which is approximately 50 km. Both these conditions are equally important.

In order to better compare the single and dual spacecraft FAC products, ISS-JUB provided filtered single spacecraft FAC products, obtained using the same the low-pass band filter adopted for the dual FAC product, which are shown in Fig. 8.3 panel (b). The agreement of the filtered single spacecraft FAC products with the GFZ dual spacecraft FAC product is now remarkable. This confirms that the current sheet is time stationary and invariant in longitude during this event since all the three spacecraft, at different times, recover the same FAC structure.

ISS-JUB also provided a FAC product based on data measured from all the three Swarm spacecraft, which can be obtained only during the conjunction events. Since this three-spacecraft ISS-JUB FAC product is based on simultaneous multipoint

**Fig. 8.3** Some of FAC products delivered to ESA for FAC-MICE, for the auroral event observed on 2014-07-11, 19:30–19:44 UT. Panel **a** shows the single and dual s/c products provided by GFZ, which are equivalent to ESA Level 2 FAC products. The filtered single S/C products, the three S/C products and the unfiltered dual S/C products provided by ISS-JUB are displayed in panels **b**, **c**, **d** respectively. In panel **e** the products from FMI based on SECS method, in panel **f** the multi-point estimates from RAL, and in panel **g** the single S/C products from UnCA. In all these panels, the dual S/C product from GFZ is shown for comparison

measurements, it does not require the filtering in order to satisfy the time-stationarity assumption. However, as for the dual spacecraft Level 2 FAC product, also in this case only current structures with dimensions larger than the separation among the three spacecraft can be recovered, otherwise spatial aliasing can occur. This threespacecraft ISS-JUB FAC product is computed both from the unfiltered 1 Hz data, and from the filtered magnetic field data, using the same low-pass filter as the dual S/C Level 2 product.

Figure 8.3 panel (c) shows the comparison of dual GFZ FAC with the three ISS-JUB FAC products. This comparison can be considered as a sort of validation for the dual S/C product, given that the three spacecraft product is based on simultaneous measurements, not requiring the stationarity assumption. The general trend of these FAC products agree quite well and the unfiltered ISS-JUB three s/c product provide more details about small scales currents features.

For the events characterized by very high correlation coefficients between the unfiltered magnetic field series (higher than 0.99, like this event), ISS-JUB delivered also a dual S/C FAC product obtained from 1 Hz unfiltered magnetic field data. The justification for this product relies on the fact that, when the correlation coefficient is very high, local FAC structures with dimension smaller than separation between Swarm A and C, should be absent.

The comparison of this ISS-JUB dual S/C unfiltered, with the GFZ dual S/C, and ISS-JUB three S/C FAC product is illustrated in Fig. 8.3, panel (d). We can note that this ISS-JUB unfiltered dual S/C product agree well with the ISS-JUB unfiltered three S/C product, and provide more details about the small scale structure of FAC with respect to the GFZ dual S/C product. This suggests that, at least for the 'stationary events' characterized by very high correlation coefficient, it is possible to remove the low-pass filter in the dual S/C product, obtaining more information about the small scale structures with respect to the Level 2 dual S/C product.

In Fig. 8.3 panel (e) are displayed the FAC products obtained with SECS by FMI, using Swarm A and Swarm C alone, or the lower pair combined. These SECS FAC products recover only the large scale features, but agree remarkably well with the dual S/C FAC product from GFZ and also with the three S/C FAC product from ISS-JUB, also displayed in Panel (e) for comparison. This agreement is particularly important, since the FAC estimates from FMI are based on the expansion in Spherical Elementary Current System (SECS, see Chaps. 2 and 3 of this book), which differ substantially with respect to the methods based on the Ampere's law, used for all the other FAC products.

In Fig. 8.3 panel (f) are shown the various RAL FAC estimates based on 4 nonplanar points using the full curlometer technique (Dunlop et al. 1988, 2002). These 4-point measurements are obtained from 3 S/C data by time-shifting the data from one S/C. The various FAC products displayed in panel (f) differ each other according to which S/C is time-shifted to obtain the fourth point measurement, and the names in the caption adopt the same convention of Dunlop et al. (2015), where p refers to positive time-shift and n to negative time-shift. The various FAC products provided by RAL show a similar behaviour, with a similar trend also with respect to the Level 2 dual S/C product. The time-shift among the various RAL products is related to

**Fig. 8.4** A comparison among the single S/C FAC products from UnCA based on residuals obtained from CHAOS-6 model, with the single S/C Level 2 FAC product, for the auroral event displayed in Fig. 8.3

the change in effective barycentre and different tetrahedral shapes formed by each configuration. All the RAL products recover a significantly larger current density than the standard Level 2 dual S/C product, also larger with respect to the other FAC products from ISS. This can be due to the fact that RAL estimated FAC from the three components of the current density vector, while the other FAC products used the radial component only, which can omit some of the actual FAC. Dunlop et al. (2015) suggested that the degree of stationarity of the current sheet can be obtained from the variance over all different FAC estimates.

The single s/c products obtained from UnCa, based on magnetic residuals from CHAOS-6 model are compared to single s/c Level 2 products provided by GFZ in Fig. 8.3, panel (g). These FAC products from UnCa are perfectly coincident with the Level 2 products, which are instead obtained from the residuals from different magnetic field models (AUX\_CORE, Lithosphere-MF7, magnetosphere-POMME-6, see Chap. 6 of this book). This perfect agreement is also confirmed in the scatter plot in Fig. 8.4, where the single s/c FAC products obtained from UnCa using the CHAOS-6 models are displayed as a function of the Level 2 single s/c FAC products. This agreement demonstrates that the choice of the model used to compute the residuals does not affect FAC calculation.

Other examples of events characterized by high correlation coefficients are the events numbered 25 and 26 in Table 8.2, observed on 20151008 in the intervals 17:10–17:24 UT and 17:24–17:35 UT.

In Fig. 8.5, the magnetic field residuals are reported, in NEC reference frame measured by the three Swarm spacecraft, with the same format as Fig. 8.2. Note that

**Fig. 8.5** The magnetic field residuals detected by the three Swarm spacecraft, in NEC reference frame, for the auroral event on 2015-10-08, 17:10–17:34 UT. Note that this is not a conjunction event among the three Swarm spacecraft. Indeed, Swarm A and C were in the northern hemisphere, and the magnetic signatures detected around 17:16 and 17:28 UT correspond to the entry and exit from the polar cap. Swarm B instead was flying in the southern hemisphere

this is not a conjunction event among the three Swarm spacecraft. Indeed, Swarm A and C were in the northern hemisphere, and the magnetic signatures detected around 17:16 and 17:28 UT correspond to the entry and exit from the polar cap. Swarm B instead, even if it detected similar magnetic signatures approximately at the same times as Swarm A and C, was flying in the southern hemisphere.

In Fig. 8.6 are reported some of the FAC estimates provided for FAC-MICE: In panel (a) are displayed the single and dual S/C products provided by GFZ, which are equivalent to ESA Level 2 FAC products. In panel (b) the filtered single S/C products provided by ISS-JUB and in panel (c) the products from FMI based on SECS method, together with the dual S/C product from GFZ for comparison.

With respect to the previous event (20140711), this event is characterized by higher geomagnetic activity, with AE ≈ 1000, and FACs are more intense. The comparison illustrated in Fig. 8.6 confirms the features highlighted in the previous event: in panel (a) the single and dual S/C products provided by GFZ, which are equivalent to the ESA Level 2 FAC products, show large discrepancies related to the high frequency fluctuations characterizing the single S/C FAC products. This discrepancy is related to the low-pass filters used for the dual S/C FAC product. Indeed, the comparison in panel (b) the dual S/C products from GFZ show a reasonably good agreement with the filtered single S/C FAC products delivered by ISS, obtained using the same low-pass filter.

**Fig. 8.6** Some of the FAC products delivered to ESA for FAC-MICE, for the auroral event observed on 2015-10-08, 17:10–17:34 UT. In panel **a** the single and dual S/C products provided by GFZ, which are equivalent to ESA Level 2 FAC products. In panel **b** the filtered single S/C products provided by ISS-JUB and in panel **c** the products from FMI based on SECS method, together with the dual S/C product from GFZ for comparison

Also in this event, the FAC products obtained from SECS method agree reasonably well with the dual S/C product from GFZ (see panel c).

#### *8.4.2 Comparison for the 'Low Correlation' Events*

The comparison among the various FAC products delivered for FAC-MICE has been performed also for a number of events characterized by lower values of the correlation parameters (Fig. 8.1), to see how the various FAC estimates compare with each other during these 'low-correlation' events.

As an example, we considered the event observed on 2014-05-17 in the interval 03:05–03:15 UT in the southern hemisphere. This is a conjunction event, during which the three Swarm spacecraft were close together, and the magnetic field residuals measured by the three spacecraft are shown in Fig. 8.7, in the NEC reference frame, with the same format as Fig. 8.2. The current sheet is observed approximately at the same time by the three spacecraft, but the differences among the magnetic signatures measured by the three spacecraft are now noticeable.

In Fig. 8.8 are reported some of the FAC estimates provided for FAC-MICE: in

**Fig. 8.7** The magnetic field residuals detected by the three Swarm spacecraft, in NEC reference frame, for the auroral event on 2014-05-17, 03:05–03:14 UT. Note that this is a conjunction event, when the three spacecraft are close together during the auroral crossing

**Fig. 8.8** Some of the FAC products delivered to ESA for FAC-MICE, for the auroral event observed on 2014-05-17, 03:05–03:14 UT. In panel **a** the single and dual S/C products provided by GFZ, which are equivalent to ESA Level 2 FAC products. In panel **b** the filtered single S/C products provided by ISS-JUB, together with the dual S/C product from GFZ for comparison

panel (a) are displayed the single and dual S/C products provided by GFZ, which are equivalent to the ESA Level 2 FAC products, and in panel (b) the filtered single S/C products provided by ISS.

Again, the single and dual S/C products provided by GFZ in panel (a) show large discrepancies, related to the high frequency fluctuations characterizing the single S/C FAC products. However, in this case, also the filtered single S/C FAC products provided by ISS-JUB differ substantially from each other, and also with respect the dual S/C L2 product. It can be noted that the larger deviation is shown by Swarm A single S/C filtered FAC product, while the filtered FAC products obtained from Swarm B and C data agree better to each other, at least until 03:09:20 UT.

This substantial difference is probably related to the presence of local FAC structures, with dimension smaller than the spacecraft separation, rather than nonstationarity of the current sheet. Indeed, only the FAC from Swarm A differ much, while FACs recovered by Swarm B and Swarm C agree quite well to each other, even if the relative time-shift between the two spacecraft is approximately 24 s. The presence of these local FAC structures affect the dual spacecraft FAC product, which is therefore not reliable for this event.

For this event, the three-S/C ISS-JUB FAC product is not available because the configuration of the Swarm constellation was too elongated. Moreover, also the FAC products based on SECS method are not available, since Swarm spacecraft were beyond 76° in magnetic latitude, the threshold beyond which the SECS analysis is not performed. Therefore, for this event, we cannot verify the large scale behaviour of the FAC structure.

As a further example of a 'low correlation' event, we considered the event observed on 2015-03-16, in the interval 22:35–23:00 UT in the northern hemisphere. This is not a conjunction event, and the magnetic field residuals measured by the lower pair are shown in Fig. 8.9, in the NEC reference frame, with the same format as Fig. 8.2.

In Fig. 8.10 are reported some of the FAC estimates provided for FAC-MICE: In panel (a) are displayed the single and dual S/C products provided by GFZ, which are equivalent to ESA Level 2 FAC products. In panel (b) the filtered single S/C products provided by ISS-JUB and in panel (c) the products from FMI based on SECS method, together with the dual S/C product from GFZ for comparison.

The single S/C products in panel (a) show high frequency fluctuations, while the dual S/C product is more stable. The comparison in panel (b) with the filtered single S/C products from ISS-JUB shows a better agreement in the first part of the event, while later on, i.e. after 22:43 UT, also the filtered single S/C FAC estimates differ noticeable from the dual S/C product. This can again be explained in terms of non-stationarity of the current sheet, and/or to strong longitudinal gradients. The signatures in the dual S/C product around 22:45, at ≈80°MLAT, could instead be due to high latitude filamentary FAC, which cannot be detected by single S/C methods because these currents are not planar, and therefore violate the orientation assumption needed for the single S/C FAC calculation (Lühr et al. 2016).

The FAC products based on SECS illustrated in panel (c) follow reasonably well the L2 FAC in the first part of the event (especially the dual spacecraft SECS). After

**Fig. 8.9** The magnetic field residuals detected by Swarm A and C, in NEC reference frame, for the auroral event on 2015-03-16, 22:35–23:00 UT

**Fig. 8.10** Some of the FAC products delivered to ESA for FAC-MICE, for the auroral event observed on 2015-03-16, 22:35–23:00 UT. In panel **a** the single and dual S/C products provided by GFZ, which are equivalent to ESA Level 2 FAC products. In panel **b** the filtered single S/C products provided by ISS-JUB and in panel **c** the products from FMI based on SECS method, together with the dual S/C product from GFZ for comparison'

22:44:40 the SECS product are not available because the magnetic latitude is larger than 76°. Therefore, it's not possible to verify the signature at 22:45 UT.

#### **8.5 FAC-MICE Comparison Summary**

Three different groups participating to FAC-MICE (ISS, UofA, UCL-MSSL) provided quality parameters based on cross-correlation of lagged magnetic field time series (or FAC single spacecraft time series) measured by Swam A and C. These parameters agree well to each other for the various FAC—MICE events, and suggest the presence of two classes of events: the 'high-correlation' events characterized by the higher correlation coefficients and the 'low-correlation' events.

For the 'high-correlation' events, the various FAC estimates provided by the various groups show a reasonably good agreement. In particular, the noticeable discrepancies between the single and dual spacecraft Level 2 FAC products is highly reduced when the same low-pass filter used for dual S/C product is adopted also for the single S/C products. This suggests that the current sheet is both stationary and approximately invariant in longitude during these events. Indeed, the three spacecraft cross the current sheet at different times and in different locations, and recover approximately the same FAC structure.

Another important validation for FAC structure is given by the three spacecraft FAC products, provided for the conjunction events when Swarm B was near the lower pair (Swarm A and C), computed both from filtered and unfiltered 1 Hz data. Indeed, the three S/C products are based on simultaneous measurements from the three spacecraft, and do not rely on the stationarity assumption or orientation assumptions which are instead needed for the single or dual S/C FAC products. However, like the other FAC estimations based on multipoint measurements, these three spacecraft FAC products can recover only FAC with spatial scales larger than the spacecraft separations. During these 'high-correlation' events, also these three S/C products agree well with the other FAC estimates.

During these 'high-correlation' events, since the large scale FAC structure is reliably recovered by these filtered products, it seems possible to look at the smaller scale features of the FAC sheet, provided by various FAC estimates based on unfiltered 1 Hz magnetic field data. However, also during these 'high correlation' events, some of the high frequency fluctuations in the single S/C Level 2 FAC products differ among the various spacecraft. In particular, the differences observed between Swarm A and C, which are close together, suggesting the presence of local structures that do not satisfy the one-dimensional current sheet assumption needed for single spacecraft FAC calculation (Lühr et al. 2015).

During the low correlation events, the comparison among the various FAC products looks instead quite different: the large discrepancies among single and dual S/C FAC estimates remain, even when considering the filtered single S/C products. This suggests that the stationarity assumption is violated, and therefore FAC densities provided by the various products can be affected by non-negligible errors. Some of the large-scale features of FAC structure given by the filtered dual S/C FAC products are confirmed by the SECS analysis also during these low-correlation events.

Other FAC estimates can be used in synergy with Swarm data to investigate the high latitude current system. For example, measurements from incoherent scatter radars can be very useful for this purpose. These radars provide various properties of the ionosphere, like electron density, ion velocity, ion and electron temperatures, ion composition, and collision frequencies, which can be used to infer the other important ionospheric parameters. However, these radars have a limited field of view, and can be used in conjunction with Swarm data only during Swarm overflights. Other methods, like AMPERE or the AMIE procedure, illustrated in Chaps. 7 and 10 of this book respectively, can provide the configuration of FACs on a global scale. Even if the resolution is not comparable with FAC obtained from Swarm, these methods can provide the global context for the more detailed Swarm measurements.

#### **8.6 FAC-MICE Round Table Discussion and Way Forward**

During the second day of the workshop in ESA-ESTEC, a Round Table discussion with all the participants debated the possible evolution of the current Level 2 FAC products, and/or potential new FAC product and FAC quality indicators.

During the discussion, it was pointed out that the various new FAC products developed for FAC-MICE could be useful to examine different features of the FAC structure:

It has been suggested to develop a toolbox for user definable FAC calculation, where the users can choose among the various algorithms to compute FAC. The suggested solution is to implement an open source platform, where the various codes for FAC calculation can be uploaded and shared among the users. In order to guarantee the reproducibility of the results, all the algorithms, filters and documentation used in this FAC toolbox need to be published.

**Acknowledgements** The FAC-MICE activity has been supported by ESA. The authors thank the International Space Science Institute in Bern, Switzerland, for supporting the ISSI Working Group: "Multi-Satellite Analysis Tools, Ionosphere", in which frame this study was performed. The Editors thanks an anonymous reviewer for his assistance in evaluating this chapter.

#### **References**

#### **The FAC-MICE Team**

K. Kauristie1, S. Käki1, Heikki Vanhamäki2, L. Juusola1, Adrian Blagau3,4, Joachim Vogt4, Octav Marghitu3, M. W. Dunlop5,6, Y.-Y. Yang7, J.-Y. Yang5, Hermann Lühr8, Guram Kervalishvili8, Jan Rauberg8, Claudia Stolle8, Ivan P. Pakhotin9, Ian R. Mann9, C. Forsyth10, I. J. Rae10, Jiashu Wu11, J. Gjerloev12, S. Ohtani12, M. Friel12

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 9 Spherical Cap Harmonic Analysis Techniques for Mapping High-Latitude Ionospheric Plasma Flow—Application to the Swarm Satellite Mission**

#### **Robyn A. D. Fiori**

**Abstract** This chapter describes spherical cap harmonic analysis (SCHA) for mapping ionospheric plasma flows measured by the Swarm satellites. In Sect. 9.1, SCHA is introduced as a tool for mapping a variety of one, two, and three-dimensional parameters. Section 9.2 provides a detailed summary of the theory pertaining to SCHA including a discussion of the spherical cap coordinate system, boundary conditions and basis functions, calculation of non-integer degree, and practical considerations. Section 9.3 provides a practical example of SCHA mapping of ionospheric plasma flow for a ground-based data set, and Sect. 9.4 focuses on two-dimensional SCHA mapping of Swarm ion drift measurements both independently and in conjunction with measurements from other instruments.

#### **9.1 Introduction**

A wealth of data exists to describe physical parameters in the near-Earth environment with data sets taken from instruments that are either discretely located or in constant motion taking measurements that range in spatial and temporal resolution. Although examination of these individual data sets provides important information on the evolution of the parameter at a single location over time, or at several discrete locations at a single time, it is often desirable to ascertain what is happening over a specific spatial region at a given time through interpolation or mapping techniques. The technique examined in this Chapter is a form of harmonic analysis, which involves the series expansion of an infinite set of appropriately chosen basis functions. The choice of which basis functions to use is dependent on the distribution of data, the area of interest, and the nature of the function being mapped. For data distributed about the entire spherical Earth, the normal approach is spherical harmonic analysis with basis functions that are comprised of the Associated Legendre polynomials in co-latitude and trigonometric functions in longitude. However, when either the data or region of

R. A. D. Fiori (B)

Geomagnetic Laboratory, Natural Resources Canada, Ottawa, ON, Canada e-mail<EMAIL_ADDRESS>
interest are confined to only a portion of the sphere, a spherical harmonic expansion cannot be applied because it is not defined across the entire natural interval (i.e. the sphere), and new basis functions must be defined for the relevant portion of the sphere. Although the new basis functions are less convenient because the function is not likely periodic over the subinterval, they are mathematically correct and allow the expansion to be applied. In the case of mapping data over a portion of a sphere, a spherical cap harmonic analysis (SCHA) technique has been developed.

Spherical cap harmonic analysis was originally developed by G. V. Haines for mapping regional magnetic fields, and has been applied to create regional maps across the globe (e.g. Haines 1985a, 1988, 2007; An 1993; An et al. 1998; De Santis et al. 1997b; Feng et al. 2015; Ji et al. 2006; Kotze, 2014; Pavón-Carrasco et al. 2008; Pothier et al. 2015; Stening et al. 2008; Tozzi et al. 2013; Weimer 2013; Waters et al. 2015). As a practical example, the Canadian Geomagnetic Reference Field (CGRF) is a regional magnetic field map that was produced every 5 years from 1985 to 2005 over Canada and adjacent regions (Haines and Newitt 1997).

Spherical cap harmonic analysis has been used for mapping a wide variety of ionospheric and magnetospheric parameters including secular variation over Canada (Haines 1985c, 1993), China (Chen et al. 2011), Europe (Korte and Haak 2000), Italy (De Santis et al. 1997a), the North Atlantic (Pavón-Carrasco et al. 2013), South Africa (Nahavo et al. 2011) and Spain (Garcia et al. 1991; Torta et al. 1992); vertical magnetic field anomalies poleward of 40° (Haines 1985b); magnetic field intensity (Kotzé 2002); Birkeland current systems (Green et al. 2006); equivalent current systems (Gaya-Pique et al. 2008; Haines and Torta 1994); ionospheric total electron content (TEC) (Liu et al. 2011, 2014; Ghoddousi-Fard et al. 2011; Otsuki et al. 2011); electric fields (Fiori et al. 2010), ionospheric conductance estimates (Green et al. 2007); and the critical frequency of the F2 layer (*f oF2* peak) (De Santis et al. 1991, 1992, 1994; Lazo et al. 2004). Notably, SCHA algorithms have been incorporated into the Assimilative Mapping of Ionospheric Electrodynamics (AMIE) mapping algorithm (Richmond and Kamide 1988). In addition to ionospheric and magnetospheric parameters, SCHA has been applied for mapping other quantities. One example is regional gravity field modelling on the moon (Han 2008), and on the Earth to which SCHA was originally applied by Jiancheng et al. (1995) with clarifications and corrections on its use by De Santis and Torta (1997) and additional work by Hwang et al. (2012). Another example is modelling sea level data (Hwang and Chen 1997).

This Chapter discusses the application of SCHA for mapping plasma flow in the high-latitude ionosphere based on data from multi-satellite missions such as is available from the Swarm mission.

#### **9.2 Theory**

The theory of SCHA was originally described in Haines (1985a) for mapping regional magnetic fields. It is assumed that a data set or region of interest is described on a limited portion of the spherical Earth described by a spherical cap. Data may be mapped in one, two or three dimensions. For the purpose of the electrostatic potential and plasma flow, two-dimensional mapping is considered. The necessary theory is summarized here with discussions on the relevant boundary conditions and recommendations on the appropriate choice of mapping parameters.

#### *9.2.1 Spherical Cap Geometry*

The spherical cap is constructed by extracting a conical section of angular radius θ<sup>c</sup> from a sphere and considering the outer slab of the section which extends from radius r1 to r2, see Fig. 9.1. A spherical cap coordinate system is defined with respect to the centre of the outer shell of the spherical cap with θ and φ representing co-latitude and longitude, respectively, in the spherical cap coordinate system. Co-latitude spans from 0° to θ<sup>c</sup> whereas ϕ spans the full 0°–360° range. For the two-dimensional mapping considered here, r1 = r2 and the spherical cap reduces to a thin spherical cap shell. In general, the parameter being mapped is given in a spherical coordinate system (i.e. geographic or geomagnetic) and it is necessary to first transform both coordinate locations and the vector-pointing directions to the spherical cap coordinate system to simplify calculations.

Figure 9.2 illustrates the transformation from an arbitrary spherical coordinate system to the spherical cap coordinate system. Consider an arbitrary point *P* in the spherical coordinate system having coordinates (θ- , φ- ) measured with respect to the North Pole (NP). Point P has corresponding coordinates of (θ, φ) in the spherical cap

**Fig. 9.1** Three-dimensional spherical cap

**Fig. 9.2** Cartoon illustrating the transformation of the north-pole centred coordinate system to the spherical-cap coordinate system

coordinate system which is centred about the spherical cap pole (SCP), which does not necessarily coincide with NP. If the location of SCP corresponds to (θSCP, ϕSCP) in the spherical coordinate system, then (θ, φ) can be determined from (θ- , φ- ) using spherical trigonometry as follows:

$$\cos\theta = \cos\theta^{\prime}\cos\theta\_{SCP} + \sin\theta^{\prime}\sin\theta\_{SCP}\cos(\varphi^{\prime} - \varphi\_{SCP})\tag{9.1}$$

$$\cos(180^\circ - \theta') = \cos\theta\cos(180^\circ - \theta\_{SCP}) + \sin\theta\sin(180^\circ - \theta\_{SCP})\cos\varphi,\quad(9.2)$$

where Eq. (9.2) can be rearranged to solve for φ

$$\cos\varphi = \frac{\cos\theta\cos\theta\_{SCP} - \cos\theta'}{\sin\theta\sin\theta\_{SCP}}.\tag{9.3}$$

Appropriate selection of the size and location of the spherical cap is necessary to ensure optimal performance of the SCHA algorithm. Both the central location of the spherical cap and θ<sup>c</sup> should be chosen to ensure adequate data coverage across the entire spherical cap to support the desired spatial resolution of the resultant mapped data (e.g. see Sect. 9.5). When data are not distributed uniformly techniques can be applied to describe the region of the spherical cap which is adequately constrained by data to produce an accurate map. Such techniques are discussed in Sect. 9.5.

#### *9.2.2 Mapping Electrostatic Potential and Ion Flow by Series Expansion*

It is often desirable to represent the electrostatic potential (*<sup>E</sup>* ) over a global or discrete region. For regional mapping, SCHA may be applied. Consider a spherical cap located on a thin spherical shell within the ionosphere on which the potential is to be mapped. Regional mapping of *<sup>E</sup>* is performed by representing *<sup>E</sup>* by series expansion (e.g. Fiori et al. 2010)

$$\Phi\_E(\theta,\varphi) = \sum\_{k=0}^{K\_{\text{max}}} \sum\_{m=0}^{M\_{\text{max}}} [A\_{km}\cos(m\varphi) + B\_{km}\sin(m\varphi)]P\_{n\_k(m)}^m(\cos\theta) \tag{9.4}$$

where θ and ϕ represent co-latitude and longitude in the spherical cap coordinate system, *P<sup>m</sup> nk* (*m*)(*cos*θ) are the Associated Legendre functions of the first kind having a non-integer degree *nk (m)* and order *m*, *k* is an integer degree-indexing term, and *Kmax* and *Mmax* are the maximum degree-index and order which truncate the series expansion. *Mmax* ≤ *Kmax* and the choice of both *Mmax* and *Kmax* is dependent on the desired latitudinal and longitudinal mapping resolution and the distribution of data within the spherical cap. In the following notation, *Mmax* is taken to be *k*. In Eq. (9.4), coefficients *Akm* and *Bkm* are constants for each combination of *k* and *m* which represent the amplitude of the harmonics. These coefficients must be determined to describe *<sup>E</sup>* at any point within the spherical cap. Equations for calculating the Associated Legendre functions with non-integer degree follow in Sect. 9.2.3.

Electrostatic potential is related to the electric field (**E**) and velocity (**v**) on the spherical shell through

$$E = -\nabla \Phi\_E,\tag{9.5}$$

and

$$\begin{aligned} \mathbf{v} &= \frac{\mathbf{E} \times \mathbf{B}}{\mathbf{B}^2} \\ &= \frac{E\_\varphi}{B}\hat{\theta} + \frac{E\_\theta}{B}\hat{\varphi} \end{aligned} \tag{9.6}$$

where **B** is magnetic field and is commonly taken to be the radial component of either the dipole magnetic field or of the International Geomagnetic Reference Field (IGRF), although measurements could be used. Breaking Eq. (9.5) into components yields

$$E\_{\varphi} = -\frac{1}{r\sin\theta} \frac{d\Phi\_E}{d\varphi} = \frac{1}{r\sin\theta} \sum\_{k=0}^{K\_{\text{max}}} \sum\_{m=0}^{k} [A\_{km}\sin(m\varphi) - B\_{km}\cos(m\varphi)] m P\_{n\_k(m)}^m(\cos\theta) \tag{9.7}$$

$$E\_{\theta} = -\frac{1}{r}\frac{d\Phi\_E}{d\theta} = -\frac{1}{r}\sum\_{k=0}^{K\_{\text{max}}}\sum\_{m=0}^{k} [A\_{km}\cos(m\varphi) - B\_{km}\sin(m\varphi)]\frac{dP\_{n\_k(m)}^m(\cos\theta)}{d\theta} \tag{9.8}$$

which can be substituted into Eq. (9.6) to describe velocity as

$$v\_{\theta} = \sum\_{k=0}^{K\_{\text{max}}} \sum\_{m=0}^{k} [A\_{km} \sin(m\varphi) - B\_{km} \cos(m\varphi)] \frac{m P\_{n\_l(m)}^m(\cos\theta)}{r \text{ B} \sin\theta} \tag{9.9}$$

$$v\_{\varphi} = \sum\_{k=0}^{K\_{\text{max}}} \sum\_{m=0}^{k} [A\_{km} \cos(m\varphi) + B\_{km} \sin(m\varphi)] \frac{1}{rB} \frac{dP\_{n\_k(m)}^m(\cos\theta)}{d\theta}. \tag{9.10}$$

In Eqs. (9.7)–(9.10) *r* represents the radius of the spherical shell on which the surface is mapped. Typically this is taken to be the radius of measurement. For representing ionospheric plasma flow this is more commonly the radius of the Fregion ionosphere (typically ~300 km altitude), and it is sometimes necessary to map measurements to this altitude.

Coefficients *Akm* and *Bkm* are determined by minimizing the difference between parameters that are measured by the instrument and represented using Eqs. (9.4), (9.7) and (9.8), or (9.9) and (9.10), depending on the parameter being measured. Minimization is performed using a method such as linear regression or singular value decomposition (SVD) (Press et al. 1992; Brandt 1998). To minimize anomalies due, for example, to a poorly constrained solution over a portion of the spherical cap, the coefficient rejection procedure described in Haines and Fiori (2013), which eliminates non-significant fitting coefficients, can be applied.

To better understand how the fitting coefficients are determined regardless of whether *<sup>E</sup>* , *E*, or *v* is known, consider observations of some arbitrary parameter *f* (*r*,θ,ϕ) scattered across a spherical cap which can be represented as

$$\mathbf{f}(r,\theta,\varphi) = \sum\_{k=0}^{K\_{\text{max}}} \sum\_{m=0}^{k} A\_{km} \mathbf{c}\_{km}(r,\theta,\varphi) + B\_{km} \mathbf{d}\_{km}(r,\theta,\varphi) \tag{9.11}$$

where *Akm* and *Bkm* are fitting coefficients as before, and c*km*(*r*,θ,ϕ) and d*km*(*r*,θ,ϕ) are basis functions. Here there are (*Kmax* + 1) <sup>2</sup> fitting coefficients. To model f(*r*,θ,ϕ) anywhere within the spherical cap, physical measurements*f <sup>i</sup>* must be substituted into Eq. (9.11) to solve for the fitting coefficients where

$$\mathbf{f}\_{i}(r\_{i},\theta\_{i},\varphi\_{i}) = \sum\_{k=0}^{K\_{\text{max}}} \sum\_{m=0}^{k} A\_{km} \mathbf{c}\_{km}(r\_{i},\theta\_{i},\varphi\_{i}) + B\_{km} \mathbf{d}\_{km}(r\_{i},\theta\_{i},\varphi\_{i}) \tag{9.12}$$

Consider *N* observations of some measured quantity f*i*(*ri*, θ*i*, ϕ*i*). An array of *N* equations may be written as follows:

$$\begin{aligned} f\_1 &= A\_{00}c\_{00\_1} + B\_{00}d\_{00\_1} + A\_{10}c\_{10\_1} + B\_{10}d\_{10\_1} + \dots + A\_{KK}c\_{KK\_1} + B\_{KK}d\_{KK\_1} \\ f\_2 &= A\_{00}c\_{00\_2} + B\_{00}d\_{00\_2} + A\_{10}c\_{10\_2} + B\_{10}d\_{10\_2} + \dots + A\_{KK}c\_{KK\_2} + B\_{KK}d\_{KK\_2} \\ &\vdots \\ f\_N &= A\_{00}c\_{00\_N} + B\_{00}d\_{00\_N} + A\_{10}c\_{10\_N} + B\_{10}d\_{10\_N} + \dots + A\_{KK}c\_{KK\_N} + B\_{KK}d\_{KK\_N} \end{aligned}$$

where, for simplification, the notation *K* = *Kmax*, c*kmi* = c*km*(*ri*, θ*i*, ϕ*i*), and d*kmi* = d*km*(*ri*, θ*i*, ϕ*i*) has been used. Such equations may be written in matrix form as

$$
\begin{pmatrix} f\_1 \\ f\_2 \\ \vdots \\ f\_N \end{pmatrix} = \begin{pmatrix} c\_{00\_1} \ d\_{00\_1} \ c\_{10\_1} \ d\_{10\_1} \cdots \ c\_{KK\_1} \ d\_{KK\_1} \\ c\_{00\_2} \ d\_{00\_2} \ c\_{10\_2} \ d\_{10\_2} \cdots \ c\_{KK\_2} \ d\_{KK\_2} \\ \vdots \\ c\_{00\_N} \ d\_{00\_N} \ c\_{10\_N} \ d\_{10\_N} \cdots \ c\_{KK\_N} \ d\_{KK\_N} \end{pmatrix}, \quad \begin{pmatrix} A\_{00} \\ B\_{00} \\ A\_{10} \\ B\_{10} \\ \vdots \\ A\_{KK} \\ B\_{KK} \end{pmatrix}
$$

or more simply as

$$F = CD \cdot AB \tag{9.13}$$

where *F* represents an *N* × 1 matrix of the measured quantities, *CD* is an *N* × (*Kmax* + 1)<sup>2</sup> matrix, and AB is a (*Kmax* + 1)<sup>2</sup> × 1 matrix composed of the *Akm* and *Bkm* fitting coefficients. In the case of electrostatic potential, electric field, and plasma flow considered in this Chapter, parameter *f <sup>i</sup>* is taken to represent measured values of *<sup>E</sup>* , or measured components of *E* or **v**, and CD is derived from Eq. (9.4) for *<sup>E</sup>* , and Eqs. (9.7) and (9.8) for *E*, or Eqs. (9.9) and (9.10) for *v*.

It is worth noting that matrix *CD* is entirely dependent on measurement coordinates. Applications of SCHA requiring the repeated mapping of data from fixed stations may take advantage of this fact as *CD* need only be calculated once. Such a consideration greatly enhances the speed of calculation.

#### *9.2.3 Associated Legendre Functions*

It is useful to define the Associated Legendre Functions and their calculation with a non-integer degree term. Haines (1985a, 1988) describe the Associated Legendre Functions of non-integer degree, *nk (m)*, and order, m, as a power series given by

$$P\_{n\_k(m)}^m(\cos \theta) = \sum\_{j=0}^{\infty} A\_j(m, n) \sin^{2j} \frac{\theta}{2} \tag{9.14}$$

with a derivative with respect to θ given by

$$\frac{d P\_{n\_k(m)}^m(\cos \theta)}{d\theta} = \begin{cases} \frac{\sin \theta}{2} \sum\_{j=1}^{\infty} j A\_j(m, n) \sin^{2(j-1)} \frac{\theta}{2} & m = 0\\ \frac{\sin \theta}{2} \sum\_{j=1}^{\infty} j A\_j(m, n) \sin^{2(j-1)} \frac{\theta}{2} + \cos \theta \Big[ \frac{m}{\sin \theta} P\_{n\_k(m)}^m(\cos \theta) \Big] m > 0 \end{cases} \tag{9.15}$$

where *Aj*(*m*, *n*) is given by

$$A\_j(m,n) = \begin{cases} K\_n^m \sin^m \theta & j=0\\ A\_{j-1}(m,n) \frac{(j+m-1)(j+m)-n(n+1)}{j(j+m)} & j>0 \end{cases} \tag{9.16}$$

and *K <sup>m</sup> <sup>n</sup>* is a normalization constant. Normalization is necessary to ensure manageable calculations in determining the solution of the expansion coefficients. Sample normalizations include the Schmidt or Neumann normalizations which take *K <sup>m</sup> <sup>n</sup>* = 1 for m = 0, and

$$K\_n^m = \frac{2^{\frac{1}{2}}}{2^m m!} \left[ \frac{(n+m)!}{(n-m)!} \right]^{\frac{1}{2}} \tag{9.17}$$

$$K\_n^m = (-1)^m \frac{1}{2^m m!} \left[ \frac{(n+m)!}{(n-m)!} \right] \tag{9.18}$$

for m = 0 where Eq. (9.17) represents the Schmidt normalization and Eq. (9.18) represents the Neumann normalization with the Condon–Shortley phase term. Ultimately, the choice of normalization coefficient is not in itself crucial, as long as the same normalization is used consistently throughout the calculation.

When dealing with a non-integer degree term, the factorials in the normalizations constants given by Eqs. (9.17) and (9.18) are replaced by

$$n! = \Gamma(n+1)\tag{9.19}$$

where the gamma function () is given by

$$\Gamma(\mathbf{x}) = \int\_0^\infty t^{x-1} e^{-t} dt. \tag{9.20}$$

In practice, the power series in Eq. (9.16) is truncated at *J* determined by evaluating each term during computation. Haines (1988) suggest a truncation limit of 60 terms.

#### *9.2.4 Boundary Conditions and Basis Functions*

The solution of Eq. (9.4) is, of course, subject to boundary conditions. The first boundary condition requires continuity in longitude such that

$$\begin{aligned} \Phi\_E(\mathbf{r}, \theta, \varphi) &= \Phi\_E(\mathbf{r}, \theta, \varphi + 2\pi) \\ \frac{\partial \Phi\_E(\mathbf{r}, \theta, \varphi)}{\partial \varphi} &= \frac{\partial \Phi\_E(\mathbf{r}, \theta, \varphi + 2\pi)}{\partial \varphi} \end{aligned} \tag{9.21}$$

restricting *m* to be both real and integer valued. The second boundary condition requires that *<sup>E</sup>* be independent of ϕ at the spherical cap pole

$$\begin{array}{c} \Phi\_E(\mathbf{r}, 0, \boldsymbol{\varphi}) = 0 \text{ m} \neq 0\\ \frac{\partial \Phi\_E(\mathbf{r}, 0, \boldsymbol{\varphi})}{\partial \boldsymbol{\theta}} = 0 \quad \mathbf{m} = \mathbf{0} \end{array} \tag{9.22}$$

thereby ensuring regularity of *<sup>E</sup>* . This condition has already been satisfied by using Associated Legendre Functions of the first kind in Eq. (9.4) and excluding those of the second kind. The final boundary condition requires *<sup>E</sup>* and ∂*<sup>E</sup>* ∂θ be arbitrary at the boundary of the spherical cap (θ = θc):

$$\Phi\_E(\mathbf{r}, \theta\_c, \varphi) = \mathbf{f}(\mathbf{r}, \varphi) \tag{9.23}$$

$$\frac{\partial \Phi\_E(\mathbf{r}, \theta\_c, \varphi)}{\partial \theta} = \mathbf{g}(r, \varphi) \tag{9.24}$$

where f(r, ϕ) and g(r, ϕ) are arbitrary functions. This boundary condition is met by choosing *nk (m)* such that

$$\left.P\_{n\_{\mathbf{k}}(m)}^{m}(\cos\theta)\right|\_{\theta=\theta\_{\mathbf{c}}} = 0 \quad \mathbf{k}-\mathbf{m} = \text{odd} \tag{9.25}$$

$$\left. \frac{d P\_{n\_k(m)}^m(\cos \theta)}{d \theta} \right|\_{\theta = \theta\_c} = 0 \quad \text{k}-\text{m} = \text{even} \tag{9.26}$$

forming two sets of basis functions having a real degree *nk (m)* which is not necessarily an integer (Haines 1985a).

Although two basis functions are defined (Eqs. (9.25) and (9.26)) it is not necessary, and sometimes not advisable, to use both; there may be a physical significance behind choosing only one basis function, or it might be desirable to simplify the solution. The choice of basis functions depends on the parameter being mapped with consideration of the significance and cost of using one or two basis functions.

There are many examples in the literature where only the k–m = odd basis function (Eq. (9.25)) is used. Green et al. (2007) apply SCHA to mapping ionospheric conductance combining magnetic data from satellites and magnetometers, and electric field data from the Super Dual Auroral Radar Network (SuperDARN) and the Defense Meteorological Satellite Programs (DMSP) satellites. SCHA techniques are applied to all data sets using only the k–m = odd basis function, which forces the mapped parameter to zero at the boundary of the spherical cap. They do this by setting the spherical cap boundary to zero equatorward of the current systems to reduce the impacts on the mapped current systems. In one instance they set θ<sup>c</sup> = 50° MLAT, but limited the plotting region to θ > 40°. Limiting the solution by calculating *nk (m)* using the k–m = 0 condition can be useful if it physically significant for the mapped parameter to be zeroed at the boundary of the spherical cap. Green et al. (2006) follow similar techniques in mapping the Birkeland currents based on Iridium and SuperDARN data. In another example, Weimer (2005) and Fiori et al. (2010) make use of this technique for mapping the ionospheric convection pattern using a spherical cap centred at the geomagnetic pole and extending to the equatorward boundary of the convection zone, where the plasma flow is expected to be zero. It should be noted that using only the k–m = odd basis function can lead to Gibb's phenomenon (often referred to as ringing) at the boundary of the spherical cap. Torta et al. (1992) propose using a spherical cap much larger than the region mapped to accommodate the larger wavelength to solve this problem. Solutions outside of the mapping region are ignored. Limiting the solution using the k–m = 0 basis function forces the mapped parameter to zero at the boundary of the spherical cap which is a useful tool when it is physically significant to do so.

Other situations have required the use of only the *k*-*m* = even basis function (Eq. (9.26)), thereby forcing the derivative of the mapped parameter to zero at the spherical cap. Examples include Weimer et al. (2010), Weimer (2013) and Pothier et al. (2015) who use the *k*-*m* = even basis functions to map geomagnetic perturbations from ground-based magnetometer data for cap-sizes ranging from 50° to 90°. These three papers are related and use the SCHA algorithm described here with only the *k*-*m* = even basis functions to reduce unnatural oscillations at the low-latitude boundary of the spherical cap where data are more sparsely distributed. A widely used application of the *k*-*m* = even basis function is in the Assimilative Mapping of Ionospheric Electrodynamics (AMIE) technique which maps various electrodynamic quantities in the high-latitude ionosphere. AMIE fits data over the hemisphere using SCHA functions with *k*-*m* = even basis functions poleward of 56° magnetic latitude followed by another tapering function from 56° magnetic latitude to the equator (Richmond and Kamide 1988).

The selection of what basis function(s) to use requires a discussion of both uniform convergence and orthogonality of the function. Haines (1990) explains that independently these basis functions are uniformly convergent and orthogonal, but their derivatives are not uniformly convergent unless of course the solution really does have a zero slope at the boundary of the spherical cap. Considering both basis functions in the solution ensures that both the parameter being mapped, and its derivative is uniformly convergent.

The consequence of the uniform convergence gained by using both sets of basis functions is sparse non-orthogonality of the system over the expansion interval. The two functions *P<sup>m</sup> <sup>n</sup> <sup>j</sup>*(*m*)(*cos*θ) and *<sup>P</sup><sup>m</sup> nk* (*m*)(*cos*θ) are orthogonal when j =k and when j-m and k-m are either both even or both odd but they are not orthogonal when one is even and one is odd. Non-orthogonality means that the expansion terms must be expressed using a specific truncation limit (i.e. not infinity) and the coefficients determined are therefore dependent on the maximum truncation level and are dependent on the choice of *Kmax*. In terms of computation of the solution, non-orthogonality causes nondiagonal terms (ill-conditioning) in the least-squares matrix affecting the solution of the fitting coefficients. Korte and Holme (2003) criticize the commonly used practice of coefficient rejection discussed in Sect. 9.2.2 as having no physical justification in SCHA due to the incomplete orthogonality of the basis functions. Problems are more pronounced for high truncation levels, particularly if the field being mapped has a large amplitude, and for small spherical caps. Although this non-orthogonality is sparse, and it is possible to minimize the impacts, it should be considered when using both sets of basis functions.

Extreme care should be taken in determining whether or not one or both basis functions should be used. In general, limiting the solution with the *k*-*m* = odd basis functions severely limits the solution and should only be applied when it is physically significant to force the solution to zero at the boundary of the spherical cap. To allow the solution to take on an arbitrary value at the boundary of the spherical cap (unless of course the solution is required to be zero at the boundary), it is, therefore, recommended that if SCHA is to be applied with only one boundary condition, then the *k*-*m* = even boundary condition should be used. Although it may be desirable to promote orthogonality by using only one set of basis functions, keep in mind that the consequence is essentially throwing out half the solution, and this is a consequence that must be considered.

#### *9.2.5 Non-integer Degree*

SCHA requires the solution of a non-integer degree term*nk (m)*to describe the Associated Legendre Functions and their derivatives along the boundary of the spherical cap (i.e. Eqs. (9.25) and (9.26)). *P<sup>m</sup> nk* (*m*)(cos θ) θ=θ*<sup>c</sup>* and *d P<sup>m</sup> nk* (*m*)(cos θ) *d*θ θ=θ*<sup>c</sup>* are oscillating functions of *nk (m)* given *m* and θc. The degree-index term *k* starts at *m* and is incremented by one every time a root to either Eqs. (9.25) or (9.26) is found. These roots represent the *nk (m)* values and occur at a zero amplitude for k-m = odd (Eq. (9.25)) and at maximum amplitude for k-m = even (Eq. (9.26)). In this way, the difference *k*-*m* fluctuates between even and odd values as *P<sup>m</sup> nk* (*m*)(cos θ) θ=θ<sup>0</sup> fluctuates between maximum and zero amplitude. Figure 9.3 shows an example of *P<sup>m</sup> nk* (*m*)(cos θ)

 θ=θ<sup>0</sup> versus *nk (m)* for *m* = 0 and θ<sup>0</sup> = 10°, 20°, and 30°. Table 9.1 provides an example of the *nk (m)* values for all combinations of *k* and *m* for θ<sup>0</sup> = 30°.

Before continuing, it is useful to take a moment to fully understand what the non-integer degree term *nk (m)* represents. In ordinary SHA, *nk (m)* is equivalent to the integer-valued degree-indexing term *k* and describes the sectioning of features in both latitude and longitude across the sphere; there are k-m divisions in latitude and

**Fig. 9.3** Associated Legendre function *P<sup>m</sup> nk* (*m*)(cos θ) θ=θ*<sup>c</sup>* versus non-integer degree *nk (m)* for *m*=0 and θ<sup>c</sup> = 10°, 20°, and 30°

2 m divisions in longitude. In SCHA, *k* and *m* also describe sectioning, but over the spherical cap opposed to the entire sphere; there are (k-m + 1)/2 divisions in latitude and 2 m divisions in longitude. Here sectioning defines the resolution of features that can be mapped using either SHA or SCHA, with a minimum resolution described for *k* = *Kmax* and *m* = 0 (for example see Fig. 9.2 of Fiori et al. 2014).

Regional mapping of ionospheric parameters using SCHA has advantages over global-scale mapping of the same parameters. Most notably, the distribution of ground-based observatories is highly dependent on the location of land and on political boundaries. It is, therefore, possible to obtain a higher resolution map over a specific region of interest (i.e. a specific country or continent) than over the globe. Both Haines (1990) and Amm and Viljanen (1999) point out that since *nk (m)* > *k* (see Table 9.1), the SCHA techniques uses higher order Associated Legendre functions than the SHA technique, and is, therefore, able to achieve a higher spectral resolution over SHA. As an example, the CGRF, which is a regional version of the IGRF with a spherical cap centred over Canada, is able to achieve a higher resolution depicting spatial variations of 833 km compared to 4000 km for the IGRF. In the case of a sparse global data set with irregularly spaced pockets of high-resolution data an SHA solution may require more basis function than data points resulting in a poor fit whereas an SCHA fit over a specific region of dense measurement offers a higher resolution fit with fewer basis functions. Practical considerations

In a theoretical setting, data are uniformly distributed about the spherical cap and the resolution of the mapped structures is a known quantity. Realistically conditions are much less ideal. This Section addresses some of the considerations necessary to apply SCHA based on a more realistic data set.

In practice, data are not uniformly distributed about the entire spherical cap as instrument locations are limited by physical parameters such as a satellite orbit or the availability of land. Mapping using a measurement-only data set has an optimal solution in the sense that the differences between the measured and calculated parameters are minimized. However, the coefficients are only well constrained over the region of measurement. Over regions not constrained by measurements, the same coefficients can cause wildly unrealistic solutions for the mapped parameter. If a smooth solution is required over the entire spherical cap, the fitting must be constrained over regions lacking measurements, using, for example, data from a statistical or empirical model.

There are numerous examples where a statistical model has been used to constrain a non-uniform data set (e.g. Richmond and Kamide 1988; Shepherd and Ruohoniemi 2000). Although the use of a model is sometimes statistically significant or necessary, care should be taken to ensure there are enough data points to overwhelm the model. Alternatively, weighting factors can be used to put more emphasis on measured data compared to the model. For an example of weighting applied to spherical harmonic mapping see Ruohoniemi and Baker (1998) who mapped ionospheric plasma flow using a modified form of spherical harmonic analysis, or Rogers and Honary (2015) who mapped absorption derived from riometer measurements.

Often it is not desirable to constrain an under sampled parameter with contributions from a statistical or empirical model, and other solutions must be explored. Walker (1989) faced a similar problem in applying SCHA techniques for mapping magnetic data based on a sparse distribution of magnetometer stations across Canada. He performed his mapping by reducing the number of longitudinal coefficients by selecting *Mmax* < *Kmax* in the spherical cap harmonic expansion to effectively reduce the longitudinal resolution of the model. He argues that this is applicable to magnetic data as the mapped features tend to have more latitudinal than longitudinal structure and reducing the number of longitudinal coefficients enhances zonal features. Such a process also serves to prevent the number of coefficients from exceeding the number of observations and provides a first order estimate of the mapped parameter. Weimer (1995) also proposes smoothing the data by selecting *Mmax* < *Kmax* for the purpose of mapping electrostatic potentials at high-latitudes, although it should be noted that their technique relies on a modified SHA algorithm opposed to SCHA. Weimer et al. (2010), Weimer (2013) and Pothier et al. (2015) handle the problem of sparse data at the boundary of a spherical cap by limiting the solution using only the k-m = 0 basis functions to reduce ringing at the boundary of the spherical cap. Although these sparse mapping techniques work well in the examples discussed, their use is not universal.

Another solution to handle an under sampled spherical cap is to apply a masking algorithm to map the mapped parameter only where constrained by data. Fiori et al. (2014) describe such an algorithm which they apply to limit mapping to a region of reliability surrounding measurement data points. They call the mapped region the 'region of constraint' or RoC. The method of determining the RoC is summarized here.

Essentially, a radius θ is defined to describe the region surrounding each data point around which the function is considered to be well constrained:

$$
\Delta\theta = \frac{4\theta\_c}{pK\_{\text{max}}}\tag{9.27}
$$

where θ<sup>c</sup> and *Kmax* are described as before, and *p* represents the number of points that must be known along a given wavelength to accurately map a given wavelength. For each data point, the number of additional data points located within θ of the data point being evaluated is counted. If there are at least *p* data points, the region defined by θ is considered valid. The sum of all valid regions constitutes the RoC.

For a given θ<sup>c</sup> and *Kmax*, accurate portrayal of the mapped parameter requires an appropriate definition of the RoC, which, from Eq. (9.27) is entirely dependent on the selection of *p*. The appropriate choice of *p* is somewhat ambiguous. Theoretically, a uniform distribution of data points would require *p* = 2 (i.e. the Nyquist frequency) for each spatial dimension being mapped. For a more realistic distribution of points, larger values of *p* are recommended. An incorrect selection of *p* can unnecessarily restrict (*p* too large) or inflate (*p* too small) the RoC. Restricting the RoC would potentially cause the loss of valuable information, resulting in a smoothed solution, whereas inflating the RoC risks the inclusion of poorly constrained portions of the mapped function, which are not representative of the true mapped parameter.

Figure 9.4 presents examples of the RoC applied to hypothetical coordinates for the Swarm satellite mission. For the purpose of this Figure, data are considered at a 2 Hz resolution over a 4-min interval along a possible satellite path where the lower satellite pair (Swarm A/C) is illustrated in blue and green and the upper satellite (Swarm B) is illustrated in red. Figure 9.4a, b are for the case where the spherical cap is large and covers the entire high-latitude convection zone with θ *<sup>c</sup>* = 32°, and the RoC is calculated for *KMAX* = 6. The RoC in Fig. 9.4a is shown for *p* = 6 and is much larger than the RoC of Fig. 9.4b where a more conservative value of *p* = 28 was chosen. In Fig. 9.4c a spherical cap centred over the area of measurement is considered and the RoC is shown for *Kmax* = 2 and *p* = 6.

Regardless of the data distribution, it is necessary to select an appropriate value for *Kmax*, and therefore *nk (m)*, to ensure an accurate representation of the mapped parameter based on the available data. Selecting too low a value of *Kmax* will result in an overly smooth map whereas too large a value will cause overfitting of the data. Consider, for example, mapping some arbitrary function on a spherical cap having a surface area *A* where

**Fig. 9.4** Grey shading indicates the region of constraint (RoC) for measurements taken at a 2 Hz resolution along the hypothetical path of the Swarm A (blue), C (green), and B (red) satellites for **a** θ *<sup>c</sup>* = 32°, *KMAX* = 6, and *p* = 6, **b** θ *<sup>c</sup>* = 32°, *KMAX* = 6, and *p* = 28, and **c** θ *<sup>c</sup>* = 10°, *KMAX* = 2, and *p* = 6. In **a** and **b** the spherical cap is centred on the north magnetic pole, and in **c** the spherical cap is indicated by a grey solid line. Illustration is in the magnetic latitude/MLT coordinate system

$$\begin{split} A &= \int\_0^{2\pi} ds = \int\_0^{2\pi} \int\_0^{\theta\_0} r^2 \sin\theta d\theta d\varphi \\ &= 2\pi r^2 (1 - \cos\theta\_0) \\ &= 4\pi r^2 \sin^2\frac{\theta\_0}{2} \\ &= 4\pi \left(\frac{180^\circ}{\pi}\right)^2 \sin^2\frac{\theta\_c}{2} \end{split} \tag{9.28}$$

If there are *n* data points uniformly distributed about *A*, then the spacing, θ, between points is given by

$$
\Delta\theta = \sqrt{\frac{A}{n}},\tag{9.29}
$$

and the minimum wavelength (λmin) that can be mapped is given by (e.g. Bullard 1967)

$$\begin{split} \lambda\_{\text{min}} &= p \Delta \theta \\ &= \frac{360^{\circ}}{n\_{k}(m)} \end{split} \tag{9.30}$$

where p is the minimum number of points that must be present on a given wavelength. By rearranging Eq. (9.30), *nk* (*m*) can be written as

$$\begin{split} n\_k(m) &= \frac{360^\circ}{p\Delta\theta} \\ &= \frac{360^\circ}{p} \sqrt{\frac{n}{A}} \\ &= \frac{2\pi r}{p} \sqrt{\frac{n}{4\pi r^2 \sin^2\frac{\theta\_0}{2}}} \\ &= \frac{\sqrt{\pi n}}{p \sin\frac{\theta\_0}{2}}. \end{split} \tag{9.31}$$

Equation (9.31) defines the maximum values of *nk (m)* appropriate for mapping a given data set assuming a uniform distribution of points.

For the case of small *m, nk (m)* can be approximated by (Haines 1988)

$$m\_k(m) \approx \frac{k90^{\circ}}{\theta\_c}.\tag{9.32}$$

Equation (9.31) can then be used to approximate the maximum degree-index as

9 Spherical Cap Harmonic Analysis Techniques for Mapping … 205

$$K\_{MAX} = \frac{4}{p} \sqrt{\frac{n}{\pi}}.\tag{9.33}$$

In reality, data are not distributed uniformly and it is necessary to select a lower value of *Kmax*, and therefore *nk (m)*. Too large a value results in the formation of unphysical flow structures in the mapped field in regions of the spherical cap where the λmin defined by *Kmax* is less than the θ characterizing measurement locations within that specific region of the spherical cap. Such problems are more pronounced along the edges of the spherical cap as the solution is only bounded on one 'side' (i.e. within the spherical cap). However, this problem can be improved by including data from outside the spherical cap where possible, even though the solution to *<sup>E</sup>* is limited to within the spherical cap. Alternatively, or additionally, a coefficient rejection procedure can be implemented to reject any insignificant coefficients which, although they may add a slight reduction to the minimization performed in the fitting procedure, do not contribute to the solution (Haines and Fiori 2013).

Green (2006) points out that the choice of *Kmax* must also necessarily depend on the number of available observations within the spherical cap. As *Kmax* increases (i.e. the spatial resolution decreases), the number of basis functions (i.e. the number of coefficients in the series expansion) increases. Care should be taken to ensure that the number of basis functions does not exceed the number of data points. Green (2006) provides an excellent example in their Sect. 9.3.

From Eq. (9.33) it is obvious that determination of the appropriate maximum value of *nk (m)* to map a given data set is dependent on an appropriate choice of *p*. As noted above, *p* = 2 is the absolute minimum, requiring 2 observations per wavelength. However, for a real data set with irregularly spaced observations most likely contaminated by noise, a larger value of *p* is more practical to ensure a good fit (Fiori et al. 2010). Bendat (1958) actually suggests that *p* should be chosen in the range of 10–20, but the author's experience suggest that such a dense distribution of measurements is an unreasonable assumption. In practice the choice of *p* is largely dependent on the typical distribution and quality of data for a given data set and is up to the discretion of the user. As an example, Fiori et al. (2010) suggest *p* = 6 for mapping two-dimensional plasma flow based on a data set smoothly distributed throughout the spherical cap. In examining satellite data confined to a narrow track transecting the spherical cap, Fiori et al. (2014) provide recommendations for selecting *p* based on evaluating the goodness-of-fit of the resultant map.

An additional consideration is the location and size of the spherical cap, given by θ*c*. Ideally, the cap should be placed either centred over the data set having a radius which fully encompasses the data, or over a specific geographic (or geomagnetic) region of interest. Errors associated with cap-size may arise when attempting to map wavelengths that are larger than the spherical cap, for example, when mapping secular variation. Such a problem is commonly encountered when mapping magnetic field data and is solved by subtracting a known reference field, such as the IGRF, to represent the residual value, and adding the reference field back after solving for the fitting coefficients. Haines (1985a) proposes this as a solution to a poor fit at the boundary of the spherical cap due to not including data outside of the spherical cap.

#### **9.3 Practical Example of Mapping Ionospheric Plasma Flow Using SHA and SCHA Approaches**

The Super Dual Auroral Radar Network (SuperDARN) (Greenwald et al. 1995; Chisham et al. 2007) samples the high-latitude plasma flow pattern every 1-2 min in both the northern and southern hemisphere. The advantage of the large spread of measurements has been modelling of the ionospheric electric potential. Restriction of measurements to only roughly 1/3 of the sphere leads to ill-conditioning of the leastsquares modelling of the potential using SHA approaches. Ruohoniemi and Baker (1998) instead take the approach of redistributing the data in co-latitude across a larger region in a method commonly referred to as either the map potential or FIT technique. The FIT technique redistributes the data located poleward of a predefined boundary (θFIT) across the entire sphere opposed to the hemisphere. The redistribution is performed by scaling co-latitude by α = 180◦/θ*FIT* to effectively 'stretch' the data across the sphere. Following redistribution, SHA is applied in the usual way by representing the electrostatic potential as a series expansion built with Associated Legendre Functions of integer degree and order. Once the electrostatic potential is modelled, the model is carefully 'un-stretched' from the entire sphere to the highlatitude region poleward of θFIT.

Redistribution of data across the entire sphere forces the boundary defined by θFIT to a single point at θ = 180° thereby forcing the electrostatic potential at θFIT to zero. Although undesirable for general mapping, the consequence to mapping the electrostatic potential does have physical significance, as the ionospheric convection pattern is expected to close at some low-latitude boundary (i.e. Hairston and Heelis 1990; Heppner and Maynard 1987; Rich and Hairson, 1994; Weimer 1995; Fiori et al. 2016). In the FIT technique, the low-latitude boundary given by θFIT is determined based on the distribution of low-velocity (<100 m/s) measured velocities. Zeroing the electrostatic potential at θFIT would be difficult to achieve using SHA without redistribution of the data and would require a very-high-order fit. In contrast, the SCHA algorithm could be used with only the k-m = odd basis functions, which has been investigated by Weimer et al. (2005) and Fiori et al. (2010, 2013).

Care should be used when applying the FIT technique in general. It is not clear what impact the changing data density caused by the redistribution of data has on the overall fit. For example, consider two pairs of observations separated by 111 km located at θ<sup>1</sup> = 29.5° and θ<sup>2</sup> = 10° for θFIT = 30°. In the transformed regime, θ<sup>1</sup> - = 177° and θ<sup>2</sup> - = 60°, and the pair of points are now separated by 12 km and >500 km, respectively. The impact of this redistribution is unclear. Mapping electrostatic potential through the FIT technique also requires the exact location of the low-latitude boundary of the ionospheric convection pattern, defined by θFIT to be known. For example, Fiori (2011) show that one consequence of a boundary located too far equatorward is the reduction of the calculated cross-polar cap potential (CPCP), which is a significant finding as SuperDARN electrostatic potential maps are shown to consistently report lower CPCP compared to values determined by other means using other instruments (i.e. Shepherd 2007).

Fiori et al. (2010) point out that constraints applied by the FIT technique to ensure a smooth fit over the high-latitude region, namely forcing a low-latitude boundary to the flow and the use of a statistical convection model (see Sect. 9.5), potentially distort the mapped convection. They apply an SCHA algorithm for mapping ionospheric convection based on SuperDARN measurements using the techniques discussed above and examine the impacts on both global-scale and regional-scale mapping. The SCHA technique is shown to perform comparably to the FIT technique over regions of good data coverage, and to provide a better solution over regions of highly variable flow, particularly near the low-latitude boundary of the convection zone near θFIT.

#### **9.4 Application of SCHA for Mapping Ionospheric Plasma Flow Based on Swarm Satellite Ion Drift Measurements**

An application of the SCHA algorithm discussed in this Chapter is the mapping of ionospheric data from multi-spacecraft missions. Satellite data from various missions have been incorporated into single or multi-instrument data sets using techniques not limited to SCHA for magnetic field mapping (e.g. An et al. 1998; AMIE; Green et al. 2007; Haines 1985b; Kotzé 2002), or for mapping statistical convection patterns based on satellite data collected over a long period of time (e.g. Heppner and Maynard 1987; Rich and Maynard 1989; Hairston and Heelis 1990; Rich and Hairston 1994; Papitashvili et al. 1999; Papitashvili and Rich 2002; Weimer 1995). However, as of the time of writing, the authors of this Chapter are not aware of the application of SCHA to the instantaneous mapping of satellite data. However, extensive work has been performed in Fiori et al. (2013, 2014) to evaluate the potential for mapping ionospheric plasma flow based on artificially generated Swarm ion drift measurements. This Section provides a summary of this work.

One objective of the Swarm mission is an investigation of electric currents flowing in the Earth's ionosphere which drive ionospheric electric fields and ion drift. Each satellite is therefore equipped with an electric field instrument (EFI) capable of measuring the three-dimensional ion drift (Knudsen et al. 2003, 2017). Ion drift measured in the along-track and cross-track directions by two orthogonal bi-dimensional Thermal Ion Imagers (TIIs) are combined to create a two-dimensional ion drift vectors describing plasma flowing in the horizontal plane. Such vectors may be processed using the SCHA mapping procedures described in this Chapter to map the twodimensional plasma flow. The algorithm is fairly straightforward and is described as follows. (1) For the interval of interest, combine along-track and cross-track ion drift measurements to create two-dimensional vectors of the horizontal ion drift. If desirable these ion drift vectors can be mapped down from the satellite altitude to the F-region peak of the ionosphere (e.g. Walker and Sofko 2015). For example, ion drifts mapped down from 470 km to 300 km, would be reduced by approximately 4%. (2) Build the matrices described by Eq. (9.13). If there are *N* measurements, then matrix *F* represents a 2 *N* × 1 matrix of the measured ion drift in the θ (*v*<sup>θ</sup> ) and ϕ (*v*ϕ) directions. Note that each ion drift component is calculated from the twodimensional ion drift and then considered separately. *CD* is a 2 *N* × (*Kmax* + 1)<sup>2</sup> matrix derived from Eqs. (9.9) and (9.10) for *v*<sup>θ</sup> and *v*ϕ, and *AB* is a (*Kmax* + 1)<sup>2</sup> × 1 matrix composed of the *Akm* and *Bkm* fitting coefficients. (3) The fitting coefficients are solved as described in Sect. 9.2.2. (4) Fitting coefficients are used to calculate the two-dimensional ion drift within the spherical cap where adequately constrained by data. This SCHA algorithm has been applied to (1) examine convection mapping on a global-scale spanning the entire convection zone in the high-latitude region based on both a Swarm-only data set and a multi-instrument data set (Fiori et al. 2013), and (2) explore the possibility of convection mapping with a Swarm-only data set over a localized region centred over the satellites' orbit, and determine parameters for evaluating the appropriate size of the localized region and the goodness-of-fit that can be obtained (Fiori et al. 2014). Note that due to the unavailability of definitive Swarm data at the time of publication, Fiori et al. (2013, 2014) rely on artificially generated ion drifts emulated from a statistical model along hypothetical Swarm satellite tracks. For a full description of the emulation technique, see either paper. The SCHA techniques described here could equivalently be applied to hypothetical or real data.

Fiori et al. (2013) examined SCHA for mapping the high-latitude ionospheric convection pattern based on a Swarm data set using a global-scale spherical cap. In order to map the entire high-latitude convection patter, the spherical cap was centred at the north magnetic pole extending to ~60° MLAT (θ<sup>c</sup> = 30°). They show that even when all three satellites cross the spherical cap at the same approximate time, the Swarm data coverage is insufficient to support such large-scale convection mapping based on a Swarm-only data set; the RoC does not exceed the measurement location. However, they do show that inclusion of the Swarm data with a complementary data set has the potential to increase the region of the spherical cap populated by data thereby improving resultant convection maps.

In the case of Fiori et al. (2013), the complementary data set under consideration was the SuperDARN (Greenwald et al. 1995; Chisham et al. 2007) which samples the high-latitude plasma flow pattern continuously every 1–2 min in both the northern and southern hemisphere. Currently, SuperDARN consists of radars with a collective field of view that covers the majority of the high-latitude region (within ~30° of either magnetic pole), extending to mid-latitudes. However, despite this superior coverage, there are frequently gaps in the data due to a variety of circumstances including geomagnetic activity, solar activity, season, and magnetic local time (MLT) of observation (Ruohoniemi and Greenwald 1996; Danskin et al. 2002; Shepherd et al. 2002; Koustov et al. 2004; Shepherd 2007). It is common practice within the SuperDARN community to supplement real data with data from a statistical convection model to produce a well-constrained and smooth map (Ruohoniemi and Baker 1998; Shepherd and Ruohoniemi 2000; Ruohoniemi and Greenwald 2005; Cousins and Shepherd 2010; Pettigrew et al. 2010). Statistical models are useful in expanding convection data so that the RoC covers the entire high-latitude region during periods of stable convection but do not always accurately represent the true plasma flow during more turbulent conditions (e.g. Fiori et al. 2010). In such instances, it is more useful to include data from other instruments, such as is available from the Swarm satellites to fill in gaps in the SuperDARN data set in order to increase the overall RoC that can be mapped.

Examples presented in Fiori et al. (2013) demonstrate how effectively data from the Swarm satellite can increase the RoC for global-scale convection mapping with SuperDARN. For instance, Fiori et al. (2013) present one case (reproduced in Fig. 9.5) where despite excellent SuperDARN coverage the RoC is insufficient to determine the cross-polar cap potential (CPCP) without adding contributions from a statistical model or other instruments. Figure 9.5a illustrates hypothetical measurement locations for the SuperDARN and Swarm instruments. In Fig. 9.5b, convection is mapped based on a global-scale spherical cap (θ<sup>c</sup> = 32°) centred at the north magnetic pole for *Kmax* = 6 and *p* = 6. Black contours indicate the electrostatic potential associate with the convection pattern, and the plus sign indicates the maximum potential. Because the minimum potential cannot be mapped, it is not possible to determine the CPCP. Addition of data from a single Swarm satellite (here data has been artificially generated based on a statistical model based on solar wind and IMF conditions), slightly enlarges the RoC. Inclusion of hypothetical data from all three Swarm satellites increases the RoC sufficiently to allow determination of the CPCP.

Fiori et al. (2013) show that the addition of Swarm ion drift data to a SuperDARN data set can cause an average relative increase in the RoC of 12.5% for periods when 1–3 Swarm satellites are located in the high-latitude region, with the greatest contributions being made when Swarm data are located at auroral latitudes.

In addition to supplementing other data sets for the purpose of global-scale convection mapping, Swarm data may also be considered independently for mapping more localized regions. For localized convection mapping data may be considered for three different satellite configurations: a single satellite track (upper satellite Swarm B), a double satellite track (lower satellite pair Swarm A/C), or for the case where all three satellites overlap. As an example, consider a regional spherical cap with θ<sup>c</sup> = 10°, which corresponds to roughly 4-min of observation processed with *Kmax* = 3. Figure 9.6 shows two statistical convection models within such a spherical cap centred at 80° magnetic latitude and 12 MLT. Figure 9.6c indicates the location of sample satellite tracks. Artificially generated Swarm ion drift vectors generated from these statistical models along the satellite paths illustrated in Fig. 9.6c can be used to map convection on the localized spherical cap. Such maps are shown in Fig. 9.7 where the RoC is described by *p* = 6. Maps in Fig. 9.7a, b can be compared to the original statistical model in Fig. 9.6a, b (and mapped using black vectors in Fig. 9.7), respectively. The visual comparison indicates good overall agreement between the mapped vectors and the statistical model, with the greatest difference occurring at the outlying region of the RoC, suggesting the RoC could be reduced by increasing the value of *p*.

Although a three-satellite map offers more data for convection mapping, and is, therefore, more likely to produce a more accurate map, it is not a frequent occurrence. In the next example, consider the case where only the lower satellite pair is available for mapping. In Fig. 9.8a the IMF Bz > 0 convection pattern illustrated in Fig. 9.6b is

**Fig. 9.5 a** Gridded coordinate locations for the Swarm A (blue), Swarm C (green), Swarm B (red), and SuperDARN (black) data sets. SuperDARN data are based on gridded measurement coordinates for 20 January 2011 11:20–11:24 UT. **b**–**d** Convection calculated for θ<sup>c</sup> = 32°, *Kmax* = 6, *p* = 6 based on artificial measurements generated from a statistical mode at possible coordinates of **b** SuperDARN, **c** SuperDARN and Swarm B, and **d** SuperDARN and Swarm A, B and C measurements shown in (**a**). Convection vectors are plotted along an evenly spaced grid of size 1° of magnetic latitude throughout the RoC. Contours are plotted with a 6 kV contour spacing. The plus and cross indicate the maximum and minimum potential and the CPCP is indicated where it can be determined. (Reproduction of Fiori et al. 2013, Fig. 9.6.)

**Fig. 9.6** Portion of two sample statistical model convection patterns for the case of **a** IMF Bz< 0 and **b** IMF Bz> 0 within a spherical cap centred at 80° MLAT and 12 MLT with ϑ*<sup>c</sup>* = 10◦. The Cousins and Shepherd (2010) models employed in Fiori et al. (2013, 2014) are used. **c** shows a 4-min segment of the Swarm A (blue), C (green), and B (red) satellite orbits

**Fig. 9.7** Coloured vectors indicate convection determined based on Swarm measurements artificially generated along the satellite tracks shown in Fig. 9.6c based on the statistical models mapped in **a** Figs. 9.6a, b for a spherical cap centred at 80° MLAT and 12 MLT with θ *<sup>c</sup>* = 10° and *K* = 3. The RoC is defined with *p* = 6

used to artificially generate Swarm ion drifts for the Swarm A/C satellite pair which follow a dawn/dusk orientation. Artificially generated ion drifts are processed to create the convection pattern shown where the RoC is limited by *p* = 12. In Fig. 9.8b a similar map is generated for a noon/midnight satellite orientation using the IMF Bz < 0 statistical convection model shown in Fig. 9.6a. Both maps indicate good agreement between the mapped convection and the statistical convection model.

Fiori et al. (2014) examined the use of SCHA for mapping convection using a localized spherical cap centred over the lower Swarm satellite pair. Convection maps

**Fig. 9.8** Coloured vectors indicate convection calculated for θ *<sup>c</sup>* = 10° and *Kmax* = 3 for a spherical cap centred at 70° MLAT and 12 MLT for Swarm A and C measurements artificially generated using the **a** IMF Bz> 0 convection pattern for a Swarm A satellite track offset 80° from the 12 MLT meridian and *p* = 12 and **b** IMF Bz<0 convection pattern for a Swarm A satellite track offset 170° from the 12 MLT meridian and *p* = 12. Black vectors in (**a**) and (**b**) indicate the direction of the statistical model convection

were created using the SCHA technique based on artificially generated Swarm measurements within a localized spherical cap having an angular radius, θ *<sup>c</sup>*, ranging from 4° to 15°. These maps used a *KMAX* that was calculated so that the SCHA mapping resolution matched that of the statistical model used to generate measurements. Within the study, the artificially generated Swarm data were used to generate 2502 maps for a variety of satellite orientations with RoC determined for varying values of *p*. Comparisons were made between the statistical model used to create the artificially generated measurements and the convection mapped using SCHA. These comparisons were used to categorize the maps as being well constrained or poorly constrained. Based on this analysis it was determined that the SCHA algorithm ideally maps Swarm data for θ *<sup>c</sup>*, = 10°, and was unreliable for θ *<sup>c</sup>*, ≥ 14°. Fiori et al. (2014) determined that Swarm-based measurements could successfully map convection over a localized region surrounding the satellite track to examine smallscale features in the flow, and that larger scale convection mapping can be achieved through combining multiple spherical caps.

#### **9.5 Summary**

This Chapter provides detailed explanations of how to map ionospheric plasma flows in the high-latitude region based on hypothetical ion drifts from the Swarm satellites using spherical cap harmonic analysis. SCHA is a useful technique for mapping data limited in space, and is closely related to spherical harmonic analysis which is commonly applied on a larger global scale. Although originally developed for mapping the three-dimensional magnetic field on a regional scale, SCHA techniques have been used to map a variety of one-, two- and three-dimensional properties. This Chapter provides a summary of the relatively new application of SCHA to two-dimensional mapping of ion drift.

Section 9.2 provides the necessary theory for applying SCHA to mapping ionospheric plasma flow, or convection. Equations for mapping electrostatic potential, electric field, and plasma flow in the ionosphere are provided along with discussion regarding the relevant basis functions and boundary conditions considered in the solution. Practical considerations such as the appropriate cap-size to use for a given data set and how to determine where the solution is adequately constrained by measurements are discussed. Section 9.3 provides an example of SHA and SCHA mapping of ionospheric convection based on the SuperDARN data. In Sect. 9.4 the theory developed in this Chapter is applied to the mapping of ionospheric ion drift measured by the Swarm satellites. Sample convection maps are presented which illustrate SCHA mapping of Swarm ion drift data (1) on a global scale in combination with measurements from additional ground-based instruments and (2) on a local scale using satellite measurements alone.

This Chapter clearly demonstrates the application of spherical-cap harmonic analysis techniques for mapping plasma flow patterns based on ion drift measurements from the Swarm satellite.

**Acknowledgements** This work was supported by the Natural Resources Canada (NRCan) Lands and Minerals Sector and Public Safety Geosciences program, and the Canadian Space Agency. The author wishes to thank G. V. Haines for communications regarding spherical cap harmonic analysis and D. Knudsen and J. Burchill for communications regarding Swarm. Useful contributions and discussion from David Boteler were appreciated. The author acknowledges 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. SuperDARN data are available from the SuperDARN website hosted by Virginia Tech (http://vt.superdarn.org). We thank the International Space Science Institute (ISSI) in Bern, Switzerland for supporting the Working Group "Multi Satellite Analysis Tools—Ionosphere" from which this chapter resulted. The Editors thank David Green for his assistance in evaluating this chapter.

#### **References**

Bendat, J.S. 1958. *Principles and applications of random noise theory*. New York: Wiley.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 10 Recent Progress on Inverse and Data Assimilation Procedure for High-Latitude Ionospheric Electrodynamics**

#### **Tomoko Matsuo**

**Abstract** Polar ionospheric electrodynamics plays an important role in the Sun– Earth connection chain, acting as one of the major driving forces of the upper atmosphere and providing us with a means to probe physical processes in the distant magnetosphere. Accurate specification of the constantly changing conditions of high-latitude ionospheric electrodynamics has long been of paramount interest to the geospace science community. The Assimilative Mapping of Ionospheric Electrodynamics procedure, developed with an emphasis on inverting ground-based magnetometer observations for historical reasons, has long been used in the geospace science community as a way to obtain complete maps of high-latitude ionospheric electrodynamics by overcoming the limitations of a given geospace monitoring system. This Chapter presents recent technical progress on inverse and data assimilation procedures motivated primarily by availability of regular monitoring of high-latitude electrodynamics by space-borne instruments. The method overview describes how electrodynamic state variables are represented with polar-cap spherical harmonics and how coefficients are estimated from the point of view of the Bayesian inferential framework. Some examples of the recent applications to analysis of SuperDARN plasma drift, Iridium, and DMSP magnetic fields, as well as DMSP auroral particle precipitation data are included to demonstrate the method.

#### **10.1 Introduction**

The most dynamic electromagnetic energy and momentum exchange processes between the upper atmosphere and the magnetosphere take place in the polar ionosphere. Physical processes producing aurora involve ionization and excitation of atmospheric constituents due to energetic charged particles precipitating into the upper atmosphere from the magnetosphere along the geomagnetic field lines, which

© The Author(s) 2020 M. W. Dunlop and H. Lühr (eds.), *Ionospheric Multi-Spacecraft Analysis Tools*, ISSI Scientific Report Series 17, https://doi.org/10.1007/978-3-030-26732-2\_10

T. Matsuo (B)

Ann and H.J. Smead Department of Aerospace Engineering Sciences, University of Colorado, Boulder, CO 80303-0429, USA e-mail<EMAIL_ADDRESS>
in turn modulates the ionosphere's ability to conduct electric currents. Polar ionospheric electrodynamics plays an important role in the Sun–Earth connection chain, acting as one of the major driving forces of the upper atmosphere and providing us with a means to probe physical processes in the distant magnetosphere. Accurate specification of the constantly changing conditions of high-latitude ionospheric electrodynamics has long been of paramount interest to the geospace science community.

Global monitoring of high-latitude geospace has dramatically improved thanks to a recent expansion of ground-based and space-based observing capability. International consortiums of ground-based instrumentation such as the Super Dual Auroral Radar Network (SuperDARN) (e.g., Greenwald et al. 1995), International Real-Time Magnetic Observatory Network (e.g., Love 2013) and SuperMAG (e.g., Gjerloev 2009) have made a large volume of quality-controlled, standardized data accessible to the public. Acquisition, processing, and distribution of engineering-grade magnetometer data from the Iridium satellite constellation for scientific purposes by the ActiveMagnetosphere and Polar Electrodynamics Response Experiment (AMPERE) program (Anderson et al. 2000) have been instrumental in making continuous, global monitoring of geomagnetic-field-aligned currents (FAC) possible. Defense Meteorological Satellite Program (DMSP) space environment instruments have long been providing valuable measurements of precipitating electron and ion particles, magnetic fields, and ultraviolet spectrographic images (e.g., Rich 1984; Hardy et al. 1984; Paxton et al. 2002). And the Swarm multi-satellite mission (Friis-Christensen et al. 2006) provides high precision measurements of magnetic fields that complement theses existing geospace observing systems.

Data assimilation techniques such as the Assimilative Mapping of Ionospheric Electrodynamics (AMIE) procedure of Richmond and Kamide (1988) have long been used in the geospace science community as a way to obtain complete maps of high-latitude ionospheric electrodynamics by overcoming the limitations of a given geospace monitoring system. The procedure combines a number of different types of space-based and ground-based observations with an empirical model of ionospheric electrodynamics to infer distributions of ionospheric electric fields and currents, FAC, associated geomagnetic perturbation fields at both ground and low-Earth-orbit altitudes, Hall and Pedersen conductance, and Joule heating. AMIE maps have yielded a number of important insights into the coupling of the magnetosphere, ionosphere, and thermosphere that takes place at high latitudes. Lu (2017) provides a comprehensive overview of AMIE applications.

This paper presents an overview of the recent technical developments of the inverse and data assimilation procedure for high-latitude electrodynamics. Some of these developments are a consequence of a reformulation of the best linear unbiased estimation problem presented in Richmond and Kamide (1988) as a Bayesian estimation problem (Matsuo et al. 2005). Under the assumption that electrodynamic variables are Gaussian distributed, these two estimation problems are equivalent. A Bayesian perspective has helped to clarify the role of the prior model (background) error covariance as a key component in the modeling of Gaussian processes, and thus guided modeling and estimation of prior covariance functions from a large volume of SuperDARN data (Cousins et al. 2013a), DMSP particle precipitation data (McGranaghan et al. 2015, 2016), and Iridium magnetic perturbation data (Cousins et al. 2015b; Shi et al. 2019). Even though ionospheric conductivity serves as a critical linkage in electromagnetic energy and momentum exchange processes, direct monitoring of this conductivity is almost nonexistent. Another notable development led by McGranaghan et al. (2016) is an assimilative mapping of the conductance using the auroral ionization derived from DMSP electron energy flux spectra with help of the GLobal airglOW (GLOW) model (Solomon et al. 1988) without the assumption of Maxwellian distribution. Since the AMIE has been developed with an emphasis on inverting ground-based magnetometer observations for historical reasons (Kamide et al. 1981; Richmond and Kamide 1988), it is not tailored to analyses of space-based magnetometer data from DMSP, Iridium, and Swarm. In order to solve the optimization problem in terms of electrostatic potential, the space-based magnetometer data first need to be converted to electrostatic potential through the application of Ohm's law and current continuity. To minimize the impact of conductance on the inversion of space-based magnetometer data for FAC, the optimization problem is now being solved in terms of both magnetic potential and electrostatic potential (Matsuo et al. 2015; Cousins et al. 2015a).

#### **10.2 Method Overview**

#### *10.2.1 Representation of Electrodynamic State Variables Using Scalar and Vector Polar-Cap Spherical Harmonic Basis Functions*

The ionosphere is treated as a thin conductive slab centered at a reference height *hr* = 110 km, and the current above the ionosphere is assumed to be strictly radial. The effect of the neutral wind dynamo is not considered. Electrodynamic variables analyzed here include the electrostatic potential -, electric fields **E**, Pedersen and Hall conductance (height-integrated conductivity) *p*, *h*, height-integrated horizontal ionospheric current density **J**⊥, toroidal magnetic potential associated with field-aligned current density **J**, and equivalent current potential associated with ground-based magnetic fields. These variables are presumed to be related to each other as follows.

$$\mathbf{E} = -\nabla \Phi \tag{10.1}$$

$$\mathbf{J}\_{\perp} = \boldsymbol{\Sigma} \cdot \mathbf{E} \tag{10.2}$$

$$\mathbf{J}\_{\parallel} = -\nabla \cdot \mathbf{J}\_{\perp} \tag{10.3}$$

$$\nabla \times \mathbf{J}\_{\perp} = -\nabla\_{\text{hor}}^{2} \Omega \tag{10.4}$$

$$\mathbf{J}\_{\parallel} = \frac{1}{\mu\_o} \nabla\_{\text{hor}}^2 \Psi \tag{10.5}$$

where *-* = *<sup>p</sup>* −*<sup>h</sup> <sup>h</sup> <sup>p</sup>* is the conductance tensor, ∇<sup>2</sup> hor is the horizontal Laplacian, and μ*<sup>o</sup>* is permeability of free space. Equation (10.4) results from the assumption of strictly vertical **J** that allows equating the curls of **J**<sup>⊥</sup> and the equivalent current (i.e., Fukushima Theorem). If the Pedersen and Hall conductances are given, the relationship among all electrodynamic variables (10.1)–(10.5) becomes linear.

In the procedure, electrodynamic variables are expressed in terms of the polarcap spherical harmonic basis functions developed by Richmond and Kamide (1988). Suppose that represents a matrix of the polar-cap spherical harmonic basis functions evaluated at discrete grid locations specified by the Modified Magnetic Apex longitude φ*<sup>m</sup>* and latitude λ*<sup>m</sup>* at the altitude of *hr* (Richmond 1995) and that **x** denotes a vector of the coefficients. is furthermore given by a set of 244 polar-cap spherical harmonic basis functions up to order *m* = 12, with non-integer degrees *n* up to a maximum of *n* = 72.6 for *m* = 0, with a polar-cap co-latitude for the functions of 40◦. Therefore, **x** is a column vector of 244 elements and is an *n* × 244 matrix, where *n* is the number of grid points. Using the Nyquist sampling rate, the effective resolution is 15◦ longitude and 2.5◦ latitude. Let's suppose that the electrostatic potential at φ*<sup>m</sup>* and λ*<sup>m</sup>* is given by

$$\Phi(\phi\_m, \lambda\_m) = \Psi \mathbf{x}\_{\mathcal{E}} + \epsilon\_l \ , \tag{10.6}$$

where *<sup>t</sup>* is the truncation error, and the electric fields **E** by

$$\mathbf{E}(\phi\_m, \lambda\_m) = -\Psi' \mathbf{x}\_{\mathcal{E}} + \boldsymbol{\epsilon}\_l \ , \tag{10.7}$$

where a (*n* × 244) matrix contains the gradients of the polar-cap spherical harmonic basis functions, which discretizes (10.1). The toroidal potential at φ*<sup>m</sup>* and λ*<sup>m</sup>* is then given by

$$
\Psi(\phi\_m, \lambda\_m) = \Psi \mathbf{x}\_{\mathcal{M}} + \epsilon\_{\iota} \ , \tag{10.8}
$$

where *<sup>t</sup>* is the truncation error, and the FAC magnitude *J* by

$$J(\phi\_m, \lambda\_m) = \Psi^\prime \mathbf{x}\_{\mathsf{M}} + \epsilon\_{\iota} \ , \tag{10.9}$$

where a (*n* × 244) matrix contains a simplified evaluation of (10.5) using the analytical expression of the horizontal Laplacian of polar-cap spherical harmonic basis functions applicable to spherical coordinates, rather than the full expression applicable to M(110) coordinates. As explained in (Matsuo et al. 2015), this computational simplification introduces errors on the order of 10%. For a given *<sup>p</sup>* and *h*, **x**<sup>E</sup> and **x**<sup>M</sup> are related linearly through the current continuity and Ohm's law (10.2)–(10.3).

#### *10.2.2 Bayesian State Estimation for Gaussian Processes*

Suppose that **y** represents a vector of *j* observations that may consist of electric field, ground-based magnetic field, and/or space-based magnetic field measurements at discrete observation locations. By evaluating the polar-cap spherical harmonics and their derivatives at observation locations, **y** can be expressed as

$$\mathbf{y} = \mathbf{H}\mathbf{x} + \boldsymbol{\epsilon}\_r \; , \tag{10.10}$$

where **H** is a (*j* × 244) matrix that contains the polar-cap spherical harmonic basis functions and their spatial derivatives with corresponding vector calculus operations as specified in (10.1)–(10.5), **x** denotes a vector of the 244 coefficients, and *<sup>r</sup>* is the sum of observational and truncation errors. The objective of the Bayesian state estimation is to infer the polar-cap spherical harmonics coefficients **x** given observations **y** according to Bayes rule: [**x**|**y**]∝[**y**|**x**][**x**] The vectors **x** and **y** are herein assumed to be distributed according to the multivariate normal distribution denoted by MN as

$$\mathbf{x} \sim \mathcal{MN}[\mathbf{x}\_b, \mathbf{C}\_b] \,, \tag{10.11}$$

$$\mathbf{y} \sim \mathcal{MN}[\mathbf{Hx}, \mathbf{C}\_r] \,, \tag{10.12}$$

where **x***<sup>b</sup>* is the prior mean, **C***<sup>b</sup>* is the prior (background) model error covariance < (**x***<sup>b</sup>* − **x**)(**x***<sup>b</sup>* − **x**)*<sup>T</sup>* >, and **C***<sup>r</sup>* is the observational error covariance < *<sup>r</sup> <sup>T</sup> <sup>r</sup>* >. **x***<sup>b</sup>* is specified by using an empirical model. **C***<sup>b</sup>* is described in the following section. The errors *<sup>r</sup>* are assumed to be uncorrelated, so **C***<sup>r</sup>* is given by a diagonal matrix of the variance of observational error. The posterior distribution or the conditional distribution of **x** given observations **y** is given by the multivariate normal distribution as

$$\mathbf{[x|y]} \sim \mathcal{MN}[\mathbf{x}\_a, \mathbf{C}\_a] \,, \tag{10.13}$$

where **x***<sup>a</sup>* is the posterior mean or the data assimilation analysis and **C***<sup>a</sup>* is the analysis error covariance < (**x***<sup>a</sup>* − **x**)(**x***<sup>a</sup>* − **x**)*<sup>T</sup>* >. In the case of normally distributed **x** and **y** and linear **H**, there are closed formulae for **x***<sup>a</sup>* and **C***<sup>a</sup>* (e.g., Jazwinski 1970; Lorenc 1986):

$$\mathbf{x}\_a = \mathbf{x}\_b + \mathbf{C}\_b \mathbf{H}^T (\mathbf{H} \mathbf{C}\_b \mathbf{H}^T + \mathbf{C}\_r)^{-1} (\mathbf{y} - \mathbf{H} \mathbf{x}\_b) \,, \tag{10.14}$$

$$\mathbf{C}\_{a} = [\mathbf{I} - \mathbf{C}\_{b}\mathbf{H}^{T}(\mathbf{H}\mathbf{C}\_{b}\mathbf{H}^{T} + \mathbf{C}\_{r})^{-1}\mathbf{H}]\mathbf{C}\_{b} \,. \tag{10.15}$$

By specifying **C***b*, **C***r*, **H**, and **x***b*, the analysis **x***<sup>a</sup>* and error covariance **C***<sup>a</sup>* can be computed for given observations **y**. The prior model error covariance **C***<sup>b</sup>* plays an important role here, not only balancing the weighting between observations and the prior model but also spreading the observation-model discrepancy information spatially according to the correlation represented in the covariance.

#### *10.2.3 Nonstationary Covariance Modeling*

Following the approach adopted in Matsuo et al. (2005) as a way to incorporate anisotropic and inhomogeneous characteristics of the prior (background) model errors into the analysis (10.14) in a computationally tractable manner, **C***<sup>b</sup>* is modeled using the empirical orthogonal functions (EOFs, i.e., principal components). EOFs and their coefficients are estimated in advance of the data assimilation, for instance, from 50 million total SuperDARN plasma drift data points over January 2011 through August 2012 for electrostatic potential (Cousins et al. 2013a), from over 60 million DMSP electron energy flux spectra during the solar cycles 22 and 24 for conductance (McGranaghan et al. 2015), and from over 300 days of Iridium magnetic perturbation data from 2010 to 2015 for field-aligned currents (Shi et al. 2019).

Since observation sampling is often irregular and incomplete, a straightforward eigenvalue decomposition of sample covariance cannot be applied to the dataset. Instead, the nonlinear regression analysis of Matsuo et al. (2002) is used, wherein *p* principal components are expressed by a linear combination of the polar-cap spherical harmonic basis functions of Richmond and Kamide (1988), and each component is estimated sequentially by a back-fitting technique along with orthonormalization of the regression coefficients for each component. Each EOF can be expressed as *β*, where *β* is a 244 × *p* matrix. Then **C***<sup>b</sup>* is given as

$$\mathbf{C}\_b \approx \mathfrak{BC}\_\mathcal{V} \mathfrak{B}^T,\tag{10.16}$$

where **C**<sup>γ</sup> is the covariance < *γ γ <sup>T</sup>* > of the EOF coefficients *γ* , where *γ* is a *p* × 1 column vector. EOFs estimated by the method of Matsuo et al. (2002) are equivalent to the eigenfunctions of a covariance matrix computed from observational data.

**Fig. 10.1** Two-dimensional correlation functions of electrostatic potential with respect to the 4 points indicated by crosses. The contour lines represent the correlation level of 0.9, 0.8, 0.7, and 0.6 (Fig. 7 of Cousins et al. 2013a)

As with other principal component analysis methods, a certain replication of data samples is required to estimate **C**<sup>γ</sup> and *β* from the observations.

Figure 10.1 shows 2-dimensional correlation maps for electrostatic potential computed from the EOF-based covariance derived from SuperDARN data (Cousins et al. 2013a) where *p* is set to 30. It is evident that the correlation structures are highly anisotropic with a larger correlation length scales in the zonal direction in comparison to the meridional direction, and correlations vary depending on reference point locations. These are features of strong nonstationary correlation, which will enable the data assimilation procedure to spatially distribute the impact of observations with consideration of realistic location-specific correlation structures of SuperDARN plasma drifts or electric fields.

#### **10.3 Analysis of Electrostatic Potential and Electric Fields**

Cousins et al. (2013b) presents an inverse and data assimilation procedure designed to specifically estimate **x**<sup>E</sup> as defined in (10.6) and (10.7) from SuperDARN data. A comprehensive cross-validation study (Cousins et al. 2013b) wherein observations are systematically set aside for validation and compared to predictions by data assimilation outperforms the standard SuperDARN mapping procedure (Ruohoniemi and Baker 1998; Shepherd and Ruohoniemi 2000). The inverse and data assimilation procedure is found to reduce median prediction errors by up to 43% as compared to the standard SuperDARN mapping procedure. The procedure is built using the prior covariance modeled with EOFs obtained by Cousins et al. (2013a) and the prior mean specified by the empirical plasma convection model of Cousins and Shepherd (2010). Figure 10.2 compares the maps of electrostatic potentials obtained by the standard SuperDARN mapping procedure (Ruohoniemi and Baker 1998; Shepherd and Ruohoniemi 2000 to the ones by Cousins et al. (2013b) along with maps of the uncertainty associated with assimilative mapping as given by the diagonal elements of **C***<sup>a</sup>* (10.15). The uncertainty reflects the observation distributions with higher uncertainty found in the area of the SuperDARN data gap. The comparison also highlights the role of the nonstationary covariance in the inverse and data assimilation procedure that help regularize assimilative mapping analysis.

#### **10.4 Analysis of Toroidal Magnetic Potential and Field-Aligned Currents**

Matsuo et al. (2015) presents an inverse and data assimilation analysis of spacebased magnetometer data that directly solves for **x**<sup>M</sup> as defined in (10.8) and (10.9) to circumvent the need to use conductance in analysis of space-based magnetometer data for FAC as has been originally done in Richmond and Kamide (1988).

**Fig. 10.2** Potential distributions for three selected times. **a**, **d**, and **g** Results from the standard SuperDARN mapping procedure. **b**, **c**, **e**, **f**, **h**, and **i** Results from the inverse and data assimilation procedure, with associated uncertainty shown as background coloring on the right side (Fig. 5 of Cousins et al. 2013b)

Figure 10.3 displays maps of the toroidal magnetic potential and FAC estimated from both AMPERE and DMSP data under four distinctive interplanetary magnetic field conditions during a magnetic cloud event on May 29, 2010, and demonstrates the Interplanetary Magnetic Field (IMF) control of high-latitude electrodynamics. Note that the uncertainty associated with the magnetic potential analysis is shown in the black-and-white contour in the background, with darker shades indicating greater errors. For comparison, the bottom row shows maps of the FAC provided by the AMPERE program. The AMPERE data product obtained from the spherical harmonic fit has an effective resolution of 3◦ latitude and 36◦ longitude (Anderson et al. 2014). As discussed in Matsuo et al. (2015), the overall distribution of FAC

**Fig. 10.3** The plots in the top and middle rows are maps of toroidal magnetic potential and FAC on May 29, 2010, estimated from the AMPERE and DMSP data over a 4min interval: **a-1**, **a-2** IMF Bz positive, **b-1**, **b-2** IMF By positive, **c-1**, **c-2** IMF Bz negative, and **d-1**, **d-2** IMF By negative. The plots in the bottom row are maps of the FAC provided by the AMPERE program, estimated from the AMPERE data over a 10min interval using Altitude Adjusted Corrected Geomagnetic Coordinates (Fig. 5 of Matsuo et al. 2015)

is similar to the one obtained by the current procedure, except for a few notable differences in the detail, such as the absence of high-frequency features and more longitudinally continuous FAC spatial structures are seen in the present analysis. Thanks to the regularization through the use of the prior model error covariance in solving the inverse problem, there is no need to fill the data gap with synthetic data to make a regression analysis stable, as is required in the AMPERE inversion.

#### **10.5 Dual Optimization Approach**

The framework for the inverse and data assimilation procedure described in Sect. 10.2 has thus far been applied to assimilative analysis of individual electromagnetic variables. In this section, the same framework is applied to the analysis of multiple

**Fig. 10.4** Distributions of **a**, **d** electrostatic potential, **b**, **e** field-aligned current density, and **c**, **f** Poynting flux for 07400750 UT, November 29, 2011. Background color indicates estimated uncertainty. **a**, **b**, **c** Results using the SuperDARN and AMPERE data, independently. **d**, **e**, **f** Results using the SuperDARN and AMPERE data together. While SuperDARN observation locations are indicated by black dots, Iridium satellite tracks (where there are AMPERE data) are indicated by black lines (Fig. 5 of Cousins et al. 2015a)

variables. The relationship among electrodynamic variables given in (10.1)–(10.5) is nonlinear, requiring a nonlinear optimization approach. As an intermediate step toward implementing a fully nonlinear solver, Cousins et al. (2015a) presents a dual optimization approach by combining the two linear optimization approaches presented in Sects. 10.3 and 10.4 but using both SuperDARN and Iridium magnetic perturbation data. For a given conductance *<sup>p</sup>* and *h*, optimal values for **x**<sup>E</sup> and **x**<sup>M</sup> are estimated independently. Specifically, the optimal interpolation (or Kalman filter update) Eqs. (10.14) and (10.15) are applied to **x**<sup>E</sup> with the prior error covariance for electrostatic potential estimated from the SuperDARN data (Cousins et al. 2013a) and with **y** being composed of SuperDARN plasma drifts and Iridium magnetic perturbation fields. For estimation of **x**M, (10.14) and (10.15) are applied with the prior error covariance for toroidal magnetic potential estimated from the Iridium magnetic perturbation data (Cousins et al. 2015b).

Figure 10.4 demonstrates the benefit of incorporating both SuperDARN and Iridium magnetic perturbation observations into the estimation of both electrostatic and magnetic potential (Cousins et al. 2015a). For example, as shown in the orangeshaded background contour in Fig. 10.4a and d, the uncertainty for electrostatic potential distributions estimated from SuperDARN data alone is higher in comparison

**Fig. 10.5** Median Absolute Differences (MADs) between the prediction by data assimilation and the validation observation over the November 26–December 2, 2011 analysis time period. MADs are computed over the polar region. **a**, **d** Borovsky coupling function (black trace, left y-axis), **b**, **e** SuperDARN plasma drifts are used to predict Iridium magnetic fields, and **c**, **f** Iridium magnetic fields are used to predict SuperDARN plasma drifts (Fig. 8 of McGranaghan et al. 2016)

to the uncertainty when both data are assimilated. This is particularly evident in the dawn cell where there is no SuperDARN data but there is Iridium data. McGranaghan et al. (2016) have examined the effects of using different conductances in this dual optimization approach to assimilative mapping. When *<sup>p</sup>* and *<sup>h</sup>* are estimated by assimilation of the DMSP electron precipitation data (blue in Fig. 10.5) rather than specified by a climatological model (red in Fig. 10.5), the prediction of SuperDARN plasma drifts by assimilative analysis of Iridium magnetic perturbation data becomes more consistent with SuperDARN plasma drifts observations, as shown in Fig. 10.5c and f. Note that SuperDARN data are not used here for prediction of Iridium magnetic perturbation data, and vice versa.

#### **10.6 Summary**

This paper demonstrates that simultaneous analysis of multiple types of space-based and ground-based global geospace observations enabled by the inverse and data assimilation procedure provides a global perspective of high-latitude ionospheric electrodynamics. The paper summarizes important technical developments that have been made in response to the expansion of high-latitude geospace observing systems. The primary areas of the methodological extension to the AMIE (Richmond and Kamide 1988) are (a) the optimization in terms of both magnetic and electrostatic potential to minimize the impact of conductance on the inversion of space-based Iridium and DMSP magnetometer data for FAC mapping (Matsuo et al. 2015; Cousins et al. 2015a); (b) the use of realistic prior error covariance estimated from a large data set of SuperDARN (Cousins et al. 2013a), DMSP (McGranaghan et al. 2015) and Iridium magnetic perturbation data (Shi et al. 2019); (c) improved assimilative conductance/conductivity mapping (McGranaghan et al. 2016).

**Acknowledgements** This study is supported by the NSF grants PLR-1443703 and ICER-1541010. We thank the International Space Science Institute (ISSI) in Bern, Switzerland for supporting the Working Group "Multi Satellite Analysis Tools—Ionosphere" from which this chapter resulted. The Editors thank Gang Lu for her assistance in evaluating this chapter.

#### **References**

Kamide, Y., A.D. Richmond, and S. Matsushita. 1981. Estimation of ionospheric electric fields, ionospheric currents, and field-aligned currents from ground magnetic records. *Journal of Geophysical Research* 86 (A2): 801–813. https://doi.org/10.1029/JA086iA02p00801.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 11 Estimating Currents and Electric Fields at Low Latitudes from Satellite Magnetic Measurements**

#### **Patrick Alken**

**Abstract** Low-latitude ionospheric electric currents produce prominent signatures in the magnetic field measurements made by low Earth-orbiting satellites. Analyzing these magnetic signatures not only provides insight into the currents themselves, but also many other important and interesting phenomena in the low-latitude ionosphere and thermosphere. The low-latitude currents are modulated by thermospheric winds, so attaining a global knowledge of the spatial structure of the currents can give insight into the neutral tidal harmonics present at ionospheric altitudes. Furthermore, the equatorial electrojet (EEJ) current is driven by an equatorial electric field which in turn is generated by a dynamo process. This electric field is additionally responsible for the vertical plasma fountain and equatorial ionization anomaly at low-latitudes. Magnetic measurements of the EEJ, therefore, allows the study of low-latitude plasma motion in the E and F regions of the ionosphere. This chapter will present techniques developed for processing magnetic measurements of the EEJ to extract information about the low-latitude currents and their driving electric fields. This chapter will present a line current approach to recover the EEJ current strengths, with an emphasis on cleaning the satellite data and minimizing magnetic fields from other internal and external sources. The electric fields will be determined using a combination of physical modeling and fitting the EEJ current strengths from the satellite measurements.

#### **11.1 Introduction**

This chapter will be concerned with the calculation of ionospheric current flow and electric fields at low-latitudes, using magnetic field measurements from low Earth orbiting (LEO) satellite missions, such as Swarm (Friis-Christensen et al. 2006). In the ionosphere, neutral particles are ionized by solar extreme ultraviolet (EUV) radiation. The resulting charged plasma then interacts with the Earth's electromagnetic

© The Author(s) 2020 M. W. Dunlop and H. Lühr (eds.), *Ionospheric Multi-Spacecraft Analysis Tools*, ISSI Scientific Report Series 17, https://doi.org/10.1007/978-3-030-26732-2\_11

P. Alken (B)

University of Colorado at Boulder, Boulder, CO, USA e-mail<EMAIL_ADDRESS>
field, neutral wind field, gravitational forces, and pressure-gradient forces. Each of these forces drives electric current flow, which has complex spatial and temporal structure. In the ionospheric E-region, which extends from about 90 to 120 km altitude, ion-neutral collisions are significant, while the electrons are mostly frozen to magnetic field lines. The frictional forces between the ions and the neutral wind field drive current which in general is not divergence free. Therefore, polarization electric fields build up globally to ensure divergence-free current flow. At low and mid-latitudes, the large-scale current system resulting from this ion-neutral coupling is called Solar-quiet (Sq). "Solar" refers to the current system's dependence on solar local time, since ionospheric conductivity at mid-latitudes peaks during daytime hours and diminishes during the night. The term "quiet" indicates geomagnetically quiet conditions, since during strong geomagnetic storms, the Sq system can experience large perturbations. At the magnetic equator, the horizontal geometry of the field lines leads to an enhanced zonal current called the equatorial electrojet (EEJ). An eastward component of the electric field at the equator, when coupled with the northward geomagnetic field, will drive vertical drift of electrons. Above about 120 km altitude, the Hall conductivity decreases substantially, since the reduced density of neutral particles result in less ion-neutral collisions, and the ions are essentially free to move with the electrons. This effect causes a nonconducting layer at the top of the E-region, and so the charge will accumulate at about 120 km altitude near the magnetic equator. This charge accumulation will cause a strong vertical polarization electric field which will drive zonal electron drift. Because the vertical polarization electric field is typically about 10 times stronger than the eastward component, the zonal **E** × **B** drift results in a strong current system, the EEJ. The EEJ has prominent magnetic signatures in both LEO satellite and ground observatory data, offering a convenient means of studying equatorial E-region dynamics. Since the EEJ is driven by the eastward component of the equatorial electric field (EEF), and is modulated by the neutral wind field, studying its magnetic signature reveals a lot of important information about the equatorial electrodynamics and also tidal features in the neutral winds.

This chapter will discuss recent methods of fitting an equivalent current model, based on line current geometry, to Swarm satellite magnetic measurements in order to recover the EEJ current strength flowing in the E-region. Line current methods have long been applied to studies of both the equatorial electrojet (Lühr et al. 2004; Alken et al. 2013a, 2014) and polar electrojets (Olsen 1996; Ritter et al. 2004; Vennerstrom and Moretto 2013; Aakjær et al. 2016). Since the electrojets typically follow lines of constant magnetic latitude, they have relatively simple spatial flow patterns, and line current models are particularly suitable for their study. For current systems which have more spatial complexity, such as Sq, pressure-gradient currents, interhemispheric field-aligned currents (IHFAC), and high-latitude field-aligned currents (FAC), it is necessary to utilize other methods to determine equivalent current flow, such as spectral methods (Fiori and Boteler 2018) or spherical elementary current systems (SECS) (Vanhamäki et al. 2018). This chapter will follow the work of Alken et al. (2013a, 2014) in order to define line currents which follow lines of constant quasi-dipole (QD) latitude in the equatorial region to estimate realistic EEJ equivalent current flow. These current estimates will then be inverted using a physics-based modeling procedure in order to recover the eastward component of the equatorial electric field.

#### **11.2 Satellite Data Preprocessing**