Distributions of Birkeland Current Density Observed by AMPERE are Heavy‐Tailed or Long‐Tailed

We analyze probability distributions of Birkeland current densities measured by the Active Magnetosphere and Planetary Electrodynamics Response Experiment. We find that the distributions are leptokurtic rather than normal and they are sometimes heavy‐tailed. We fit q‐exponential functions to the distributions and use these to estimate where the largest currents are likely to occur. The shape and scale parameters of the fitted q‐exponential distribution vary with location: The scale parameter maximizes for current densities with the same polarity and in the same location as the average Region 1 current, whereas the shape parameter maximizes for current densities with the same polarity and in the same location as the average Region 2 current. We find that current densities |J| ≥ 0.2 μA m−2 are most likely to occur in the average Region 1 current region, and second most likely to occur in the average Region 2 current region. However, for extreme currents (|J| ≥ 4.0 μA m−2), we find that the most likely location is colocated with the average Region 2 current region on the dayside, at a colatitude of 18°−22°.

found that Northern Hemisphere currents were higher than those in the Southern Hemisphere during the Northern summer as well as spring and autumn; later studies found the same effect using Swarm (Workayehu et al., 2020) and DMSP (Xiong et al., 2020).
It was outlined in Cowley (2000) that the region 1 (R1) Birkeland current should be colocated with the open/ closed field line boundary (OCB). This has been confirmed experimentally, showing that Birkeland current morphology is intimately related to the location of the polar cap Clausen et al., 2012;Coxon et al., 2014b). Studies have found that accounting for the resulting latitudinal variability can be very important (Coxon et al., 2017;S. E. Milan et al., 2015) and recently, a method has been made available to the community in order to allow current densities to be plotted in a coordinate system defined by the OCB (Burrell et al., 2020;S. Milan, 2019); this coordinate system shifts the OCB to 20° colatitude such that all current densities poleward of that are R1 current densities and equatorward of that are R2 current densities.
The probabilities of extreme GICs have been considered by measuring the rate of change R of the magnetic field that induces GICs at the Earth's surface and hence estimating where the largest GICs will occur (Rogers et al., 2020;Thomson et al., 2011) and what causes them Smith et al., 2019). Substorms are phenomena that cause large currents to flow into and out of the ionosphere in the substorm current wedge (e.g., Coxon et al., 2014b;. However, Freeman et al. (2019) showed that only ∼14% of the largest values of R correspond to magnetic perturbations observed during the substorm current wedge (termed DP-1), indicating that most of the highest R values are not directly attributable to substorms. Smith et al. (2019) performed a similar study with sudden commencements (SCs) of geomagnetic storms, finding that over 90% of the R values 3σ above the mean were observed within 72 hr of a SC, thus indicating that geomagnetic storm arrival predictions are key to predicting the impacts of space weather on Earth's surface. One method of examining the behavior of GICs is to model how they are distributed using a q-exponential function (Barbosa et al., 2017), which has also been investigated in some detail in ionospheric vorticity using SuperDARN data Chisham & Freeman, 2010. Chisham et al. (2009) showed that ionospheric vorticity measured with the SuperDARN radars  is closely related to Birkeland current behavior, as expected theoretically Sofko et al., 1995). Chisham et al. (2009) utilized SuperDARN measurements from 2000 to 2005 inclusive in order to obtain enough vorticities to evaluate the distributions of the ionospheric vorticities. Chisham and Freeman (2010) showed that radar-observed ionospheric vorticities demonstrated heavy-tailed distributions, finding that the best fit to the distributions was a Weibull or a q-exponential (Tsallis) distribution rather than a Gaussian. These distributions have higher kurtosis than a Gaussian distribution (meaning they are leptokurtic) and often have higher kurtosis than an exponential (meaning they are heavy-tailed); this implies that extreme ionospheric vorticities are more common than would be expected if the data were distributed normally or exponentially. Although Chisham et al. (2009) demonstrated that ionospheric vorticity could be used as a proxy for Birkeland current, the work done to quantify the vorticity distributions found using SuperDARN data (Chisham & Freeman, 2010 has not been done using Birkeland current data directly, because until now a suitable data set was not available for this type of analysis. We address this by using the Active Magnetosphere and Planetary Electrodynamics Response Experiment (AMPERE, Anderson et al., 2000Anderson et al., , 2021Waters et al., 2001Waters et al., , 2020. Furthermore, SuperDARN studies of ionospheric vorticity have been constrained by the geographical locations of the radars, whereas AMPERE is truly global. Recently, Chisham and Freeman (2021) showed how the parameters of a q-exponential fit to the ionospheric vorticity distribution vary spatially across the high-latitude ionosphere, and used these parameters to calculate the probabilities of seeing extreme vorticities in different parts of the polar cap and auroral region. The q-exponential is based on two parameters: q and κ. (The q-exponential is described in detail in Section 3, below.) q here is a free parameter which is known to be associated with intermittency; extreme values of a variable are much more likely than would be expected from a Gaussian distribution if that variable has an intermittent distribution. Higher values of q indicate more intermittency. Intermittency is one property of turbulence, and examinations of ionospheric vorticity in the paradigm of turbulence have shown that intermittency is similar between the ionosphere and solar wind (Abel et al., 2007) and that this only holds true when the two regimes are coupled by magnetic reconnection during Southward IMF . Freeman (2010, 2021) argued that the fact that ionospheric vorticity was well-modeled by non-Gaussian distributions such as the q-exponential provided evidence that there is a spatial variation in intermittency and therefore a spatial variation in the likelihood of

Birkeland Current Densities From AMPERE
AMPERE is a data set that combines measurements of magnetic perturbations from Iridium satellites with spherical harmonic fitting to estimate the spatial variation of Birkeland current densities. The data processing method was first outlined by Anderson et al. (2000); Waters et al. (2001) and a more detailed consideration is presented in Anderson et al. (2021); Waters et al. (2020). This current density j is presented on a grid of 1° latitude by 1 hr of magnetic local time (MLT) in altitude-adjusted corrected geomagnetic coordinates (AACGM), with a convention that upward currents are j > 0 and downward currents are j < 0. j is estimated from measurements in a ten-minute long sliding window, which is evaluated every 2 min. AMPERE has been used widely to study Birkeland current densities (Coxon et al., 2018). In this study, we use AMPERE data from 2010 to 2012; this is because there are some temporal gaps in AMPERE data in subsequent years which introduce uneven sampling of seasons.
Figure 1 (left) shows the spatial variation of the resulting mean current density ̄ at each latitude-MLT measurement location, which looks very similar to previous images of the mean Birkeland current densities. The morphology of these currents is as follows (Iijima & Potemra, 1976a, 1976b, 1978: The means in both hemispheres comprise two rings of current, one poleward (Region 1, or R1) and one equatorward (Region 2, or R2). The R1 currents are downward on the dawn side and upward on the dusk side, and the R2 currents are oppositely directed on the dawn and dusk sides. The regions are statistically latitudinally broad, spanning up to 10° of latitude from their poleward to equatorward extent. The mean j in the Northern Hemisphere is stronger than in the Southern Hemisphere (Coxon et al., 2016;Workayehu et al., 2020;Xiong et al., 2020). The center of the two rings of current is displaced toward the nightside. There are also two regions of mean current within the polar cap, which are some combination of Region 0 and NBz currents. The mean currents are in excess of 0.1 μA m −2 , consistent with the quiet time averages reported by Christiansen et al. (2002); Xiong et al. (2020) but lower than the averages reported by Juusola et al. (2014); Workayehu et al. (2019); this may be because 2010-2012 was in general less active than the periods reported by those authors, or may be an indication that AMPERE is underestimating the current densities.
Figure 1 (right) shows the mode current density j m . We again see rings of current which are positive/negative on different sides of the polar cap, but these are poleward of the ̄ of the same polarity. j m is smaller than ̄ , and the regions of current are latitudinally narrower, extending ∼5° in latitude. As with ̄ , the rings are displaced toward the nightside. In cases where j m and ̄ are of different signs, this tells us about the current distributions-for example, at 6 hr of MLT and 18° colatitude, where j m > 0 but ̄0 , this implies that the most common current density is positive, but the negative current densities are typically larger.
Previous work has shown that the currents reported by AMPERE in the polar cap are subject to large uncertainties (D. J. Knipp et al., 2014). To examine this, we present Figure 2 which shows the standard error on the mean for the same spatial regime as Figure 1. The standard error is higher in the Northern Hemisphere than in the Southern Hemisphere, which we infer is a signature of the generally higher Northern current densities. The standard error is highest colocated with the boundary between the R1 and R2 current regions and the NBz region shown in Figure 1, which is consistent with the work of D. J. Knipp et al. (2014). However, since we have approximately 7 × 10 5 data per bin, the standard errors are of the order 10 −7 , which is substantially lower than the means and modes reported in Figure 1; as such, we conclude that the error on the data presented herein is small.
The q-exponential is positive definite with a maximum at the origin. This means that it can be fit to any data distribution in which all the data are of the same sign (Chisham & Freeman, 2010). As AMPERE data can be positive or negative in any given bin, we take the mode in each bin (plotted in Figure 1, right) and subtract it from the distribution of current densities in each bin prior to performing the fitting. This shifts the distribution such that the mode is at zero, and means that we can fit two q-exponentials with a maximum at the origin (zero) and current densities of the same sign on each side of the origin. In this way we fit the q-exponential separately to the distribution of the modulus of the measured "positive" current density (j > j m ) and to the modulus of the "negative" current density (j < j m ). We refer to j > j m as positive current density and j < j m as negative current density hereafter to aid readability. We plot the number of AMPERE estimates of j in each bin across the polar cap for the The regions in which the number of data is lower or higher than the average correspond well to the ̄ plotted in Figure 1 (left). We note that this does not show the number of satellites contributing to the fits, but rather the number of fits in each bin, such that the number of fits is equal at all locations for total positive plus negative measurements. We also note that there are over 0.2 million data per bin, which is typically much larger than the number of data when using SuperDARN (Chisham & Freeman, 2021).

Distributions of Birkeland Current Densities
In this study, we fit q-exponential distributions to the Birkeland current densities from AMPERE, following the logic of previous work on the underlying ionospheric vorticities Chisham & Freeman, 2010. The q-exponential probability density function (PDF) is given by where f q,κ (j) is the probability density, κ is the scale parameter, q is the shape parameter, and j was defined earlier as current density. Values of q < 1 have less kurtosis than an exponential fit, indicating that the distribution decays more quickly than an exponential would. Values of q > 1 have more kurtosis, and so decay more slowly; they are heavy-tailed. The κ parameter governs the scale of the distribution; the decay rate is constant with respect to variations in κ but higher κ implies that the vorticity varies more from the mode magnitude of the current density and low κ means the variability is lower. Therefore, both higher q and higher κ imply that more extreme current densities will occur, but the underlying reason for this are higher intermittency and higher variability respectively. Figure 4 shows the way in which the shape of a q-exponential distribution changes with the two key parameters q (left) and κ (middle, right). In these schematics the value of the other parameter is kept constant (κ = 0.1 and q = 1.1 respectively). Where q = 1, the q-exponential tends to the exponential distribution; where q < 1 the distribution is less heavy-tailed than the exponential (it decays more quickly), whereas at larger q the distribution becomes heavier-tailed (it decays more slowly). Increases in κ make the distribution longer-tailed or wider, which does not affect the decay and therefore does not affect the shape of the distribution. This is shown by Figure 4 (right) which shows the value of j at which the value of f(j) has decreased to a value of f(j) = f(j = 0) × q −35 (chosen to be visible on the schematic); it can be seen that as kappa increases monotonically so does j, showing that the distribution is being scaled upward (and thus becoming longer-tailed) with increasing κ. Figure 5 shows how current densities are distributed in two of the coordinates underlying Figure 1. We use maximum likelihood estimation (MLE) to determine the best fit PDF to the distributions of current density (Chisham & Freeman, 2010. Figure 5 suggests that the q-exponential fit is better than the exponential or Gaussian fits to the distribution for both positive and negative current densities; this is confirmed by comparing the Akaike Information Criteria (AIC) of the three fits.
Figure 5 (top) shows the distribution for a colatitude of 15° and an MLT of 18 hr. Both the mode and mean are positive, and the mean is larger than the mode. This coordinate is chosen because it is in the region of average R1 current on the dusk side of the polar cap, and demonstrates the way in which the asymmetry and the tails of the underlying distributions dictate whether ̄ or ̄ in a given cell. The tails of the R1 current distribution decay more slowly than a Gaussian (they are leptokurtic) but with approximately the same decay as an exponential (they are not heavy-tailed). The top of the distribution is wider than an exponential, especially for positive current densities. These inferences are consistent with the values of the fit, which shows that κ is high (implying a wide or long-tailed distribution) and that q < 1 (implying that the distribution decays more quickly than an exponential).  density on average (̄) , but we see that there are some large negative current densities within the distribution. There are fewer negative current densities than positive for any given magnitude, and the logarithmic scale makes it obvious that the vast majority of the negative current densities in this coordinate are in the range −0.4 < j < 0 μA m −2 : that is to say the R2 distribution is narrower than the distribution for R1, but clearly leptokurtic and heavy-tailed compared to the exponential. These inferences are also consistent with the q-exponential fit, which shows that q > 1 (implying that the distribution decays more slowly than an exponential) but shows that κ is small (implying a narrow distribution). Figures 6 and 7 show the values of q and κ respectively for the distributions in each bin for the period of interest. Light gray bins correspond to where q = 1. Because neither q nor κ are negative, we multiply them by the sign of the mean of the currents ̄ when plotting them to make it easier to compare them to the large-scale current systems shown in Figure 1 (left). Figure 6 shows the values of q. The locations of q > 1.2 are mirrored for positive and negative current densities, which are the mirror image of one another across the noon-midnight meridian, but which are in almost identical locations in the Northern and Southern Hemispheres; as such, we focus on the Northern Hemisphere henceforth.
Since high values of q are colocated with negative ̄ on the dusk side and positive ̄ on the dawn side, we infer that high values of q are predominantly associated with the average R2 current system. For positive current densities (left), the region of high q on the dawn side is further from the pole than the region on the dusk side; the dusk side region encroaches onto the equatorward edge of the R1 current system. The distribution at a given location will be a mixture of the distributions of the processes which occur in that location (Chisham & Freeman, 2021), and so we interpret this as an effect of the mixing of the R1 and R2 distributions due to the expansion and contraction of the current ovals in response to flux transfer into and out of the polar cap (Cowley & Lockwood, 1992;Lockwood et al., 1990;Lockwood & Freeman, 1989;Siscoe & Huang, 1985). We see the mirror image of this for negative current densities (right). This implies that the R2 currents are characterized by high intermittency and therefore the distributions are heavy-tailed, with a decay rate lower than an exponential distribution. The R1 currents do not display high intermittency and instead the distributions are close to exponential.
At higher latitudes in the Northern Hemisphere, for positive current densities (left) there are regions of q ∼ 1.1 on the dayside poleward of 10°, with negative ̄ on the dusk side and positive ̄ on the dawn side. These appear to be associated with NBz currents which we associate with northward IMF. For negative current densities (right) there is a region of q ∼ 1.2 poleward of 10° pre-noon and equatorward of 10° post-noon. The same is seen with slightly lower q in the Southern Hemisphere. We associate these four high q regions with the NBz current system caused by northward IMF B Z ; these current systems do not have intermittency as high as the R2 current systems but are heavy-tailed. Figure 7 shows the values of κ. The locations of κ > 0.1 are mirrored for positive and negative current densities and are very similar between the two hemispheres; as such, we focus again on the Northern Hemisphere. For positive current densities (left), κ is highest (κ > 0.2) for the positive ̄ on the dusk side, corresponding to the R1 current system, but encroaches very slightly into the R2 current system. For negative current densities (right) we see similar morphology of the opposite polarity. There is also a region of κ ∼ 0.1 corresponding to the poleward Figure 6. Maps showing the value of q in the best q-exponential fit to the underlying distribution per coordinate, in the same format as Figure 3. White bins indicate bins in which the fitting routine to identify the best parameters for the q-exponential was unsuccessful (note that the color bar is gray at 0, rather than white). The values are multiplied by the sign of the mean current ̄ , such that positive values correspond to positive ̄ and vice versa. Note that these maps are not maps of current density.
edge of positive currents on the dawn side, corresponding to the edge of the R2 current system. Again, we interpret this encroachment into R2 as evidence of a mixture distribution due to the motion of the current ovals.
The high values of κ imply that although the R1 current distribution is not heavy-tailed, it is wide or long-tailed. κ > 0.2 on either side of noon at high latitudes (<10°) corresponding to the NBz current system, from which we infer that NBz currents have a long-tailed or wide distribution. Signatures associated with NBz are more difficult to discern in the Southern Hemisphere (bottom).
As such, we have concluded that R1 current is characterized by long-tailed or wide (high-κ) but not heavy-tailed (q ∼ 1) distributions and R2 current is characterized by narrower distributions (low κ) which are heavier-tailed (high q). NBz currents sit in the middle in both respects; heavier-tailed than R1 but less than R2, and longertailed than R2 but less than R1. This implies that current densities closer to the center of the distribution are more likely in R1 than in R2, but that very extreme currents are more likely to be seen in R2 current. It further implies that moderate current densities are less likely in NBz but extreme current densities are also less likely, which is consistent with Figure 1 showing that the average NBz currents are smaller than R1 and R2. We examine this in more detail in the next section. Maps showing the value of κ in the best q-exponential fit to the underlying distribution per coordinate, in the same format as Figure 6. The values are multiplied by the sign of the mean current ̄ , such that positive values correspond to positive ̄ and vice versa. Note that these maps are not maps of current density.

Probabilities of Birkeland Current Densities
Having calculated f(j) per coordinate in Section 3, we can now calculate the probability of extreme values of current density by using the survival function (Chisham & Freeman, 2021). For j > j m the survival function is given by and for negative currents (j < j m ), the survival function is given by where P is the probability of current flowing and n is the number of measurements. P + , n + , q + , κ + are those quantities for positive current densities and P − , n − , q − , κ − are those quantities for negative current densities.
We assess the probabilities of seeing current densities above 0.2, 1.0, 2.0, and 4.0 μA m −2 in the following figures, and discuss the probabilities we find and our interpretations of these probabilities in the context of the scientific literature in Section 5. Figure 8 shows the probability (| | ≥ 0.2 A m −2 ), which is sometimes used as the bottom of the color scale when plotting current densities with AMPERE to highlight large-scale structure (e.g., Coxon et al., 2014a). Qualitatively, the probabilities presented here show that the highest probabilities are co-located with the NBz, R1, and R2 current regions seen in Figure 1 (left) and are also co-located with the regions of high q or κ described in Section 3, which is as we would expect. Quantitatively, (| | ≥ 0.2 A m −2 ) is almost 30% for Region 1 current in both hemispheres, and is closer to 15% for Region 2 and NBz current. Thinking back to the discussion of q and κ in Section 3, the fact that these current densities are more probable in the R1 current region is consistent with the fact that R1 is characterized by wider or longer-tailed (higher-κ) distributions than R2 current; moderate current densities are more probable in R1 as a result.
Notably, the chance of seeing current (| | ≥ 0.2 A m −2 ) is higher in the Northern Hemisphere than in the Southern Hemisphere, especially toward the dayside. In the Northern Hemisphere the probabilities for R1 current decrease on the nightside compared to the dayside; this does not occur in the Southern Hemisphere, for which the probabilities for R1 current stay close to 30% at all dusk local times for positive j and all dawn local times for negative j. Figure 9 shows (| | ≥ 1.0 A m −2 ). The difference between this and Figure 8 is striking. We see two main zones of likelihood P > 0; for positive current densities these zones are located between noon and 22 MLT (on the dusk side) and between midnight and 6 MLT (on the dawn side). We refer to these zones as Zone A and Zone B respectively hereafter, and these zones are highlighted in Figure 12 for the convenience of the reader. For negative current densities Zones A and B are mirrored, such that Zone A is on the dawn side (2 MLT to noon) and Zone B is on the dusk side (15-21 MLT). The probability (| | ≥ 1.0 A m −2 ) in Zone A is ∼1% (∼30× lower than in the previous figure) and in Zone B is ∼0.1%. The hemispheric difference here remains the same as in Figure 8, in that Northern Hemisphere likelihoods are higher than those in the Southern Hemisphere. Figure 10 shows (| | ≥ 2.0 A m −2 ), which looks qualitatively similar to Figure 9. The key difference is that the probability in Zone A is now ∼0.1% and in Zone B is P ∼ 0.01%, which is an order of magnitude smaller than in the previous figure. The Northern Hemisphere probabilities are 2× higher than the Southern Hemisphere. The probabilities in both regions again decrease by an order of magnitude in Figure 11, which shows (| | ≥ 4.0 A m −2 ) and which looks very similar to Figure 10.

Discussion
The mean currents presented in Figure 1 (left) imply that R1 currents are typically larger than R2 currents, consistent with previous averages (Anderson et al., 2008;Weimer, 2001) but potentially attributable to spatial smearing at more equatorward colatitudes (Coxon et al., 2017). Studies have found that R1 currents are typically stronger than R2 when the amount of current flowing is high, but have found that R2 currents dominate when the amount of current is lower (Coxon et al., 2014a); Figure 8 shows that the probability of seeing current densities (| | ≥ 0.2 A m −2 ) is higher for R1 (10°-15° colatitude) than it is for R2 (18°-25°), which is consistent with this. Chisham et al. (2009) reported that the large-scale morphology of ionospheric vorticity is similar to the expected large-scale morphology of Regions 1 and 2 Birkeland current. Plotting the mean ionospheric vorticity on a grid identical to the ones in this manuscript (but with an equatorward cutoff at 24° colatitude) showed two large regions between 10°−20° colatitude on either side of the noon-midnight line corresponding to R1 current and two regions of the opposite sign corresponding to R2 current between 20°−24° colatitude (the R2 vorticity regions extend to the equatorward edge of the grid and so may extend further equatorward than that range implies). Additionally, Chisham et al. (2009) observed regions of opposite polarity to the poleward edge of the R1 regions; they subdivided by IMF clock angle to differentiate these into NBz and R0 current systems. Figure 1 shows that the mean Birkeland current measured by AMPERE is broadly similar to that reported by Chisham et al. (2009), with the R1 and R2 currents approximately colocated with their corresponding vorticity regions. The one area in which this was not the case was on the dayside, where our R1 and R2 currents are slightly closer to the pole than the vorticity regions; however, this could be due to the different time domains of Figure 8. Maps of (| | ≥ 0.2 A m −2 ), in the same format as Figure 6. Bins in which the probability could not be computed were set to zero. Note that these maps are maps of probability and not maps of current density. the two datasets. We see mean current densities in the polar cap consistent with the NBz vorticities observed by SuperDARN. Chisham and Freeman (2021) presented the values of the mode for the SuperDARN vorticities and reproduced the means shown in Chisham et al. (2009). The mode values reported were approximately half the mean values, which is consistent with the values reported from AMPERE in Figure 1. However, the morphology of the mode vorticities is quite different from AMPERE, with the mode values associated with R1 current smaller than those associated with R2 current, as opposed to the values reported by AMPERE in which the modes in the R1 current are larger than those in R2. In the AMPERE data we also see mode currents which are of the same polarity as the mean NBz currents, whereas there is no signature of the NBz vortices in the SuperDARN modes reported by Chisham and Freeman (2021).

Correspondence to Ionospheric Vorticity
In terms of data coverage, Chisham et al. (2009) reported that there were between 100 and 100,000 data per bin in the colatitude range 4°−24° and zero data per bin outside that range . The data coverage was uneven due to the uneven nature of SuperDARN observations: there is a region in which most of the data is found which is a circle offset from the pole toward the nightside, such that most of the data is found in the colatitude range 17°-23° at midnight; 15°-20° at dawn and dusk; and at 10°-15° at noon. The region with the most data is this region on the nightside; the region with the least is the equatorward edge on the dayside and the polar edge at all magnetic local times. Contrasting this to Figure 3 herein, we can see that there are ∼700, 000 data per bin (the sum of the two panels), and the number of positive or negative current densities per bin varies between ∼200, 000 and ∼550, 000, which is between twice the number and 55,000 times the number of data per Figure 9. Maps of (| | ≥ 1.0 A m −2 ), in the same format as Figure 8. Note that these maps are maps of probability and not maps of current density. bin. Additionally, the figures in this paper extend to 40° instead of cutting off at 24°, and in this paper we are able to examine the Birkeland currents in the Southern Hemisphere whereas this has not yet been possible with SuperDARN Chisham & Freeman, 2010.
The vorticities observed using SuperDARN data have been shown to be generally well-described by fitting a q-exponential distribution to each side of the distribution from the mode (Chisham & Freeman, 2010). The best values for the distribution in each bin were examined by Chisham and Freeman (2021), corresponding to the top rows of our Figures 6 and 7. Interestingly, their results differ somewhat from the ones presented in our figures. They found, like we have, that the Region 2 current system is characterized by high values of q but low values of κ, and found that dayside R1 currents are close to exponential (q ∼ 1 and high κ), similar to our findings. However, they found that the characteristic distribution for R1 was different on the nightside, which showed higher values of both parameters; we do not observe this. They also found that the polar cap was characterized by intermediate values of q (which we observe herein) but by low values of κ (which we do not observe). There will be some differences between the analyses because the current is measured in space in the inertial frame, whereas the vorticity is measured from the ground in the co-rotation frame. However, the differences between the current and vorticity distributions must mainly be related to spatial changes in conductance, and this is an area for future research.
In terms of the latitude ranges, our finding is similar to Chisham and Freeman (2021), who found that large vorticities were most likely between colatitudes of 15°−20° but that extreme vorticities were most likely at colatitudes equatorward of 20°. They found that the largest probabilities were between 10 and 18 MLT, which is similar to Figure 10. Maps of (| | ≥ 2.0 A m −2 ), in the same format as Figure 8. Note that these maps are maps of probability and not maps of current density. our result in that it is biased toward dusk, but they did not see a bimodal probability distribution in MLT sector as we do. Figures 8-11 show that the largest probabilities of observing a current above a given threshold occur at 05-09 MLT (for negative current density) and 14-20 MLT (for positive current density) and that the peak probability moves equatorward with increasing threshold. This movement is shown more clearly in Figures 13 and 14. Figure 13 shows cuts of the probability through the meridian of MLT = 15 and 3h for positive current densities and through MLT = 21 and 9 hr for negative current densities. Zone A corresponds to the 15 and 9 hr meridians, whereas Zone B corresponds to the 3 and 21 hr meridians. Figure 14 then plots how the colatitude of peak probability varies with current threshold for zones A and B. In both figures, we can see that the colatitude of peak probability moves equatorward as the threshold current density is increased. Figure 12 also shows that the peak current is most likely to be seen in the Northern Hemisphere, as has been previously found for total current (Coxon et al., 2016) and Poynting flux (D. Knipp et al., 2021;Pakhotin et al., 2021). At the lowest thresholds (Figures 13a and 13b) the Northern Hemisphere currents are ∼50% larger than those in the Southern Hemisphere, and Northern Hemisphere dominance increases with current threshold (∼7× larger at the highest thresholds, shown in Figures 13g and 13h).

Interpreting the Probabilities in Terms of R1 and R2 Current
Zone A occurs in the same sense as Region 1 current; that is to say, for positive current densities it is on the dusk side, and for negative current densities it is on the dawn side. R1 current is positive on the dusk side and negative on the dawn side. It is tempting to interpret this as evidence for this Zone comprising R1 current. However, as we have discussed previously in this study, the large-scale sense of R1 and R2 current does not preclude large current densities on the "wrong" side of the distribution. For example, in Figure 5, the R1 current (top) and R2 current (bottom) are taken from bins which are plotted as positive current density in Figure 1 because the mean of both distributions is positive. However, there is clearly a lot of negative current densities contributing to both distributions, and some of them are large (especially in Figure 5, bottom). Therefore, it is not implausible that there are (at times) very large negative current densities in a spatial region which is characterized on average by positive current.
As such, we can explain Zone A in one of two paradigms, both of which have interesting repercussions for the behavior of the system.

Paradigm 1
At any given location, we will measure either R1 or R2 and so the negative and positive sides of the distribution are driven by R1 and R2 dominating at a location at different times. It has been shown that the current ovals expand and contract with the polar cap Clausen et al., 2012; as part of the expanding/contracting polar cap (ECPC) paradigm (Cowley & Lockwood, 1992; Lockwood et al., 1990;Siscoe & Huang, 1985), and thus the current ovals are at their most expanded when the amount of magnetic flux (and therefore energy) in the system is at its highest.
It has further been shown that more energy is deposited by substorms when more magnetic flux is present in the polar cap (Coxon et al., 2014b; and that more current flows during and immediately after magnetic reconnection Coxon et al., 2014aCoxon et al., , 2017. As such, if the first paradigm is the explanation behind our observations, we can interpret this in terms of more energy and more current flow during periods of magnetic reconnection. In this paradigm, the peak probability of extreme current results from R1 current which has moved equatorward as magnetic flux is opened by dayside reconnection, which is consistent with the sense of the currents. However, this model implies that the intermittency of R1 current increases (i.e., q increases) depending on the size of the polar cap. We would not expect the underlying distribution of each current region to change in this way, and this result would be unintuitive.

Paradigm 2
In any given current system, the large-scale current system will comprise multiscale currents which have fine structure. There can be both intense upward and downward currents, but these average to look like the R1 and R2 current systems. In this model, R1 distributions are characteristically wide or long-tailed (high κ) and R2 distributions are characterized by high intermittency and therefore heavy-tailed (high q); the probabilities at high thresholds are related to distributions with high intermittency.
We find that the strongest currents are seen in the R2 system, and we find that the strongest currents are seen at dawn and dusk, biased toward the dayside. Coxon et al. (2014b) showed that the ratio of R1 to R2 current is larger during substorms than during other times, and Coxon et al. (2017) showed that although R2 contributes to the substorm current wedge it is to a smaller extent than R1. Additionally, if these signatures are due to the substorm current wedge they would occur on either side of midnight; therefore, we infer that these very strong currents are not seen during substorms.
It has long been thought that R2 currents close through the partial ring current on the nightside of Earth (Iijima & Potemra, 1976a, 1978, and it is well-known that the ring current is enhanced during geomagnetic storms (Chapman & Ferraro, 1933). As such, if paradigm 2 is underpinning our observations, we infer that highest probabilities  Figure 13. The left-hand column shows the Northern Hemisphere and the right-hand column shows the Southern Hemisphere. Panels (a and b) show the location of peak probability of the Zone A plotted in Figure 13 and panels (c and d) show the Zone B meridians. The values for the positive current densities are plotted in red and have 0.1° added; for negative current densities, blue is used and 0.1° is subtracted. The additions and subtractions are to make it easier to observe the trend for the reader. occur when R2 is closing through the enhanced ring current during geomagnetic storms and that these closures happen around the dawn and dusk sectors.
However, in this paradigm it is unintuitive that the largest R2 currents occur in the opposite sense to the average R2 current system. We would expect the strongest currents to occur on the side of the distribution that dominates the mean current of a given region, although it may be that our observations are consistent with the embedded Birkeland currents recently reported by Liu et al. (2021). In a future study, we will repeat this analysis based on times at which we see ring current intensifications in order to verify this inference.

Why Not Both?
It is possible that both paradigms 1 and 2 are true, and that our observations are the result of an interplay between the two. We also need to note that AMPERE does not have fine-scale spatial resolution, and so it may be that both are true but AMPERE observations bias the observer toward interpreting the data in paradigm 1. Figure 15 shows the data shifted into a coordinate system defined by the location of the open-closed field line boundary (OCB) (Burrell et al., 2020;S. Milan, 2019;S. E. Milan et al., 2015) in order to try to address this ambiguity. In this plot the OCB is located at 20°, such that R1 currents should be poleward and R2 currents should be equatorward. For positive current densities, the region of highest probability is equatorward of the circle, which may suggest that the largest currents observed here are R2 currents (paradigm 2). However, for negative current densities, the region of highest probability is poleward of the circle, which may suggest that the largest currents  (Burrell et al., 2020). Note that these maps are maps of probability and not maps of current density. observed here are R1 currents (paradigm 1). This confirms the need for further work which examines the distributions of the currents during different types of driving.

Conclusion
We take maps of AMPERE data and analyze them statistically. We find that the mean current densities resemble the Regions 1 and 2 current systems which have been investigated throughout the history of the field (Birkeland, 1908;Iijima & Potemra, 1978;Liu et al., 2021;Weimer, 2001). However, we find that the modes of the current density distributions are displaced toward the pole relative to the means, and that the mean current densities traditionally interpreted as R1 and R2 are explained mainly by asymmetries in the distribution of the currents rather than by the mode current.
We analyze the current density distributions and find that they are well-described by a q-exponential distribution rather than a Gaussian distribution. We infer that the current densities exhibit complex behavior due to their heavy-tailed distributions, and we fit q-exponential distributions to describe the behavior of the distribution.
From the fits we derive the probability of current densities above certain thresholds and find that R1 currents are relatively more likely to be of clearly measurable magnitude than R2 currents, and that currents are most likely to occur at colatitudes of 10°-15°. However, the most extreme currents are seen at 18°-22° and are potentially linked to one of two paradigms: 1. The most extreme currents occur in the R1 currents at the point when the current ovals are most expanded. 2. The most extreme currents occur in the R2 currents due to closure through an intensified ring current during geomagnetic storms.
Adjusting the coordinate system to one defined by the OCB implies that these paradigms both contribute to the observations, and further work is needed to understand the contributions of each.