Decadal period external magnetic field variations determined via eigenanalysis

We perform a reanalysis of hourly mean magnetic data from ground‐based observatories spanning 1997–2009 inclusive, in order to isolate (after removal of core and crustal field estimates) the spatiotemporal morphology of the external fields important to mantle induction, on (long) periods of months to a full solar cycle. Our analysis focuses on geomagnetically quiet days and middle to low latitudes. We use the climatological eigenanalysis technique called empirical orthogonal functions (EOFs), which allows us to identify discrete spatiotemporal patterns with no a priori specification of their geometry—the form of the decomposition is controlled by the data. We apply a spherical harmonic analysis to the EOF outputs in a joint inversion for internal and external coefficients. The results justify our assumption that the EOF procedure responds primarily to the long‐period external inducing field contributions. Though we cannot determine uniquely the contributory source regions of these inducing fields, we find that they have distinct temporal characteristics which enable some inference of sources. An identified annual‐period pattern appears to stem from a north‐south seasonal motion of the background mean external field distribution. Separate patterns of semiannual and solar‐cycle‐length periods appear to stem from the amplitude modulations of spatially fixed background fields.


Introduction
Studies of the electrical conductivity of the Earth's mantle are crucial to the accurate separation of the magnetic fields produced internal to the Earth's surface, in the Earth's liquid iron core, in the lithosphere, and induced in the mantle. In turn, accurate estimates of mantle conductivities depend upon a good knowledge of the inducing fields which originate in the Earth's magnetosphere and ionosphere.
As discussed by Olsen [1999] and more recently by Kuvshinov [2012] there are multiple different mathematical approaches to probing the mantle at depths greater than approximately 1000 km, but each requires a description of the external magnetic field variations of long periods, i.e., months to years. This overlaps with the periods on which the core field varies significantly; hence, the estimates of long-period external fields are subject to substantial uncertainty. Typically, [e.g., Banks, 1969;Constable and Constable, 2004], the assumed geometry for the inducing source on periods greater than a few months is the spherical harmonic P 0 1 representing the symmetric magnetospheric ring current field. It has been known for some time that the ring current exhibits asymmetry [e.g., Daglis et al., 1999;Balasis et al., 2004;Balasis and Egbert, 2006]. Recent studies [e.g., Maus and Lühr, 2005;Lühr and Maus, 2010] have shown that the nonaxisymmetric magnetospheric source dominates the ring current at quiet times. More recently, Kuvshinov [2013, 2014] used a fuller set of harmonics for the magnetospheric terms based on the study of Sabaka et al. [2013] (though as discussed later, this is still possibly a simplification of the true case since it does not include the ionospheric contribution). Püthe et al. [2015] have introduced a method allowing the resolution of the induced response from an inducing field of effectively unlimited complexity.
In this study, we describe the spatial and temporal morphology of the long-period external magnetic fields for quiet days over a full solar cycle. We use the method of empirical orthogonal functions (EOFs) in a reanalysis of 13 years (over a solar cycle) of ground-based observatory magnetic vector hourly means. EOFs are used to decompose a large, highly variable data set into a small number of independent spatiotemporal subgroups, called "modes," which describe the patterns of maximum variance in the data set. The modes arise

10.1002/2015JA022066
Key Points: • Spatiotemporal patterns of external inducing fields resolved over a full solar cycle • Relative importance of the ring current to the total inducing signal quantified over solar cycle • Resolved spatial patterns of current systems with annual and semiannual periods of variability from the intrinsic dynamical properties of the data; hence, we need apply no simplifying assumptions of the long-period field morphology prior to analysis.
In section 2 we describe the data and the processing applied to them. In section 3 we outline the theory governing the EOF decomposition and describe our numerical implementation of the EOF method. In section 4 we present the results obtained from combining the outputs of several EOF analyses via spherical harmonic analysis and discuss these further in section 5. In section 6 we summarize our findings.

Data Processing
The magnetic field data are hourly mean vector values from a ground-based network of 149 observatories spanning the period 1997-2010 (inclusive), described in full in the ESA Swarm Level 2 processing report of Olsen et al. [2011] and by Macmillan and Olsen [2013]. Due to higher proportions of missing data in 2010, our analysis will span 1997.0-2010.0. Seven of the observatories had unexplained baseline drifts, which are unlikely to have a physical explanation, and are therefore ignored in the subsequent analysis. To resolve local time (LT) structures in the external fields for a full solar cycle, we must perform various data selection and processing procedures prior to the EOF analysis of the data.
To isolate the external field component, core field and secular variation estimates are removed using version 4 of the CHAOS-4 model of Olsen et al. [2014]. This model does not describe the local lithospheric field, so this contribution is removed by subtracting the mean value for each component (at each observatory) at geomagnetic quiet times-defined as when Kp < 2 + and |dDst∕dt| < 2 nT/h-over the 13 year time series [Olsen et al., 2011].

Data Rotation
The EOF decomposition does not take the positions of the input data into account, so we aim to make the signals of interest as spatially static as possible in the frame of reference of our analysis. Given the source fields under consideration, we achieve this aim by choosing a reference frame based on the Sun. We rotate all data (coordinates and components) to the inertial frame of the SM (solar magnetic) system [e.g., Hapgood, 1992], in which the (r, , ) components are, respectively, radius (positive up), dipole colatitude, and dipole local time (DLT) in a geocentric spherical polar sense. The horizontal plane of SM is the geomagnetic dipole equatorial plane, and the SM axis is perpendicular to a plane containing the Earth-Sun line and the geomagnetic dipole axis. After this rotation we force the observatory distribution to be (relatively) static in LT throughout the timespan of the analysis, by selecting data at only one UT (universal time) hour per day. The distribution of the observatories in LT is then dependent on which UT is chosen to form the data subsets-this is termed the "start-UT" of the data arrangement. We perform 24 EOF analyses, each with a data subset from a different start-UT (and later combine the results of these analyses).

Data Selection
We need each of the 24 EOF decompositions of the data to be similar. We have found that processing data from periods of very high activity makes it less likely that the decomposition for each start-UT will lead to a consistent ranking (in terms of data variance explained) of the characteristic temporal oscillations. Thus, we select data from only the internationally defined five magnetically quietest days per month.

Data Weighting
Since the EOF analysis responds only to correlated data variance, any spatial clustering of observatories (e.g., Europe) will have a disproportionate contribution to the EOF solution, an effect we must mitigate. We apply weighting to counteract the tendency for high concentrations of observatories to have highly correlated signal-weighting procedures are discussed in detail by Baldwin et al. [2009] andJolliffe [2002, p. 382]. The amount of independent information contained in each observatory's data time series is approximated as follows. First, the observatories' spatial distribution is triangulated (via calculating the convex hull) [Weisstein, 2003], as shown in Figure 1. We then sum the surface areas of each spherical triangle which has a given observatory as a vertex, which provides the weighting value for that observatory. The weights have the effect of decreasing the contribution of dense clusters of observatories to the decomposition. The vector of weights has one value for each observatory. To avoid the auroral regions dominating the decomposition, we multiply this vector by a polar mask (not shown) to enforce zero contribution to the EOF solution from observatories within 30 ∘ of the geomagnetic poles (this reduces the effective observatory count from 142 to 112). The resulting weights w are applied to the data prior to EOF analysis. The magnetic field signal which remains after these SHORE ET AL. DECADAL PERIOD EXTERNAL FIELDS 5173 corrections has contributions from multiple current systems which affect the middle-to-low latitude region. Since we aim to discover the relative contribution of each of these to the long-period external inducing field, we do not correct the data for any further estimated contributions.

Data Infill
The EOF method requires that missing data in each observatory's time series are infilled prior to analysis. This entails making assumptions about the nature of the unmeasured signal, so we choose to replicate only the basic spatial signal of the daily variation in our interpolations. To perform the infill, we apply an hour-by-hour spherical harmonic analysis (SHA) to the observatory data. The proportion of missing data in each year ranges from 16 to 33%. The SHA has harmonics up to degree 2 and order 1 assuming an external field only, which we consider a suitable representation of the basic signal. The SHA model is used to predict missing data in each hourly epoch.
For each magnetic component we format the processed observatory data into matrices whose columns are time series of 780 hourly means (comprising one UT from each of the five quietest days per month for 13 years) recorded at each of the 142 observatories, for a given single start-UT. We then concatenate these matrices from each magnetic component to form the full "data matrix" X, which has dimensions of n × p, where n is the number of hourly means (780) and p is the number of spatial parameters (three times the number of observatories).
r 1,1 r 2,1 ⋮ r n,1 r 1,2 r 2,2 ⋮ r n,2 · · · · · · r 1,p∕3 r 2,p∕3 ⋮ r n,p∕3 where (r, , ) indicate the magnetic components in the respective directions, not the coordinates themselves. We will decompose the variance of X, so prior to the EOF analysis of these data we remove the time mean of each column known as "centering" the data matrix where X are the columnar means of X, averaging over time. The removed means are stored and will be used in later processing. The variance fieldX is now ready for EOF analysis.

Method
The following description of the EOF process is based on Bjornsson andVenegas [1997, p. 12], von Storch andZwiers [2002, pp. 294-295], and Jolliffe [2002, p. 5]. Given the centered variance fieldX, the EOF analysis defines p independent spatial patterns v (represented by , where the v j are column vectors of length p) and p independent temporal patterns t (represented by T = ( t 1 , t 2 , ..., t p ) , where the t j are column vectors of length n) which, when combined, reconstructX The V are the eigenvectors of R, the covariance matrix ofX, defined as R =X TX . The T are given by projecting the spatial eigenvectors onto the original data their projections onto the original data as the Principal Components (PCs), or just "temporal oscillations." Each pair of a spatial pattern and its temporal oscillation are termed an EOF "mode." The modes (which are ranked in order of decreasing eigenvalue of R) are the orthogonal patterns which describe the maximum possible variance of the data. The proof that the basis vectors which diagonalize R also maximize the variance of their projection ontoX is given by Hannachi et al. [2007] and von Storch and Navarra [1999]. The spatial and temporal nature of V and T arises from how we have defined the data matrix- Richman [1986] discussed several methods of forming X from the data, but we only use the layout given in equation (1) here.
The application of the (essentially scalar) EOF method to the vector data matrixX follows from the approach termed Common (or Combined) Principal Component Analysis by Jolliffe [2002]. We have treated the three magnetic components as additional parameters in the spatial dimension. Hence, a given spatial pattern v i describes the relative variation in amplitude between the components, and the three components are each subject to the same (associated) temporal oscillation t i . This approach simplifies the physical interpretation of the modes and avoids differences in mode ranking between magnetic components of characteristically different variability.
The value of the EOF spatiotemporal decomposition is that the vast majority of the variance of the data can be described with some small number k ≪ p of modes. The spatial patterns of the external fields from different source regions are not orthogonal (i.e., mutually physically exclusive). However, certain amounts of the fields from each source region will be increasingly temporally correlated on seasonal (and longer) timescales. Hence, it is reasonable to expect that the modes will identify distinct spatiotemporal patterns corresponding to (for instance) the semiannual, annual, and solar-cycle-length trends in the external magnetic field. It is for this reason that we find the decomposition based on successively and orthogonally maximizing variance to be useful. Note that we should not expect the patterns to stem individually from a single source region.
Recall that we must weight the EOF decomposition to account for the observatories' spatial distribution. The weighting matrix W has dimension p × p, composed of the vector of weights w repeated three times on the diagonal and zeros elsewhere. The weighting matrix is applied to the data matrix prior to the EOF analysis viã where ′ indicates that weighting has been applied toF and the exponent of 1/2 applied to W accounts for the squaring implicit when the covariance matrix is calculated. The change in terminology fromX toF reflects the fact that some information is lost in the application of the polar mask applied to the weighting matrix; hence, the deweighted data matrix (described below) is not equal to the unweighted data matrix (i.e., the original data). Now we form the weighted covariance matrix R ′ =F ′TF′ and define B ′ to be the weighted EOFs, which diagonalize R ′ . The temporal coefficients (Y ′ ) of the weighted EOFs are given by The deweighted EOFs are computed via B = W −1∕2 B ′ [Jolliffe, 2002]; they are used to calculate a deweighted version of Y ′ for display purposes, denoted Y, via where y j is the jth columnar vector of Y, j runs from 1 to p (there are p EOFs and p PCs), and the scalar B j is the maximum-amplitude element of b j , itself the jth column vector of B.
For display purposes and to interpret the modes in the units of the original data, we extract a "deweighted" partial version of the data matrix for a single mode (denoted here as j) as follows [Jolliffe, 2002] All temporal oscillations shown are deweighted versions, whereas the spatial reconstructions may be of either the original or deweighted data. Unless otherwise stated, we depict the component in all spatial maps. The weighted EOF analysis method described in this section is applied to each of the 24 start-UTs. We describe the results of these analyses in the next section.
SHORE ET AL. DECADAL PERIOD EXTERNAL FIELDS 5175 Figure 2. Proportion of variance accounted for by each of the first 10 modes: a normalized eigenvalue spectrum. This mean eigenspectrum is computed from the eigenspectra of each start-UT-the envelope of these contributory spectra is shown in the form of error bars on the mean spectrum.

Eigenspectrum and Mode Order
In Figure 2 we show the eigenspectrum, taken as a mean across each of the eigenspectra from the 24 EOF analyses for each start-UT. Modes 1-4 explain the majority (67.4%) of the variance of the data. The error bars in Figure 2 denote the differences between the proportions of variance explained (for a given mode) between the different start-UTs. We will later combine the different start-UTs to assess the nature of the modes in more detail-here we assess the extent to which the temporal oscillations of modes 1-4 are the same between the different start-UTs.
In Figure 3 we show the (absolute values of the) correlations between the temporal oscillations of a given mode from a reference start-UT of 18:00 and the same mode from all other start-UTs. This method of autocorrelating the modes reveals the consistency of their temporal patterns between different start-UTs. We expect that a mode will be most different when its start-UT is about 12 h away from 18:00, though there will be some complexity associated with the structure of the mode in dipole local time (DLT) and the irregular observatory distribution. For mode 1 (the dark blue line in Figure 3), we see that the correlation values start at 1 at 18:00 and dip to slightly below 0.9 at start-UTs near 06:00, indicating that this mode is consistent across all start-UTs. The same is true of mode 2 (red line), though the correlations drop off faster away from 18:00 and exhibit more variability. The correlations for mode 3 (green line) exhibit a similar pattern, except for start-UTs between 20:00 and 23:00, for which the autocorrelations are minimal. The autocorrelation values for mode 4 are not shown since they exhibit high variability and drop to low values (below 0.4). We perform a correlation (Auto)correlations between the temporal oscillations of a given mode from a reference start-UT of 18:00 and the same mode from all other start-UTs. Mode 1 is shown in dark blue, 2 in red, and 3 in green. A consistently high correlation indicates that a consistent pattern is being represented for a given mode across the EOF analyses of each start-UT. This is important since we later combine the signals from each start-UT. The light blue line shows the correlation between mode 3 at start-UT 18:00 and mode 4 for each other start-UT-see text for further details.
between mode 3 at a start-UT of 18:00 and mode 4 at all start-UTs (light blue line) and note that these correlations improve markedly between start-UTs between 20:00 and 23:00. This indicates that the spatiotemporal pattern typically described by mode 3 is described instead by mode 4 between these start-UTs. Henceforth, the term "mode 3" refers to the collection of mode 3 from start-UTs of 01:00-19:00 and 24:00 and mode 4 from start-UTs between 20:00 and 23:00. In the next section we show that with this substitution, we can associate a physical meaning to these first three modes, and since they describe the majority (62%) of the variance, we focus on modes 1-3 for the remainder of our analysis. For our purposes, all other modes are considered to be noise.

Identification of Temporal Oscillations
In Figure 4a we show the mean time series (based on all temporal oscillations from each start-UT) for modes 1-3 (respectively, SHORE ET AL. DECADAL PERIOD EXTERNAL FIELDS 5176 blue, red, and green). These quiet time daily averages are irregularly distributed and are shown with a piecewise linear interpolation-this is for visual purposes only and does not affect the following analysis. We show the power spectra of these mean series in Figure 5. From Figures 4 and 5 we can begin to identify the physical meaning behind each mode.
For comparison with mode 1, we show (in Figure 4a) the daily means of the summed external and induced components of the (negative) P 0 1 term of the RC index (representing the symmetric ring current described by Olsen et al. [2014]) for the appropriate quiet days. RC is calculated in the same coordinate system of magnetic latitudes that we use in our analysis. The time series of RC and mode 1 have a Pearson correlation value of 0.98 and are very similar in amplitude-from this we infer that mode 1 describes the ring current. Note that correlations of the RC index with modes 2 and 3 do not exceed an (absolute) value of 0.06. The residual of mode 1 and the RC index is shown in Figure 4c-while the residual amplitude is small, it exhibits an annual oscillation which modulates with the solar cycle. From the power spectra shown in Figure 5, we see that RC has a smaller annual-period spectral peak than mode 1, indicating that mode 1 is the origin of the annual signal in Figure 4c. Our data do not enable us to determine whether this annual signal is truly representative of the ring current or if it is due to an imperfect separation of the modes. The power spectra of mode 1 and RC show that both modulate most strongly with the 11 year solar cycle.

SHORE ET AL.
DECADAL PERIOD EXTERNAL FIELDS 5177 The time series of mode 2 shown in Figure 4a exhibits a clear annual oscillation which modulates (weakly) with the 11 year solar cycle-both these aspects are reflected in its power spectrum shown in Figure 5. Mode 2 cannot be meaningfully compared with indices of solar-terrestrial coupling since it is an apparent oscillation, resulting from the motion of the Sun (and thus, the latitude of the peak current flow) in the SM system. The physical nature of this mode will be discussed in section 4.4.
The time series of mode 3 in Figure 4a has a clear modulation with the solar cycle and a correlation of 0.73 with the F 10.7 index shown in Figure 4b. It is therefore likely that this mode represents an ionospheric signal. From the power spectrum of mode 3 (shown in Figure 5) we can confirm its solar cycle modulation and note that the next most powerful spectral peak is at a period of 6 months. This semiannual oscillation present in mode 3 is clearest between the years 2006 and 2007 (Figure 4a), when the mean external field activity is lower. As for mode 2, it is likely that the semiannual component of the mode 3 pattern is due to apparent motions within the SM frame, perhaps according to the Russell-McPherron effect which arises due to the tilt of the Earth with respect to the Sun at the equinoxes [Russell and McPherron, 1973].
We have identified three characteristic oscillations: mode 1 represents the ring current, mode 2 is an annual oscillation (as we show in section 4.4, it is more correctly the annual modulation of the background mean field), and mode 3 is a semiannual oscillation (likely of ionospheric current fields), each of which modulates with the solar cycle. When averaged over each start-UT, the proportion of the variance accounted for by these patterns is, respectively, 40.3%, 14.3%, and 7.4%. In the following two sections, we assess the spatial patterns of these modes.

Combining EOF Analyses
We combine the results from separate start-UTs so that our spatiotemporal description of the inducing fields uses all the information we have available. We compute the data matrix reconstructions (given by equation (7)) for the three modes identified in the previous section. For each different start-UT, the associated reconstruction will have a specific distribution of the observatories in LT. We are able to recombine the data predictions from all start-UTs in the SM frame to obtain a more complete (but still irregular) global coverage of discrete data predictions from the EOF method. We take advantage of this improvement in coverage to perform a SHA of the isolated modes to form continuous models in space with a daily temporal resolution.
Since ∇ ⋅ B = 0, we can relate the magnetic field measurements B to the scalar potential V by B = −∇V. Since V satisfies Laplace's equation it can be expressed in spherical polar coordinates as [e.g., Langel, 1987] V(r, , )=a where a = 6371.2 km is a reference radius, r is radius, and are SM colatitude and longitude as defined in section 2.2, P m l (cos ) are the associated Schmidt seminormalized Legendre functions [Langel, 1987], l and m are the degree and order of the expansion up to a maximum degree of N, g m l and h m l are the internal (induced) field Gauss coefficients, and q m l and s m l are the external field coefficients. The EOF analysis will respond to both inducing (external) and induced (internal) fields-the use of a joint internal-external expansion allows us to isolate just the (external) inducing fields of interest.
We will form the symmetric ring current SH model (which we term "SR") from mode 1, the annual oscillation model (which we term "AO") from mode 2, and the semiannual model ("SA") from mode 3. These three models do not describe patterns which are stationary in the SM frame over the period of analysis. Hence, it is also desirable to produce a model of the "ground state" that these three models are describing the variance of. For this reason, we will also produce a model of the removed means ("RM") described by equation (2).
As an example, to construct the data vector for the AO model SHA, we first extract the EOFs and PCs for mode 2 and a given start-UT. We use a notation such that y ′ 2(01) is the second column vector (i.e., eigenvector of the second mode) of Y ′ and for a start-UT of 01:00. Likewise, we also extract b ′ 2(01) from B ′ . We calculate the data matrix reconstruction for all analysis times and spatial positions, for mode 2, and for a start-UT of 01:00 as Note thatF 2(01) contains contributions from all three components and has the same structure as X in equation (1). We perform SHA on each day of data-due to the combination of 24 start-UTs within each of these periods, the epoch of the SHA is at a time of 12:00, for an analysis covering all of that day. Hence, we must extract a row vector (corresponding to all spatial locations and one time) from each reconstructed data matrix, for each component. For example,f i2r (01) is the ith row (epoch)-each of its j elements (locations) is from a different column ofF 2(01) which corresponds to the r component. We repeat this process for all start-UTs, appending 24 transposed row-vector reconstructions for a given component (here, r) into a single column vectorf The result of this concatenation is the improved coverage of observatories in the SM frame. The improved distribution is shown in section 4.4.
We concatenate the start-UT series in a similar manner as shown above for the and components, following which these single-component series are appended into a data vector for the SHA of the ith epoch Using these data we solve for an estimate of the vectorm of spherical harmonic coefficients of equation (8) via iteratively reweighted least squares [Olsen, 2002;Olsen et al., 2009], utilizing Huber weights with a tuning constant of 1.5 (but without spatial regularization). We performed the joint internal-external SHA with an initially low value of N and investigated stability as it was increased. Our final data predictions are from harmonics up to degree 5 of an analysis with N = 16, since these coefficients were the most stable. The analysis involves forming the equations of condition matrix G [Menke, 1989] which relates the vector of spherical harmonic expansion coefficients from equation (8) to the magnetic components in equation (11) and calculating the data predictiond i (e.g., of the AO model for the ith epoch) bŷ Unless otherwise stated, all data predictions presented here are calculated from just the external field coefficients. The residual of the SHA data prediction to the input data (the "SHA residual") is given by d i −d i .

Spatial Patterns of the SHA Data Predictions
The full description of the solar cycle long-period oscillation of the external inducing fields provided by our analysis is given by the sum of the four SH models (SR, AO, SA, and RM). However, identifying the apparent physical meaning of each individual model will help place the relative importance of their spatial patterns in a wider context. As stated in section 3, it is not guaranteed that each mode will correspond to the field of a single current system. Where we are able to assign a probable source for a given pattern, we refer to the current system which appears to dominate that pattern. In section 4.2 we have assessed the modes' temporal oscillations, both visually and via correlation with independent indices of geomagnetic activity, and shown that modes 1 and 3 are likely dominated, respectively, by magnetospheric and ionospheric signals. Here we compound this evidence with an assessment of the spatial patterns of the external component of each SH model shown in Figure 6. The component is shown because it highlights the trends we discuss more clearly than the r or components. Note that in each panel of Figure 6 the data prediction in the polar regions merely results from the global nature of the SH expansion and is not indicative of the true external inducing field geometry there. We display the polar regions since this conveys the extent to which the coefficients of a given SH model will depend upon regions without data coverage.
The component of the background means (given by SH model RM) is shown in Figure 6a, where a positive perturbation indicates a westward directed current located above the Earth's surface. It is reasonable to expect that the mean field will have contributions from several current systems. Indeed, we see that a pattern SHORE ET AL. DECADAL PERIOD EXTERNAL FIELDS 5179 consistent with the magnetospheric ring current dominates its nightside, while the dayside is dominated by a pattern consistent with the ionospheric Sq (solar quiet) current system. Note that the polar mask means that fields in the polar regions are poorly resolved and not indicative of the true polar inducing fields.
The spatial pattern of the annual oscillation (AO) model is shown in Figure 6b. When the north-south derivative (not shown) is taken of the RM model pattern shown in Figure 6a, the result bears a good similarity to Figure 6b. We show (as highlighted in Figure S1 in the supporting information) that the AO model describes the north-south oscillation (on an annual timescale) of the background field represented by the RM model. Addition of the AO model (at a given phase of its temporal oscillation) to the RM model produces the appropriate seasonal motion of the background field. Since the AO model is an oscillation resulting from the motion of the Sun in the SM system, it is expected to entail a mixture of current contributions. From our analysis, it appears that the AO and RM models likely contain contributions from the same set of source currents.
The distribution of the SR model is shown in Figure 6c. It shows a high level of zonal (i.e., LT) invariance across the nightside, consistent with the pattern expected of the magnetospheric ring current. Thus, our earlier interpretation of the meaning of this pattern (based on its temporal oscillation) appears to be justified. The dayside part of the SR model is less zonally invariant than the nightside-we cannot confirm the origin of this aspect of the SR pattern, but we offer two possibilities. First, the SR pattern on the dayside could be due to a weakening (as measured at the Earth's surface) of the magnetic signal of the ring current according to the (eastward) Chapman-Ferraro currents [e.g., Cowley, 2007] at the dayside magnetopause. Second, since the mode 1 oscillation was shown (in Figure 5) to exhibit a weak annual periodicity, the dayside part of the SR model might reflect an incomplete separation of the annual oscillation pattern from the ring current by the EOF analysis.
The SA model pattern is shown in Figure 6d. The pattern has a good resemblance to the magnetic perturbations of the Sq current system, with minimal nightside signal. This ties in with the correlation of mode 3 (on which the SA model is based) with the F 10.7 index, and it therefore appears that Sq dominates the SA model. This interpretation is consistent with the semiannual portion of mode 3 arising from the modulation of the Sq system due to the Russell-McPherron effect, as discussed in section 4.2.
To summarize the spatial interpretation of the models, it appears that SR and SA are dominated, respectively, by the day-to-day amplitude modulation of the ring current and Sq systems, while both of these systems contribute approximately equally to the background mean field (RM) and its apparent annual north-south motion in the SM frame (AO). The significance of these results in light of the findings of previous studies is discussed in section 5.
SHORE ET AL. DECADAL PERIOD EXTERNAL FIELDS 5180  Figure 6 are based on these models. The full time series of internal and external coefficients up to degree and order 16 are given in the supporting information.
The external SH coefficients (from which the data predictions in Figure 6 are made) are shown for each model in Table 1, in units of nanotesla (nT) at their respective maximum-amplitude epochs, up to degree and order 5. The q 0 1 harmonic of the SR model is the strongest term, indicating that the symmetric ring current is (as expected) the dominant source in the inducing signal. We also resolve substantial (cumulatively strong) complexity in both zonal and nonzonal higher-degree coefficients, which are subject to long-period oscillations. The SR model's next strongest harmonics are the q 0 3 and q 0 5 terms, which are of similar (or lesser) SHORE ET AL. DECADAL PERIOD EXTERNAL FIELDS 5181  Figure 6, though a different color scale is used. Note that the data prediction used to compute the residual is the sum of the internal and external parts of the AO model. amplitude than the harmonics q 1 1 in the AO model and q 1 2 in the RM and SA models. Our results justify the use of a more complete description for the inducing field than just the symmetric ring current. An assessment of the likely improvement in mantle conductivity estimates as a result of the inclusion of this information is beyond the scope of our study.
We provide the internal Gauss coefficients in the supporting information. These are of smaller amplitude than the corresponding external coefficients, justifying our assumption that the EOF approach responds primarily to the inducing fields of interest.
An example of the SHA residuals-those of the AO model at its strongest epoch-is shown in Figure 7. The absolute values are the highest near observatory clusters (particularly on the dayside) but are generally low in amplitude and are not globally coherent; hence, we believe that the models are unbiased in their representation of SM frame patterns. The residuals for the other models and components are similarly distributed, though the polar mask does increase their amplitudes in high-latitude regions for those models which have most power in their low-degree coefficients.

Discussion
Based on the amplitudes of the coefficients shown in Table 1, it appears that the ring current is the most important pattern, both on short timescales and on the scale of a full solar cycle. The annual and semiannual oscillation model coefficients are smaller in amplitude but still important in describing the long-period external field variations. Although the application of a SHA in recombining the 24 start-UTs resolves each mode's oscillation into discrete spherical harmonics (each with rigid geometries), the SHA is sufficiently high degree (and the replicated data distribution dense enough) that this step does not impose significant geometrical simplifications on the solutions.
The compactness of the EOF description of the external fields is due to the correlation between the observatories' data in the chosen frame of reference. Thus, the coordinate system we use is an implicit prior assumption of our decomposition. This choice does not affect the completeness of the representation of the data's variance by the EOFs, but it does affect our physical interpretation of a given mode (remembering that these are data modes, which do not necessarily describe individual physical phenomena). Ultimately, our interpretation of the meaning of each mode is subjective but (we hope) intuitive. On the timescales covered by our analysis, our system of latitudes is effectively geographically fixed, creating an apparent motion of the Sun (and with it, the current systems). This apparent motion will be the same as that "perceived" by the conducting interior of the Earth; thus, we consider our choice of coordinate system to be reasonable.
It is worth commenting on the longest periods of the obtained modes in more detail. Recall that each identified mode has both a periodic oscillation (annual, semiannual, or the length of the solar cycle) and an aperiodic day-to-day variance, which vary in their relative importance to the composition of the mode's signal. For instance, we have stated above that the AO model describes principally a source region motion, while the SR and SA models describe principally a source region amplitude modulation (of certain parts of the static fields, represented by the RM model). It is likely that the apparent motion of the source regions on the scale of a solar cycle is minimal in comparison to their apparent seasonal motion (in our chosen coordinate frame). Hence, we expect that the solar cycle-related trends will be primarily amplitude modulations, contributing to both the day-to-day and periodic variances. This is probably why we see the solar cycle trends split across several modes. While it is beyond the scope of our study, if we were to take one "snapshot" of the data per year in the same manner as the start-UTs are selected once per day, we could likely compress the solar cycle variations into fewer modes, with minimal interference from the seasonal motion of their source regions.

SHORE ET AL. DECADAL PERIOD EXTERNAL FIELDS 5182
We have referred to the possible source regions of the magnetic signals multiple times and shown evidence to back up our physical interpretation, but we reiterate that the EOF method cannot be used to distinguish reliably between current systems. This is unfortunate, since the relative importance of the ionosphere to the inducing fields has been recognized since at least the studies of Banks [1969] and Malin and Işikara [1976]. More recently, Balasis et al. [2004] have analyzed CHAMP satellite data to show that the long-period inducing source has significant nonaxisymmetric (i.e., nonzonal) structure, with evidence to suggest that this nonzonal signal is magnetospheric rather than ionospheric. In response to this finding, Egbert [2006, hereafter BE2006] applied EOF analysis to 4 years of midlatitude ground observatory hourly mean magnetic data to search for evidence of nonzonal source fields. The BE2006 study is perhaps closest to ours, though the results are not directly comparable since BE2006 removed an estimate of the Sq and ring current signals and filtered their data to retain only periods of between 5 and 50 days. The dominant pattern resolved by BE2006 exhibited cross-nightside asymmetry in LT and had a semiannual temporal periodicity (additionally well correlated with the Dst index). We also resolve a semiannual pattern-indeed, we resolve nonzonal structure in each of the models we present. However, we find that the cross-nightside nonzonal structure (i.e., within the 12 h range centered on local dipole midnight) is always dominated by the contributions from the dayside. Furthermore, the dayside contributions appear as important (if not more so) to the inducing signal on periods of months and longer than the nightside contributions.
The more recent use of low-Earth orbit (LEO) satellite data in induction studies means that there is a tendency to parameterize only the magnetospheric inducing terms and to consider the other contributions as correlated noise. This is primarily because the ionospheric component is internal to LEO altitudes and thus is not distinguishable from the induced fields. Our findings confirm the prevailing view [cf. Püthe and Kuvshinov, 2014] that the symmetric ring current is the dominant inducing field. Furthermore, we show which additional inducing field harmonics are important to ground-based data at all LTs and have provided information of the behavior of the baselines of these inducing fields over the course of a solar cycle.

Conclusions
We have applied the EOF method to 13 years of ground-based observatory hourly mean data and extracted the dominant spatiotemporal patterns of the quiet time long-period external fields, without making explicit assumptions about their geometry. While the EOF method cannot distinguish between source regions, we find that the majority of the inducing signal is consistent with the pattern of the symmetric ring current, which is commonly assumed to be the sole inducing source. We also find significant contributions from patterns representative of other (nonzonal) current systems, which are important to the total signal on annual, semiannual, and solar cycle periods. These patterns collectively describe the variance of the data over the 13 year analysis span-we also describe the pattern of the ground state of the data, in the form of the mean field over the same period. We hope that these new estimates of the inducing fields will aid future studies of mantle induction and highlight the continued utility of ground-based observatory data.