As mentioned before, extending B N measurements to regions without magnetometers can be a difficult task due to the influence of localized phenomena throughout the magnetosphere and the lithosphere. This study seeks to overcome such difficulties and precisely reconstruct geomagnetic field variations by interpolating with the Spherical Elementary Current Systems (SECS) technique (Amm and Viljanen, 1999). SECS is commonly used to estimate geoelectric fields for the purpose of calculating GICs over regional-scale infrastructure (e.g., Wik et al. (2008), Kelbert and Lucas (2020)). A summary of the technique is presented in this section; readers looking for a more in-depth explanation are encouraged to refer to Vanhamäki and Juusola (2020).

SECS interpolation operates by representing the ionospheric current distribution with a parametric model. Empirical determination of the model’s parameters at each moment in time is carried out based on magnetometer observations, which allows for the assumption that the model current distribution causes all observed ground-level geomagnetic activity. Once the model is fitted using this approach, one can apply the Biot-Savart law to the model current distribution in order to calculate the magnetic field deflection at any point within a region of interest.

The model current distribution is a superposition of curl-free and divergence-free vector functions defined on a common spherical surface, each with a different location for its pole. These basis functions are called elementary current systems. One is free to choose the placement of the poles of elementary current systems, but typically they are arranged in a latitude-longitude grid for convenience. By Helmholtz’s Theorem, any smooth, differentiable vector function can be exactly decomposed into such basis functions; using a finite number of basis functions at worst only approximates the original vector field. The divergence-free part ( J ⃗ df ) and curl-free part ( J ⃗ cf ) of each elementary current system are defined as (Equation 1): J ⃗ df ( ϕ ′ , θ ′ ) = I 0 , df 4 π R I cot θ ′ 2 ϕ ̂ ′   , J ⃗ cf ( ϕ ′ , θ ′ ) = I 0 , cf 4 π R I cot θ ′ 2 θ ̂ ′ (1)

One would typically consider the coordinate system in which r is the radial distance from the center of Earth, and ϕ and θ are geographic longitude and colatitude, respectively. The “primed” spherical coordinate system ( r ′ , ϕ ′ , θ ′ ) differs from the standard “unprimed” ( r , ϕ , θ ) system in that its north pole, where θ ′ = 0 , is defined to be at the location of that current system. The azimuthal angle ϕ ′ is defined accordingly. As opposed to the geographic unprimed coordinate system, we will hereafter call a primed coordinate system an egocentric coordinate system, since the pole of each current system has its own unique primed coordinate system centered on itself.

Here the currents are assumed to lie in an infinitely thin current shell of radius R I (which typically takes the value R E + 110 km, where R E is the radius of Earth, so as to approximately correspond with the heights of the maximum Pedersen and Hall conductivities). The factors I 0 , df and I 0 , cf are empirically determined at every moment in time in order to fit this model ionosphere to magnetometer data. The total vector field representing ionospheric currents may not actually equal the currents in the real ionosphere, especially if the associated field-aligned currents are not approximately radial. However, the elementary currents still form a basis from which the magnetic field can be later reconstructed.

Under the assumptions of this current-shell model, it can be shown that the ground magnetic effect of the curl-free current systems vanishes within the current surface (Fukushima, 1976). This allows I 0 , df to be rewritten more simply as I 0 . Thus, the horizontal component of the magnetic field due to a single current system at a point r ′ = ( R E , ϕ ′ , θ ′ ) is solely due to J ⃗ df and can be shown to be (Equation 2): b ⃗ H ( θ ′ ) = − μ 0 I 0 4 π R E ⁡ sin θ ′ ⋅ R E R I − cos θ ′ 1 − 2 R E ⁡ cos θ ′ R I + R E R I 2 + cos θ ′ θ ̂ ′ (2)

The constant μ 0 in this equation denotes the magnetic permeability of free space. There is also a radial component to this magnetic field, but it will not be used in this analysis because we are only concerned with horizontal field components (specifically, this study eventually only examines the northward component). Then, the total horizontal vector magnetic field B ⃗ H ( θ ′ ) is the sum of the magnetic contributions b ⃗ H ( θ ′ ) from each current system.

To determine the contributions from each current system, one can define a vector Z ⃗ containing the ground-truth magnetometer measurements at stations 1 through n , a vector I ⃗ of scaling factors for current systems 1 through m , and a 2 n × m transfer matrix T ̄ (Equation 3): Z ⃗ = z 1 , θ z 1 , ϕ z 2 , θ z 2 , ϕ ⋮ z n , θ z n , ϕ , I ⃗ = I 0,1 I 0,2 ⋮ I 0 , m , T ̄ = T 1,1 , θ T 1,2 , θ … T 1 , m , θ T 1,1 , ϕ T 1,2 , ϕ … T 1 , m , ϕ T 2,1 , θ T 2,2 , θ … T 2 , m , θ T 2,1 , ϕ T 2,2 , ϕ … T 2 , m , ϕ ⋮ ⋮ ⋮ T n , 1 , θ T n , 2 , θ … T n , m , θ T n , 1 , ϕ T n , 2 , ϕ … T n , m , ϕ (3)

Each element T i , j , ( θ , ϕ ) represents the θ or ϕ component of the magnetic effect that current system j has on magnetometer i . One must then solve the matrix equation Z ⃗ = T ̄ I ⃗ for I ⃗ and plug the resulting scaling factors into Equation 2 to obtain interpolated magnetic field values. This is typically done using singular value decomposition and setting small singular values to 0 to properly condition the matrix. One usually considers “small singular values” to mean those that are less than the maximum singular value multiplied by some cutoff coefficient ϵ . The best value of ϵ depends on the specifics of the SECS model in question; refer to Section 3.2 for information on how this value was determined for the present case.

When using Equation 2 to find the contribution of each current system to the magnetic field, it is important to notice again the difference between the geographic and the egocentric coordinate systems. These contributions are most easily added together when recognizing that a contribution vector b ⃗ ego = ( b θ ′ , b ϕ ′ ) at geographic colatitude and longitude ( θ i , ϕ i ) due to a current system at geographic colatitude and longitude ( θ j , ϕ j ) may be transformed from the egocentric to the geographic coordinate system with Equations 4, 5: θ ego = arccos cos θ i ⁡ cos θ j + sin θ i ⁡ sin θ j c o s ( ϕ i − ϕ j ) , (4) b ⃗ H = b θ b ϕ = − b N b E = b θ ′ ⁡ cos cos θ j − cos θ i ⁡ cos θ ego sin θ i ⁡ sin θ ego b θ ′ ⁡ sin sin θ j ⁡ sin ( ϕ j − ϕ i ) sin θ ego (5)

The total instantaneous magnetic field in the geographic coordinate system is then the sum of each contribution vector b ⃗ H , which we calculate component-wise to find the northward component, B N , and the eastward component, B E . Since B N is this study’s variable of interest, we focus solely on B N without using B E .

We define our model current distribution as the superposition of a set of current systems whose poles are placed on a quadrangular latitude-longitude grid extending from 40.5 ° N to 60.5 ° N and from 70.5 ° W to 140.5 ° W (the green quadrangle in Figure 1) in geographic coordinates. The current system array is offset half a degree from whole-number placement to avoid colocating the current systems with points at which the magnetic field is evaluated, as this would cause a singularity to arise in Equation 2. There is also the danger of a singularity if a current system is colocated with a magnetometer; however, this becomes obvious with a cursory inspection of the interpolation outputs. If one places the current systems such that this is not a problem, one can rest assured that no singularities will arise later because the locations of current systems and magnetometers do not change over time.

Poles are placed every 2°, making this a 10-by-35 array of current systems. The performance of SECS interpolation is not especially sensitive to small changes to the number and spacing of the current systems, so the 2-degree separation was chosen based on its use in McLay and Beggan (2010). We further validated this choice by performing joint Bayesian optimization of the number of rows and columns of current systems (while keeping the boundaries of the array fixed) using the Optuna package for Python (Akiba et al., 2019). The objective of this optimization was to minimize the absolute error between the model outputs at the location of each magnetometer and the measurements of the magnetometers themselves, averaged across magnetometers in each timestep, and then averaged across all timesteps in the dataset. The errors at each magnetometer were calculated using a leave-one-out scheme in which the model was fit using data from all magnetometers except the one in question, meaning a new model was fit for each magnetometer, all identical except for which magnetometer was omitted. The optimization found that the error was minimized for an 8-by-39 array of current systems. This agrees with our similar choice of array dimensions, whose error is less than 1% higher than the optimum. As such, we kept the 10-by-35 array so as to maintain a rounded 2-degree separation, as well as for possible future comparison to studies such as McLay and Beggan (2010).

Together with the number of rows and columns of current systems, the above procedure also jointly optimized a third quantity: the cutoff coefficient ϵ . We use the best value determined through this optimization, ϵ = 0.0932 . This value is comparable in magnitude to what has been used in previous related studies (McLay and Beggan, 2010; Weygand et al., 2011).

When creating a SECS-based model, it is sometimes helpful to introduce an additional current sheet located at a negative altitude (that is, internal to the planet). Doing so may allow one to delineate which part of geomagnetic perturbations arise from ionosphere-associated phenomena and which perturbations arise from currents induced in the subsurface (Pulkkinen et al., 2003). However, the goal of the SECS model as it is used here is only to reproduce geomagnetic perturbations, not to explain them in terms of their contributing equivalent currents. As a result, the additional current sheet has not been employed in this study.

Measurements of d B N are interpolated to give estimates of this field component at a set of points within the region of interest. The points could have any spatial density since they are sampled from the output of the SECS method, which is a continuous function. We limit ourselves to making estimations on a regularly spaced 14-by-32 grid of points to keep computation time under control. Each of these “heatmaps” (e.g., Figure 2) is associated with a set of simultaneous magnetometer measurements, so exactly one heatmap is created per timestep in our dataset.

Once the heatmaps have been created, one can visually inspect them to view the estimated extent of LGMDs. However, we need to define a method to locate LGMD in a given heatmap. We do so here by granting a quantitative aspect to the definition of LGMD given in Section 1; that is, by defining a threshold value d B threshold . Any point in a heatmap whose absolute value exceeds that of the threshold (i.e., | d B N | > d B threshold ), is considered part of an LGMD. Points identified this way that are adjacent, including diagonally, are considered to be part of the same LGMD. We then draw a contour line around each LGMD in a heatmap using the marching squares contour-finding algorithm implemented in the measure.find_contours function of the scikit-image Python package (Van Der Walt et al., 2014).

To identify LGMDs in this study, we define our perturbation threshold as d B threshold = 25.95 nT, which corresponds to the absolute value of the median of the distribution of the minimum values from each heatmap. It is not a particularly strict threshold, but that is because we are interested in all LGMDs, not only anomalously intense ones. LGMDs that extend past the edge of the region of interest, and thus are not associated with a closed contour, cannot be characterized with confidence, and are thus completely omitted from our analysis. Examples of both closed and open contours drawn with our method can be seen in Figure 2.

All our results in Section 4 are based on the heatmaps created using the procedure outlined in this section. Some are, more specifically, based on the perimeters of LGMDs identified in the heatmaps, each of which is calculated by adding up the distances between successive vertices of the contour that surrounds it. One must respect the spheroidal shape of Earth when computing these distances. Given two contour vertices with (latitude, longitude) coordinates ( λ 1 , ϕ 1 ) and ( λ 2 , ϕ 2 ) in radians and the radius of Earth R E = 6,378 km, we find the haversine distance d between vertices (Equation 6) rather than the Euclidean distance as d = 2 R E ⁡ arcsin sin 2 λ 2 − λ 1 2 + cos λ 1 ⁡ cos λ 2 ⁡ sin 2 ϕ 2 − ϕ 1 2   . (6)

In this work, we focus on perimeter as a measure of LGMD size. While there is no fundamental problem with investigating the behavior of LGMDs with different measures – for example, diameter – concentrating on the perimeter of LGMDs removes a degree of arbitrariness from the study. In the example of studying LGMD diameter, one would be forced to confront the question of along which direction to calculate that diameter, especially considering that LGMDs can vary widely in shape (an effect which is displayed in Figure 4). This is not an insurmountable issue; we hope that future work will extend the simpler baseline focus of this study to include other size metrics.

We also wish to study the shape of LGMDs, and pursue analysis of their aspect ratios as a way of doing so. We define the aspect ratio A for a given LGMD as (Equation 7): A = ( ϕ E − ϕ W ) cos λ E + λ W 2 λ N − λ S   , (7)

with λ denoting the latitude and ϕ denoting the longitude of the easternmost, westernmost, northernmost, and southernmost points (subscripts E, W, N, and S, respectively) on the LGMD’s perimeter. This quantity describes the ratio of the LGMD’s longitudinal extent to its latitudinal extent, with the longitudinal extent being measured along the line of latitude halfway between the easternmost and westernmost points. Section 4 contains the results of the application of this definition to our set of computed contours.

After applying Spherical Elementary Current Systems and our contour-finding algorithm, we can study the behavior of LGMDs by examining relevant statistics of the resultant LGMD perimeter dataset. For that, we first consider the distribution of the number of LGMDs identified in each of our 395,463 heatmaps shown in Figure 3a. The distribution has a significant skew to the right, and no heatmap was ever observed to contain more than six LGMDs. A large proportion – about 58% – of heatmaps exhibit no LGMDs, though the exact amount is somewhat sensitive to the d B N threshold used to discriminate LGMDs from background geomagnetic activity.

We can also consider the perimeters of the LGMDs, as shown in Figure 3b. This figure describes the overall distribution of perimeter values. It is agnostic toward simultaneity and simply includes all LGMDs; that is, any two LGMDs contribute their perimeters toward this distribution even if they reside in the same heatmap. This distribution has a significant positive skew and shows that the mean, median, and mode LGMD perimeters are 1,620 km, 1,372 km, and 550 km respectively. Also of note is that the largest LGMD was observed to be 12,024 km in perimeter. As a point of comparison, the region of interest measures 9,344 km in perimeter, but it is possible for an LGMD perimeter to exceed this value if its associated contour forms a concave shape. However, most LGMDs are smaller than 9,000 km in perimeter, with only a handful of cases in the entire dataset close to the size of the region of interest.

Importantly, the typical perimeter of LGMDs as per these results is of the order that one would expect from estimates of the few-hundred-kilometer spatial scale derived using other methods, such as those mentioned in Section 1 (e.g., Dimmock et al. (2020); Dimitrakoudis et al. (2022)). This approximate correspondence between perimeter and diameter holds unless the LGMD in question is extremely oblong, which is quite rare according to the results contained in Figure 4.

Next, we examine the distribution of the aspect ratios of observed LGMDs. The distribution shown in Figure 4 is roughly symmetric about its central value at log ( A ) = 0.29 ( A = 1.95). This value is markedly positive, meaning that we observe LGMDs with a longitudinal extent greater than their latitudinal extent more frequently than those with the opposite orientation. In fact, only 5.7% of the mass of the log aspect ratio distribution lies below zero. This means that over 94% of observed LGMDs are at least slightly oblong along the east-west direction. The cause of this asymmetry around log ( A ) = 0 is not entirely clear, but it may be related to the fact that current systems which induce perturbations to B N , such as the auroral electrojets, are also extended longitudinally.

One must note that our choice of the region of interest could be biasing this distribution in the positive direction by “cutting off” some data points on the negative end of the distribution. This is because the region of interest is wider in the longitudinal direction than it is in the latitudinal direction. To illustrate this, an LGMD which has been found to have a longitudinal extent of 1,400 km might be counted among those in Figure 4, but one with a latitudinal extent of 1,400 km would surpass the boundaries of the region of interest; its perimeter contours would not be closed, and thus it would not be represented in Figure 4. If this bias does exist, it could be corrected by extending the region of interest to have the same latitudinal extent as its longitudinal extent. However, for our purposes, this is not practical, because there are not sufficient magnetometers with available data to support such an expansion without reducing interpolation quality. Instead, we examined the distribution of aspect ratios of only LGMDs less than 1,000 km in perimeter. The LGMDs in this subset are too small to be cut off by the boundaries of the region of interest, even if they are oblong in the north-south direction. When normalized, the distribution of this subset is almost identical to the full distribution. We take this to be sufficient evidence that the dimensions of the chosen region of interest do not significantly impact the distribution of LGMD aspect ratios, at least for statistical purposes.

We have inspected the distributions of three attributes: the number of LGMDs per heatmap, the perimeter of LGMDs, and the aspect ratio of LGMDs. For readers interested in additional details concerning these LGMD attributes, plots of the distributions of perimeter and aspect ratio, separated by the number of simultaneous LGMDs, are available in the Supplementary Material. Additional statistics, whose significance is not analyzed in this work, are also available in the Supplementary Figure S3–S8.

We are also interested in whether the phase of the solar cycle has any relationship to the statistics explored in the previous sections. Therefore, we now proceed with an analysis of two different subsets of our dataset: the set of timesteps surrounding solar maximum (14 October 2012 through 31 January 2015) and those surrounding solar minimum (1 January 2009 through 31 October 2009 and 4 March 2018 through 31 December 2019) for Solar Cycle 24. The dates defining these intervals are set following Reyes et al. (2021), who precisely defined the phase boundaries for that solar cycle.

In this section, we reproduce each histogram from the previous section to include the solar maximum and solar minimum data subsets alongside the full dataset. These subsets each contain different numbers of data points because each represents a different range of instances in time. Therefore, the distributions shown have been normalized such that the area under each is 1. This is in contrast to the distributions in Section 4.1; those essentially represent the relevant attributes’ probability distributions, whereas the plots in this section show probability density functions (PDFs) because they have been normalized in this way. We show PDFs so that the data subsets may be directly compared. The distribution of the number of LGMDs per heatmap, though it represents a discrete variable as opposed to a continuous variable, can be represented by this normalization as a PDF because the “bin” width is 1.

The effects of separating solar minimum and solar maximum from the overall dataset are shown on the number of LGMDs (Figure 5), their perimeter (Figure 6), and their longitudinal-latitudinal aspect ratios (Figure 7). For each of these three statistics, the trends of the solar maximum and solar minimum subsets follow each other very closely. This motivates the question: are the differences that do exist between the subsets simply due to aleatoric variability (that is, do they only exist by chance?), or are they caused by some higher-order physical phenomena? To answer this question, we can compute a statistical divergence between the distributions corresponding to the solar maximum and solar minimum subsets. If this divergence is sufficiently large, that will be evidence against the hypothesis that the distributions differ only because of aleatoric effects.

We first quantify the distance between a PDF associated with the solar maximum dataset, X , and that associated with the solar minimum dataset, N , using J KL , the symmetrized form of the Kullback-Leibler divergence D KL (Kullback and Leibler, 1951) (Equation 8): J KL ( X , N ) = D KL ( X , N ) + D KL ( N , X ) = ∑ i w X i ln X i N i + ∑ i w N i ln N i X i (8)

Here, the sums are over the “bins” of the discretized PDFs. The width of each bin is w , which is written without subscript because, for this study, each bin within the same PDF has the same width. We use the symmetrized form of the divergence to avoid confusion that may arise due to D KL ( X , N ) generally not being equal to D KL ( N , X ) .

From an information theory perspective, J KL measures how much information would be absent if one PDF were to be used to approximate the other. As the distance between X and N increases, J KL ( X , N ) will also increase. However, there is some ambiguity as to the threshold above which J KL can be said to indicate that the differences between X and N are due to more than just aleatoric variability. Some clarification appears as one takes J KL to various limits. First, the symmetrized divergence between a distribution and itself is zero. Second, while there is no upper bound on J KL even for PDFs, one can create a “benchmark” PDF U known a priori to be sufficiently distant from X . If J KL ( X , N ) is much closer to J KL ( X , U ) than it is to 0, that may point to a phenomenological difference between X and N .

For this comparison, we simply define U as the PDF describing a normalized uniform distribution, which is fundamentally different from any of the PDFs shown in Figures 5–7. The J KL ( X , N ) values for each of the statistics presented in this section are shown in Table 1. Also shown for comparison are the symmetrized divergences between X and U and between N and U for LGMD number, perimeter, and the base-10 logarithm of aspect ratio, log 10 ( A ) . Because any conclusion drawn about log 10 ( A ) has a linear equivalent that can be drawn about A , we will refer to the “base-10 log aspect ratio” just as the “aspect ratio” for the sake of simplicity.

It is clear from the reported values that the PDFs of each statistic shows much less divergence between solar maximum and solar minimum than they do between a uniform PDF and either solar maximum and solar minimum. The values of J KL ( X , N ) are considerably larger for perimeter and aspect ratio than they are for LGMD number, but this likely due to the larger number of bins in those two quantities’ PDFs accentuating the divergences between PDFs; J KL ( X , N ) for perimeter and aspect ratio is in each case more than an order of magnitude smaller than either J KL ( X , U ) or J KL ( N , U ) . The fact that J KL ( X , N ) is so small in comparison does not suggest than these three LGMD attributes exhibit any dependence on the phase of the solar cycle.

The information theory perspective makes the J KL criterion interpretable. However, there are some disadvantages to this approach. The effect of the logarithms in the definition of J KL (Equation 8) emphasizes the relative importance of the tails of the PDF to the symmetrized divergence. This could raise the concern that outliers in each of the attributes’ distributions may be influencing the conclusions drawn from Table 1. Additionally, bins in which one of the PDFs is equal to zero must be ignored to keep J KL well defined. This means that some small number of bins in the perimeter and aspect ratio PDFs are not taken into account by this criterion. These drawbacks do not invalidate the use of J KL in this context, but they do suggest another criterion may be needed to verify the preliminary conclusions about the effect of the solar cycle. We have used the total variation distance D TV (Equation 9) to supplement the use of J KL : D TV ( X , N ) = 1 2 ∑ i w X i − N i (9)

The total variation is computed outside of “logarithm space”, which avoids both drawbacks mentioned regarding J KL . The D TV values for each of the statistics presented in this section are shown in Table 2. The divergence between X and N remains the smallest divergence by over one order of magnitude for the number of LGMDs. For their perimeter and their aspect ratio, it is smaller by about a factor of eight. Using the total variation distance as a criterion reinforces the results obtained with the initial Kullback-Leibler analysis, which indicate insufficient evidence to suggest a phenomenological difference between solar maximum and solar minimum LGMD aspect ratios.

One might expect the solar cycle to have an effect on the PDFs of these three LGMD attributes. After all, the solar cycle has an effect on the general geoeffectiveness of the solar wind–strong geomagnetic storms, storms due to coronal mass ejections, and strong ground-level geomagnetic activity are all more common during solar maximum, for instance (Webb, 1991; Richardson and Cane, 2012; Borovsky and Denton, 2006). However, as we move away from global size scales toward the smaller scales of LGMDs, we can not reject the hypothesis that the solar maximum and solar minimum distributions differ only because of aleatoric effects.

This does not necessarily mean that the solar cycle has no influence on the properties of LGMDs, only that any influence that does exist could not be detected with the methods used here. This is perhaps due to the fairly strict choice of the uniform distribution’s PDF as a benchmark PDF. One distribution being more similar to a uniform distribution than to the other solar cycle phase’s distribution – the criterion necessary to reject the baseline hypothesis – is a high bar to meet. There still may be some much smaller solar cycle effect that could require a more sensitive test to extract. We have used the uniform distribution’s PDF as a benchmark even still because other tests would include further amounts of arbitrariness. For instance, one could repeat this analysis on LGMD aspect ratio, replacing U with the PDF of a log-normal distribution, but it would be difficult to justify a choice of mean and variance for this new benchmark. Because of this, first investigating the influence that other phenomena have on LGMD attributes would likely be more fruitful than identifying a more sensitive test.

Though this numerical analysis suggests that the distribution of aspect ratios does not strongly depend on the phase of the solar cycle, it does not reveal the physical processes that do determine the shape of the PDFs shown. This motivates a cause-and-effect analysis of LGMD aspect ratio, which likely deserves a study all to itself. Such a study could explore phenomena inside the magnetosphere, where the effect of the solar wind is modulated by a multitude of current systems, and it could connect these phenomena to the LGMD attributes lain out here. For instance, it could investigate the shape of localized structures embedded within the auroral electrojet and how they correspond to LGMDs. Similarly, heeding the results of Vandegriff et al. (2024) wherein points on Earth’s surface were magnetically mapped to spatially dispersed locations in the magnetotail, it may be of interest to examine the influence of large-scale features in the magnetosphere on LGMD aspect ratio.

Better establishing the empirical relationship between the solar wind, magnetospheric and ionospheric processes, and the characteristics of LGMDs will provide knowledge about the Sun-to-Earth phenomenological chain, and help to describe the link between LGMDs’ attributes and their drivers.