AE, DST, and Their SuperMAG Counterparts: The Effect of Improved Spatial Resolution in Geomagnetic Indices

For decades, geomagnetic indices have been used extensively to parameterize space weather events, as input to various models and as space weather specifications. The auroral electrojet (AE) index and disturbance storm time index (DST) are two such indices that span multiple solar cycles and have been widely studied. The production of improved spatial coverage analogs to AE and DST is now possible using the SuperMAG collaboration of ground‐based magnetometers. SME is an electrojet index that shares methodology with AE. SMR is a ring current index that shares methodology with DST. As the number of magnetometer stations in the SuperMAG network increases over time, so does the spatial resolution of SME and SMR. Our statistical comparison between the established indices and their new SuperMAG counterparts finds that, for large excursions in geomagnetic activity, AE systematically underestimates SME for later cycles. The difference between distributions of recorded AE and SME values for a single solar maximum can be of the same order as changes in activity seen from one solar cycle to the next. We demonstrate that DST and SMR track each other but are subject to an approximate linear shift as a result of the procedure used to map stations to the magnetic equator. We explain the observed differences between AE and SME with the assistance of a simple model, based on the construction methodology of the electrojet indices. We show that in the case of AE and SME, it is not possible to simply translate between the two indices.


Introduction
Long-term geomagnetic indices have proven to be key parameters for the space weather community. AE and D ST are geomagnetic indices based on ground-based magnetometer observations. Both AE (Davis & Sugiura, 1966) and D ST (Sugiura, 1964) have been almost continuously recorded over multiple solar cycles. AE indices have been widely used to characterize magnetospheric substorm events (e.g., Hajkowicz, 1991;Wik et al., 2009) and to study physical processes in the magnetosphere-ionosphere system. For example, Echer et al. (2006) used the AE index in a comparison of the geoeffectiveness of different solar wind interplanetary structures, and Weimer (1994) used the index to study the time evolution of substorm disturbances. Many studies have identified statistical behavior in the distributions of AE (Consolini & De Michelis, 1998;Hnat et al., 2004;Hush et al., 2015;March et al., 2005;Takalo et al., 1993;Uritsky et al., 2001). The D ST index is an effective index as it responds to the strength of the equatorial ring current. D ST is also a function of magnetopause currents and tail currents (Burton et al., 1975). Acting as a monitor for geomagnetic storms, it is commonly used to classify storms (Bisi et al., 2010;Gonzalez et al., 1994;Gopalswamy et al., 2005) 10.1029/2020JA027828 and to track the evolution of the storm activity over multiple cycles (Banerjee et al., 2012;Love et al., 2015;Shen et al., 2017;Yermolaev et al., 2013). Extreme value theory and risk analysis have been applied to D ST in order to predict the likelihood of extreme geomagnetic activity (Acero et al., 2018;Riley, 2012;Silrergleit, 1996;Tsubouchi & Omura, 2007). Both AE and D ST are regularly considered as benchmarks for the specification of the overall magnetosphere-ionosphere system as they capture overall activity (Chapman et al., 2018;Lockwood, 2013).
There has long been concern about the limitation of the indices (Rostoker, 1972). In the original derivation of the AE index, Davis and Sugiura (1966) posed the question of whether the network of 12 stations was of high enough spatial resolution to capture the dynamic behavior of interest. D ST relies upon an assumed symmetry of the ring current system, and this is no longer explicitly accepted (Campbell, 2004;Greenspan & Hamilton, 2000).
Using a global network of magnetometers, the SuperMAG (Gjerloev, 2012) collaboration compiles the SME (Newell & Gjerloev, 2011) and SMR (Newell & Gjerloev, 2012) indices with as many stations as are available from the appropriate region. SME and SMR have been introduced as high spatial resolution counterparts to AE and D ST , respectively. They share methodology, but, for recent solar cycles, the number of magnetometer stations used is typically 10 times larger than those used in the original indices. The use of all available stations leads to improved timing, intensity, and location specification of events. This is simply due to the higher probability of having a ground station at the right location at the right time to monitor the overhead ionospheric current.
In order to elucidate the effect of increased spatial resolution in these indices, we examine the behavior of the indices over time. Geomagnetic activity is modulated by the solar cycle (Campbell, 1979), and each of these cycles varies in activity and duration (Hathaway, 2015). In this paper, we explore the influence of spatial resolution on the AE and D ST geomagnetic indices with a particular focus on how well they capture large-to-extreme events. In section 2, we describe the indices and highlight their differences. We characterize the solar cycle variation of the indices and find that the distributions of electrojet indices, AE and SME, diverge for severe events. However, the distributions of ring current indices, D ST and SMR, have the same distribution functional form and can be mapped onto each other by a linear shift in the mean and variance. In section 3, we demonstrate that the overall differences between AE and SME can be captured with a simple model that incorporates the spatial density of stations that changes with epoch. Finally, we present our conclusions in section 5.

Data Sets and Preparation
The official IAGA-approved AE index (Nose et al., 2015a) is produced at 1-min cadence using data from up to 12 magnetometer stations at latitudes that correspond to the average location of the auroral oval. We use the final AE index for the years 1975-1988 and the provisional AE index for the years 1989-2018. The monthly baseline at each station is defined as the average horizontal (H) component recorded for the internationally defined five quietest days of the month. Once a baseline has been removed, following procedure outlined by Allen (1972), the H components of each station are compared. The auroral upper index, AU, is defined as the record from the station with the largest positive H-component disturbance. The auroral lower index, AL, is the record from the station with the largest negative H-component disturbance. The AE index is then the difference between the upper and lower values, AE = AU − AL.
SuperMAG now produces SME, an equivalent to AE, at 1-min cadence. SME is the difference between upper (SMU) and lower (SML) indices, that is, SME = SMU − SML. SMU and SML are based on the H-component measured at stations in the latitudes of the auroral oval, with baseline removal carried out. The SuperMAG baseline does not depend on quiet days but is instead determined from the data itself. The baseline takes the form of a daily baseline, a yearly trend, and a residual offset (Gjerloev, 2012). The key difference between AE and SME is the number of stations used in their derivation. While AE uses 12 stations, the number of stations used to derive SME increases with time. The index has been retrospectively produced from 1975. As seen in Figure 1, in the earliest years available, the number of stations used to derive SME is of the same order as the number used to produce AE, in some cases, even fewer stations are used. However, for the last solar maximum, the SME index is based on data from over 110 stations. As the AE and SME indices are based on the maximum and minimum values measured by the network, the location of the station that contributes records to AU/SMU and AL/SML varies with time of day and geomagnetic activity.
The official IAGA-approved D ST index (Nose et al., 2015b) is available at 1-hr cadence and uses data from four magnetometer stations located around the equator. We use the final D ST index for the years 1967-2014, the provisional D ST index for the years 2015-2016, and the real-time D ST index for 2017. Secular variation and solar quiet daily variation are removed via a baseline removal procedure based on the five quietest days of the month, as described by Sugiura and Kamei (1991). The average of all four H-component disturbances is calculated. To correct for the effect of geomagnetic latitudes of the stations, this average H-component disturbance is scaled by a factor of 1/cos( ), where is the average geomagnetic latitude of the four stations.
The SuperMAG equivalent to D ST is SMR. It is produced at 1-min cadence. SMR is calculated using H-component disturbances from as many equatorial stations as possible, with SuperMAG baseline removal applied (Gjerloev, 2012). The stations are separated into four magnetic local time (MLT) zones. A partial ring current index is defined for each MLT zone; these are SMR 00 , SMR 06 , SMR 12 , and SMR 18 . The subscripts indicate the central MLT of the corresponding zone. Each partial ring current index is calculated using the same procedure. For example, SMR 00 uses each available station in a 6-hr MLT zone centered at midnight. SMR 00 is the average H-component disturbance measured by these stations, scaled by a factor of 1/cos( ), where is the average latitude of the stations. As stations change MLT zone with the Earth's rotation, the number of stations used to calculate each partial ring current index varies with time of day. The SMR index is the average of the four partial ring current indices, that is, SMR = SMR 00 + SMR 06 + SMR 12 + SMR 18 ∕4. Averaging each sector rather than the available stations is done to eliminate the UT dependence of the SMR index. Overall, the number of stations used in the derivation of the SMR index increases with time ( Figure 1); it is below 10 until 1985 and increases to a typical value of 100 stations for the last solar maximum. We produce a reduced SMR data set for comparison with the hourly D ST index. As in Wanliss and Showalter (2006), this consists of one average SMR value per hour. Solar cycles vary in duration and intensity. This is reflected in the distributions of geomagnetic indices. We present results for periods of solar maximum selected using the definition in Chapman et al. (2018); the solar cycle maximum is 3.5 years in duration and begins 2.5 years after the previous minimum. The dates of the solar minima used for the last four solar cycles are determined from the smoothed Monthly International Sunspot Number provided by SILSO World Data Center (1976); they are the first of the following months: March 1976, September 1986, May 1996, and December 2008. These periods are designed so that they are long enough to capture a statistically significant number of extreme events while being short enough that  the data set may be treated as quasitime stationary. We use samples centered on each solar maximum as this sample captures the dynamic range of the geomagnetic indices. Differences between samples of the indices at periods of weaker solar activity may be more difficult to distinguish.

Characterization of Solar Cycle Variation in Indices
To give an example of how the indices perform, in Figure 2a, we plot the AE and SME data series as recorded over the course of the Bastille Day event of July 2000. Overall, SME can be seen to record higher values than AE; however, the difference is not uniform. Figure 2b plots the D ST and SMR indices recorded for the same time period. D ST and SMR can be seen to track each other closely. We quantify the similarities and differences through the rest of this section and in section 3.
In order to investigate the effect of increased station numbers in AE and D ST , we compare the statistical distribution of observations of the established indices and their new SuperMAG counterparts. We compare the empirical survival distribution, S(x), of the indices, where x may indicate any of the four indices to be studied.
emphasizes the large-to-extreme values in the data set and is convenient in that it does not require any binning of the data. For each data set, the empirical survival distribution plots the rank-ordered observations on the x axis, that is, x 1 , … , x k , … , x N , where N is the number of observations, k is the rank order, and x N is the largest value. The y axis plots the corresponding values of S(x) = 1 − k∕N. We use the formula by Greenwood (1926) to calculate uncertainties in the survival functions, and these uncertainties are indicated by shading in the plots.
We use samples from the last four solar maxima to form survival distributions. Figure 3 shows the variation of the individual indices from one solar cycle to the next. Comparing AE and SME (Figures 3b and 3c), we see that in both cases, Cycles 21 and 24 display weak activity and have similar distributions. The comparison shows a larger difference between midrange values (1,000-2,000 nT) of Cycles 21 and 24 in the SME observations than in AE. These differences are not as significant as those seen for Cycles 22 and 23. While AE records more activity during Cycle 22, for the majority of the observations, SME indicates a similar level of electrojet activity during Cycles 22 and 23. Here, we also see a difference in the case of the most active cycle recorded by the indices. AE records the most extreme events during Cycle 22, while SME exhibits the largest readings during Cycle 23. Thus, the use of either index in a space climatological study would yield a different result. We examine the difference between the two indices quantitatively in section 3. Considering Figures 3d and 3e, we see the survival distributions of −D ST and −SMR track each other in how the indices compare from one cycle to the next. We quantify the relationship between D ST and SMR later in this section. Figure 4 overlays the distribution of AE and SME for each individual solar cycle maximum. The left-hand panels plot probability density function (PDF). The main panels show the PDF on a linear plot with a truncated y axis, and insets are a full semilog plot. The linear PDF emphasizes the response of the index to moderate disturbances, which are most common in the data set. The full data set is seen in the semilog PDF. The right-hand panels plot the survival distributions, S(x) = 1 − C(x). The survival distributions emphasize the tail of the distribution that consists of responses to large-to-extreme disturbances. In Figure 4, we compare the distributions of the electrojet indices, AE and SME. We see that the difference in the response of AE and SME is not consistent between solar cycles. For the early Cycles 21 and 22, we see that there is some crossover in the PDFs (Figures 4a and 4c) and the survival distributions (Figures 4b and 4d). In Figure 4a (inset), we see that in the region of the plot <1,500 nT, the observed AE index values are larger than the observed SME index values. This is not seen in any of the other solar cycles. For later Cycles 23 and 24, we see a significant shift between the PDFs (Figures 4e and 4g) and the tail of the SME survival distributions extend further than the tail of the AE distributions (Figures 4f and 4h). For these cycles, the AE index systematically undersamples when compared to SME.
We compare the distributions of −D ST and −SMR in Figure 5. For all four solar cycles, we see similar behavior. The PDFs in Figures 5a, 5c, 5e, and 5g show a shift between −D ST and −SMR, which is consistent across solar cycles. −SMR has a higher peak at the center of the distribution. The −D ST distribution then crosses over to display a slightly longer-tailed distribution. The shift is more clearly seen in the survival distributions in Figures 5b, 5d, 5f, and 5h. The distributions track each other remarkably well; however, they do not lie precisely over one another. Note also that the plot is on a semilog axis, so the shift is more pronounced than one may initially imagine. Overall, the distributions track each other with a consistent systematic shift.
The data quantile-quantile (QQ) plot can be used to determine whether or not two data sets share the same underlying distribution (Embrechts et al., 2013). Where two data sets contain the same number of observations, the data may be arranged in ascending order, and these ordered data sets may be plotted against one another, one on the x axis and one on the y axis. Where observations share the same underlying distribution, one expects the QQ plot to be a straight line along x = y. In the case where the data sets do not contain the same number of observations, a constant width Gaussian kernel density estimator may be used to evaluate a particular number of quantiles of the data set (Braun & Murdoch, 2016). As the array of quantile values are all of the same size, these may now be plotted against one another to generate a QQ plot. Where observations share the same underlying distribution, one expects the QQ plot of the data quantiles to be a straight line along x = y. If one of the distributions has undergone a linear transformation, the plot will be transformed by the same transformation. A difference in slope or y intercept of the QQ plot when compared to the x = y diagonal is indicative of such a linear transformation. The QQ plot can also be used to identify different regimes of behavior in a distribution. For example, where two linear regions exist in a plot, one may be indicative of a core component of the distribution where the second component corresponds to the tail in  Figure 4. The full range of data is shown in Figure 3. the distribution. This kind of dual-distribution QQ plot is common in space weather variables where the core of the distribution corresponds to small-scale turbulent behavior and the tail of the distribution identifies large-scale driven behavior (Tindale & Chapman, 2017). Here, we will compare the distribution of −SMR to −D ST using the QQ plot. We will return to the comparison of AE and SME index distributions in section 3.
In Figure 6, we show a QQ plot of −SMR versus −D ST ; 1,000 quantiles have been calculated using kernel density estimation. For all four solar cycles, the QQ plot indicates that −D ST and −SMR share the same functional form as the plots are almost straight lines. Note that there is some deviation from the single linear behavior in some cases, particularly noticeable in the small values (−D ST < 25 nT) of Cycle 22 in   Figure 6e, we see that a change in slope relative to the diagonal indicates that one of the distributions has undergone a linear transformation at larger values. Where such differences in moments of a distribution of observations occur, the distributions rescale onto one another by the mean and standard deviation of each individual data set. Such an operation can be used to confirm that −SMR and −D ST are truly the same subject to only a linear shift.
where z q is the baseline-corrected magnetic disturbance at station q. SMR p , the partial ring current index for MLT zone p, is constructed as follows: where N p is the number of stations in MLT zone p and z i,p refers to the baseline-corrected magnetic disturbance at station i in p. p is the average geomagnetic latitude of the magnetometer stations in p. Note that Different stations are used to generate D ST and SMR so that the value of m is also different, D ST ≠ SMR . This can introduce a linear shift to the observations, that is, if a factor applied to D ST is z and that to SMR is z ′ , the shift will be z ′ = Az + B. represented. On the right-hand side, the spatial current profile is shown. Fixed magnetometer stations are randomly distributed in the x-y plane. The number of stations used to calculate a model AE index is varied, for example, 10 stations (blue, square) followed by 100 stations (red, circle). The center sketch shows the y dependence of the current density, indicating the parameters of its Gaussian profile. The parameters w i and c i specify the current density profile width and the y coordinate of the peak. A i is the area integrated Gaussian profile that scales the total current.
For each cycle, the lowest moments of the distributions of −D ST and −SMR, that is, mean, , and variance, , are used to rescale the observations, x, to (x − )∕ . In Figure 7, the rescaled PDFs and survival distributions are plotted. We see that the simple linear shift accounts for much of the systematic shift between the two indices. Considering the PDFs (Figures 7a, 7c, 7e, and 7g), the distributions show the same peak behavior, within uncertainties. The survival distributions (Figures 7b, 7d, 7f, and 7h) lie more precisely over one another, within uncertainties. The small remaining difference between the two distributions, seen in the crossing of the survival functions, is therefore due to a change in the higher moments of the distributions. Sources of such a nonlinear shift may be the different baseline removal procedures used in the calculation of the indices or the variation of stations within the MLT zones defined for SMR. While the results of our rescaling procedure support a linear relationship between the combined SMR, which sums over the partial SMR p values, and D ST , we note that the individual SMR p values can differ from each other during active times, consistent with an asymmetric ring current (Newell & Gjerloev, 2012).

Simplified Model of Increased Station Number in AE
The QQ plot analysis suggests that the distributions of D ST and SMR are the same subject to a linear transformation. On the other hand, when we compare the distributions of AE and SME using the QQ plot in Figure 8, we find that there is no such simple linear translation. As these QQ plots are not linear, the distributions of observations do not share the same functional form. There are a number of factors that could contribute to the discrepancy between AE and SME, such as the number of stations, their geographical location, their calibration, relative sensitivity and availability, and the method for baseline removal. In particular, the number of stations used to form SME increases monotonically with time, and this may introduce a time (epoch)-dependent variation in the observed space climate that will bias comparisons between solar cycles as shown in Figure 3. To test this hypothesis, we will now develop a simple model to explore the impact of increased station number when we compare the AE and SME indices.
We construct a simplified model of a current system and magnetometer stations. A schematic of the model is shown in Figure 9. A square domain is defined of size 1 × 1 units. Magnetometers are placed in the grid; the x and y coordinates of the stations are drawn from a uniform random distribution. These magnetometer locations remain fixed. We use a simple configuration for a single current sheet over the square plane. Each sheet realization, i, varies in space, y, and time t.
The current sheet variation in space is as follows. The sheet current is infinitely thin in the z direction and is at uniform height above the ground, z 0 . The current, J, is constant in the x direction and follows a Gaussian profile in the y direction. This reflects that the auroral electrojet can be reasonably represented by a Gaussian profile in the latitudinal direction (e.g., Ahn et al., 1984Ahn et al., , 1986Knudsen et al., 2001;Sun et al., 1993).
For each current sheet realization, i, the current varies in space as We will generate an ensemble of many realizations of this current sheet model by drawing each of the model parameters from a distribution at random. We will choose distributions for the random variables that are input to the model that are well known, are few parameter, and are the simplest available that support the relevant properties of the model. For each realization of the current sheet, i, there are three model parameters. w i determines the extent of the current sheet in the y direction; it is positive definite and its distribution has finite moments. We present results where w i is drawn from a log-normal distribution with mean = 0.1 and variance = 1. The average current width is 10% of the domain width. A i determines the height of the Gaussian profile and thus controls the current sheet intensity. It is again positive definite, and its distribution has finite moments. We present results where A i is drawn from an exponential distribution with mean = 0.1. c i defines the location of the center of the Gaussian profile on the y axis; it is drawn from a uniform random distribution such that all locations on the grid are equally probable.
Many studies have examined the evolution of substorm activity in time, for example, Weimer (1994) compare the decay of substorm activity to an exponential curve and Baker et al. (1986) describe the response of the magnetosphere as being on two separate timescales, a relatively short-driven response and a slower "unloading" response. Bargatze et al. (1985) associate the fast-driven response with strong activity levels and the slow unloading process with more moderate activity levels. These findings inform the time evolution of the current sheets in this model. Each current sheet realization lasts a duration of 6 hr. The model assumes a simple pulse behavior for each sheet, that is, a single pulse that rises instantaneously and decays exponentially. J i,0 is the initial current of a sheet realization i at time t = 0, that is, J i,0 = J i (t = 0). The current sheet decays in time as J i (t) = J i,0 exp(−t∕ i ). i is a decay constant defined for a given storm. We initially examine a case where i is the same for every sheet realization. Later in this section, we will examine the case where i is dependent on J i,0 . The current sheet in space is then Taking ∇ × B = 0 J where z∕ t is negligible, the H component of the magnetic field, B H , as it varies in the y direction may be defined as The B H value at all magnetometer stations is calculated. The modeled index value for a given current sheet realization at a given time is then defined as the maximum B H value obtained at any station in the domain. We generate 3 × 10 3 current sheet realizations. Taking 6 hr of 1-min measurements, we obtain a distribution of 1.08 × 10 6 modeled index values. The model distribution sample size is then comparable with the sample size of observations during one solar cycle maximum.
We compare the survival distributions of modeled index values for different numbers of stations in Figure 10. Figures 10a and 10b show that as the number of stations in the grid increases (i.e., average station spacing decreases relative to the average current sheet width, <w i >), the modeled index resolves larger magnetic √ 60 = 0.13. This is of order the average width of the model current sheet (w i = 0.1). As we would expect, once the average station spacing is smaller than the average current sheet width, the spatial resolution of the network is sufficient to resolve the current system. Figure 10c shows this more clearly. As the number of stations increases, so too does S(x = 0.15), the fraction of observations larger than x = 0.15, determined from Figure 10a, where = 30 min. The increase in the median values of the box plot is exponential, saturating as spatial resolution increases.
Comparing Figures 10a and 10b, we demonstrate the effect of decay constant on the model AE index. Each current sheet lasts for 6 hr but decays with a decay constant of 30 min (Figure 10a) or 180 min (Figure 10b). In both cases, the distributions of the model indices have a long tail shape and are close in appearance to those observed for the real AE and SME indices (Figure 4). It has been shown that substorm behavior exhibits two characteristic timescales (Baker et al., 1986). Large and extreme storms will decay much more quickly than average substorms. We can incorporate this behavior into our model with a dependence of i on J i,0 . A threshold J T is defined at a quantile, q T , of the distribution of initial currents. Figure 11 shows the behavior of the modeled indices where i depends on the initial current J i,0 . A threshold J T is determined by the quantile q T . For the majority of the current sheet configurations, with J i,0 < J T , i = 3 hr. Where J i,0 > J T , the exponential decay constant is faster, with i = 0.5 hr. We see that this introduces a "kink" in the distribution. As the quantile of the threshold is increased and fewer storms decay at the faster rate, the kink moves to larger values and becomes more pronounced. This feature of the distribution is comparable to the distributions of AE and SME observations. Figure 12 is the QQ plot of the modeled indices with increasing station numbers against the modeled index with the highest station number, AE(110). We take AE(110) to approximate the ideal modeled AE index distribution, based on the convergence of the distribution tails in Figure 10. Similar to the results of the analysis of the AE and SME index observations (Figure 8), we see that the QQ plots for AE(10) or AE(35) do not follow the straight line diagonal that is seen for larger station numbers. Where the modeled AE network does not contain enough stations to resolve extreme values, there is no linear relationship between the modeled index and other modeled indices with a higher number of stations in the network.

Discussion
Figures 5-7 show that there is a simple scaling relationship between the D ST and SMR indices so that studies based on one index may be quantitatively related to studies based on another. There is already a considerable literature in which D ST is routinely used to parameterize the overall intensity of geomagnetic storms, but on the other hand, SMR offers the prospect of improved resolution. Our results identify a procedure to obtain the read-across between D ST and SMR.
We have found that there is no simple relationship between AE and SME; indeed, our results suggest that the driving factors in how they differ depend on the epoch (solar cycle) of interest. The AE index has been consistently based on data from 8-12 ground-based magnetometer stations. The SME index calculated for Cycle 21, the solar cycle maximum that began in 1978, was based on less than 20 magnetometer stations. When we compare the AE and SME index distributions for Cycle 21 maximum (Figures 4a and 4b), we find that the AE distribution samples exceed those of SME except at the largest values. Comparing Figures 4a  and 4b to the distributions for more recent solar cycle maxima in Figure 4, we see that Cycle 21 is the only solar cycle maximum for which the AE index significantly and systematically exceeds the SME index. This is qualitatively consistent with our simple model, which shows that as the station number is increased, the survival distribution moves to the right ( Figure 10) as the typical "observed" values increase in magnitude. In terms of our modeled indices, for Cycle 21, both AE and SME would be comparable to the AE(10) index in Figure 10; they should share the same distribution. This suggests that the observed differences between the AE and SME index during Cycle 21 may be attributed to their different baseline removal procedure and/or that different stations are used to derive the index, the geographical locations of the stations, their calibration, relative sensitivity, and/or availability. For more recent solar cycles, where SME is composed of significantly more stations than AE, our model suggests that we should see that SME is systematically returning higher values compared to AE, and this is indeed what we see in Figure 4. We can then attribute the difference between AE and SME to the improved spatial sampling that SME can achieve in more recent solar cycles.
Our results do not straightforwardly suggest the exclusive use of either AE or SME. For studies over a time interval within a given solar cycle, the higher-resolution SME index is a more appropriate auroral electrojet index to use. The higher spatial resolution of the magnetometer network used in the index construction will be more sensitive to the extremes of the auroral electrojet behavior. The AE index may underestimate the true level of activity. In contrast to this, for space climatological studies over multiple solar cycles, there will be a trade-off between the consistency of station number in AE compared to the improved spatial sensitivity of SME, which increases from one solar cycle to the next as the number of SME stations increases.

Conclusions
We have conducted a statistical comparison of AE and D ST with their SuperMAG counterparts across all data for the last four solar maxima. To find the difference between the traditional and updated indices, we identified how they vary across multiple solar cycles and compared the data sets for individual solar cycles. We found that AE systematically undersamples when compared to SME. The difference between the two indices is of the order of difference seen within a single index between different solar cycles. D ST and SMR were shown to track each other with a small systematic shift in the observations. By linearly rescaling D ST and SMR to their cycle specific sample mean and standard deviation, we have shown that the systematic shift is linear. This shift may arise since the operation that is used to map the baseline-corrected H-component disturbances at D ST and SMR stations to the magnetic equator is linear. The remainder of the difference between the indices may be a result of the different baseline removal procedures in the indices.
We have found, using a simplified model, that we can recreate the observed impact of higher spatial resolution in the SME index. We initially demonstrated that the undersampling is a result of the average station spacing in the AE index being less than the average spatial scale of the magnetic disturbances below the electrojet currents. We found that once station spacing in our model is reduced sufficiently, the extreme values in the modeled current systems are effectively resolved. By including time dependence in our model, we successfully qualitatively recreated other behavior in the electrojet indices such as the long-tailed shape of the distribution of the indices and the kink in the distribution seen at the transition from moderate to extreme values.

10.1029/2020JA027828
Our results provide insight into how the considerable body of work obtained using D ST can be directly related to new information arising from SMR. We found that the relationship between AE and SME is more complicated. For work regarding specific space weather events, SME has the advantage as it resolves the large and extreme events more effectively. However, for long-term studies regarding the behavior of the indices over several solar cycles, care should be taken as the increase of station number with epoch introduces subtle effects to the data set, which can be of the order of the space climatological effects being studied.