Linear Versus Nonlinear Methods for Detecting Magnetospheric and Ionospheric Current Systems Patterns

There is a growing interest in the development of models and methods of analysis aimed to recognize in the geomagnetic field signals the different contributions coming from the various sources both internal and external to the Earth. Many models describing the geomagnetic field of internal and external origin have been developed. Here, we investigate the possibility to recognize in the magnetic field of external origin the different contributions coming from external sources. We consider the measurements of the vertical component of the geomagnetic field recorded by the European Space Agency (ESA) Swarm A and B satellites at low and mid latitudes during a geomagnetically quiet period. We apply two different methods of analysis: a linear method, that is, the empirical orthogonal function (EOF), and a nonlinear one, that is, the multivariate empirical mode decomposition (MEMD). Due to the high nonlinear behavior of the different external contributions to the magnetic signal the MEMD seems to recognize better than EOF the main intrinsic modes capable of describing the different magnetic spatial structures embedded in the analyzed signal. By applying the MEMD only five modes and a residue are necessary to recognize the different contributions coming from the external sources in the magnetic signal against the 26 modes that are necessary in the case of the EOF. This study is an example of the potential of the MEMD to give new insights into the analysis of the geomagnetic field of external origin and to separate the ionospheric signal from the magnetospheric one in a simple and rapid way.


Background
The Earth's magnetic field results from different sources, both internal and external with respect to the solid Earth. The largest part of the magnetic field is of internal origin (the so-called main field), being mainly due to a self-sustaining hydrodynamic dynamo operating in the Earth's fluid outer core and only for a small part to the magnetized material in the crust. In addition to the internal field, there is the magnetic field generated by electric currents flowing in the ionosphere and the magnetosphere, called external field, whose strength ranges from less than 1 nT to some thousands of nT, according to different geomagnetic activity levels and latitudes. Lastly, in order to have an overall view of the different sources of the Earth's magnetic field we have to consider the magnetic fields generated by the electric currents in the crust and mantle, which are induced by the time-varying main and external fields. Similar induced currents can be also found within the salty waters of the oceans, which produce weak magnetic fields of the order of a few nanotesla at ground level (Baumjohann & Nakamura, 2009).
Of course when we make a measurement of the Earth's magnetic field on the ground or from a satellite in low Earth orbit it will collect the contributions from all the different examined sources, both internal and external to the solid Earth. For this reason, the recognition of individual contributions to the overall geomagnetic field is quite challenging. In recent years, there has been an increasing interest in the development of geomagnetic field models of increasing complexity and accuracy based on the combined analysis of both ground-based observatory magnetic measurements and data derived from several satellite missions. Among these models we mention GRIMM (it is an acronym for the GFZ Reference Internal Magnetic Model) (e.g., Lesur et al., 2010), POMME (POtsdam Magnetic Model of the Earth) (e.g., Maus et al., 2006), CHAOS (CHamp, Ørsted and Sac-C data) (e.g., Finlay et al., 2016;Olsen et al., 2014), and the well-known series of "Comprehensive Models" (CMs) (e.g., Sabaka et al., 2002Sabaka et al., , 2004Sabaka et al., , 2015. They are capable of adequately representing the different (internal and external) sources. In principle, these models were born with the goal of providing an accurate representation of the internal field, but very quickly it was clear that to push them to higher spatial and temporal resolution it was necessary to constrain at best also the magnetic field of external origin. Thus, the study of the external field is of cross-interest to the scientific community. For scientists working on the core and crustal fields the contribution of the external field is unwanted and represents essentially a source of noise which is useful to characterize (see, e.g., Finlay et al., 2016;Kunagu et al., 2013;Maus & Lürh, 2005). At the same time, for scientists working on ionosphere and magnetosphere, the external field is of central interest and permits the investigation of processes involving small magnetic strengths but fast timescales with respect to the dominant contribution represented by the internal field. Different methods have been developed and used to study the spatial and temporal structure of the ionospheric and magnetospheric current systems at various latitudes, which are the sources of external fields. Standard methods, such as spherical harmonic analysis (SHA) or spherical elementary current systems (SECs) (Amm, 1997;Amm & Viljanen, 1999), have been introduced to reconstruct the complex spatial and temporal features of these currents, but they have not often been capable of reproducing realistic current systems due to a priori constraints, the use of fixed basis functions, and intrinsic limitations caused by the unavailability of data.
In this paper we investigate the capabilities of two different methods of analysis to recognize and characterize the various sources responsible for the generation of the magnetic field of external origin recorded at low and mid magnetic latitudes. To this aim, we analyzed the magnetic data acquired by two of the satellites of the Swarm constellation (see, e.g., Friis-Christensen et al., 2006) spanning 2 years at 1 Hz cadence. We used the CHAOS-6 geomagnetic field model (Finlay et al., 2016;Olsen et al., 2014) to remove from the observed data the main field and its secular variation to obtain in the residual signal the geomagnetic field of external (magnetospheric and ionospheric) origin. We applied to the obtained external magnetic field both the empirical orthogonal function (EOF) analysis (Ghil et al., 2002) and the multivariate empirical mode decomposition (MEMD) method (Rehman & Mandic, 2010). The aim is to extract from the analyzed signal the main intrinsic modes describing the different magnetic spatial features inside it. We recognize in the various intrinsic modes the different ionospheric and magnetospheric contributions and compare the results from the two different methods in order to find the method that is capable of recognizing better the structures present in the analyzed signal.
The paper is organized as follows. Section 2 is dedicated to the description of the analyzed data set, while in section 3 we illustrate the two different chosen methods (EOF and MEMD) and their applications. Finally, in the last section we summarize the main findings and discuss the obtained results comparing the two different methods.

Data Description
We used Level-1b low-resolution (1 Hz) vector magnetic field data recorded on board two of the three satellites of the Swarm constellation (see, e.g., Friis-Christensen et al., 2006). In detail, we considered data recorded by Swarm A satellite during a period of 2 years from 1 April 2014 to 31 March 2016 and data recorded by Swarm B satellite for comparison. During this time interval the Swarm A (B) satellite flew around the Earth at an altitude of about 460 (510) km thus exploring the F region of the ionosphere. Data are freely available at ftp://swarm-diss.eo.esa.int upon registration.
We analyzed the vertical component of the geomagnetic field (B z , being measured inward to the Earth's surface) at low and mid latitudes (within ±65 • magnetic latitude) recorded during periods characterized by very low geomagnetic activity levels, which were selected using simultaneously two different geomagnetic indices: AE (Davis, 1966) and SYM-H (Iyemori, 1990). In particular, we considered the following simultaneous conditions: AE < 80 nT and −10nT < SYM-H < 5nT that permitted us to select periods where the magnetic disturbances due to storm and substorm events were excluded. AE and SYM-H data with 1 min time resolution were downloaded from the OMNI website (www.cdaweb.gsfc.nasa.gov/istp-public/).
As the main target of this work is to characterize the geomagnetic field of external origin and its spatial structure, we removed the internal geomagnetic field from the original data recorded by Swarm A (B) by using the CHAOS-6 model (Finlay et al., 2016). It is the latest generation of the CHAOS series of global geomagnetic field models introduced by Olsen et al. (2006Olsen et al. ( , 2010Olsen et al. ( , 2014. It is derived from Swarm, CHAMP, Østered, and SAC-C satellite magnetic data and ground observatory data, respectively. It is able to estimate the internal geomagnetic field with high resolution in time and space. It includes a parametrization of the quiet time, near-Earth magnetospheric field due to ring current, magnetotail, and magnetopause currents, but it does not take into account the contribution coming from the ionospheric currents. In other words, CHAOS-6 does not model all the sources of external origin in representing the geomagnetic field potential but only the magnetospheric ones. To remove from our data the internal field we have used the CHAOS-6 geomagnetic model up to the spherical harmonic degree N = 110. We binned data into 5 × 5 • -sized square bins across the Earth's surface after conversion to quasi-dipole (QD) latitude ( qd ) and local time (LT). We used the QD coordinates reference system (Richmond, 1995) mainly for two reasons: (i) with respect to orthogonal systems it captures the features (and the distortions) at all latitudes and is well defined everywhere (Emmert et al., 2010), and (ii) with respect to other nonorthogonal systems, due to its dependence on the geodetic altitude, it is very useful for magnetically localized phenomena with a specific height distribution, such as the current systems confined in the conducting layer of the ionosphere (Laundal & Richmond, 2016). Moreover, we considered the LT to better visualize the effects on the geomagnetic field due to the dynamical processes affecting the magnetosphere-ionosphere system. Figure 1 shows the qd versus LT global map of the geomagnetic field of external origin along theẑ (vertical) component computed from Swarm A observations. The mapped values are the average values falling within each bin (5 × 5 • -sized square bin). The minimum bin population is 3,009, the maximum is 10,025, and less than 5% of all the bins is populated with less than 4,000 data points. Thus, each bin of our 2-year-long observations is adequately populated, and the statistics is robust enough to make the average as representative of each data bin, allowing us to describe the mean geometry of the currents in the near-Earth space; that is, these patterns are clearly invariant with time, although seasonal variations are present, which will be reported in a forthcoming paper. As shown in Figure 1, the bin average external vertical field ranges between −20 and 20 nT, and a two-lobe structure is clearly visible. It is consistent with the solar quiet (S q ) daily variation of the geomagnetic field, a regular variation due to electric currents flowing in the ionosphere (e.g., Campbell, 2003). The basic pattern of the equivalent S q current system consists in a near-two-dimensional current circuit centered around noon at ∼110 km altitude fixed with respect to the Earth-Sun line and flowing in counterclockwise direction in the Northern Hemisphere and clockwise direction in the Southern Hemisphere. This current system generates an induced magnetic field alongẑ directed outward in the Northern Hemisphere and inward in the Southern Hemisphere, in both cases opposite to the main geomagnetic field vertical component, and thus it is revealed by Swarm observations as a decrease of the geomagnetic field in theẑ direction in the Northern Hemisphere and an increase in the Southern Hemisphere (e.g., Campbell, 2003). The regular magnetic variation associated with this ionospheric system is visible mainly when solar-wind-driven disturbances are absent. During geomagnetically disturbed periods, associated with the occurrence of storms and substorms, the S q signal tends to be easily masked. At low and mid latitudes other magnetic signatures can be detectable such as the magnetospheric ring current, magnetotail, and magnetopause currents. All these currents become stronger during times of enhanced geomagnetic activity, and for this reason their magnetic signatures become visible during geomagnetic disturbed periods. Nevertheless, a certain amount of ring current, which is the nearest magnetospheric current to the Earth, is always flowing even during quiet times. This current, centered in the magnetic equatorial plane, provides at Earth a uniform magnetic field which is aligned with the magnetic dipole axis and pointing southward. Thus, on our global map of the geomagnetic field of external origin along the vertical component, the field associated with the ring current appears as a positive contribution to B z in the Northern Hemisphere and a negative in the Southern Hemisphere.

Methods and Analysis
Usually, both univariate and multivariate analysis methods are based on a priori fixed decomposition basis, obtained by exploiting linearity and stationarity conditions (Chatfield, 2016). The above requirements, 10.1029/2019EA000559 strictly assumed to satisfy mathematical properties, are not generally verified when natural signals are analyzed, requiring adaptive analysis methods (Huang et al., 1998). In the following, we describe two different decomposition methods based on clearly different requirements: a linear method, that is, the EOF analysis (see, e.g., Chatfield, 2016;Ghil et al., 2002;Lorenz, 1956), and a nonlinear one, that is, the empirical mode decomposition (EMD) and its extensions (Huang & Wu, 2008;Rehman & Mandic, 2010).

EOF analysis
The EOF analysis, often called principal component analysis (PCA) in Earth sciences (see, e.g., Chatfield, 2016;Ghil et al., 2002), is a decomposition technique for both univariate and multivariate data. Generally, the univariate method is used for decomposing data into a sum of (orthogonal) components obtained by the diagonalization of the covariance matrix of the data based on embedding a given series of discrete data x(n) (of length N) in a matrix M of dimension m × N, m being the embedding dimension (see, e.g., Chatfield, 2016;Ghil et al., 2002;Takens, 1981). In the multivariate case, the data set is described by a data matrix {s(n)}| n∈N = {s 1 (n), s 2 (n), … , s k (n)}, assumed to be related to k observations for a given length N. Then, the set of observations is converted into a set of values of linearly uncorrelated variables, that is, the PCs Φ l (n), as l being the transpose of the lth eigenvector of the covariance matrix of the data obtained as C = {s} T {s}. Both EOFs and PCs can be also retrieved by applying the singular value decomposition (SVD) on the matrix of the data {s(n)}| n∈N = S under investigation as S = UΣV T , where U and V are orthonormal matrices, and Σ is a diagonal matrix. The columns of U are called left singular vectors, the rows of V T contain the elements of the right singular vectors, and the elements of Σ are called the singular values (Ghil et al., 2002). The right singular vectors are equivalent to the eigenvectors of the covariance matrix C, while the singular values l are equal to the square root of the eigenvalues l of C (Chatfield, 2016;Ghil et al., 2002). Thus, the decomposition is complete and orthogonal (by construction), the normalized eigenvalue l captures the partial variance (i.e., the energy content) of the lth principal component, and their sum exploits the total energy content (Ghil et al., 2002). In summary, the main steps of the EOF method are as follows: 1. to organize data as a matrix (by using the embedding theorem for univariate data (see, e.g., Takens, 1981)); 2. to evaluate the covariance matrix of data (embedded data for univariate data); 3. to diagonalize the covariance matrix to find eigenvectors and eigenvalues; 4. to project data on eigenvector directions to find the uncorrelated variables, that is, the principal components.
This method has been applied to different fields such as solar physics (see, e.g., Consolini et al., 2009;Vecchio et al., 2005), geomagnetic variations (see, e.g., Balasis & Egbert, 2006;De Michelis et al., 2010;Rotanova et al., 1982;Shore & Whaler, 2016;Xu & Kamide, 2004), and extensively in climate research (see, e.g., Ghil et al., 2002;Lorenz, 1956;Lovejoy & Schertzer, 2013). Here, we apply it to our data set. Having binned data into 5 × 5 • bins across the Earth's surface, the data matrix has a dimension (m × T) = (26 × 72), and consequently the method extracts a set of m = 26 components (L l ). However, to correctly deal with boundary effects, we show our results within ±60 • , without considering the boundary latitudinal bins. Since our data set consists of spatial measurements, we obtain eigenfunctions (i.e., EOFs and PCs) that depend on geomagnetic latitude and longitude, the latter expressed in terms of local time variations. Thus, we are investigating spatial variations at different scales by exploiting the local properties of the covariance matrix of the external geomagnetic field measurements. This means that we are able to detect the different spatial structures of the external components of the geomagnetic field. Figure 2 reports the results obtained by applying the EOF method to our data. The partial variance of each eigenvalue is shown in the upper panel, while some components (L l ) resulting from the analysis are reported in the other panels of the figure. From the values of the variance we notice that L 1 captures the most variance of the signal ( 1 ∼ 90%) and contributes with L 2 and L 3 to the reconstruction of the ∼98% of the total variance. L 4 − L 6 captures ∼1% of the variance, and the remaining components are below the noise level (Ghil et al., 2002).
The first three components (from L 1 to L 3 ), shown in the left column of Figure 2, are characterized by large-scale spatial patterns. Interestingly, the most energetic contribution given by L 1 does not reproduce the ALBERTI ET AL. main spatial pattern that is visible in the original data associated with the Sq daily variation. This structure is captured by L 2 . Indeed, L 1 is characterized by a symmetric spatial pattern both in latitude and in LT, which remembers the magnetic signature of the ring current. Conversely, L 2 is characterized by a two vortex-like structure centered around noon and symmetric with respect to the geomagnetic equator, in agreement with the Sq main pattern structure. On the contrary, L 3 seems to be characterized by a symmetric pattern in qd , with no evidence of LT symmetry. The right column panels of Figure 2 present some of the main characteristics of components L 4 -L 6 which show striped patterns, characterized by latitudinal ribbons of alternate positive and negative amplitudes. Finally, the remaining components (not shown) can be attributed to the noise, due to the low variance they account for (see, e.g., Ghil et al., 2002).

EMD
The EMD, differently from traditional data analysis techniques (like Fourier analysis or wavelets) (see, e.g., Chatfield, 2016), works directly in the data domain rather than in a conjugate one to extract the so-called intrinsic mode functions (IMFs) which satisfy two requirements: (i) the number of extrema and the number of zero crossings must be either equal or differ at most by 1, and (ii) at any data point, the mean value of the envelope defined using the local maxima and that obtained from the local minima is 0 (Huang et al., 1998). They are derived through a direct and adaptive process, called sifting process (Huang et al., 1998), which acts on a series x(t) as follows: 1. The local extrema are identified (i.e., local maxima and minima, corresponding to data points where abrupt changes are observed); 2. both local maxima and minima are separately interpolated by using a cubic spline, in order to have continuous (and smoothed) functions with smaller error than other polynomial interpolation, also avoiding the Runge's phenomenon (see, e.g., Prenter, 1975); 3. the spline interpolation produces the so-called upper u(t) and lower (t) envelopes; 4. the mean envelope m(t) is obtained as m(t) = u(t)+ (t) 2 ; 5. the so-called detail or candidate IMF is evaluated as The previous steps are iterated n times until the obtained detail h(t) can be identified as an IMF (often called empirical mode) (Huang et al., 1998), while the complete sifting process stops when no more empirical 10.1029/2019EA000559 modes, e.g., IMFs c i (t), can be extracted from data such that where r(t) is the residue of the decomposition, which can be a constant function, a monotonic function, or a function with only one extremum not containing an oscillatory component physically meaningful (Huang et al., 1998).
Analytically, the mathematical requirements for detecting an IMF are satisfied only when n → ∞; numerically, the sifting process is stopped after n * iterations according to a defined stopping criterion (Huang & Wu, 2008). The first criterion has been proposed by Huang et al. (1998) such that, being the sifting algorithm stops at the step n * when n * < 0 , 0 being between 0.2 and 0.3 (Huang et al., 1998). Another stopping criterion, e.g., the so-called threshold method proposed by Rilling et al. (2003), sets two thresholds, that is, 1 and 2 , to guarantee globally small fluctuations (as in Huang et al., 1998) and, in the meanwhile, to take into account locally large excursions (see, e.g., Flandrin et al., 2004;Rilling et al., 2003, for more details).
The decomposition procedure is completely adaptive, exclusively based on the local characteristic of the data, and highly efficient for processing nonlinear and/or nonstationary data (Huang & Wu, 2008). From a mathematical point of view, convergence is assured by construction, while orthogonality of the basis is satisfied in all practical senses, unless it is not theoretically guaranteed. However, by construction all empirical modes are locally orthogonal, since they are obtained by local maxima and minima properties (i.e., by the zeros of the first derivative), and also a posteriori globally orthogonal (e.g., Huang & Wu, 2008).
One of the novelties introduced by the EMD, beyond its adaptive character, is the concept of instantaneous amplitude and instantaneous phase (Huang et al., 1998). Indeed, once the decomposition is completed, by applying the Hilbert transform to each empirical mode it is possible to construct a complex analytical signal described by an amplitude wave modulation model. In this way, assuming to consider a time series, each empirical mode can be seen as an oscillating function with both time-dependent amplitude a i (t) and phase Both a i (t) and i (t) can be obtained by the Hilbert transform of the ith empirical mode, which is defined as P being the Cauchy principal value, such that from the complex analytical signal From the above concepts of instantaneous amplitude and phase, the mean energy content of each empirical mode can be simply derived as e i = 1 dt , the mean frequency as ⟨ i (t)⟩ t = 1 T ∑ T =1 i (t ), and the mean timescale as i = 2 ⟨ i (t)⟩ t (see, e.g., Alberti et al., 2017).
Based on numerical experiments on white noise data, Wu and Huang (2004) found that the EMD acts a dyadic filter, the empirical modes being all normally distributed and covering the same area on a semilogarithmic scale (see also Flandrin et al., 2004). This means that the product between the energy density of the ith empirical mode, defined as E i = 1 N ∑ N =1 [c i ( )] 2 with N being the length of the data, and its corresponding mean timescale i is constant such that the energy density is chi-square distributed (Wu & Huang, 2004).

10.1029/2019EA000559
This method can be used to assess the significance of empirical modes with respect to those derived from purely white noise processes, giving us theoretical spread function values on different confidence levels.
Being direct and intuitive, the EMD method is one of the most used adaptive methods, which is able to carefully analyze all those data resulting from nonlinear and/or nonstationary processes (see, for example, Consolini et al., 2017;Guhathakurta et al., 2008;Piersanti et al., 2017). It is capable of overcoming some limitations of different decomposition techniques (as for example a required fixed decomposition basis), also avoiding misleading results (as for fixed eigenfunction analysis) when complex and chaotic time series are analyzed (see, e.g., Alberti et al., 2019;Consolini et al., 2018). However, some outstanding problems, mostly dealing with end effects and/or stopping criteria, need to be outlined (see, e.g., Alberti et al., 2018;Huang & Wu, 2008;Wu & Huang, 2009), although various methods have proposed to avoid and/or mitigate these effects, as mirror and data extending methods (see, e.g., Huang & Wu, 2008;Yang et al., 2014).

MEMD
Although the EMD allows us to overcome some limitations when univariate signals are analyzed, it cannot be directly applied to multivariate data. The problem is that local extrema cannot be well defined on a n-dimensional space, and, consequently, the computation of the local mean is not possible and the concept of empirical mode is rather unknown (Rehman & Mandic, 2010). First attempts to approach multivariate signals by using EMD were based on channel-wise processing by applying univariate EMD to each channel (Huang & Wu, 2008). The algorithm idea was to generate a pseudomultivariate EMD by translating the univariate algorithm on n directions, grouping modes on similar scale by processing ensemble EMD over each direction (Huang & Wu, 2008).
To extend the concept of local extrema on k-dimensional space and to produce more suitable multivariate decompositions, Rehman and Mandic (2010) proposed to consider the k variate signal as formed by k-dimensional data sets, each of which was projected to appropriate directions over the k-dimensional space. In this way for each projected signal the envelops can be calculated for each direction, and, by averaging over the k-dimensional space, the local mean of the multivariate signal can be obtained using two different methods able to create a suitable set of direction vectors in the k-dimensional space. They are (i) the uniform angular sampling coordinates method and (ii) quasi-Monte Carlo-based low-discrepancy sequences. These methods provide an uniform distribution of direction vectors and more accurate local mean estimates in k-dimensional spaces (see, e.g., Rehman & Mandic, 2010, for more details).
Then, the usual steps (e.g., multivariate spline interpolation and IMF properties check) of the standard EMD are used to evaluate the multivariate IMFs such that a k variate signal {s(n)}| n∈N = {s 1 (n), s 2 (n), … , s k (n)} can be written as where the set of k-dimensional embedded patterns {c i (n)}| n∈N is affine to the univariate decomposition basis formed by the IMFs and {r(n)}| n∈N is affine to the univariate residue. This process decomposes a multivariate signal in several local monocomponent k-dimensional functions, each of which containing the same frequency distribution.
A characteristic scale for each MEMD mode can be obtained as  averaging over the k directions, we obtain the mean energy associated with each MEMD mode, through which the relative contribution can be derived as Finally, as for EMD modes (Huang et al., 1998), also MEMD modes empirically and locally satisfy orthogonal and completeness properties (Rehman & Mandic, 2010) in the k-dimensional space such that partial sums of equation (8) can be obtained.
When spatiotemporal signals are analyzed, MEMD is able to extract intrinsic spatiotemporal components with different characteristic spatial and temporal scales that can be used to investigate spatial patterns evolving in time without any a priori fixed assumption on linearity and stationarity of the signal. This means that MEMD is able to describe local (in terms of space) nonstationary (in terms of time) variations due to nonlinear components (in terms of amplitude variations in space and time). In our case, we applied the MEMD to spatial measurements such that the MEMD modes depend only on spatial coordinates (i.e., geomagnetic latitude and local time). In this way, we are able to detect variations of the external components of the geomagnetic field measurements at different spatial scales, which can be used to investigate the different spatial patterns of both ionospheric and magnetospheric current systems. We chose the threshold method proposed by Rilling et al. (2003) to stop the sifting process and we used the improved characteristic wave algorithm to prolong the data series at the boundaries to deal with the edge effect (see, e.g., Huang et al., 1998;Huang & Wu, 2008). However, the results are not significantly sensitive to the chosen threshold parameters and/or boundary algorithms. Figure 3 reports the results of the MEMD decomposition of B z for the Swarm A satellite observations. In the top panel of the same figure we report the percentage energy, calculated from equation (10), associated with each IMF as a function of the corresponding number. The first three modes contain less than 3% of the total energy of the signal (brown dots); conversely, the modes with i = 4, 5 contain ∼97% of total energy and consequently the signal obtained from the superposition of these modes represents the main part of the original one. In the other panels the IMFs (c i ), obtained applying the MEMD technique, are shown and sorted in an increasing-scale order from 1 (the smallest spatial scale) to 5 (the largest spatial scale). At last, the residue is shown in the bottom right panel of Figure 3. The number of detected IMFs and their 10.1029/2019EA000559 characteristic spatial scales are automatically found by the algorithm according to the criteria described above, the procedure being completely adaptive (in contrast with EOF analysis, where the number of components depends on data matrix dimension). Moreover, due to the nonlinear behavior and spatial dependence of the different components, the adaptive nature of the MEMD method can be really helpful in detecting the different spatial features and variations of both magnetospheric and ionospheric source processes and currents. The left column panels of Figure 3 illustrate some of the main characteristics of the first three IMFs. These IMFs (c 1 , c 2 and c 3 ) are characterized by an amplitude in the range ±5 nT and their spatial structures are similar to latitudinal ribbons alternating positive and negative amplitudes. In the right column of Figure 3 the IMFs 4-5 are shown (c 4 and c 5 ). The large scale patterns in the maps have strengths spanning the range from ∼ ±5 to ∼ ±10 nT and represent the main structure originated by the S q current in quietness, c 5 being the main component and c 4 its spatial harmonics. In fact, the component with the largest spatial scale (c 5 ) contains patterns which have the right characteristics in order to represent the main contribution to the S q : they are centered at noon, have a negative (positive) field variation in the Northern (Southern) Hemisphere in a background of opposite sign, and extend for about 12 hr, which is the time period marking the transition from the dayside to the nightside and vice versa. On the other hand, the features appearing in c 4 may be considered as harmonics embedded in the main variation.
The MEMD technique provides also the residual of the original map (referred to in equation (8) as r(t)), that is, the part of the original signal that cannot be decomposed into IMFs, as shown in the bottom right panel of Figure 3. It ranges within ∼ ±20 nT, is positive in the Northern Hemisphere, and is negative in the Southern Hemisphere. This implies that the MEMD residual of B z is inward in the Northern Hemisphere and outward in the Southern Hemisphere. At qd between ∼ ±20 • the residual assumes very small values, which increases at increasing qd . We also note that the increase of the external field at higher latitudes is a feature common to all longitudes, and no localized patterns appear, unlike it happens in all the detected IMFs (and also at high latitudes for the most energetic component L 1 detected by EOF analysis, see Figure 2). Similar results have been found for the B x and B components (see, e.g., Alberti, 2018).

Results and Conclusions
We applied two different methods of analysis to our data set consisting of the spatial measurements of the geomagnetic field vertical component at low and mid latitudes during a geomagnetically quiet period. The aim is to compare the results coming from the two methods, one linear (EOF) and one nonlinear (MEMD), in order to understand which one is the best to recognize the magnetic spatial structures of external origin embedded in the data. Figure 4 reports a comparison between the results obtained from the two different methods of analysis. In detail, we report the results obtained from EOF decomposition method in the panels on the left of Figure 4, while the results obtained from MEMD are shown in the panels on the right of the same figure. From Figure 4 we notice that by applying the MEMD method we are capable of separating the different modes that contribute to the magnetic field of external origin during quiet periods. We find that our patterns can be represented as a linear combination of five empirical modes and a residue. The first three modes, that is, those characterized by the smallest spatial scales in LT, appear in form of spurious north-south patterns. The other two modes, that is, those with the largest spatial scales, seem to describe the effects on the geomagnetic field of the electric currents flowing in the ionosphere, that is, mainly the S q ionospheric current pattern. Lastly the residual, which represents the long-term trend of the analyzed data, seems to be due to the electric currents flowing in the magnetosphere and describes the effect on the geomagnetic field of the magnetospheric ring current. In fact, when considering only theẑ component of the magnetic field, the presence of the magnetospheric ring current should add a contribution to the magnetic field which is basically null at and nearby the magnetic equator, and should increase with the latitude, like what can be observed looking at the residual map. It is important to notice indeed that the ring current, which is known to lead to a global-scale reduction in the horizontal component of the geomagnetic field during the geomagnetic storm, is a magnetospheric current which always exists, also during quiet periods (Shore & Whaler, 2016). Only its intensity and distance from the Earth change during the disturbed periods (De Michelis et al., 1997), together with the partial ring current (Milan et al., 2017). By applying the EOF analysis we are able, also in this case, to separate the different modes, which contribute to the magnetic field of external origin. However, in this case, the different magnetic spatial structures embedded in the analyzed signal are more difficult to recognize. We can recognize the magnetic field due to the ring current in the first EOF (L 1 ) and the magnetic (From top to bottom) L 1 and the residue of the MEMD method can be attributed to the ring current contribution, L 2 and c 5 patterns can be related to the main S q pattern, L 3 and c 4 can be attributed to a subharmonic structure of the S q current, while short-scale reconstructions L 4−26 and C 1−3 could be related to different source mechanisms (external driver, magnetopause current).
field due to the S q ionospheric current pattern in the second EOF (L 2 ). Conversely, the third EOF (L 3 ) does not seem to describe the effect on the magnetic field produced by a particular current system but it could be a subharmonic of the EOF L 2 and consequently to partially describe the effects on the magnetic field of the S q ionospheric current pattern. However, all these three modes are contaminated by the solar quiet daily variation. Thus, the method does not seem to be capable of completely separating the different spatial structures probably due to the nonlinear nature of the analyzed signal. Moreover, other 23 EOFs are necessary to completely reproduce the original data. To confirm our interpretation about the origin of the different contributions (ionospheric Sq or magnetospheric ring current) we have repeated our analysis on magnetic data recorded by Swarm B satellite, which flows at an higher altitude than Swarm A (about 50 km). By analyzing the difference between the results obtained by the two satellites (data not shown here) we found that the residual magnetic field increases with the altitude, as it is expected in the case of a contribution due to the magnetospheric current systems, while the contribution due to the Sq current system decreases with the altitude. Furthermore, by analyzing the ionospheric field, obtained by removing from the original data the internal magnetic field and the magnetospheric one modeled by CHAOS-6, the contribution due to the ring current cannot be revealed (data not shown).
In order to show more clearly the differences between the two methods, we compare the longitudinal (i.e., local time) behavior of the ionospheric contribution obtained from the original data by using CHAOS-6 model at fixed latitudes with the signals describing the magnetic field due to sources localized in the ionosphere obtained from the two methods. The results are reported in Figure 5. First, we notice that the behavior of B iono z (red asterisks) is that expected in quiet conditions, being a few nT from dusk to dawn, with a negative bump up to ≃10-15 nT in the Northern Hemisphere and a positive bump in the Southern Hemisphere around noon. The comparison among the three signals shows that the MEMD analysis is able to reconstruct the magnetic signal of ionospheric origin better than the EOF analysis. This is clearly visible at midlatitude where the trend reproduced by the combination of the IMFs c 4 and c 5 (green line) well describes the effect of S q ionospheric pattern on the magnetic field. Conversely, the EOF reconstruction of the magnetic field of ionospheric origin (blue line, L 2 + L 3 ) is not very good as can be realized comparing it with the original data at midlatitude, due to an incorrect estimation of the nonlinear residue (note that nor L 1 neither the residue of the MEMD have been included in reconstructions of EOFs and IMFs). To quantify the different fits to the B iono z data we have estimated the correlation coefficients between B iono z and both MEMD and EOF reconstructions of the Sq variability in the local time interval between 06:00 LT and 18:00 LT, where the Sq current . Longitutinal (LT) behavior of B z from Swarm A observations at different latitudes, from 37.5 ± 2.5 (top panel) down to −37.5 ± 2.5 (bottom panel), respectively. Red asterisks mark the ionospheric contribution derived by CHAOS-6 (B iono z ); the blue solid line represents the summed EOFs L 2 + L 3 ; the solid green line represents the summed IMFs c 4 + c 5 . The error on QD latitudes is computed as half of the bin size. r EOF and r MEMD refer to the values of correlation coefficient between B iono z and Sq reconstructions by using EOF (blue text) and MEMD (green text), respectively. systems are localized. The results, reported in Figure 5, confirm that a higher correlation is found between B iono z and MEMD reconstructions. Moreover, it is interesting to note that similar large-scale structures have been found by using both EOF and MEMD which is an indication of the robustness and significance of the detected spatial variability on these scales.
In general, therefore, it seems that MEMD method can help in the interpretation of the external magnetic field signals better than EOF method. Using MEMD analysis a few modes are necessary to recognize in the magnetic signal the different contributions coming from external sources. They are not the result of a model but can be directly extracted from the original signals with no a priori assumption on the nature of data. These modes, each associated with a characteristic spatial scale, describe the basis representing the data and 10.1029/2019EA000559 are able to identify various dynamical components of the analyzed signals that can be related to different physical scales and sources. This study is an example of the potential of the MEMD method to give new insights into the analysis of the different sources responsible for the geomagnetic field of external origin; and at the same time, it can be used as a good filter in the analysis of the geomagnetic field of external origin, permitting to separate the ionospheric signal from the magnetospheric one. The results presented in this paper rely on data collected by one of the three satellites of the Swarm constellation which are freely available at ftp:// swarm-diss.eo.esa.int upon registration. We thank the European Space Agency that supports the Swarm mission. Geomagnetic indices data were downloaded from the OMNI website (www.cdaweb.gsfc.nasa.gov/ istp-public/). The authors kindly acknowledge N. Papitashvili and J. King at the National Space Science Data Center of the Goddard Space Flight Center for the permission to use 1 min OMNI data and the NASA CDAWeb team for making these data available. F. G., P. D. M., and G. C. acknowledge the support by ESA under Contract 4000125663/18/I-NB (INTENS).