Principal component analysis of Birkeland currents determined by the Active Magnetosphere and Planetary Electrodynamics Response Experiment

Principal component analysis is performed on Birkeland or field‐aligned current (FAC) measurements from the Active Magnetosphere and Planetary Electrodynamics Response Experiment. Principal component analysis (PCA) identifies the patterns in the FACs that respond coherently to different aspects of geomagnetic activity. The regions 1 and 2 current system is shown to be the most reproducible feature of the currents, followed by cusp currents associated with magnetic tension forces on newly reconnected field lines. The cusp currents are strongly modulated by season, indicating that their strength is regulated by the ionospheric conductance at the foot of the field lines. PCA does not identify a pattern that is clearly characteristic of a substorm current wedge. Rather, a superposed epoch analysis of the currents associated with substorms demonstrates that there is not a single mode of response, but a complicated and subtle mixture of different patterns.


Introduction
The Active Magnetosphere and Planetary Electrodynamics Response Experiment (AMPERE) [Anderson et al., 2000;Waters et al., 2001;Anderson et al., 2002Anderson et al., , 2008 has provided measurements of the Birkeland currents or field-aligned currents (FACs) in the Northern and Southern Hemispheres at 10 min cadence from 2010 to 2013, using magnetometer observations from the Iridium constellation of close to 70 satellites. A number of studies have employed AMPERE observations to study the structure of the current systems and their response to solar wind-magnetosphere coupling. For instance, Clausen et al. [2012Clausen et al. [ , 2013aClausen et al. [ , 2013b and Coxon et al. [2014aCoxon et al. [ , 2014b have demonstrated that the region 1 and 2 current systems [Iijima and Potemra, 1976] observed by AMPERE undergo cycles of expansion to lower latitudes and contractions to higher latitudes in tune with the substorm cycle, and Anderson et al. [2014] showed dayside followed by nightside activations of the current systems in response to southward turnings of the interplanetary magnetic field (IMF). Both sets of observations are consistent with the expanding/contracting polar cap (ECPC) model of the Dungey cycle of magnetospheric and ionospheric convection [Dungey, 1961;Siscoe and Huang, 1985;Cowley and Lockwood, 1992;Lockwood and Cowley, 1992;Milan et al., 2007;Milan, 2015]. In addition, Murphy et al. [2013] and Sergeev et al. [2014] have used AMPERE observations to explore the structure of nightside FACs during substorms, specifically attempting to elucidate the structure of the substorm current wedge [McPherron et al., 1973], while Wilder et al. [2013] studied cusp currents associated with magnetic tension forces during periods of nonzero IMF B Y .
The aim of the current paper is to use principal component analysis [Jolliffe, 2002] to determine the repeatable patterns of currents that make up the constantly varying observations. Principal component analysis (PCA) decomposes a data set into a series of basis functions that reveal the structure of the underlying correlations within the data. For instance, PCA has also been used for facial recognition [Sirovich and Kirby, 1987;Turk and Pentland, 1991] and it is potential for recognizing otherwise hidden structure within AMPERE current maps that motivates our study. Recently, PCA has seen wide adoption by the space and geophysics community [e.g., Ohme et al., 2013;Natali and Meza, 2010]. He et al. [2012] used an empirical orthogonal function decomposition (closely related to PCA) to investigate FAC patterns observed by the CHAMP satellite. Kim et al. [2012] applied PCA to polar cap ionospheric convection observations and found two dominant modes of response: uniform Red and blue indicate upward and downward FACs, between ±1 μA m −2 , respectively. The cross displaced toward midnight from the geomagnetic pole is the assumed centroid of the region 1/2 current system. The black circle of radius Λ fit is fitted (see below) for normalization purposes. (b) The summed currents around the circumferences of circles of radii Λ, centered on the cross in Figure 1a. Positive/negative bipolar signatures are identified (open circles), and the zero crossings (closed circles) are found. The largest bipolar signature is identified with the R1/R2 current system, and the radius of the corresponding zero crossing, Λ fit , is used for normalization. (c) The normalized current map, centered on the geomagnetic pole and stretched to a normalization radius of 20 ∘ , is indicated by the black circle.
antisunward flow, correlated with the B Z component of the IMF, and a mode which introduced a dawn-dusk asymmetry in the flows, correlated with IMF B Y . Cousins et al. [2015] performed PCA on AMPERE current maps and found two dominant components related to IMF B Z and B Y , the FAC counterparts of the convection modes of Kim et al. [2012], and a third related to expansions and contractions of the polar cap.
In this paper we perform PCA on AMPERE current maps that have been preprocessed to remove variations in the radius and center of the current ovals (related to the ECPC), so as to reduce smearing and to better resolve small-scale features in the patterns. We then investigate the response of the principal components to solar wind parameters and the occurrence of substorms.

Observations and Discussion
This study employs AMPERE observations of the Birkeland currents from the Northern Hemisphere only. Each AMPERE map covers the region poleward of 50 ∘ geomagnetic colatitude, on a 24 × 50 grid in magnetic local time and colatitude, centered on the northern or southern geomagnetic pole. Although AMPERE provides current maps with 2 min cadence, these are produced by a sliding 10 min average of the Iridium observations, so 144 independent maps are available per day. An example map is shown in Figure 1a. The region 1 and 2 (R1/R2) current system is clearly visible as upward/downward (red/blue) currents near a latitude of 70 ∘ . The polarity of the currents is opposite in the dawn and dusk sectors.
The size of the FAC pattern changes continuously due to expansions and contractions of the polar cap [Milan et al., 2003[Milan et al., , 2007Clausen et al., 2012]. Before undertaking PCA, we normalize the current patterns to a consistent size to remove this effect. Using a technique similar to that described in Milan [2009], we assume that the R1/R2 current regions are approximately circles centered on a point displaced a few degrees antisunward from the geomagnetic pole. For instance, in Figure 1a it is assumed that the center of the current systems (marked by a cross) is offset by 3 ∘ of latitude, Λ 0 , along the midnight meridian.
We determine the current strength integrated around circles of differing radii, centered on this point. Consider a circle of radius Λ. We find the mean of the current density, j, at 48 equally spaced points around the circumference of the circle, first multiplying currents in the dusk sector (12-24 MLT) by −1. In this way, if the circle coincides with region 1 currents, a positive value is measured, and a negative value is measured if the circle coincides with region 2 currents. The variation of this sum, Σj, with Λ is shown in Figure 1b, clearly showing a bipolar signature associated with the R1/R2 currents.
The largest bipolar signature in Σj is identified, and the latitude of the zero crossing, Λ fit , is used as the boundary between R1 and R2 currents. This procedure is repeated with values of Λ 0 between 0 ∘ and 5 ∘ , and the combination of Λ 0 and Λ fit which gives the greatest R1/R2 peak-to-peak variation is used as the fit. This is indicated in Figure 1a by the black circle. If the currents are too weak, such that the bipolar signature represents an average peak-to-peak current density of less than 0.15 μA m −1 , then the fit is considered unreliable and the map is discarded. Seventy-nine percent of maps were successfully fitted.
The current map is then stretched onto a 24 × 40 grid (with n = 960 elements), such that the reference circle is centered on the pole and has a radius of 20 ∘ . The resulting pattern is shown in Figure 1c. In all, 123,340 maps (m = 123, 340) from 2010 to 2012 were normalized in this fashion.
We performed principal component analysis on these normalized maps. Each map is represented by an n−dimensional vector, J. Each J is mean centered-that is, the mean of all the elements, 1 n ∑ n i=1 J i , is subtracted from each element-but otherwise no scaling is applied. The n × m matrix X represents all the observations, being a concatenation of all the Js. From this, the covariance matrix is computed as = 1 m X T X, where X T is the transpose of X. The symmetric n × n covariance matrix represents the correlations between variations in all n(n + 1)∕2 pairs of elements in the normalized maps. Eigendecomposition of is performed, resulting in n eigenvectors, F i , (each with n elements), and corresponding eigenvalues, i . The eigenvectors are the principal components of the data set, encoding the directions in n-space along which correlated variations in the data set are best described. The eigenvectors with the largest eigenvalues represent the majority of the variation in the data set. In the field of facial recognition, these eigenvectors are called eigenfaces. We will refer to our eigenvectors as eigenFACs.
These eigenFACs can be interpreted as maps indicating features which are consistently present in the data set. The 12 most significant eigenFACs are presented in Figure 2. The percentage of the variance in the data set described by eigenFAC F i is given by 100 i ∕ ∑ n j=1 j . Figure 2 (bottom) shows this contribution for the first 20 eigenFACs; the i = 1 level is shown by the horizontal dashed line. The cumulative variance explained by the first 300 eigenFACs is indicated in Figure 2 (bottom). For reasons that will be discussed below, we also computed the eigenFACs for just the months June and July 2010-2012 (m = 24251, 93% success rate), presented in Figure 3. The eigenFACs calculated from the observations for all months are labeled F i A (Figure 2), those from summer months are labeled F i S (Figure 3). The PCA technique decomposes the data set into as many eigenvectors as there are elements in each data vector, in our case n = 960. Which of these is significant in describing the data set? One of the aims of PCA is dimension-reduction: there are 960 elements in each FAC map, but can the gross structure of the data be described by just a few variables, and how many variables are necessary? The eigenvalue associated with each eigenvector is a measure of this significance. However, there is no hard-and-fast rule for determining the threshold for significance, though there are several commonly applied heuristics. Kaiser's criterion states that eigenvectors for which i ≥ 1 are significant [Kaiser, 1960]. In our case, this includes the first five or so eigen-FACs. The Scree test [Cattell, 1966] orders the eigenvalues, as in Figures 2 and 3, and then fits a straight line through the lower eigenvalues and retains eigenvectors whose eigenvalues rise above this line. Application of the Scree test to our data indicates that the first 20 or so eigenvectors should be considered significant. On the other hand, another measure of significance is the number of eigenvectors necessary to explain the majority, say 95%, of the variance in the data set [Hair et al., 1995]. In our case, approximately 70 eigenFACs are required to describe 80% of the variance and near 150 for 90% of the variance.
We first discuss the eigenFACs presented in Figure 2. The principal eigenFAC, F 1 A , corresponding to approximately 25% of the variance in the data set, has the form of a concentric pair of upward/downward rings of Birkeland current, of opposite polarity in the dawn and dusk sectors, rotated such that the polarity changes along the 11/23 MLT meridians. This clearly represents the region 1 and 2 current system described by Iijima and Potemra [1976]. The rotation of the current rings is consistent with the rotation often seen in the ionospheric convection pattern [e.g., Ruohoniemi and Greenwald, 1996]. The magnitude of the contribution of this eigenFAC, 1 A , to a given current map, J, will represent the strength of magnetospheric convection. It might be expected that 1 A will become large when the interplanetary magnetic field is negative, IMF B Z < 0, and low-latitude magnetopause reconnection is ongoing. We note that F 1 A is consistent with the first principal component found by Cousins et al. [2015] in AMPERE data and by Kim et al. [2012] in their decomposition of polar cap ionospheric flow observations. EigenFAC F 2 A , representing approximately 10% of the variance, is a bipolar pair of latitudinally separated currents, strongest in the noon sector. The presence of F 2 A , for which the coefficient 2 A could be positive or negative, will modify the F 1 A pattern such that the R1 currents in the noon sector become upward or downward, with opposite polarity currents at higher latitude. This is consistent with currents expected to be associated with magnetic tension forces on newly reconnected field lines for nonzero IMF B Y [Cowley et al., 1991]. F 2 A also contains a region of nightside R1 currents which possibly contribute to a dawn-dusk asymmetry in the convection flows on the nightside associated with IMF B Y (see, for instance, the discussion in Milan [2015]). Again, we note that the dayside portion of F 2 A is consistent with the second principal component found in the study of Kim et al. [2012]. It is also consistent with the third principle component found by [Cousins et al., 2015]; their second principal component was associated with expansions and contractions of the current ovals, an effect that we have removed by preprocessing the current maps. F 3 A has a similar form to F 2 A , though suggesting a more complicated, three-current structure on the dayside associated with IMF B Y . EigenFACs F 3 A and beyond contain significant currents on the nightside, in the main consisting of east-west aligned current sheets, which are most likely associated with substorm processes. The substorm current wedge is formed from upward and downward currents at the western and eastern edges of the substorm auroral bulge [McPherron et al., 1973]. More than one eigenFAC contains upward/downward current pairs straddling the midnight meridian, e.g., F 5 A , F 7 A , and F 8 A . As discussed above, these eigenFACs have low significance in comparison to F 1 A and F 2 A , suggesting that there is not one eigenFAC that represents substorm current systems, but that these current systems are made from a variable superposition of a number of basis functions. As will be demonstrated below, although these eigenFACs have relatively low significance, they respond in a reproducible way to substorm onset. Turning to Figure 3, the first three eigenFACs are very similar to their counterparts in Figure 2. However, eigenFAC F 4 S represents a quadrupolar region of currents near noon, consistent with reverse convection cells associated with lobe reconnection occurring for positive IMF B Z . This eigenFAC also has currents in the dawn and dusk sectors of opposite polarity to the usual R1/R2 sense. This indicates that when the IMF is northward and lobe cells appear near noon, the R1/R2 currents at dawn and dusk are reduced. F 6 S in Figure 3 seems likely to be associated with IMF B Y changes in the location of the lobe reconnection site and distortions of the lobe cells [see, e.g., Milan et al., 2000]. We note that although the R1/R2 current system of F 1 S is rotated from the noon-midnight meridian, most cusp-related features in Figure 3 are symmetric about the noon meridian. Other differences between the eigenFACs of Figures 2 and 3 will be discussed further below.
The eigenvectors produced by PCA form an orthonormal basis set, such that linear combinations of eigenFACs can be used to represent any current map. The contribution of a given eigenFAC, F i , to an individual map, J, can be found from the inner product of F i and J: i = ∑ n j=1 F i j J j , where J j are the elements of J and F i j are the elements of F i . Then That is, i is the projection of J along a particular basis vector F i . This procedure can be used to find the contribution of each eigenFAC to the observed currents in a particular map, and how these change with time or geomagnetic activity.
As discussed above, F 1 A is expected to be associated with large-scale convection, so 1 A should be a measure of convection strength, and hence dependent on the activity of the substorm cycle. Figure 4a shows the contribution, 1 A , of F 1 A to each current map J from 2010, as a function of dayside reconnection rate Φ D , determined using the parameterization of Milan et al. [2012] and 1 min OMNI solar wind data [King and Papitashvili, 2005]. Figure 4b shows 1 A as a function of the absolute value of the AL electrojet index, which we use as a proxy for the nightside reconnection rate Φ N [see also Holzer et al., 1986and Coxon et al., 2014a, 2014b. 1 A increases with both Φ D and Φ N , though the dependence is rather weak for Φ D . The expanding/contracting polar cap model Milan et al., 2007] suggests that the cross-polar cap potential should be dependent on both dayside and nightside reconnection. Figure 4c shows the mean value of 1 A as a function of both Φ D and |AL|, confirming that 1 A increases for both, consistent with the findings of Coxon et al. [2014a]. We expect 2 A to be associated with tension forces due to the B Y component of the IMF. Figure 5 presents this dependence for 2010. In this case we find a strong seasonal variation in 2 A , so we have separated the data by month. At all times there is a clear correlation, but the magnitude of 2 A is larger in summer months than in winter months. This suggests that the conductance of the ionosphere at the footprint of newly reconnected field lines significantly modulates the magnitude of the currents associated with tension forces. It is this seasonal dependence that motivated us to produce the eigenFACs for June and July shown in Figure 3. Comparing the F 2 A to F 5 A (Figure 2) with F 2 S to F 5 S (Figure 3), it is clear that noon sector currents are more pronounced with respect to nightside currents in summer months (Figure 3) than in the eigenFACs computed from the whole data set (Figure 2).
EigenFACs with significant nightside currents are expected to respond to the occurrence of substorms. As already discussed, it seems that no one eigenFAC captures the substorm dynamic, otherwise such an eigen-FAC would be expected to have a relatively large eigenvalue. Rather, although Kaiser's criterion suggests that eigenFACs beyond F 5 A have low significance, they contribute in a collective manner to describe the variance observed in different substorms. Here we test if the eigenFACs show some consistency in their response to substorm onset. We do this by performing a superposed epoch analysis of the i A s, in which we use approximately 3000 substorm expansion phase onsets as determined by SuperMAG [Newell and Gjerloev, 2011] as the zero epoch. The results are presented in Figure 6.
Figure 6 (first panel) shows the variation of Λ fit , the radius of the current ovals used in the normalization procedure. The ovals expand to lower latitude prior to onset, continue to expand until 20 min after onset, before contracting to higher latitudes again, consistent with the ECPC. Figure 6 (second to sixth panels) presents the i A s which show the largest variation; vertical bars indicate the standard error on the mean.

1
A shows a small increase prior to onset (convection associated with dayside reconnection and the growth phase) and a large increase after onset (associated with convection excited by nightside reconnection). Both Λ fit and 1 A results are consistent with the findings of Coxon et al. [2014b].

Projections 2
A to 4 A are not shown in Figure 6 as their variation with substorms is small. This is expected as they largely describe the effects of the dayside tension forces associated with IMF B Y , which should not vary consistently with substorm cycle.
The other i A s all show variations associated with the growth and/or expansion phases. For instance, 5 A is close to zero before the growth phase but becomes increasingly negative during the expansion phase. The associated eigenFAC, F 5 A , has currents on the nightside that are of opposite sense to the R1/R2 current system represented by F 1 A , suggesting that F 5 A when subtracted from F 1 A represents an enhancement in nightside R1/R2 currents during substorms. This is consistent with an enhancement of nightside convection associated with magnetotail reconnection and the enhancements in currents due to conductance changes in the nightside auroral zone. Projections 6 A and beyond also show significant substorm-correlated changes. The corresponding eigenFACs have increasingly structured currents on the nightside which will contribute to the form of the substorm current wedge. Figure 6 suggests that the substorm FAC structure evolves rapidly after substorm onset, especially with respect to the 10 min cadence of the AMPERE observations: for this reason, a single eigenFAC cannot represent the substorm current wedge at all substorm phases. Moreover, despite the low significance attributed to principal components beyond F 5 A , their consistent response to substorm onset suggests that they do represent repeatable structures within the nightside current systems. In part, their low significance will arise as a consequence of the small proportion of time they are active, whereas the R1/R2 currents and dayside tension forces (F 1 A and F 2 A ) will be present much of the time.
The ECPC predicts that magnetospheric convection is associated with the action of both dayside and nightside reconnection (in fact, that the cross-polar cap potential is the average of Φ D and Φ N ) [Lockwood, 1991]. However, two of our observations seem to somewhat contradict this: in Figure 4a the dependence of 1 A on Φ D is weaker than for |AL|, and in Figure 6 the increase in 1 A is much more significant after substorm onset than during the growth phase. We interpret this as the control of current magnitudes by the nightside ionospheric conductance, which will be enhanced during substorms, but not necessarily during the growth phase.
Finally, we discuss some similarities and differences between the eigenFACs calculated for the whole period 2010-2012, F i A (Figure 2), and for the months June and July alone, F i S (Figure 3). F 1 , F 2 , and F 3 are very similar in both sets of eigenFACs, though we note that the variance represented by F 2 S and F 3 S is relatively larger than F 2 A and F 3 A , recognizing the enhanced cusp currents when the dayside ionospheric conductance is high. The nightside currents in F 3 S are weaker than in F 3 A , again representing the dominance and consistency of noon currents in summer.
It becomes more difficult to see direct relationships between other pairs of eigenFACs in Figures 2 and 3 as the eigenvalues decrease. Indeed, many of these eigenFACs appear to share common features of dayside/nightside currents, but the contributions differ between the two data sets. This reinforces the point discussed above that substorm FACs do not have one dominant mode or spatial pattern but change in a subtle but complicated way from substorm to substorm. Hence, no single mode of response is identified by the PCA technique. Despite this, it is anticipated that a more in-depth study of these eigenFACs could lead to a greater understanding of substorm processes.
There are two additional factors to consider in future studies: the dependence of substorm behavior on the latitude of substorm onset and the IMF B Y dependence of magnetotail dynamics. Milan et al. [2009] demonstrated that the auroral signature of substorms is enhanced when the substorms occur on an expanded auroral oval (greater proportion of open magnetic flux at onset), and Grocott et al. [2009] showed corresponding changes in the ionospheric convection response to substorms. A similar dependence may be expected to occur in the contributions of different eigenFACs during substorms of differing onset latitudes. It is also known that the local time of substorm onset can be influenced by the B Y component of the IMF some time prior to onset [e.g., Østgaard et al., 2005;Milan et al., 2010], and this might be reflected in the eigenFACs.

Conclusions
We have presented a principal component analysis (PCA) of Birkeland current measurements from the AMPERE experiment. The technique reveals that the region 1/2 current system first identified by Iijima and Potemra [1976] is the most reproducible feature of the data set. Dayside currents associated with magnetic tension forces on newly reconnected field lines are the second most consistent feature. The strength of these cusp currents is significantly modulated by season, being strongest in summer months when the ionospheric conductance at the footprint of the cusp region is greatest. This will introduce a significant interhemispheric asymmetry in the currents flowing on the dayside magnetopause. Cusp currents associated with lobe reconnection and reverse convection are also apparent in the decomposition during summer months. In contrast, nightside current patterns are complicated and do not have a consistent response to substorm onset-there is no easily identifiable substorm current wedge-though the magnitude of the R1/R2 currents is enhanced.