Time-scale dependence of solar wind-based regression models of ionospheric electrodynamics

The solar wind influence on geospace can be described as the sum of a directly driven component, or dayside reconnection, and an unloading component, associated with the release of magnetic energy via nightside reconnection. The two processes are poorly correlated on short time scales, but exactly equal when averaged over long time windows. Because of this peculiar property, regression models of ionospheric electrodynamics that are based on solar wind data are time scale specific: Models derived from 1 min resolution data will be different from models derived from hourly, daily, or monthly data. We explain and quantify this effect on simple linear regression models of various geomagnetic indices. We also derive a time scale-dependent correction factor that can be used with the Average Magnetic field and Polar current System model. Finally, we show how absolute estimates of the nightside reconnection rate can be calculated from solar wind measurements and geomagnetic indices.

Scientific RepoRtS | (2020) 10:16406 | https://doi.org/10.1038/s41598-020-73532-z www.nature.com/scientificreports/ The parameters in Eq. (2) are functions of τ because of the two-step response to solar wind driving discussed above. When τ is small (i.e., of the order of minutes or less), d(τ ) is correspondingly small since D and N are not closely related on such short time scales. On the other hand, in the limit τ → ∞ , c(∞) = 0 and d(∞) = 1 since dayside and nightside reconnection rates must balance (i.e., ∞ D = ∞ N ). Because of the time scale dependent relationship between D and N , statistical models of some measurable quantity related to magnetosphere/ionosphere convection that use D are also dependent on time scales. To quantify this dependence, assume that such a measurable quantity, y, depends on V as follows: y is here a generic term for any parameter that follows Eq. (4), and will later be replaced by specific geomagnetic indices. The physical justification for Eq. (4) is discussed further in the next section. This equation is assumed to be valid on all time scales, and α and β do not depend on τ . Since we do not know V, Eq. (4) is not very useful for making empirical models of y. Instead, it is common to use a solar wind magnetosphere coupling function [9][10][11] , ε . Such coupling functions combine solar wind measurements in different ways to maximize the correlation with the transfer of energy or magnetic flux from the solar wind to the magnetosphere. We will assume that they correlate with the transfer of mganetic flux, i.e., the dayside reconnection rate D : where k is a proportionality constant that scales ε such that its amplitude matches that of the dayside reconnection rate. The unit of Eq. (5) is magnetic flux per second (Wb/s), which is equal to volt (V). The Milan coupling function 10 is a special case, where the function includes an empirically determined scale factor so that k = 1.
Using the above equations, we can derive a model of y τ that depends only on ε τ . The model will be time scaledependent (hence the superscripts τ ), since we make use of Eq. (2): Starting with the time-averaged Eq. (4), we express V τ in terms of ε τ using Eqs. (1), (2), and (5): where Equation (6) is a linear model of y τ in terms of ε τ , with model parameters a and b that depend on τ since c(τ ) and d(τ ) from Eq. (2) are time-scale dependent. The implication of this time scale dependence is that a regression model on this form which is derived using data with time resolution τ in general should not be compared with data having a different time resolution. In the next section we demonstrate this by estimating a and b for a range of time scales by replacing y with specific geomagnetic indices. Later we discuss how the time-scale dependence of solar-wind based regression models affects the interpretation of more complex empirical models, using the Average Magnetic field and Polar current System (AMPS) model 3 as an example. The AMPS model is an empirical global model of ionospheric currents and associated magnetic field based on magnetic field measurements in low Earth orbit from the CHAMP and Swarm satellites.
In the limit τ → ∞ , Eq. (7) reduces to a(∞) = α since c(∞) = 0 , and Eq. (8) reduces to b(∞) = βk since d(∞) = 1 . Using the above results, we can derive an expression for N in terms of y and D that is independent of time scales: By solving Eq. (1) for N , and using (4) to replace V, we get where in the last step, we used Eqs. (7) and (8) in the limits τ → ∞ to replace α and β . As will be demonstrated in the next section, a(∞) and b(∞) can be estimated from data, by averaging over several days. If we use a coupling function that is scaled so that k = � D /ε = 1 , such as the function developed by 10 , N can be calculated from observations of y and of the solar wind. Equation (9) is in principle valid on any time scale of N , D , and y. We will later put this equation to use and discuss its limitations. time-scale dependence of solar wind based linear regression models of geomagnetic indices. Here we apply the equations from the previous section to a set of geomagnetic indices: AL, AU, PC, and ASY-H. These are four geomagnetic indices with 1 min time resolution, which are believed to describe various aspects of geomagnetic activity. The AL and AU indices monitor the westward and eastward auroral electrojets, respectively. The PC (polar cap) index is based on a magnetometer that is typically poleward of the auroral oval, and it is ostensibly an indirect measurement of the magnitude of the solar wind-magnetosphere coupling. The ASY-H index is a measure of the longitudinal variability of the magnetic perturbations at mid latitudes, and it tends to correlate with the auroral electrojet indices. For definitions and details about interpretation, see the review by Kauristie et al. 12 , and references therein. We will show that these parameters, when modeled as functions of a solar wind coupling function ε , depend on the time scale τ.
(4) y = α + βV . www.nature.com/scientificreports/ For the sake of argument, we assume that each index is a linear function of dayside and nightside reconnection rates, and nothing else. In other words, they obey (4), but with different values for α and β . This is of course an approximation, since they would otherwise be exactly correlated. For example, we ignore conductivity effects for all indices except AU, which we scale by (1 + sin ψ) −1 , where ψ is the dipole tilt angle. This scaling is meant to roughly account for seasonal variations in solar EUV induced conductivity. The scaling improves the results with AU, but not with the other indices which may be more dependent on precipitation induced conductivity. In this analysis ε is the coupling function defined by Newell et al. 9 , but similar results are obtained when using other functions 10,11 . The Newell coupling function depends on the solar wind speed and the component of the interplanetary mganetic field (IMF) that is perpendicular to the Sun-Earth line. The solar wind data and indices are taken from the 1 min resolution OMNI database, between 1995 and 2018. Figure 1A shows, for each of the four indices, the ratio b(τ )/b(∞) where b(τ ) is defined in Eq. (6) and estimated using ordinary least squares. b(∞) was calculated with τ = 8 days. The ratio b(τ )/b(∞) is related to the slope d(τ ) in Eq. (2), which describes the relationship between dayside and nightside reconnection rates on time scale τ . Specifically, d(τ ) = 2b(τ )/b(∞) − 1 , which can be found by solving Eq. (8) for d(τ ) and using that βk = b(∞) . All curves in Fig. 1A follow a similar variation with τ , with a steep increase at τ < 3 h ( τ = 3 h is indicated with a dashed bar), and then a more gradual increase towards 1. The transition at about 3 h is expected from previous studies 13,14 indicating that this is an average substorm cycle period (i.e., the time between bursts of tail reconnection associated with individual substorms). Figure 1B shows the squared of the Pearson correlation coefficient r 2 between model and data. This coefficient is a measure of the fraction of variation in data that is explained by the model. For each index, r 2 is rapidly increasing for larger τ . For the AL and PC indices, r 2 > 0.8 at large τ . This indicates that (4) is indeed a good representation of these indices, and it reflects that ε = k� D = k� N as τ → ∞ . Figure 1C shows the root mean square error (RMSE) of the model compared to the data, relative to the RMSE at τ = 1 min. The basic result shown in Fig. 1 is that the models improve as τ increases; the physical reason for this trend is that ε becomes representative of both dayside and nightside reconnection.
The indices studied here follow a skewed distribution. The peak of the distribution is near zero for each of the four indices that we consider. This represents a quiet day baseline. Active events tend to give negative excursions in the AL index and positive excursions in the other three indices. Since the indices do not follow a normal distribution, a quantitative interpretation of the Pearson correlation coefficient is questionable. However, in this study, the main purpose of the correlation coefficient is to qualitatively indicate how the association between our simple model and the observation that it describes vary with time scales. Given the very large number of measurements ( > 20 years of measurements at 1-min resolution), we believe this use is appropriate.

Effect of time-scale dependence on climatological models
Several climatological models of ionospheric electrodynamics depend on solar wind parameters, and lack parameters that represent the high-frequency portion of N 2,3 . The AMPS model 3 describes the ionospheric disturbance magnetic field in space, and the associated ionospheric currents. It is formulated in a way that makes it possible to write the model magnetic field as where B 0 and B ε are two different functions of space and of the external parameters used in the AMPS model (the IMF, solar wind velocity, dipole tilt angle and F 10.7 solar flux index). In the sum in Eq. (10), B ε is multiplied by Newell's solar wind-magnetosphere coupling function 9 ε , and B 0 is not. This is conceptually similar to the much simpler Eq. (6), and it implies that the AMPS model may also be time-scale dependent.
Relating magnetic field disturbances in space and ground. To address the time scale dependence of the AMPS model, we compare model estimates with measurements of the magnetic field disturbances from ground magnetometers, instead of satellites. This is because time averaging satellite measurements mixes spatial and temporal effects. However, the AMPS model is based only on magnetic field measurements in space, above the conducting layer of the ionosphere where horizontal currents flow. This region is the only place where the full current system can be estimated using only magnetic field measurements 15 . The reason for this is that below the ionosphere, the magnetic field is partially canceled. In this section we derive mathematical expressions which relate ground magnetic field perturbations to the AMPS model coefficients, thus allowing for calculation of model predictions of the AL index, which we will use in comparisons with measurements in the next section.
We model the ionospheric current system as a two-dimensional sheet current on a sphere at some height h R , set to 110 km here. The sheet current is connected to the magnetosphere by a vertical volume current that extends out from the sphere. At polar latitudes, where the Earth's magnetic field lines are almost vertical, such currents are approximately equal to the Birkeland currents (magnetic field-aligned electric currents). The sheet current can be decomposed further into divergence-free ( J df ) and curl-free ( J cf ) components. The latter is related to the radial current density J u by the current continuity equation ∇ · J cf = −J u . The magnetic fields of J cf and J u cancel below the ionosphere, according to the Fukushima theorem 16 . Consequently, only the magnetic field of J df must be considered. From now on we call that magnetic field B , disregarding contributions from Birkeland/curl-free currents. We also assume that that contributions from the Earth's core, crust, and large-scale magnetosphere are accounted for by subtracting CHAOS-6 model predictions 17 . The CHAOS-6 model is a high resolution, regularly updated, geomagnetic field model derived from magnetometer measurements from ground observatories and several satellites in low Earth orbit.
Above and below the current sheet, the magnetic field of J df can be expressed as a gradient, B = −∇V . In the AMPS model, V is expressed as www.nature.com/scientificreports/  www.nature.com/scientificreports/ where q is quasi-dipole latitude 18 , φ mlt is the magnetic local time 19 , h (km) is the altitude, R E = 6371.2 km is the Earth radius, P m n (sin q ) are Schmidt semi-normalized associated Legendre functions of degree n and order m, and the spherical harmonic coefficients g m n and h m n are specific functions of solar wind speed, IMF B y and B z , dipole tilt angle, and the F 10.7 solar flux index 3 . In the AMPS model, the truncation level of the double sum is set to n ≤ 45, m ≤ 3 . The radial dependence of Eq. (11) corresponds to a magnetic field of internal origin 20 . It can be modeled by a sheet current J V at some height h R < h , where J V = k × ∇ , k is an upward unit vector, and where µ 0 is the vacuum permeability ( 4π · 10 −7 H/m). Since Eq. (11) refers to a magnetic field of internal origin, relative to altitudes above the ionospheric horizontal currents, it is not suitable for use below the ionosphere, where the ionospheric currents are external. Below the ionospheric currents, the magnetic potential field is often split in two parts, B = −∇V e − ∇V i , one ( V e ) that corresponds to external sources (ionospheric currents) and another ( V i ) that corresponds to currents below the Earth's surface, that are induced by the external currents. A spherical harmonic expansion of V i is exactly analogous to the expansion of V in Eq. (11), while the expansion of V e has a different radial dependence: Just like the internal potential of Eq. (11), the external potential can be modeled by a sheet current J e at height h R > h , where J e = k × ∇ e and We now make the assumption that the currents J V = J e = J df , so that = e + c , where c is an arbitrary constant. This equation allows us to express the spherical harmonic coefficients in Eq. (14), a m n and b m n , in terms of the AMPS model coefficients, g m n and h m n : These expressions can be inserted in Eq. (13), and thus we can calculate the magnetic field component of the external magnetic field below the ionosphere, B e = −∇V e , in terms of AMPS model coefficients: The magnetic field components B φ , B q , and B r , must be regarded as quasi-dipole components. For conversion to geographic east, west, and upward components, we use the quasi-dipole base vectors g i and f i , which vary with the geometry of the main magnetic field 19 : At least two critical approximations are involved in the approach described above: (1) The choice of 110 km altitude is somewhat arbitrary; lowering it would increase the ground magnetic field; (2) We ignore effects of ground-induced currents. Induced effects are in theory present in the AMPS field, but they will be stronger on ground, and will also be time-scale dependent 21 .
In our analysis we derive a synthetic AL index from the AMPS model, and compare with measured AL values. We chose the AL index since it measures the westward electrojet, which typically peaks on the nightside in response to increases in magnetotail activity, most notably substorms. It is also comparatively straightforward to calculate: The synthetic AL values are calculated as the lower envelope of the H component of the AMPS model magnetic field at the AL station locations. The H component is the magnetic field along the horizontal part of Earth's main magnetic field (here we use the CHAOS model 17 ). An updated version (version 0103) of the AMPS model was used, defined and calculated as described in the original AMPS model paper 3 , but with a 43% larger dataset which includes more recent Swarm data. In the updated version we have also corrected a programming error that damped the longitudinal variation in the toroidal magnetic field and associated magnetic field-aligned electric currents (FACs) of the original model. The correction leads to peak FACs that are significantly increased  www.nature.com/scientificreports/ (by typically ∼ 20% ) compared to the original model, however the integrated FACs differ by typically less than ∼ 10% . The changes in poloidal magnetic field and associated horizontal divergence-free currents and ground magnetic field disturbances are small (typically ∼ 1% ). The new version is now integrated in the latest update (version 1.4.0) of the publicly available AMPS forward code, pyAMPS 22 , which was used for this study. The model is also published as an ESA Swarm Level 2 data product. Figure 2 summarizes the results of the comparison between AMPS model predictions and ground magnetic field measurements. To the left, we show the r 2 and RMSE between model and measured AL, in black. The RMSE is calculated using iterative reweighting with Huber weights, to reduce the effect of outliers. The coefficient of determination ( r 2 ) increases with time scale, as expected. The increase in RMSE with time scale is more surprising. The model increasingly underestimates the data points as time scale increases. The reason for this is probably the time-scale dependent relationship between ε and N : The AMPS model is based on relatively high time resolution data, and is thus scaled to fit the directly driven part of the magnetic disturbances. When compared with data on longer time scales, it correlates with both the loading and unloading parts of the disturbances, but the amplitudes still represent only one part. The time-scale dependence discussed in the previous sections suggest that the AMPS model estimates could be improved by including a time scale-dependent scale factor for the εB ε term in Eq. (10). We estimate such a scale factor using the model and measured AL time series from the period 1995 to 2018. For each time scale considered in the scatter plots in Fig. 2, we find the scale factor that minimizes the Huber weighted root mean square model data mismatch. The resulting points, plotted against time scale τ , are approximately sigmoidal, and we therefore estimate the τ-dependent scaling factor as a logistic function plus a constant:  Fig. 2A and B. We see that the correlation is almost unaffected. This is expected, because the scaled version does not include any more information about the high-frequency component of nightside processes. However, the RMSE is significantly improved: At time scales of the order 24 h the RMSE is only 20 nT, and about 90% of the variation is explained. Example time series are shown in Fig. 2C,D. Figure 2C shows the AL index measurements averaged over 5 h in grey, and corresponding AMPS estimates in black and orange. The black curve is unscaled, and clearly underestimates the fluctuations in AL. The orange curve is scaled using C(τ ) in Eq. (20), and is much closer to the measurements. Figure 2D shows measurements of the ground magnetic field at the Thule station in Northern Greenland together with AMPS model estimates. The plot represents perturbations in the southward direction averaged over 5 h. The black curve represents unscaled model estimates, and the orange curve represents model estimates scaled using C(τ ) in Eq. (20). The close fit indicates that the scale factor derived using AL index comparisons can be applied to other ground magnetometers as well.

Results
The scale factor C(τ ) given by Eq. (20) is 1.51 at small τ . This may imply that the AMPS model ground perturbations underestimate the measurements at short time scales, even though the model was derived using high time resolution data: 1 Hz satellite measurements and τ = 20 min solar wind data. It is possible that this is because the AMPS model lacks the spatial resolution to see localized instantaneous features. It is also possible that this discrepancy is due to the neglect of induced effects in our calculations of ground magnetic field disturbances. The induction effect, which is likewise time scale dependent, also influences the scale factor C(τ ) in Eq. (20), such that C(τ ) can not be interpreted as a pure effect of the time-scale dependence between dayside and nightside reconnection rates. estimating nightside reconnection rate using geomagnetic indices and solar wind data. We have shown how measurements of y and D can be used to estimate the nightside reconnection rate N on any time scale via Eq. (9).  23 . We see that both substorm onsets are preceded by a period of strong dayside reconnection (growth phase), and followed by a period of strong tail reconnection (substorm expansion). Figure 3 shows only a small extract of a time series calculated for the years 1995 through 2014, containing 8.3 million 1 min values. During these years the Pearson correlation coefficient between N and D is 0.2, which means that only 4% of the variation in N can be explained by D for a time scale τ = 1 min. This tells us that N , estimated from Eq. (9), contains information which is not in D . Since the calculation can be done whenever solar wind data is available to calculate D , and since geomagnetic indicies typically have very few data gaps, this may be a useful parameter to include in future empirical models of ionospheric electrodynamics.
The correlation r ( r 2 ) between D and N at τ = 8 days is 0.85 (0.72). If Eq. (9) and Milan's coupling function 10 gave perfect estimates of N and D , we would instead find r ≈ 1 at this time scale. The difference indicates that there is a need for closer scrutiny of the N estimates. In future work we plan to carry out more detailed calibration and validation of our estimates of N by comparing with other independent estimates. For example,   [24][25][26] . It can in principle also be done by inferring the polar cap boundary position from magnetometer measurements at low Earth orbit 27 .

Discussion
The two-step response of the magnetosphere-ionosphere system to changes in the solar wind has been known for decades, and it is an integral part of the expanding contracting polar cap paradigm 7,8 . The contribution of the present study is to quantify the time-scale dependent statistical relationship between nightside and dayside reconnection rates, as well as its effect on solar wind based regression models of ionospheric electrodynamic parameters; and to use this relationship to derive Eq. (9), which following more detailed calibration and validation may prove useful as a monitor of the nightside reconnection rate. We have shown that solar wind-based regression models of ionospheric electrodynamics are time scale dependent: Model parameters depend on the time resolution of the input data. We have demonstrated this with simple linear regression models of the AL, AU, PC, and ASY-H indices. We have also shown that this applies to the more complex AMPS model. By comparing AMPS model predictions with measured values of the AL index, we have derived the time-scale dependent correction factor C(τ ) given in Eq. (20), which can be used with the AMPS model when estimating time-averaged magnetic field disturbances. The Python AMPS forward code, pyAMPS 22 , includes support for this scaling in the get_B_ground function, via the epsilon_multiplier keyword argument.
We have also presented a simple model for the nightside reconnection rate, Eq. (9), which depends on geomagnetic indices and on the dayside reconnection rate. The latter can be estimated from solar wind parameters 10 . This derivation assumes that the geomagnetic indices we have considered are proportional to the cross polar cap potential (i.e., they obey Eq. (4)), and takes into account that nightside and dayside reconnection rates balance over time. The high correlation between linear models and corresponding indices on long time scales, seen in Fig. 1, indicate that this assumption is reasonable. However, it is not perfect, since that would imply that all the indices are perfectly correlated with each other. In this study we have not attempted to assess which index, when used in Eq. (9), gives the most accurate estimate of nightside reconnection rate. By comparing different realizations of this equation, using combinations of indices with independent estimates of nightside reconnection, it should be possible to derive an empirical function that describes the "unloading" response, similar to how empirical solar wind-magnetosphere coupling functions describe "loading", or direct driving. Such a function will be the topic of a future study. The inclusion of an unloading parameter in physics-based empirical models, like the AMPS model, may help to determine the role of individual driving processes in controlling the ionospheric electrodynamics.