On optimum solar wind – magnetosphere coupling functions for transpolar voltage and planetary geomagnetic activity

.


Introduction
Coupling functions are combinations of interplanetary parameters that are used to quantitatively predict terrestrial space weather indicators and indices.They should have a linear relationship with the index or measured parameter that they aim to predict.There are a great many combinations that have been proposed and tested since correlations between interplanetary parameters measured by spacecraft and terrestrial disturbance indices became possible (Arnoldy, 1971).The concept of a combination of parameters capturing their net influence (i.e., a coupling function) grew out of the PhD studies of Perreault (1974).An excellent review of the development of coupling functions, the theories behind them and the empirical fits, has been given by McPherron et al. (2015).

10.1029/2021JA029946
2 of 41 testing the coupling function in unusual regions of parameter space; however, as always with models, the validity of the results depends on the assumptions, parameterizations, and resolutions used in setting up the model.
Coupling functions have generally, but not exclusively, taken the basic mathematical form of the product of measured parameters, each to the power of an exponent.Parameters used have been the interplanetary magnetic field (IMF), B = |B| or its transverse component perpendicular to the Sun-Earth line, B ⊥ ; the solar wind speed, V SW ; the solar wind number density N SW or its mass density ρ SW = m SW N SW (where m SW is the mean ion mass); and (for timescales shorter than about 1 yr), a factor to allow for the orientation of the IMF in the Geocentric Solar Magnetospheric (GSM) frame of reference, such as a function of the clock angle in GSM, θ.Note that for coupling functions θ is defined as tan −1 (|B Y |/B Z ) where B Y and B Z are the Y and Z components of the IMF in GSM: the use of the modulus of B Y means that any effects on coupling associated with its polarity would have to be added separately but avoids complications caused by discontinuities in clock angle between 0 and 2π that would otherwise arise.We here denote magnetic field exponents by a, mass density or number density exponents by b, solar wind speed exponents by c, and IMF orientation factor exponents by d.Some improvements to this basic multiplicative form have been suggested in the form of additive terms.For example, Newell et al. (2008) proposed the addition of a term designed to predict the dayside magnetopause reconnection voltage with a smaller term to predict the voltage generated by non-reconnection "viscous-like" interaction.Lockwood (2019) proposed a development to energy-transfer coupling functions whereby, in addition to the energy extracted from the dominant energy flux in the solar wind (namely the kinetic energy flux of the particles), the smaller one due to the solar wind Poynting flux is added.Given that the Poynting flux in the solar wind is two orders of magnitude smaller than the particle kinetic energy flux, this appears an unnecessary complication: however, the Poynting flux enters the magnetosphere without the relative inefficiency with which kinetic energy of the solar wind is converted into Poynting flux by currents flowing in the bow shock, magnetosheath and magnetopause (Cowley, 1991;Ebihara et al., 2019;Lockwood, 2004).
Other, more complex forms, with combinations of additive and multiplicative terms have been proposed (e.g., Borovsky, 2013;Luo et al., 2013).The formulation of Luo et al. (2013) aims to take account of daily and seasonal variations in the terrestrial space weather index predicted (that are due to station locations and orientation of the Earth's dipole) and non-linearities caused by the expansion and contraction of the polar cap as solar wind driving varies.It also removes rapid fluctuations using low-pass filters.The result is that it is highly complex and, as noted by McPherron et al. (2015), it is unclear how many free parameters are present in this coupling function, but they estimate that it is at least 35.Because these more complex formulations add to the number of free fit parameters, this greatly increases the problem of statistical "overfitting" (Chicco, 2017).Overfitting occurs when a fit has too many degrees of freedom and it can start to fit to the noise in the training data, which is not the same as the noise in the test or operational data.As a result, the fit has reduced predictive accuracy.This is a recognized pitfall when signal-to-noise ratio in the data is low, as is usually the case in disciplines such as climate science (Knutti et al., 2006) or population growth (Knape & de Valpine, 2011), but has not often been considered in space physics in the past.However, this is now changing with the advent of systems analysis of the magnetosphere and the application of machine-learning techniques to space weather data (e.g., Camporeale, 2019;Stephens et al., 2020).Overfitting is a problem for the generation of coupling functions because there are a great many sources of noise, not all of which have been recognized and some of which we cannot do much about when we take note of the need to have large datasets to cover all potential regions of solar wind/magnetosphere parameter space.The noise sources in correlative solar wind-magnetosphere studies include: instrumental observation errors in interplanetary measurements and in the terrestrial disturbance index or indicator; propagation errors between the spacecraft observing the solar wind conditions and near-Earth space (these include using the correct time lag but, more importantly, spatial structure in interplanetary space that means the solar wind sampled by the spacecraft is not always the same as that which impinges on Earth's magnetosphere); the effects of the bow shock in changing near-magnetopause characteristics of the magnetosheath, relative to the undisturbed solar wind outside the bow shock; gaps in data sequences; effects of averaging and timescale; non-linear responses of the magnetosphere; pre-conditioning of the magnetosphere, and the effects of prior solar wind/magnetosphere coupling history; dipole tilt effects on ionospheric conductivities, magnetospheric structure, and current sheets.

10.1029/2021JA029946
3 of 41 statistically significant improvement (at over the 3σ level), the improvement for annual or Carrington rotation means was not statistically significant: hence in the latter cases no statistically significant improvement was achieved, despite the number of free fit variables being doubled from 1 to 2 and the additional term being based on theory.It should also be noted that the branching ratios used with additive terms can become inappropriate if conditions move outside the region of parameter space that was used to derive them.Of course, the same is true for all coefficients and exponents in any coupling function, but effects are particularly serious for coupling functions using additive terms because an incorrect branching ratio can cause one of the two terms to become spuriously dominant.A common example is averaging timescale which, in general, has different effects on different terms and so the ratio of the two that is appropriate to one timescale may not apply on another.
Table 1 gives the timescale τ on which each coupling function was derived and/or has been tested and/or deployed.It is noticeable that at larger τ, simpler coupling functions have been very successful in yielding very high correlations (Finch & Lockwood, 2007).These high correlations are achieved because averaging over long intervals gives cancellation of noise.The averaging timescale of the interplanetary and the terrestrial data that are compared is a crucial consideration because solar wind parameters have a variety of autocorrelation times which means that their distributions of values change with τ in different ways (Lockwood, Bentley, Owens, Barnard, Scott, Watt, & Allanson, 2019;Lockwood, Bentley, Owens, Barnard, Scott, Watt, et al., 2019).However, this is not often considered when compiling a coupling function and τ is not even explicitly defined in several publications (in several cases in Table 1, τ could only be inferred from the data plots presented).
One idea that has been proposed is that there is a "universal coupling function" that best predicts all terrestrial space weather indices and indicators (Newell et al., 2007(Newell et al., , 2008)).This idea runs counter to the method now routinely used to reconstruct interplanetary parameters from historic observations of geomagnetic activity.These reconstructions exploit the finding that different geomagnetic indices have different responses to interplanetary parameters and so combinations of historic index observations can be used to infer the separate interplanetary parameters.This was inherent in the reconstruction of open solar flux from historic observations of geomagnetic activity by Lockwood et al. (1999) but first used to extract more than one parameter by Svalgaard et al. (2003), who noted that on annual timescales the IMF B and solar wind speed V SW could both be derived from any combination of geomagnetic indices that had different dependencies on these two parameters (i.e., different optimum coupling functions).This has been exploited by Svalgaard and Cliver (2007), Rouillard et al. (2007), Lockwood et al. (2009), Lockwood andOwens (2011), andLockwood et al. (2014).These methods and results have developed from simple single fits to large ensembles of fits allowing for uncertainties and been reviewed by Lockwood (2013).If different indicators of geomagnetic activity have different optimum coupling functions, it means that other space weather activity indicators, such as transpolar voltage, cannot share the same optimum coupling as all, if any, of the geomagnetic activity indices.
We here investigate the differences between the optimum coupling functions for transpolar voltage Φ PC , the global am geomagnetic index, and the nightside northern hemisphere auroral oval index, SML, using simultaneous data.The am index (Mayaud, 1980) has been shown to have the most uniform response to solar wind forcing with Universal Time and time of year by virtue of the relative uniformity of the observing network and its use of area-based weighting functions (Lockwood, Chambodut, et al., 2019).However, it has the disadvantage of a time resolution of only 3 hr.The SuperMAG project's SML index (Newell & Gjerloev, 2011) has the great advantage of being available at 1 min resolution.It is based of the same concepts as the auroral electrojet AL index (Davis & Sugiura, 1966) and is, like AL, a sensitive indicator of the substorm current wedge.We use it here because it employs data from magnetometers at middle and low geomagnetic latitudes as well as at auroral latitudes whereas AL uses only a ring of auroral stations.This means that SML avoids the non-linearity of the AL response during

References
Power input to the magnetosphere stands for "half-wave rectified").Column 8 gives the time resolution of the data on which the function was mainly developed and used.The last column is a reference to a paper using or proposing the formulation.In some cases the formulation is not proposed as a viable coupling function and has only used to make comparisons with proposed coupling functions, some are physical properties of the interplanetary medium and given here only to record the exponents a, b, and c that they yield.
6 of 41 very disturbed times when the auroral oval moves a long way equatorward of the AL stations.That having been said, we have carried out all the studies presented here using both SML and AL and the results are very similar indeed.Neither SML nor AL are global indices in the same way as am, because both use only data from northern hemisphere stations.This means that they have an additional noise source caused by the seasonal effect on the electrical conductivities of the northern hemisphere ionosphere.This noise is only averaged out in studies that use averages taken over multiples of whole years.The dataset of transpolar voltage Φ PC measurements used here are the hourly means of 2-min integrations of SuperDARN radar data from 1995 to 2019 presented by Lockwood and McWilliams (2021).Like the SML and AL data they are from the northern hemisphere only as radar coverage of the southern hemisphere polar and auroral regions is not yet adequate.Even in the northern hemisphere there are gaps in the radar coverage and, furthermore, within the radars' fields of view there is often a lack of radar echoes (see, e.g., Shepherd, 2007).To deal with this, the Doppler shifts of radar echoes are fitted to a convection model in the "map-potential technique."Lockwood and McWilliams (2021) found good agreement between hourly means of 2-min integrations of the derived Φ PC values with those from over-passes by Low-Earth Orbit (LEO) satellites, provided that the mean number of radar echoes for the 30 2-min integrations exceeded 255.This yielded a dataset of 65,133 valid hourly Φ PC values from the 210,384 potential values in the interval 1995-2020, a ratio of roughly 30%.We here use hourly means of SML from the same hours as the valid Φ PC values and interpolate the three-hourly am values to the same times as the valid Φ PC estimates.
Table 1 shows that many of the proposed coupling functions predict a role of solar wind number density N SW or mass density ρ SW = m SW N SW (where m SW is the mean ion mass) as contributing to solar wind energy coupling and/ or to the driving of magnetospheric convection.For energy considerations, this is mainly because ρ SW and N SW control the dominant (kinetic) energy flux in the solar wind, ρ SW V SW 3 /2, but it has been shown that solar wind dynamic pressure (P SW = ρ SW V SW 2 ) also has an independent effect (Lockwood, McWilliams, et al., 2020;Lockwood, Owens, Barnard, Haines, et al., 2020;Lockwood, Owens, Barnard, Watt, et al., 2020).This is partly through altering the cross-sectional area that the magnetosphere presents to the solar wind flow (Vasyliunas et al., 1982) and also via the compression of the near-Earth tail, which enhances the magnetic energy density stored there for a given open magnetospheric flux, thereby enhancing the current in the auroral electrojet of the substorm current wedge when that stored energy is released during a substorm expansion phase (see review by Lockwood, 2013).Such a dependence of geomagnetic disturbance in the substorm current wedge region was isolated and identified by Finch et al. (2008).This would be in addition to the dependence on ρ SW and V SW due to the energy flux in the solar wind and/or any effect on the magnetic reconnection at the magnetopause which generates the open flux.In addition, the squeezing of the near-Earth tail by P SW would elevate the magnetic shear across the cross-tail current sheet, and hence the total current in that sheet.This could enhance the nightside reconnection voltage Φ N that closes open field lines.The expanding contracting polar cap (ECPC) model predicts that this would elevate the transpolar voltage Φ PC which is influenced at any one instant by the reconnection voltages in both the dayside magnetopause Φ D and the cross-tail current sheet Φ N (Cowley & Lockwood, 1992;Lockwood, 1991;Lockwood & McWilliams, 2021).However, we need to consider the averaging timescale used, τ.If τ is short compared to the substorm cycle duration we would expect Φ PC to reflect the enhanced Φ N , and so show some dependence on P SW from this effect of squeezing the tail.On the other hand, if τ is long compared to the substorm cycle duration, the average Φ N tends to Φ D and we would therefore expect Φ PC to show only any dependence that Φ D has on P SW which appears to be considerably smaller (Lockwood & McWilliams, 2021).However, we note that it has long been proposed that P SW has an effect on Φ D by increasing the magnetic field strength in the near-magnetopause sheath and hence the magnetic shear across the dayside magnetopause during southward IMF (e.g., Scurry & Russell, 1991).Such an effect was inferred from cusp latitude shifts by Newell and Meng (1994).
This discussion of the role of solar wind dynamic pressure is just one example of an important general pointnamely that there are a great many processes simultaneously at play in driving the terrestrial space weather response.To allow for these, solar wind coupling functions have evolved away from having theoretically derived exponents a, b, c, and d (which were often integers or ratios of integers) to empirically fitted non-integer values.Hence for the example of P SW effects on the near-Earth tail we do not complicate the coupling function with an additional term or weighting branching ratio, rather we allow the exponents b and c (in the terms ρ SW b and V SW c ) to vary to allow for such an effect and we would expect such an effect of P SW to raise the exponent b and raise c by twice as much.Hence combinations of mechanisms can be allowed for as long as their effects are multiplicative.To bring theoretical and empirical approaches together, Borovsky (2013) used the approach of making a complex 7 of 41 theoretical derivation and then reducing to a simple multiplicative form with approximations to derive exponents; however, the uncertainties introduced by any one approximation are not always apparent.
There is one last important point to note about coupling functions that is discussed further in the final section of the present paper.None of the forms listed in Table 1 allow for the pre-existing state of the magnetosphere.There are many reasons to expect non-linear magnetospheric responses.For example, the response to a given solar wind forcing quantified by a coupling function will depend on how much open magnetospheric flux already exists at the time but in addition is very likely to also depend on how enhanced the ring current is at the time and/or the state of the mid-tail plasma sheet and cross-tail current sheet.These effects all depend upon the prior history of solar wind-magnetosphere coupling.There are also regular diurnal and annual effects to consider such as dipole tilt effects and seasonal effects in the ionosphere.If they are neglected, all these factors are a source of noise for correlation studies between interplanetary coupling functions and terrestrial disturbance indices.
In this article, we do not attempt to compare the performance of the large number of proposed coupling functions.Such test have been carried out in the past, often as part of an evaluation of a newly proposed function.Detailed tests against model output were carried out for three coupling functions by Spencer et al. (2009) and the performance of seven coupling functions in predicting mid-latitude geomagnetic range indices was compared for a range of timescales τ between 1 day and 1 yr by Finch and Lockwood (2007).Newell et al. (2007) compared 20 coupling functions against 10 terrestrial indices at hourly resolution.Rather than compare individual coupling functions, we here establish some general principles using a few examples and apply a generalized common form of coupling function to an unprecedently large dataset containing two different indicators of terrestrial space weather disturbance (the transpolar voltage and two geomagnetic indices) to see if they are significantly different or can be predicted by a common "universal" coupling function.

Coupling Functions Based on Energy Considerations
Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson (2019) and Lockwood, Bentley, Owens, Barnard, Scott, Watt, et al. (2019) have shown that the am, AL, and SML geomagnetic indices, which all respond primarily to the substorm current wedge, are all well predicted over a range of timescales by the estimated power input to the magnetosphere, P α (Vasyliunas et al., 1982).This coupling function is given by the product of the dominant energy flux in the solar wind (due to the kinetic energy of the particles), the cross-sectional area of the dayside magnetosphere that it is incident upon, and a dimensionless transfer function (t r the fraction of the incident power that crosses the magnetopause into the magnetosphere).
where L o is the effective radius of cross-section of the magnetosphere presented to the solar wind flow.
The dayside magnetosphere is assumed to be constant in shape so that L o = cL s where c = L o /L s is the dayside magnetopause shape factor (assumed constant) and L s is the stand-off distance of the nose of the magnetosphere which is derived from pressure balance between the geomagnetic field and dynamic pressure of the solar wind, P SW (Farrugia et al., 1989): where k 1 is the pressure factor for shocked supersonic flow around a blunt nose object, M E is the magnetic moment of the Earth, and μ o is the permeability of free space (a.k.a., the magnetic constant).Vasyliunas et al. (1982) use a dimensionless transfer function of the form: where the solar wind Alfvén Mach number is M A = V SW (μ o ρ SW ) 1/2 /B, and k 2 is a constant and α is called the "coupling exponent" that arises from the unknown dependence of t r on M A and is the one free fit parameter.θ is the IMF clock angle in the GSM frame of reference.The dependence of t r on M A arises from the fact that the dominant energy flux in the undisturbed solar wind, the kinetic energy flux of the particles, is converted into the Poynting flux that enters the magnetosphere by the currents that flow in the bow shock and magnetosheath (Cowley, 1991;Ebihara et al., 2019;Lockwood, 2004Lockwood, , 2019)).From Equations 1-3 10.1029/2021JA029946 8 of 41 where {M E 2/3 )} are rolled into the constant k.However, note that the secular variation in M E , and hence k, can be allowed for in long-term reconstructions of space weather conditions using models of the intrinsic geomagnetic field (Lockwood et al., 2017).Despite allowing for B, ρ SW , V SW , and θ, the coupling function P α has only the one free fit parameter, the coupling exponent α that arises from an unknown dependence of the transfer function on the solar wind Mach number.This means that P α is much less prone to overfitting than functions that have separate exponents for the parameters.(Essentially, the exponents of B, ρ SW , V SW are related by the theory, and all are determined by just α).
The IMF orientation factor sin d (θ/2) was not treated as an independent variable by Vasyliunas et al. (1982).However, these authors did outline a test which was used to find that d = 2 was the required factor for the optimum (best-fit) α.The same test for other applications of the formulation by Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson, (2019) and Lockwood, Bentley, Owens, Barnard, Scott, Watt, et al. (2019) found a slightly different α (and that it varies with timescale) and this made d = 4 marginally better.Table 1 shows that sin d (θ/2) is a commonly used IMF orientation factor for low τ, particularly with d = 4.However, a range of d between 1 and 6 has been proposed in the literature.We here note that the test by Vasyliunas et al. (1982) has the very important implication that the optimum d is not independent of the other parameters in the coupling function.Vasyliunas et al. (1982) are somewhat uncertain as to whether they should employ the transverse component of the IMF, B ⊥ (the magnitude in the GSM YZ plane), or the full IMF magnitude B = (B X 2 + B ⊥

2
) 1/2 .They found it made only a minor difference in practice but opted to use B ⊥ in their text and equations.Their argument was that B X is not relevant because the field was draped over the nose in the magnetosheath.However, this choice is somewhat inconsistent theoretically because the IMF enters into their coupling function only through the Alfvén Mach number M A in the interplanetary (unshocked) field and that depends on B and not on B ⊥ .On the other hand, B ⊥ sin d (θ/2) is physically meaningful as a way of quantifying the southward component if the IMF in GSM coordinates.

Coupling Functions Based on Voltage Considerations
In addition to planetary geomagnetic activity, we here aim to predict transpolar voltage Φ PC for which we might expect a coupling function based on the IMF to be more appropriate.Many studies (e.g., Cowley, 1984;Reiff & Luhmann, 1986), suggest that the transpolar voltage Φ PC is well predicted by the dawn-to-dusk interplanetary electric field Because the voltage applied by the solar wind across the diameter of the magnetosphere is 2L o E SW , we can define the reconnection efficiency (the fraction of incident interplanetary field lines captured by magnetopause reconnection) η as We can then make the same assumption about the dayside magnetopause as was used to generate P α and again use pressure equilibrium with the solar wind dynamic pressure (Siscoe et al., 2002) where From Equations 5-7 we have a theoretical prediction of Φ PC , which we denote by Φ SW (the predicted value of Φ PC from solar wind parameters) Note that the reconnection efficiency η is very unlikely to be a constant.For example, increased solar wind dynamic pressure may increase the magnetic shear across the relevant current shear by compressing the field and various factors may vary the fraction of the dayside magnetopause covered by the magnetopause reconnection Borovsky and Birn (2014) argue that η is determined by the local Alfvén speeds on the two sides of the magnetopause to the extent that the interplanetary electric field is irrelevant.That being the case any similarity of an empirical coupling function to predict Φ PC and Equation 8 would be a coincidence.From reconnection rate theory and by making approximations Borovsky and Birn (2014) arrive at two distinct coupling functions for predicting dayside reconnection voltage here termed Φ BB .The sharp transition point between the two regimes where these apply is solar wind Alfvén Mach number, M A ≈ 6.For M A < 6 they find the approximate form B 0.51 N sw 0.24 sin 2 (θ/2) and for M A > 6 they find the approximate form B 1.38 N sw −0.19 V sw 0.62 sin 2 (θ/2).

Coupling Functions From Empirical Fits
Like many of the papers listed in Table 1, we here make empirical fits using a general multiplicative form of coupling function C f , given by This general form which can reproduce P α for any α (a = 2α, b = 2/3−α, and c = 7/3− 2α), E SW (for a = 1, b = 0 and c = 1),Φ SW (for a = 1, b = −1/6, and c = 2/3) as well as Φ BB (for M A < 6, a = 0.51, b = 0.24, and c = 1.49 and for M A > 6, a = 1.38, b = −0.19,and c = 0.62).As shown by Table 1, this form also encompasses a wide variety of the proposed empirical coupling functions.Note that this form could also reproduce the often-used "epsilon" factor, ε (for which a = 2, b = 0, and c = 1) but that is not considered further in this paper because ε is based on the incorrect assumption that the relevant energy flux in the solar wind is the Poynting flux (see Lockwood, 2013Lockwood, , 2019) ) and, although this can be made consistent with other energy coupling functions such as P α (i.e., correctly based on the dominant solar wind kinetic energy flux) this is only achieved using an extreme value of unity for the coupling exponent α, and this does not agree at all with experimental estimates.This theoretical flaw is the reason why ε performs considerably less well than P α on all averaging timescales (see Finch & Lockwood, 2007).
It should be noted that not all proposed coupling functions, not even all the simple ones, fit the general formulation given in Equation 9, particularly those that employ additive terms.An important example is that by Boyle et al. (1997) who propose the use of 10 −4 V SW 2 + 11.7 B sin 3 (θ/2) to predict Φ PC , which it does exceptionally well: the reasons for its success will be analyzed later in this paper.In general, the problem with additive terms is that, unless each term is describing a distinct physical mechanism, they are purely numerical fits to the available data.Adding terms until a fit is achieved without a theoretical basis does make the risk of overfitting greater: essentially one can fit any time series with combinations of other time series if one is free to select enough of them until a good fit is obtained.Physics-based coupling functions are usually fundamentally multiplicative in form although some factors can be broken down into the sums of additive terms for theoretical reasons (e.g., Borovsky, 2013;Lockwood, 2019;Newell et al., 2008).
The following sections describe how there are a number of procedural issues to resolve for studies using even the relatively simple form of coupling function generalized by Equation 9.For this reason, in the present paper we do not extend the present study to formulations involving additive terms.

Frequently Neglected Factors in Deriving Coupling Functions
There are a number of factors that have often been neglected when deriving coupling functions, the most important being: (a) the effect of data gaps; (b) the effects of data averaging; (c) the effect of the number of datapoints available; (d) the differences between the various terrestrial space weather indicators; (e) overfitting; (f) non-linearity and pre-conditioning of the magnetosphere; (g) measurement errors; (h) propagation lags; (i) spatial structure in interplanetary space; (j) effects of the bow shock; (k) seasonal effects; and (l) dipole tilt effects.We address just some of these in this article.The effect of data gaps was studied by Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson (2019) who introduced synthetic gaps at random (but to give the same distribution of durations as had occurred for early interplanetary observations) into continuous and near-continuous data and studied the errors introduced.These errors were not only in the greater uncertainty of one individual fit,

The Effect of Averaging Procedure
The magnetosphere responds to integrated forcing (Lockwood et al., 2016).For example, if we have a terrestrial indicator that responds to the energy input into the magnetosphere and a coupling function that quantifies that energy input, over a period τ we require the total of that energy input.Similarly, for a predictor of transpolar voltage we want the integral of the open flux generated in the interval.Hence for any empirical coupling function C f (Equation 9) we want the integrated solar wind forcing over the time interval.By the definition of the arithmetic mean, this means we need a coupling function for the interval τ given by where the values C f , B ⊥ , ρ SW , V SW , and θ are all values from high-time resolution measurements.However, this has usually in the past been approximated using the seemingly similar formula In addition, in many cases the average clock angle has been computed from the means of the IMF Y and Z components, that is, [θ] τ is used for θ and <B ⊥ > τ is replaced by [ B ⊥ ] τ , where and This generates major errors because both B Z and B Y can have both polarities within the averaging interval τ.We here use the notation that using means of the IMF components generates a coupling function [C f *] τ that has two separate problems.The problems caused by [θ] τ and [B ⊥ ] τ were addressed by the averaging procedure adopted by McPherron et al. (2015) who evaluated both at high time resolution before averaging and so avoided using either [θ] τ and [B ⊥ ] τ (this is hereafter referred to as the MEA15 procedure and is what we will use in later sections).In Figure 1 we highlight its importance but also deconvolve it from a second effect.Note that exactly the same operations are used in generating <C f > τ , [C f ] τ , and [C f *] τ and the difference between them is purely the order in which they are carried out: <C f > τ can be characterized as the parameters being "combined-then-averaged," whereas for [C f ] τ and [C f *] τ they are "averaged-then-combined." (The difference between [C f ] τ and [C f *] τ is that for the latter "averaged-then-combined" is even applied to the derivations of clock angle θ and transverse magnetic field, B ⊥ from the IMF B Z and B Y components).
Figure 1a demonstrates that it is far from a valid assumption to take <C f > τ and [C f *] τ to be the same, using the example of the Vasyliunas et al. (1982)   Figure 1a compares <P α > τ and [P α* ] τ and the linear correlation coefficient between the two is very poor indeed, being just 0.26.Note that in Figure 1a both <P α > τ and [P α* ] τ have been normalized by dividing by P o , the overall mean of P α : this has the advantage of canceling out all the constants in the theoretical derivation of P α .Rather than presenting scatter plots with massively overplotted points, Figure 1 employs "data density" plots (2-dimensional normalised histograms) with the fraction of samples, n/Σn, color-coded on a logarithmic scale with n being the number of sample pairs in small bins.In Figure 1a there are 100 bins of width 0.08 for both axes.Note that in all the data density plots presented in this paper, the lower limit of the color scale is chosen to be below the "onecount level", log 10 (1/Σn) as this ensures that even just one sample in a bin shows up as a blue pixel.Figure 1b identifies why the agreement in Figure 1a is so poor: it is for G, which is C f (in this case is P α ) without the IMF orientation term, that is, This is a factor that we will use again later in deriving optimum values for the exponent d.  13)-in other words we have moved to the MEA15 procedure in order to remove the IMF component-averaging effect on B ⊥ .The agreement is here is very good indeed, with values close to the diagonal line.
However, the agreement in Figure 1b is still not quite perfect.Small differences remain because of the difference between "Hölder means" (or a "power means") [<X p > τ ] 1/p of a general variable X and the corresponding arithmetic means <X> τ and hence between <X p > τ and <X> τ p . Figure 1b shows these differences are very small indeed for the variables X, the exponents p and the timescales τ involved in G in the example shown in Figure 1 and can be neglected.However, in general, arithmetic and Hölder means are related by what is called the "Hölder path" which results in the Hölder mean increasing with p (the arithmetic mean being the Hölder mean for the special case of p = 1).From comparison of Figures 1a and 1b, we know that the poor correlation in Figure 1a must be arising from the IMF orientation term, F(θ) = sin 4 (θ/2) and/or not using the MEA15 procedure to averaging of B ⊥ .Figure 1c compares the combine-then-average values of the clock angle θ, <θ> τ with the average-then-combine value [θ] τ , given by Equation 12, in the same format as Figure 1a (for bins of 2° × 2°) and although a great many points line up along the diagonal, there is considerable spread, especially at θ near zero or 180° (strongly northward and strongly southward IMF, respectively).Figure 1d makes the same comparison for the transverse field estimate, B ⊥ .Note that if we use the IMF magnitude B instead of B ⊥ in the coupling function, this effect does not arise; however, as found by Vasyliunas et al. (1982), tests show that using B ⊥ usually results in somewhat higher correlations.Figure 1e is for the same comparison for sin 4 (θ/2) and the spread is greatest at the southward IMF end of the range.
Figure 1f demonstrates that the MEA15 averaging essentially removes all problems associated with B ⊥ by avoiding [B ⊥ ] τ .However, Figure 1g shows that a problem still remains with the clock angle term sin 4 (θ/2).This is because Figure 1.Comparison of combine-then-average, average-then-combine and our compromise hybrid procedure for averaging 1-min data into 1-hr data (τ = 1 hr).In all panels, the horizontal axis gives the result of the combine-then-average approach which is what we ideally would wish to use to mimic solar wind forcing of the magnetosphere.The vertical axes in (a-e) give the result of an average-then-combine procedure.In each case the fraction of samples n/Σn is color-coded on a logarithmic scale, where n is the number of samples in small bins.The raw data used are 9,930,183 valid 1-min integrations of estimated power input to the magnetosphere, P α , and 11,646,678 valid 1-min values of the IMF clock angle θ and tangential component B ⊥ observed between 1995 and 2020 (inclusive).(a) The coupling function P α for α = 1/3 and d = 4 (the normalizing factor P o is the arithmetic mean of P α for all datapoints) in bins of P α /P o of size 0.08.The x-axis shows the means of 1-min values of P α , <P α > 1hr and the y axis the values [P α *] 1hr computed from 1-hr averages (including computation of the clock angle [θ] 1g and many of the points line up along the ideal diagonal: hence it is tempting to say this is just one more (small) source of noise and so it is valid to use <sin(θ/2)> d instead of <sin d (θ/2)>.However, there is a subtle point here: the spread shown in Figure 1g increases with d because the difference between arithmetic and power means increases with exponent.Hence using <sin(θ/2)> d discriminates against higher d by introducing more noise and so such studies will tend to derive a value for d that is too low, unless the larger-d fits are artificially enhanced by overfitting.
We can understand why the IMF orientation term is so different from the other three by looking the variability of the various factors within the averaging period.Figure 1 of Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson (2019) showed that the autocorrelation time of the IMF orientation is considerably shorter than for the other parameters and so most of the variability of P α on sub-hour timescales originates from the IMF orientation term.This is true for all coupling functions.If a parameter X is constant over the averaging time, then both the Hölder mean [<X p > τ ] 1/p and the arithmetic mean are equal to that constant value of X and <X p > τ = <X> τ p .On the other hand, if X varies a great deal during the averaging interval, then the Hölder mean is greater/smaller than the arithmetic mean for p greater/smaller than unity.Hence the much greater variability in the IMF orientation is the reason why the term in coupling functions accounting for its effect behaves so differently from the other terms.(However, note that if we increase the averaging timescale τ, the other parameters will also start to suffer from the same problem as the clock angle term.) We conclude that the often-used average-then-combine procedure generates large errors for the IMF orientation terms in deriving an empirical coupling function C f , even for τ = 1 hr.The MEA15 averaging procedure removes a great deal of the problem (at last for τ = 1 hr), but a second error (due to the difference between Hölder means and arithmetic means) remains for the clock angle term.This generates a problem when using an iterative procedure, such as the Nelder-Mead simplex search method used here (Lagarias et al., 1998;Nelder & Mead, 1965) to fit the exponents a, b, c, or d.This is because of the need to compute the mean of the combination of the samples (and in the dataset used in Figure 1 there are 9,930,183 valid 1-min samples of P α ) at the start of every round of the iteration.We have achieved this in some cases, but it takes enormous amounts of computer time and sometimes fails to converge.Fortunately, Figure 1 points to a compromise.It suggests we can use a hybrid approach of using <B ⊥ > a , <ρ SW > b , and <V SW > c , but must use <sin d (θ/2)> and not <sin(θ/2)> d for the IMF orientation term.This yields a mean coupling function estimate for averaging time τ of Figure 1h compares <P α > τ and [P′ α ] τ and it shows that agreement is very good with all points lying close to the diagonal line and the correlation coefficient is 0.997.We have repeated this test for all permutations of the maximum and minimum estimates of the exponents a, b, c, and d derived here and it is always valid to this level for τ = 1 hr Equation 15 is practical for use in an iterative fit procedure because for a given d we can compute <B ⊥ > τ , <ρ SW > τ ,<V SW > τ , and <sin d (θ/2)> τ just once before each iteration and then readily iterate a, b, and c to the optimum fit using the Nelder-Mead simplex search.This can then be repeated for different values of d.We have carried out some sample tests of our analysis that compared the results of fits using the ideal mean <C f > τ and our pragmatic hybrid solution [C′ f ] τ and the results were almost identical.However, we were limited in the number of these tests that we could carry out by the extremely large compute time caused by the need to average the whole dataset at each iteration step to define the exponents when using <C f > τ .We have repeated all calculations using the average-than-combine procedure [C f ] τ (but using the MEA11 procedure for B ⊥ and θ to avoid [θ] τ and [B ⊥ ] τ ) and, as described later, the fits obtained were always poorer because of the effect highlighted in Figure 1g.

Data Employed
We use the dataset of hourly mean transpolar voltage Φ PC observed over the years 1995-2020 (inclusive) by the northern-hemisphere SuperDARN array of coherent-scatter HF radars, as described by Lockwood and McWilliams (2021).These hourly data are means of 30, 2-min integrations.We adopt the requirement that the hourly mean of the number of radar echoes available, n e , exceeds a minimum value n lim = 255.This threshold was derived by Lockwood and McWilliams (2021) as the optimum compromise between having enough echoes that the influence of the model used in the "map-potential" data-assimilation technique is small, but not so large that the distribution of Φ PC values is greatly distorted by the loss of low-flow, low-n e samples.Lockwood and McWilliams (2021) also found that this threshold gave peak correlation between the radar Φ PC estimates and those from nearby passes of low altitude polar-orbiting spacecraft.The condition that n e > n lim = 255 yields a total of 65,133 Φ PC samples in the dataset.
We wish to compare the optimum coupling function for the global parameter Φ PC with that for global geomagnetic activity.We here use the am geomagnetic index (Mayaud, 1980) and the SuperMAG SML index (Newell & Gjerloev, 2011).Results are presented here for the SML index but were also generated for AL and the results were very similar.The am index has the most uniform network, in both hemispheres, of observing stations and uses weighting functions to yield the most uniform response possible to solar wind forcing with Universal Time and time of year (Lockwood, Chambodut, et al., 2019).It is based on the range of variation of the horizontal field component in 3-hr windows.To get a data series that is simultaneous with the Φ PC data, we here linearly interpolate the three-hourly am values to the midpoints of the hours used to generate the Φ PC data.This is only done for the Φ PC samples that meet the n e > n lim = 255 criterion and so we end up with a dataset of 65,133 interpolated am samples that are simultaneous with the Φ PC data.The advantage of using am is that it is the geomagnetic index that is by far the most free of seasonal and hemispheric effects which introduce noise in correlation studies, and it is genuinely global.The disadvantage is that it is three-hourly and the interpolated values will reflect this timescale.
We also compare with simultaneous hourly means of the SML index, by averaging 1-min values over the same hour as used to average the radar data.Note that SML comes from northern hemisphere stations and so contains an annual variation caused by seasonal changes in ionospheric conductivities: this is an additional noise factor for correlative studies that could potentially be reduced using a model of the effect of the conductivities.
Figure 2 compares these hourly datasets of Φ PC , am, and SML by presenting data density plots of the normalized geomagnetic indices (the am index in Figure 2a, am/<am> and the SML index in Figure 2b, SML/<SML> where the means are taken over the whole dataset) as a function of the simultaneous normalized transpolar voltage, Φ PC /<Φ PC >.In both cases, means of the normalized geomagnetic index (with error bars between the 1−σsigma points of the distribution) are also plotted for coarser bins of Φ PC /<Φ PC >. Figure 2a shows that the am index is, on average, close to proportional to Φ PC , but with considerable scatter.This proportionality of mid-latitude range indices and transpolar voltage, such as am and kp, has been discussed by Thomsen (2004).The variation of SML with Φ PC is a bit more complex with only a small increase at Φ PC /<Φ PC > below about 0.5 (i.e., Φ PC below about 20 kV), above which SML increases in magnitude more rapidly with Φ PC than does am.The scatter is somewhat higher for SML because, coming from one hemisphere, it contains noise associated with the seasonal variation in ionospheric conductivities.In contrast, am has very little such noise, being compiled from matching rings of stations in both hemispheres (and using weighting functions to account for any inhomogeneity) and has been shown to have an extremely flat response in both UT and time-of-year to solar wind forcing as a result (Lockwood, Chambodut, et al., 2019).
To derive the coupling functions, we use 1-min resolution averages of the Omni dataset of near-Earth measurements of interplanetary space made from craft in halo orbits around the L1 Lagrange point (King & Papitashvili, 2005).From this we generate running means using 1-hr (61-point) boxcar averages of B ⊥ , ρ SW , V SW , and sin d (θ/2) for the value of d we are investigating (the using the MEA15 averaging procedure).Mean values are only considered valid when the number of samples is large enough to make the error in the mean less than 5%, thresholds that were determined by Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson (2019) for each parameter by the random removal of 1-min samples from hourly intervals for which all 60 samples were available: because of its much lower auto-correlation function (acf), the most stringent requirement is set by the IMF orientation factor which requires 82% of samples (i.e., 43 out of the 60).The averaging generates a sequence of hourly running means that are 1 min apart.We combine these into mean coupling function [C′ f ] 1hr using our hybrid averaging formula (Equation 15).For test purposes only we also generate [C f ] 1hr using the average-then-combine procedure (Equation 11, with MEA15 averaging to generate hourly means of θ and B ⊥ ).We then select the value at each time of the transpolar voltage and am dataset, allowing for the appropriate propagation lag, δt p .
To determine the required propagation lags we make the initial assumption that the IMF orientation factor is sin 3 (θ/2) (i.e., d = 3), although this is refined in Section 3 of this article.We have carried out a sensitivity test to show that this choice does not influence the optimum derived lags.The Omni data have been propagated from the point of observation to the nose of the magnetosphere (King & Papitashvili, 2005): any variable error in that propagation will be a source of noise in our correlation studies.We then add a lag δt to allow for propagation across the magnetosheath to the dayside magnetopause and then to the relevant part of the ionosphere.We then vary δt between −60 min and +120 min and for each lag evaluate the linear correlation coefficients between Φ PC and am and the optimum coupling function, C f (for the assumed value for d of 3).Note that here and hereafter we refer to the hourly coupling function generated by our hybrid averaging procedure [C′ f ] 1hr as just C f , unless we are making a comparison with the results of the often-used average-then-combine procedure, in which case we distinguish between [C′ f ] 1hr , [C f ] 1hr , and [Cf*] 1hr .We want C f to be linearly related to the terrestrial activity indicator and so we maximize the linear correlation coefficient, r.(In later sections we will investigate the use of this metric).The exponents a, b, and c at each δt were determined using the Nelder-Mead simplex method to minimize (1−r) (Lagarias et al., 1998;Nelder & Mead, 1965).From this the optimum exponents a, b, and c (for the assumed d = 3) and the correlation coefficient r were determined at each lag δt.
The lag correlograms, r(δt) obtained this way are shown in the top panel of Figure 3: mauve is for Φ PC , the blue is for the interpolated am and the green is for SML.The vertical dashed lines mark the lags δt p giving peak correlation.The bottom panel shows the best-fit exponents a, b, and c as a function of lag δt: it can be seen that they do vary somewhat with δt but only to a small extent around the optimum lags.δt p .From Figure 3, we determine the optimum lags are δt p = 18.5 min for Φ PC , δt p = 31.0min for the interpolated am, and δt p = 45.0 min for SML.Note that the much greater persistence in the plot for am, because of it is interpolated from three-hourly data, and this makes the peak for am lower and broader.The survey of the Φ PC dataset by Lockwood and McWilliams (2021) demonstrates how Φ PC responds to both the reconnection rate at the dayside magnetopause Φ D and reconnection in the cross-tail current sheet tail Φ N (good proxies for which is the AL and SML auroral electrojet indices), as predicted by the ECPC model (Cowley & Lockwood, 1992;Lockwood, 1991).Indeed, in the approximation that the polar cap remains circular at all times, Φ PC is the average of Φ D and Φ N (Lockwood, 1991).Lockwood and McWilliams (2021) show that for low −AL, the lag of Φ PC after solar wind forcing is about 5 min, which is consistent with the expected response delay of Φ D , but the lag of the AL response (and hence inferred Φ N ) is 35 min, similar to the lag for am that is derived here.Hence we would expect the average lag for Φ PC , which is generated by a combination of Φ D and Φ N , to be around 20 min, as is indeed found to be the case in Figure 3.However, we note that there is considerable variability in the lags connected with Φ N , partly because of the variability in substorm growth phase duration (Freeman & Morley, 2004;Li et al., 2013) but also because, depending on the onset location, the precipitation in the initial part of the expansion phase can suppress ionospheric flow by enhancing conductivity, giving an addition delay in the appearance of the full voltage due to Φ N (Grocott et al., 2009).2, along with their estimated uncertainties.The gray areas in Figure 3 define the 1σ, 2σ, and 3σ uncertainties in the δt p estimates.These are evaluated by looking at the significance S of the difference between the correlation at a general lag r(δt) and its peak value at the optimum lag δt p (where r = r p ) where S = 1−p, and p is the probability of the null hypothesis that r and r p are actually the same.S is computed using the Meng-Z test (Meng et al., 1992) for the significance of the difference between correlation r AB (between two variables A and B) and r AC (between A and C) allowing for the fact that B and C may be correlated (|r BC | ≠ 0).S is, by definition, zero at the optimum lag δt p , and the 1σ, 2σ, and 3σ uncertainties are the lags at which S has risen to 0.68, 0.95, and 0.997, respectively.For Φ PC the 2σ uncertainty band is ±1.3 min; for am it is ±4.0 min; and for SML it is ±6.1 min.Note that these uncertainties are smaller than in many studies because the number of samples is so large.Because Figure 3 was generated using an assumed value of d = 3, it was repeated for a range of selected values of d between 1 and 7 (which Section 3.2 shows covers the range of interest), the differences between the derived optimum lags were always considerably smaller than the above 2σ uncertainties.
For completeness and to enable reproducibility, the precise and full details of the procedure for determining exponents and their uncertainties is defined in advance here for the case of Φ PC (and is implemented in subsequent sections for SML and am as well as for Φ PC ).The optimum lag δt between Φ PC and C f of 18.5 min derived in Figure 3 is employed.The value of the exponents are then found by varying d between 1 and 7 (a range larger than the largest suggested in the literature) in steps of 0.001 (a value found by iteration to be about 1% of the uncertainty in d, δd).For each d, the values of a, b, and c that give the peak correlation r between Φ PC and C f are found using the Nelder-Mead simplex search algorithm (specifically the algorithm minimizes 1−r).The optimum value of d is then determined by searching for the value for which the combination of best-fit exponents a, b, c with d that give a C f that has a fully linear relationship with Φ PC .This is done using a quadratic polynomial fit of suitable functions of C f and Φ PC , as given by Vasyluinas et al. (1982) who devised the method.The polynomial uses least-squares fitting and the uncertainties in the derived coefficients are computed using the standard technique employing the reduced chi-square statistic.Linearity is achieved when the coefficient for the squared term in the quadratic fit falls to zero and this defines the optimum value of d.The distribution of p-values of d (and hence the uncertainty δd at a given probability) is found from the p-values of the squared term coefficient when it equals zero.The optimum d defines the optimum values of the other three exponents, a, b, and c (i.e., that give the maximum correlation between C f and Φ PC for that d).To determine the uncertainty in, for example, the exponent a For interpolated one hourly data.b For all three-hourly data.c For simultaneous one-hourly data.d For all one-hourly data.

Table 2
The Best Fit Exponents a, b, c, and d and the Resulting Peak Correlation Coefficient r p for the Terrestrial Parameters Φ PC , am, and SML From Fits Using the Data From the Range of Dates Given in the Footnotes a, we vary it around its optimum value and study how the optimum correlation coefficient varies (found using the Nelder-Mead method by fitting b and c for each a using the optimum d and δt).As before, we use the Meng-Z test (Meng et al., 1992) for the significance of the difference between two correlations to find the values of a for which the correlation is lower than the optimum value at given p-value of the null hypothesis that the correlation is the same as the peak.Finding the values of a that give these p-values of 68%, 95%, and 97.7% gives us the uncertainty bands in a at the 1σ, 2σ, and 3σ levels.The same procedure was then repeated for exponents b and c.

The IMF Orientation Factor
As discussed by Vasyliunas et al. (1982), the optimum IMF orientation factor is not independent of the other fit exponents.In addition, Section 1.4 has described how, because its much greater rapid variability, we have to deal with it differently when generating average coupling functions.Section 3.1 discusses the distributions of IMF orientation factors before in Section 3.2 we evaluate the optimum values of d for Φ PC , am, and SML.

Occurrence Distributions of IMF Orientation Factors and the Effect of Averaging Timescale
Figure 4 shows the distributions of various parameters relevant to the IMF orientation factor, all panels being for 1-min integrations of data and in the GSM frame of reference.For panels (e-j) the same vertical scale is used but an increasingly larger fraction of samples occur in bin 1 and exceed this scale: hence the value of n/Σn for bin one is given by the legend for panels (f-j).
the various coupling functions come about from a near-Gaussian distribution of the IMF B Z component, which is very close to symmetric around zero, and a double-peaked distribution of the IMF B Y component, which is also very close to symmetric around zero.As discussed above, the most commonly adopted form of the IMF orientation factor has been sin d (θ/2) with d = 4 although a range of d from 1 to 6 has been proposed.Figure 4f shows that d = 2 yields a symmetric distribution around an average of 0.5 with dominant isolated peaks in the bins closest to 0 and 1.On the other hand, Figure 4g shows that d = 4 yields a highly asymmetric distribution with an even-larger isolated peak in the bin nearest 0 and only a very small one in the bin nearest 1.The peak in the lowest bin is even larger for d = 6, shown in Figure 4h and larger again for two other commonly used "half-wave rectified" IMF orientation factors B S in Figure 4i (where B S = −B Z for B Z < 0 and B S = 0 for B Z ≥ 0) and U(θ)cos(θ) in Figure 4j (where U(θ) = 0 for θ < 90° and U(θ) = −1 for θ ≥ 90°).The distributions for B S and U(θ)cos(θ) are very similar because U(θ)cos(θ) = B S /B and the factor 4.5 is used to display B S on the same scale because it makes the mean value the same as for U(θ)cos(θ) and very similar to that for sin 6 (θ/2).
These strange distributions of IMF orientation factors have great significance for statistical studies of the performance of a proposed coupling function because they determine the weighting given to a given clock angle θ in a correlation study.This means that when we alter d, we are not just investigating the how the IMF orientation influences solar wind-magnetosphere coupling, we are also changing the statistical weighting given to certain IMF orientations in our correlation studies.For B S and U(θ)cos(θ) the value is zero for 50% of the dataset (for B Z > 0) and so the coupling function is more weighted to accurate prediction of quiet times, which is probably not what is wanted in many applications.Figure 4h shows the distribution is not quite so extreme for sin 6 (θ/2), but it has the same basic form.As we reduce d, that weighting shifts until for d = 2 the distribution is dominated by two equal peaks close to due northward and close to due southward IMF.For d = 1 (Figure 4e) it is dominated by close to purely southward IMF.The key point is that the choice of the IMF orientation factor is also setting the weighting given to certain data in the statistical fit of the coupling function if we use an overall fit-quality metric such as correlation coefficient or root-mean-square deviation.
Figure 1 of Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson (2019) shows why the IMF orientation factor has a key role in setting the variability of a coupling function.It is because its autocorrelation function (acf) falls much more rapidly with time lag for any other solar wind parameter.For a lag of 1 hr, the autocorrelation function for sin 4 (θ/2) in near-Earth space is 0.45, whereas for the solar wind number density N SW it is 0.88, for the IMF B it is 0.93, and for the solar wind speed V SW is 0.99.Hence short-term variability of a coupling function is set by that in the IMF orientation factor whereas, as shown below, this factor essentially becomes constant at timescales of a year or more.This exemplifies the general fact that the IMF orientation factor distribution depends critically on averaging timescale which is here illustrated by Figure 5 for the commonly adopted sin 4 (θ/2) factor.We take running boxcar (running) means of the 1-min data over intervals τ and deal with data gaps by only retaining averages that are made up of a fraction of the potential maximum number samples that exceeds f(τ), the minimum needed to keep errors due to data gaps below 5%.The minimum fractions f(τ) needed were computed by introducing random synthetic data gaps into continuous IMF data, computing the error caused and repeating 10 times for each hourly mean, as carried out for τ = 1 hr by Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson (2019).For example, Figure 1b of Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson (2019) shows that we require f(τ) > 0.82 to keep errors in the hourly mean IMF orientation factor to below 5%.As τ increases, it becomes increasingly hard to find intervals with no data gaps and hence to compute the threshold; however, at the same time f(τ) falls with increased τ and for τ > 1 day we set the f(τ) requirement threshold value to be the same as for 1 day.
As τ is increased, the central limit theorem (Fischer, 2010) applies and the distribution of any parameter narrows toward a delta function at the overall mean (i.e., the value derived for a τ equal to the duration of the whole dataset).However, because of the unusual form of the distribution at τ = 1 min, the distribution for sin 4 (θ/2) evolves through a series of forms and how it does so is determined by the timescales of the variability in the IMF orientation.For τ = 15 min the distribution is quite similar to that for τ = 1 min, but the peak at sin 4 (θ/2) = 0 has diminished and more samples occur at larger values.For τ = 1 hr (the timescale used in this article), this results in a near-linear distribution, but still with a pronounced peak at 0. By τ = 6 hr the distribution has evolved into very close to a lognormal form and by τ = 1 day it is close to a Gaussian form that is symmetrical about the overall mean value (the mauve vertical dashed line).Further increases in τ cause the width of the distribution about the overall mean to decrease.For τ = 1 yr, the distribution is narrow and hence the IMF orientation factor can, 10.1029/2021JA029946 19 of 41 to within a reasonably small error, be taken to be constant.This is why successful coupling functions at annual timescales usually do not contain a factor that allows for IMF orientation.Note that all parameters in a coupling function, not just the IMF orientation, follow the central limit theorem, but the other factors tend to start (for 1-min observations) from a log-normal form and then evolve into the narrowing Gaussian and do not start from the unusual distributions for the IMF orientation factors (Lockwood et al., 1999).
The averaging timescale τ has significance on two levels.Here, we study it purely in the context of averaging data and the changes of the distribution that are associated with the reduction in noise brought about by the averaging.However, it should be noted that τ also has significance on a physical level.This is because the IMF orientation in the upstream solar wind will be influenced by the passage of the solar wind through the bow shock and magnetosheath and there will be time constants for changes in the coupling of energy, mass and momentum from the near-magnetopause sheath into the magnetosphere (e.g., changes in the reconnection rate and in the X-line latitude and orientation).These will almost certainly act as a low-pass filter on the IMF orientation variations, but it is not yet clear what averaging timescale τ will best mimic the effects of this low-pass filter and how it might vary with solar wind conditions.The optimum τ will depend on the terrestrial parameter considered.For example, studies using ground-based radars show rapid responses in ionospheric flows and the location of the inferred open/closed boundary in the cusp region (almost immediately after the arrival of the Alfvén wave down the field line from the magnetopause to the ionosphere).However, flows over the polar cap (quantified by the transpolar voltage) evolve more slowly and do not fully respond until 15-20 min later (Lockwood & McWilliams, 2021), consistent with the Expanding-Contracting Polar Cap model (Morley & Lockwood, 2005, 2006)-although we note that quasi instantaneous responses are also possible if the magnetosphere has been pre-conditioned by prior (g) a solar rotation period of 27 days, and (h) 1 yr.The numbers of samples, n, as a fraction of the total number Σn, in bins 0.01 wide are shown in each case and the dataset used is the same as in Figure 4.The vertical mauve dashed lines are for the overall average of all samples.The vertical green line is at θ = 90° for which the IMF lies the GSM equatorial plane.Note.That the lowest bin in sin 4 (θ/2), which is 0-0.01,corresponds to a range in θ of 0°-36.9°whereas the highest bin (0.99-1) corresponds to 171.9°-180°.
10.1029/2021JA029946 20 of 41 magnetopause reconnection (Morley & Lockwood, 2005).Hence determining the timescale that is relevant to a given response is a multi-faceted and complex problem.
Figure 6 is the same as in Figure 5, but for another value for d that has been proposed in the literature, namely d = 2 (e.g., Borovsky, 2013;Kan & Lee, 1979;Lyatsky et al., 2007).This reveals the sin 2 (θ/2) has very different behavior to sin 4 (θ/2).At all τ, the distribution is symmetric about 0.5 and the mean value (vertical dashed line) and the value for in-equatorial field (vertical green line) are both always at 0.5.For τ up to about 15 min, this yields a uniform distribution with sin 2 (θ/2) with just small peaks at zero and unity that decay as τ is increased.This even distribution makes sin 2 (θ/2) an attractive choice if studying timescales up to about 15 min.However, for τ = 1 hr and above the distribution takes on some undesirable characteristics, with most samples coming from near-in-equatorial field and fewer from the extremes near 0 and 1.As discussed below this has some consequences In the literature values for d between 1 (Borovsky, 2008;Fedder et al., 1991) and 6 (Balikhin et al., 2010;Temerin & Li, 2006) have been proposed and used.From the above, the choice of IMF orientation factor and of the averaging timescale both have a subtle effect on the coupling function fitting by changing the weighting given to the data samples.The central limit theorem means that the same effect applies to other factors in the coupling function, but the effects are less marked because they do not start from as extreme a distribution for 1-min values as does the IMF orientation factor.One key insight here is that we should not expect a coupling function that works well at one timescale to be equally effective at another.Hence some of the differences between the coupling functions proposed in Table 1 will have arisen from the different averaging timescales used.
The behavior in Figures 5 and 6 is very different to that obtained by an average-then-combine procedure given by Equation 12(not shown).In these cases, the distribution tends to maintain its high-resolution form up to τ of about 1 day when it starts to narrow under the central limit theorem.However, as τ is further increased it gets noisy and broadens again as the means of both the Y and Z components of the IMF tend to zero.The key point is that this behavior is purely an artifact of the average-then-combine procedure, and the combine-then average is what mimics the physics of the magnetosphere.The distributions of the other parameters in the coupling function are largely log-normal and also influence the net distribution of C f , but it the IMF orientation factor that has the most marked effect, and the imprint of its strange distributions is clearly seen in C f (Lockwood, Bentley, Owens, Barnard, Scott, Watt, et al., 2019;Lockwood, Chambodut, et al., 2019).

Optimum Exponent d of the IMF Orientation Factor
In Section 2 we defined the optimum lags for the interplanetary data, δt p , and found that they were not significantly influenced by the choice of the exponent d in the sin d (θ/2) IMF orientation factor.In this section, we define the optimum d using those lags.We vary d over the full proposed range (we used values from 1 to 7.5 in steps of 0.001) and using the optimum lags δt p , we optimized a, b, and c to maximize the correlation r at each d.The results are shown for Figure 7, using the same format as Figure 3.
The top panel of Figure 7 shows that for Φ PC , am, and SML, the correlation has a peak at quite low d, specifically d = 2.1 for Φ PC (in mauve) and d = 1.3 for am (in blue) whereas for SML (in green) the peak correlation is at d = 3.7, very close to the value found by MEA17 for AL.The bottom panel shows how the other exponents (a, b, and c) depend slightly on d.Note that we have also used the MEA15 averaging methods to generate hourly coupling functions C f , [C f ] 1hr using Equation 11(not shown): as expected from Figure 1g, the correlations for [C f ] 1hr were systematically lower than for [C′ f ] 1hr .For a few sample values of d (specifically 2, 3, 4, and 6) we also repeated the computation using <C f > 1hr (Equation 10): in each case, iteration took over a 1,000 times longer than the corresponding fit using [C′ f ] 1hr , but the results for a, b, c, and r were all the same for <C f > 1hr and [C′ f ] 1hr to within the estimated uncertainties.From Figure 7a, it appears that the sin 2 (θ/2) IMF orientation factor performs best for Φ PC and that an even lower d is best for am because they yield higher correlation coefficients.
However, as discussed in the previous section, some of this is the favorable distribution of samples that averaging brings about and the subsequent weighting of IMF orientations in deriving the correlation coefficient.This is demonstrated by Figure 8 for fits to the Φ PC data.Figures 8a and 8b show that for a d value that is too low or too high the relationship between C f and Φ PC is not linear (with curvature in the opposite sense in the two cases).Figure 8c is for the peak correlation (d = 2.2) and it can be seen that the variation is not linear, but d is slightly too small, giving the same curvature as seen in Figure 8a. Figure 8d shows that it requires a slightly larger d (=2.5) to give a linear variation, even though the correlation is slightly lower and the rms deviation is slightly larger than for d = 2.2 that yields peak correlation.The reason lies in the effect of the distribution of C f values on the fits.The color contours reflect the point made in relation to Figure 4, namely that higher d causes a greater density of points at low C f and so biases the fits to lower values of Φ PC and hence northward IMF.This can be seen by comparing the color contours in the Figures 8a-8d.
An interesting point to note is that the variation in Figure 8c could be interpreted as a saturation effect at work, whereas it is in reality the application of a value of d that is too high.Saturation is identified when the observed Φ PC is not as high as we would expect for a given coupling function for the prevailing interplanetary conditions (Hairston et al., 2005;Shepherd, 2007).Such an empirical identification and quantification of a saturation affect assumes that the coupling function had been made to have a linear variation with Φ PC and Figure 8 demonstrates that deriving the coupling function using correlation coefficient can give a non-linear variation of C f with Φ PC .It seems likely that saturation is a real phenomenon-for example, it is generated by MHD simulations Kubota et al. (2017) and we note that saturation the maximum Φ PC /<Φ PC > in Figure 8 is near 2.7 which corresponds to 100 kV (<Φ PC > = 37 kV) and saturation has generally been reported at larger Φ PC , typically 150-200 kV and certainly at a level above 100 kV.In addition, the curvature caused by excessively large d extends throughout

10.1029/2021JA029946
22 of 41 all values of Φ PC -unlike saturation effects.But we conclude most of the data in Figure 8 are not influenced by saturation.Furthermore, the variation that looks like saturation in Figure 8d is generated by an exceptionally large d (=6.5) whereas the effect of statistical weighting is to tend to underestimate d when using correlation.However, we must remain aware non-linearity introduced into the coupling function, caused by statistical biasing toward certain IMF clock angles, can cause us to underestimate or overestimate the true saturation effect.
There is second way to derive d that avoids the possibility of statistical bias, and this is presented in the Section 3.3.Vasyliunas et al. (1982) provide a test for the optimum form of the IMF orientation factor F(θ), such as sin d (θ/2).This is based on the fact that we want the coupling function C f to be linearly related to the terrestrial response at all activity levels and not be biased in the way illustrated by Figure 8.To evaluate this, we use the function G (i.e., C f without the F(θ) factor, defined by Equation 14).We want C f to vary linearly with the terrestrial index T (which is either Φ PC , SML, or am in the current article).Hence we want

Test of the IMF Orientation Factor and Linear Regression Coefficients
where s T and i T are the best-fit linear regression coefficients.This yields a requirement that   10c, and 11a-11c are the corresponding plots of <(am−i am )/G> and of <(SML−i SML )/ G>, respectively, as a function of <F(θ)>.In all cases we use the derived optimum G for the value of d in question (i.e., using the coefficients a, b, and c given in Figure 7b).Averaging is carried out over 25 bins of F(θ) of width 0.04, covering the full range of 0-1.In Figures 9-11, the Parts (a), (b) and (c) are for, respectively, d below, equal to, and above the optimum value which is derived below: they show that the best fit quadratic polynomial (the red line) and this is not linear in the Parts (a) or (c) (the green line gives the best linear regression which will be the same as the red line for a linear dependence).Figures 9a, 10a, and 11a, the coefficient of the power-2 term in the Under the mauve line in three shades of gray are the 1σ, 2σ, and 3σ uncertainty bands in a Φ , the limits to which define the corresponding uncertainty bands in the optimum d, giving a 2σ uncertainty in the optimum d of ±0.08.Note.That in this case for Φ PC the differences between the uncertainty bands are often so small that they cannot be discerned; they are more clearly seen in Figure 10 for am and Figure 11 for SML.(b) Confirms this proportional relation at this optimum d = 2.50 for which the exponents are given in Table 2.
best fit quadratic polynomial is positive, whereas for the Figures 9c, 10c, and 11c it is negative-that is, the curvature of the best fit of the polynomial is in the opposite sense to in the corresponding Figures 9a, 10a, and 11a.
For the Figures 9b, 10b, and 11b, the fit is linear, and this is what makes the d used in these cases the optimum value as it means the coupling function is linearly related to the terrestrial index.
The derivation of the optimum value of d is shown in the Figures 9d, 10d, and 11d which plot the power-2 term coefficient in the best fit-quadratic (a Φ for Φ PC , a am for am, and a SML for SML) as a function of the exponent d over the full range of values proposed in the literature.The uncertainty band of this coefficient, at the 1σ, 2σ, and 3σ levels (derived for the polynomial fit using the reduced chi-square statistic), are shown in shades of gray in all Figures 9-11 (but more easily discerned in Figures 10 and 11).The optimum d for Φ PC , am, and SML are the values that make, a Φ , a am , and a SML (respectively) equal to zero-that is, for which the variation is linear.The 1σ, 2σ, and 3σ uncertainties in d are where the edges of the uncertainty bands in a Φ , a am , and a SML fall to zero and this yields the vertical uncertainty bands around the optimum d that are shown.
Figure 9 shows that the required d is 2.50 ± 0.08 (at the 2σ uncertainty level) for Φ PC , Figure 10 shows that it is 3.97 ± 0.18 for am and Figure 11 shows that it is 5.20 ± 0.41 for SML.Hence the optimum IMF orientation factors for Φ PC , am, and SML are not the same within 2σ uncertainties and in all three cases are larger than the value derived by correlation.Essentially SML requires a function, that is, most like a half-wave rectified function and Φ PC requires a function that is least like one.The optimum d and their uncertainty bands for Φ PC , am, and SML are also shown in Figure 7 which reveals that the uncertainties do not overlap even at the 3σ uncertainty level.
Note that the commonly used value of d = 4 is too large for Φ PC and am but too small for SML.Some agreement between the behavior of am and SML is to be expected because both are dominated, at high activity at least, by the effect of the substorm current wedge and so do show considerable agreement (Adebesin, 2016; supplementary information to Lockwood, Bentley, Owens, Barnard, Scott, Watt, & Allanson, 2019).However, they are different indices and, as indicated by Figure 2  The question then arises as to why the correlations r at these optimum d are slightly lower than the peak correlations that are always found at slightly lower d, as can be seen in Figure 7a.The answer can be found by referring back to the analysis of the d = 2 case and the F(θ) = sin 2 (θ/2) factor presented in Figure 5.This series of distributions shows that the dataset becomes weighted toward the middle of the range of sin 2 (θ/2) values as the timescale is increased and there are fewer data constraining the large and low values.This is clearly demonstrated by the distribution for these data with τ = 1 hr in Figure 5c.Hence although sin 2 (θ/2) gives very slightly higher r p , it is only because the dataset becomes weighted toward the center of the distribution with weaker weighting given to the extremes of low and high F(θ).To test this conclusion, we carried out correlations where the data were divided into 25 bins of F(θ) and for each bin, samples were selected at random such that all the F(θ) bins contained the same number of samples (the number that were in the least-populated bin), thereby removing the sampling bias at the expense of losing data.The peak correlations were indeed shifted to larger d and closely matched the values derived in Figure 7.These correlation tests are still not bias-free because reducing the samples to the minimum number in any one bin means that fits for some d have systematically higher sample numbers than others.Nevertheless, this test is enough to confirm that the choice of d does influence the correlation coefficients by preferentially weighting certain clock angles.
In contrast, in fitting the quadratic polynomial to the bins in Figures 9a-9c, 10a-10c, and 11a-11c, equal weight is given to the data points for the different F(θ) bins, despite the fact that there are different numbers of samples in those bins.Hence, unlike the correlation coefficient r, these fits are not influenced by the distribution of samples.Hence, they provide a better test of the optimum form of F(θ) that best describes the physics of solar-wind magnetosphere coupling than do the correlation coefficients.
It can be seen from the bottom panel of Figure 7 that, in general, the uncertainty in d introduces only small changes in the best-fit exponents a, b, and c.However, the changes across the uncertainty bands are not zero.Hence when we compute the uncertainties in a, b, and c we need to fold in both the fit uncertainties at the optimum d and effect of the uncertainty in that optimum d.With all four exponents and the linear regression coefficients now defined, the predicted terrestrial index can then be determined from:

First-Order Check for Overfitting
We here fit with three free fit parameters (a, b, and c), and are pre-determining two others (d and the optimum lag, δt p ) which can influence the results and hence, even for such a large dataset, overfitting could be a problem.An initial test is to check that correlations are not unrealistically high.We carried out tests for the effect of the noise introduced into our correlations by the use of interplanetary data from spacecraft in a halo orbit around the L1 Lagrange point: the point being that the solar wind that is sampled by the spacecraft is not, in general, the same as hits Earth because of spatial structure in the interplanetary medium.We computed our generalized coupling function, covering the full range of a, b, c, and d indicated by Figure 7b, using data from both the ACE in a "halo" orbit around the L1 point and the THEMIS-B spacecraft for 2010-2019 (inclusive) when the latter was outside the bow shock in the near-Earth solar wind.THEMIS stands for "Time History of Events and Macroscale Interactions during Substorms" and for the time interval studied the THEMIS-B craft was in Geocentric orbits and between about 55 R E and 65 R E from Earth (where R E is a mean Earth radius, 6,370 km) which resulted in it being in the undisturbed solar wind ∼70% of the time, in the shocked solar wind of the magnetosheath about 15% of the time and inside the magnetosphere for the remaining 15%.For both ACE and THEMIS-B, several sample coupling functions for d = [2:1:6] were computed at 1-min resolution and then averaged with a 60-point running mean into hourly values with 1-min cadence.The coupling functions from the two craft were correlated over intervals 2 days long, giving 48 fully independent samples in each period.The optimum lag in each case was determined as a function of time and the peak correlation evaluated from the lag correlograms.The results did vary a little with the coupling function exponents used and, in particular, correlations were lower for higher d, indicating that spatial structure in the IMF orientation was one of the larger causes of noise.For the 70% of time when THEMIS-B was in the solar wind the distribution of correlation coefficients between the hourly means of the same coupling function at the two craft typically had a mode value of 0.96 and a mean value of 0.84.As expected, the lower correlation coefficients were all associated with lower correlations of the IMF orientation term.It was found craft locations made little difference with only a very small reduction in correlation seen when the L1 craft was at its greatest distance from the L1 point.However, the mean correlation was increased to 0.92 for data when the am planetary index exceeded 70 nT, indicating spatial structure has less effect during space weather events, presumably as they are driven by large interplanetary structures such as Coronal Mass Ejections and Corotating Interaction regions.The same study was repeated for the 15% of time that THEMIS-B was in the magnetosheath.Correlation coefficients were considerably lower, with a mode value of 0.81 and a mean value of 0.75.This agrees well with studies by Šafránková et al. (2009), Borovsky andYakymenko (2017), andWalsh et al. (2019) who show that passage through the bow shock is a major source of noise for coupling functions measured in the undisturbed solar wind.From these and our own studies correlations above 0.9 for hourly data in the undisturbed solar wind are a warning that the data have most likely been overfitted.Note also we have only considered one potential source of noise and we should regard 0.8-0.9 as about the best correlation that we can achieve for hourly means using upstream data from the undisturbed solar wind.Note that correlations of up to about 0.98 are possible using annual means of data because so much of the noise is canceled out at this timescale.
When correlating with terrestrial indices such as AL, SML, am, and the transpolar voltage Φ PC , we were surprised to find that correlations were systematically slightly higher for L1 craft than for THEMIS-B as we expected the latter to have less spatial structure error: this hints that instrument calibration issues generate larger differences than spacecraft location errors.
We here also test for overfitting in a straightforward way by dividing the data into just two "folds" (while noting that machine-learning techniques often use several more folds for different tasks) of roughly equal numbers of samples and then fitting to the one half and the testing against the independent second half.Note also that testing also raises another set of complications with a variety of performance metrics available for consideration (Liemohn et al., 2018), and the most appropriate one (or ones) for the application in question should be deployed, especially in the context of forecasting (Owens, 2018).
We here use the optimum lags δt p and d exponents derived above and consider only linear correlation coefficient and root mean square (rms) error as test metrics.The results are demonstrated in Figures 12 and 13.The fit dataset used to define exponents a, b, c (for the predetermined d for the parameter in question) was for 2012-2019, inclusive and the resulting values are given in the legend to Figure 12.The same exponents and regression coefficients were then applied to generate the predicted values for both the fit and the test subsets (1995-2011) using Equation 18.Because there are so many datapoints, information is lost in a scatter plot because so many points are overplotted: Figures 12 and 13 are therefore presented as datapoint density plots.Comparing Figures 12 and 13 there are no obvious differences in behavior, which is quantified by the correlation coefficients r and the rms deviations Δ between observed and predicted values.For the predicted and observed Φ PC , r is 0.853 and 0.886 for the fit and test sets, respectively, and Δ is 12.9 and 10.4 kV.Hence, by both metrics, the test set is actually performing slightly better than the fit set.For the predicted and observed am, r is 0.813 and 0.822 for the fit and test sets, respectively, and Δ is 10.1 and 10.7 nT.Hence in this case the correlation is very slightly better for the test set, but the rms deviation is slightly better for the fit set.For the predicted and observed SML, r is 0.808 and  LOCKWOOD AND MCWILLIAMS 10.1029/2021JA029946 28 of 41 0.764 for the fit and test sets, respectively, and Δ is 84.4 and 83.8 nT.Hence in this case the situation is the opposite to that for am, but differences are again very small.In all cases, the performance of the fits on the test set is essentially the same as for the fitting set and there is no doubt that the coupling functions have predictive power.
Note from the plots presented in Figures 12 and 13 the influence that the d value has on where data are in parameter space.For Φ PC (which requires d = 2.5) there is a high density of samples over a large segment of the best-fit diagonal line.For am (which requires a higher d = 3.0) the highest density of data is more closely confined to near the origin and this effect is even more marked for SML (which requires a yet higher d = 5.23).The key point is that the influence of northward IMF conditions on the derived general coupling function is greater for SML than it is for am and Φ PC which needs to be remembered when we evaluate its performance.

Estimation of Uncertainties and the Influence of the Number of Samples
Figure 14 presents distributions of fitted values of the exponents a, b, and c for three subsets of the transpolar voltage data and compares them to the value for the full set of N = 65,133 samples (given by the vertical dashed line in each case).The distributions are generated by taking 1,000 random selections of   samples (from the total of N T = 65,133 samples with n e > n min = 255 available): the values of N used were N T /25 = 2,606 (on average, corresponding to 1 yr of data); N T /10 = 6,513 (on average, corresponding to 2.5 yr of data) and N T /2.5 = 26,503 (on average, corresponding to 10 yr of data).The fraction of samples n/Σn are plotted in bins of width (1/30) of the maximum range of each exponent shown.In each case, three histograms are shown: the light gray histogram bounded by the mauve line is for N T /25 samples, the darker gray bounded by the blue line is for N T /10 and the darkest gray bounded by the black line is for N T /2.5.The distributions are generally symmetric about the optimum value for the whole dataset, but not always so for the smallest N and, as expected, they narrow down toward the value for the full dataset as N is increased.The standard deviations of the distributions are given in each case on the plot.This analysis was repeated for the geomagnetic indices and the results were very similar (not shown).Distributions are broader and peaks lower for am and SML than for Φ PC , which is expected because all plots presented thus far have had greater noise and larger uncertainties for the fits to the geomagnetic data.Figure 14 stresses how much in error an individual fitted value can be if smaller datasets are used.However, that both the mean and the mode of some of the distributions are shifted from the value for the whole dataset when N is low, meaning that there are systematic errors as well as random errors when sample numbers are low.
To determine the uncertainties in exponents a, b, and c from our full dataset we assigned one of the three exponents a fixed value that was then varied round its optimum value and the other two were fitted using the same Nelder-Mead simplex search procedure that was used to fit all three exponents in previous plots (again, we are using the optimum d and lag δt p defined previously).The significances S of the difference between the correlation at a general value of the exponent and its peak value for the optimum exponent was then evaluated.As before, we evaluate S = 1−p (where p is the probability of the null hypothesis that the correlations are the same) using the Meng-Z test and find the exponents value at the 1σ, 2σ, and 3σ uncertainty levels.This yields the uncertainty associated with the fit at the optimum d, which was added in quadrature with the uncertainty caused by the uncertainty range in that optimum d.The resulting 2σ uncertainties are given with the optimum values in Table 2.

Significance of the Differences Between Fits for Transpolar Voltage and Geomagnetic Activity
A notable feature established earlier is that the optimum d for Φ PC , am, and SML are not the same: the shaded areas of Figure 7 show that the uncertainties do not overlap for even the 3σ level.Form Table 2 we can see that the exponents a, b, and c (of B, ρ SW , and V SW , respectively), are also, in general, different.We conclude that there is no such thing as a universal coupling function and optimum coupling functions must be tailored to the index or indicator that they aim to predict.We have carried out a number of experiments of the kind illustrated in Figure 14 using randomly sampled subsets of the data and found that some exponents that appeared to be the same, within predicted uncertainties, are found to be different, to very high significance, when we use the full dataset.

Matching the Distributions of Coupling Functions and Observed Parameters
Thus far in this article, we have noted that the optimum coupling function depends upon averaging timescale and on which terrestrial activity indicator is predicted.In this section, we make the point that it also depends on the application for which it is intended to be used.If we require a coupling function to predict the largest events it is most important to match the large-value tail of the distribution of activity indicators.However, we may require a climatology that also predicts integrated values accurately.Examples of such applications are in integrated radiation doses for humans and spacecraft electronics and transformer degradation levels cause by geomagnetically induced currents.In such cases we need the coupling function to match the full distribution of values and not just the large event tail.This point is here illustrated by Figure 15 which compares the distributions of three empirical coupling functions to the observed distribution of valid hourly transpolar estimates from the survey of Lockwood and McWilliams (2021).The observed distribution of normalized transpolar voltage, Φ PC /<Φ PC > is given by the black line that bounds the area shaded gray.The origin of this distribution is analyzed by the yellow and cyan lines that show the component distributions for northward and southward IMF, respectively.These show that above-average transpolar voltages are dominated by, but are not exclusively for, data for southward IMF.Conversely, Φ PC below about 0.8<Φ PC > are dominated by northward IMF samples.The mauve and blue lines are the distributions of the coupling functions C f defined by Equation 9 with optimum fits (peak correlation within the requirement of linearity) for Φ PC (in mauve) and for SML (in blue), the exponents a, b, c, and d for which are given in Table 2.The green line is the distribution for the coupling function C BEA predicted by the Boyle et al. (1997) formula.The plot shows that the empirical function coupling function C f matches the large tail of the observed distribution well (Φ PC /<Φ PC > above about 1.4 in the case of the mauve line and above about 2.2 in the case of the blue line) but gives an extremely poor match at low and even mode values of the distribution.On the other hand, C BEA is a very good match to the whole distribution.
A quantitative assessment of the similarity of two distributions can be made using quantile-quantile (q-q) plots.These are presented in Figure 16, overlaid on data density plots which give an indication of the scatter in the agreement between the coupling function and Φ PC .The q-q plots contain 1,000 points so the largest point is for the top 0.1% of the distributions.The two distributions are identical if the points lie along the line of equality (the dashed cyan lines in Figure 16).Figure 16c shows that the two distributions are very similar for C BEA , and only begin to diverge at Φ PC above about 2.75<Φ PC > (the 98.4 percentile, i.e., for the largest 1.6% of samples).The deviation shows that C BEA is a bit "thick-tailed" (also called fat tailed or heavy tailed) in its distribution at large values compared to the observed distribution.The q-q plot for C f with d = 2.5 (the optimum fit to Φ PC ) does not diverge from the observed distribution quite as much as C BEA at the very largest values but is very slightly thicktailed at most values above the mode.The q-q plot for C f with d = 5.2 (the optimum fit for SML) is very far from the ideal match given by the cyan dashed line and so the distributions are very different.It should be remembered that the SML data used in Figures 15 and 16 are simultaneous with the Φ PC data and so these figures stress how different the optimum coupling functions for different parameters (in this case Φ PC and SML) can be.
It is important to note that none of these differences between the observations and coupling function distributions are captured by correlation coefficient r which is largest for C f with d = 2.50 (the optimum fit to Φ PC , shown in Figure 16a) and is somewhat lower for C BEA (Figure 16c), despite the fact that it reproduces the core of the distribution exceedingly well, and considerably lower for C f with d = 5.20 (the optimum fit to SML, shown in Figure 16b).2).

Discussion and Conclusions
We have analyzed the optimum coupling functions for a dataset of 65,133 hourly mean transpolar voltage estimates Φ PC observed between 1995 and 2020 by the northern-hemisphere SuperDARN radar network and matching sets of fully simultaneous am and SML index values, in the case of am these were linearly interpolated to the center times of the radar data hours from the three-hourly index.We have fitted using a generalized mathematical function that encompasses many proposed coupling functions and have carried out only a twofold test for overfitting (i.e., dividing the data into a fitting and a test data set roughly equal sample sizes).
Our aim in this article has been to establish some important principles concerning how the data can be averaged and fitted to ensure the IMF orientation term used does not bias the data in a way that does not match the physics of solar wind-magnetosphere coupling and also to ensure that the coupling functions derived are linear predictors of Φ PC , am, and SML.
Table 2 gives optimum values and the 2σ uncertainties derived here.Also given are the correlation coefficients r obtained and the fraction of the variance explained, r 2 .Note that correlations for SML and am here are for all the available data from 1995 to 2020 (but using the exponents derived here from the data subsets that are simultaneous with the radar data that meet the n e > 255 criterion (roughly a third of the full data).In addition, for am the raw three-hourly data are used to evaluate r and r 2 rather than the interpolated hourly values.The correlations for Φ PC are for only the n e > 255 data.It should be remembered that the noise introduced by spatial structure in the solar wind and the bow shock, on its own, limits r to about 0.9 (r 2 to about 0.81) and there are other noise sources (propagation lag uncertainty, instrumental errors in both the interplanetary data and the terrestrial disturbance indicator, seasonal, and/or UT effects on terrestrial data, data gaps, effects of averaging, nonlinearity of response, and dipole tilt effects).The values in Table 2 are slightly higher than previously proposed coupling functions, Figure 16.Data density plots with overlaid quantile-quantile (q-q) plots (white points joined by a red line) normalized coupling function predictions (C f /<C f >) against normalized transpolar voltage (Φ PC /<Φ PC >) for: (a) The optimum empirical coupling function, C f , given by Equation 9for the hourly mean Φ PC data (exponents given in Table 2 and distribution given by the mauve line in Figure 15).(b) The best-fit empirical coupling function, C f , for the hourly mean SML data (exponents also given in Table 2 and distribution given by the blue line in Figure 15).(c) The Boyle et al. (1997) formula, C BEA (distribution given by the green line in Figure 15).The q-q plots use 1,000 quantiles and so the largest point is for the top 0.1% of the two distributions compared.The correlation coefficient is given in each panel and is slightly higher for C f than for C BEA , even though Figure 15 shows that it does not give as good a match to the observed distribution-a fact confirmed by the less linear form of the q-q plot shown here.(b) Shows that the optimum coupling function for SML does not work at all well for predicting Φ PC .
10.1029/2021JA029946 32 of 41 but the gains in r 2 are marginal.It appears that coupling functions are achieving correlations almost as high as is possible for interplanetary observations made at L1 and the terrestrial disturbance data that we have available.
Table 2 also gives the performance of some theoretical coupling functions.For Φ PC these are simple prediction based on interplanetary electric field given by Equation 8 Our empirical fits exceed all these theoretical values, as indeed they should as they have three free fit variables.The results are quite similar in r 2 achieved to other empirical studies: for example, McPherron et al. (2015) explained 43.7%, 61.2%, 65.6%, and 68.3% of the variance in the hourly AL index using, respectively, epsilon ε (Perreault & Akasofu, 1978), V SW B s , the universal coupling function (Newell et al., 2007) and the optimum coupling function that they had derived which was B ⊥ 0.79 N SW 0.10 V SW 1.92 sin 3.67 (θ/2) (i.e., a = 0.79, b = 0.10, c = 1.92, and d = 3.67).Unfortunately, Newell et al. (2007) did not test the 20 coupling functions they considered against the am index.The closest they used to am was the kp index for which the main coupling functions correlation gave 100r 2 that ranged from 30% for ε to 58% for their universal coupling function.However, we note that there is a ±20% peak-to-peak "McIntosh" pattern in am caused by dipole tilt effects (Lockwood, Owens, Barnard, Haines, et al., 2020) which our optimum coupling function does not attempt to allow for with a dipole tilt term.This makes predicting 66.3% of the variation in am without it very encouraging.
The correlation for our transpolar voltage coupling function is r = 0.865 which means we are predicting 100r 2 = 75% of the variation in Φ PC .This is as high as has any that has been reported previously and is for a much larger dataset.An early study by Wygant et al. (1983) from a limited number of satellite passes explained 55% of the variation in Φ PC with the coupling function BV SW sin 4 (θ/2) (i.e., a = 1, b = 0, c = 1, and d = 4).Applying this to our 25-yr SuperDARN dataset of 65,133 samples with n e > 255, and using all best practice (i.e., computing the coupling function at 1-min resolution, averaging and the determining optimum lag) we find the Wygant et al. (1983) formulation explains 66% on the variance.Mori and Koustov (2013) surveyed the effectiveness of different coupling functions in predicting a Φ PC values from 1 yr of SuperDARN radar data.They found percentages of the variance explained ranging from 13% for ε in equinox up to 61% (for B ⊥ 1/2 V SW 1/2 sin 2 (θ/2); that is, a = 0.5, b = 0, c = 0.5, and d = 2), the latter is close to the optimum found here and testing against our data set we find it explains 73.5% of the variance in Φ PC , only very slightly lower than the value for our fit.
However, the benchmark test in transpolar voltage prediction is set by the coupling function of Boyle et al. (1997) who reported correlations of up to 0.87, explaining 75% of the variance of Φ PC , from observations from a number of LEO satellites over a 3-yr interval.The coupling function they derived was the addition of two terms: 10 −4 V SW 2 + 11.7Bsin 3 (θ/2) (where V SW is in kms −1 and B is in nT).A concern of any additive fit of this kind is that it may be open to overfitting and may not apply on all timescales.However, we can now check for overfitting by testing it against the fully independent SuperDARN Φ PC data used here.The correlation we obtain is r = 0.830, and so 68.8% of the variance in our Φ PC data is explained.This is not quite as high as Boyle et al. (1997) reported for their fit dataset, nor quite as high as the correlation we have found here; however, neither is it much lower.However, if we take the two terms in the Boyle function separately, we find the correlation with V SW 2 is very low with r = 0.2 (100r 2 = 4%) but that with Bsin 3 (θ/2) is 0.831 (100r 2 = 69.0%),and actually very slightly better than for the combination of terms.Hence, the key part of the Boyle et al. (1997) function with regard correlation coefficient has exponent a = 1, b = 0, c = 0, and d = 3 and its quite similar to the empirical fit derived in this article.However, we have shown that the Boyle function matches the distributions of transpolar voltages much more closely, even though its correlation coefficient is lower.This highlights correlation coefficient is not always the best performance metric to use when deriving and/or testing a coupling function.
We have studied the effect of different procedures in deriving the hourly means.In addition to the best practice combine-the average, <C f > 1hr , we computed all proposed coupling function [C f ] 1hr using the procedure of MEA15 (with averaging of 1-min values of θ and B ⊥ ) and also [C f *] 1hr for which θ and B ⊥ are both computed using hourly means of the B Y and B Z IMF components.Using [C f ] 1hr instead of <C f > 1hr typically lowers the variance explained by between 5% and 3%, whereas using [C f *] 1hr instead of <C f > 1hr typically lowers it by about 20%-40%.For the Boyle et al. (1997) parameter the behavior is strange in that for [C f ] 1hr the value is reduced from 68.8% to 68.0% but using [C f *] 1hr it plummets to 4%.The reason is the first term has become the larger 10.1029/2021JA029946 33 of 41 of the two because the coefficients of the two additive terms are no longer appropriate.Hence the first term of the Boyle equation has actually lowered the variance explained slightly but also made it unstable to the precise implementation.This is a general risk with additive terms.

The IMF Orientation Factor
As shown in Table 1, exponents d of an IMF orientation factor sin d (θ/2) of between 2 and 6 have been suggested from empirical studies and simulations with numerical global MHD models have suggested d as low as 1.5 (Hu et al., 2009) or even 1 (Borovsky, 2008;Fedder et al., 1991).For both the transpolar voltage Φ PC and the am geomagnetic index, we find that the IMF orientation factors in the coupling function for all suggested d between 1 and 6 all perform reasonably well in terms of the correlation coefficient.We find that marginally higher correlations for hourly averages for the low d exponents, the best correlations being for Φ PC at d = 2.1, for am at d = 1.3.However, we have shown that the distributions mean that these low d values are favored mainly because they weight the statistics toward near θ = 90° and against data for strongly northward IMF (θ approaching 0) and strongly southward (θ approaching 180°).The latter bias is, of course, particularly undesirable because periods of large θ drive the strong space weather which is often what we want the coupling function to predict and quantify.
As shown by Table 1 a great many studies have used sin d (θ/2) with d = 4 and this exponent has also been found for energy transfer across the magnetopause in MHD simulations of global energy transfer across the magnetopause (e.g., Laitinen et al., 2007).From the requirement of linearity across all clock angles we find the optimum exponents d are 2.50 ± 0.07 for Φ PC , 3.00 ± 0.22 for am, and 5.23 ± 0.48 for SML.

Other Coupling Function Exponents
The values of the other exponents a, b, and c (of B, ρ SW , and V SW , respectively), do, in general, depend on the exponent d used in sin d (θ/2).Some empirical fit studies have derived values for d that are not within the optimum range derived here, and the concern is that the associated a, b, or c have also been shifted from optimum values to compensate.
Table 2 shows our best fit exponents for Φ PC are somewhat different to the values of a = 1, b = −0.167,and c = 0.667 expected for the theoretical coupling function Φ SW based on the interplanetary electric field (Equation 8) and the differences imply that the reconnection efficiency η has quite considerable dependencies on all three parameters.Specifically, from our results and Equation 8 η appears to vary as B −0.358 , ρ SW 0.185 , and V SW −0.117 .Work is needed to see if these inferred external influences are consistent with the analysis of Borovsky and Birn (2014) who concluded that the reconnection voltage is not a function of the interplanetary electric field at all.One surprising value is the relatively large c (the exponent of V SW ) for the am geomagnetic index.Table 2 shows that the estimated power input into the magnetosphere P α fitted to the am index (for the 3-hr timescale) gives d = 2 and a coupling exponent α = 0.34 ± 0.04.From Equation 4 this predicts a = 0.68 ± 0.08, b = 0.32 ± 0.04, and c = 1.65 ± 0.08.Table 2 shows that although the values of a and b close to those expected for P α , c is much larger than predicted by P α .
From energy coupling into the magnetosphere from numerical MHD simulations Wang et al. (2014) derive a = 0.86, b = 0.24, and c = 1.47 (with a d of 2.7, similar to the 3.0 found here) which is extremely close to the above exponents for P α with α = 0.44 found by Lockwood, Bentley, Owens, Barnard, Scott, Watt, and Allanson (2019).From Table 2, the difference between the exponents b for am and SML for the solar wind mass density is Δb = 0.360-0.061= 0.299 (±0.074) and the difference between the exponents c for am and SML for the solar wind speed is Δc = 2.562-1.746= 0.816 (±0.114).Hence 2Δb = 0.598 (±0.148) which almost overlaps with Δc to within the 2σ uncertainties (and does overlap to within 3σ uncertainties).This suggests that there is a dependence of am on solar wind dynamic pressure of order P sw 0.36 = ρ sw 0.36 V sw 0.72 that is absent from SML. Lockwood, McWilliams, et al. (2020) find that 75% of the variation in am is explained by the estimated power input to the magnetosphere and that some of the remaining variance is associated with the solar wind dynamic pressure P sw , combined with the dipole tilt.They argue this is the effect of squeezing the near-Earth tail, an effect Lockwood, McWilliams, et al. (2020) show is found in both global MHD simulations and in the inference of an empirical model of the magnetopause location as a function of dynamic pressure.

10.1029/2021JA029946
34 of 41 On the other hand, our results for Φ PC and SML show almost no dependence on ρ SW .The SML result is particularly surprising as SML depends on the substorm current wedge which we imagine should also be influenced by the squeezing of the tail.Figure 11 of Lockwood and McWilliams (2021) shows influence of P sw (and hence ρ SW ) on Φ PC , am, and AL; it is complex and behavior depends on the IMF B Z component, but it is stronger at all B Z for am.  and McWilliams (2021).On the other hand, for SML there is no additional dependence beyond that of Φ PC (e ≈ 0-in fact the best fit is e = −0.1)for Φ PC ≤ 20 kV and e = 0.25 for Φ PC > 20 kV case.The difference for the Φ PC > 20 kV data is Δe = 0.61-0.25 = 0.36 which is the additional P sw dependence for am (compared to that for SML) derived above from the exponents in Table 2.
Hence, it is clear that am has a dependence on P sw that is not present in Φ PC and this is reflected in the coupling function we have derived for am.The reasons why the SML coupling function does not show the same P sw effect are not yet clear.Comparing Figures 17b and 17d we can see that the effect of P sw on am during southward IMF, and consequently enhanced Φ PC , is greater than for SML.This implies range geomagnetic indices from mid-latitude stations, such as am, are responding to a factor that does not greatly influence SML in addition to the substorm current wedge (which dominates SML).Matzka et al. (2021) note that the k-index (range) variation at mid-latitude stations (and hence increases in am and kp) arises from large-scale ionosphere-magnetosphere current systems and they are sensitive to a much broader longitudinal sector of the auroral oval than is detected by auroral stations.Hence mid-latitude positive bays reflect larger scale currents as well as the more localized substorm current wedge (McPherron & Chu, 2017).We note that Thomsen (2004) attributes the proportionality of mid-latitude range indices and transpolar voltage to the effect of polar cap expansion and that is indeed a factor; however, our results indicate that a parallel factor is that they are responding to the ionosphere-magnetosphere current circuits facilitated by the region 1 and region 2 field aligned currents and not just the substorm current wedge.It seems likely that this is the cause of the greater dependence of am of P sw than SML and AL.

Universality of Coupling Functions
We have found that that although the coupling functions for Φ PC and am could appear to have the same exponents if we use small datasets, when we use a very large one, as in this paper, the differences are shown to be highly significant and real.This implies that there is no such thing as a universal coupling function that can optimally predict both voltage disturbances in the magnetosphere and all geomagnetic disturbances and the coupling function needs to be tailored to the terrestrial disturbance indicator of interest in each case.This opens up new areas of systems analysis of the magnetosphere, namely combining the different responses of the various magnetospheric state indicators to different solar wind driving coupling functions (Borovsky & Osmane, 2019).It also has implications for how we might allow for "preconditioning" of the magnetosphere which is discussed in the next section.

Preconditioning
One major limitation of all the coupling functions discussed in this paper is that they assume that the terrestrial space weather index predicted is determined by the prevailing near-Earth interplanetary conditions only (allowing for the required propagation lag).This means that any preconditioning of the magnetosphere-ionosphere system is neglected and will contribute to the noise in the fits.To start to make allowance for preconditioning we have to make a distinction between two types: (a) preconditioning caused by the Earth's dipole tilt and (b) preconditioning that depends on the prior history of the solar wind.

Preconditioning by Dipole Tilt
Preconditioning by the dipole tilt can change the response of the magnetosphere, giving a larger or smaller response to a given solar wind forcing.This is an external factor depending on Earth's orbital characteristics which means it should be highly predictable.Studies show that genuinely global geomagnetic activity indices show a pronounced "equinoctial" (a.k.a."Mcintosh") pattern with time-of-year and Universal Time, associated with the tilt of Earth's magnetic dipole axis (see reviews by Lockwood et al., 2021;Lockwood, Owens, Barnard, Haines, et al., 2020).Attempts to expand the coupling function with a factor to allow for the effect of the dipole tilt were made by Svalgaard (1977), Murayama et al. (1980), Li et al. (2007), and Luo et al. (2013) and dipole tilt effects have been included in the filters used in the linear prediction filter technique (McPherron et al., 2013).
However, the detail of how this should best be done does depend on the mechanism that is responsible and there are a large number of postulated mechanisms aimed at explaining the Mcintosh (a.k.a.equinoctial) pattern.
One invokes the dipole tilt influence on ionospheric conductivities within the nightside auroral oval and postulates that the electrojet currents are weaker when conductivities caused by solar EUV radiation are low in midnight-sector auroral ovals of both hemispheres (Lyatsky et al., 2001;Newell et al., 2002).Other proposals invoke tilt influences on the dayside magnetopause reconnection voltage (Crooker & Siscoe, 1986;Russell et al., 2003) or the effect of tilt on the proximity of the ring current and auroral electrojet (Alexeev et al., 1996) or tilt effects on the stability of the cross-tail current sheet through its curvature (Danilov et al., 2013;Kivelson & Hughes, 1990;Kubyshkina et al., 2015).All of these effects have the potential to reproduce the McIntosh dipole tilt pattern, but which if any, are effective remains a matter of debate.Recently, strong observational (Lockwood, McWilliams, et al., 2020) and modeling (Lockwood, Owens, Barnard, Watt, et al., 2020) evidence argues that the amplitude of the McIntosh pattern increases with solar wind dynamic pressure, suggesting that the dipole tilt influences the degree of squeezing of the near-Earth tail by solar wind dynamic pressure.Given that dynamic pressure effects are included in most coupling functions via the ρ SW , and V SW terms, and that the effect is reasonably simultaneous with other solar wind effects, we might expect this effect to influence best-fit coupling exponents by raising b and c for geomagnetic activity but not for transpolar voltage.Thus, this mechanism has some relevance to understanding why the coupling function for transpolar voltage may be so different from that for the am index, as discussed in the previous section.

Preconditioning Related to Prior Solar Wind History
The storage-release system manifest in substorms shows that the response of the magnetosphere is inherently non-linear: the effect of a given burst of southward-pointing IMF, for example, is different at the start of the growth phase (when the open magnetospheric flux is flow) compared to at the end of the growth phase when it is high (Milan et al., 2021).Hence, the response that depends on the state of the magnetosphere is in at the time, and that is set by the prior history of solar wind magnetosphere voltage coupling.One technique to allow for the non-linearity of response caused by this type of preconditioning is local linear prediction (Vassiliadis, 2006;Vassiliadis et al., 1995).In this technique, moving average filters are continually calculated as the system evolves and these are used to compute the output of the system for this filter.The filter used is derived or selected according to the state of the system.Another way of dealing with this non-linearity is by using neural networks (e.g., Gleisner & Lundstedt, 1999).Our finding that the coupling function is significantly different for transpolar voltage and geomagnetic activity is significant in this respect.It means that if, for example, we wanted to allow for preconditioning due to the open flux in the magnetosphere, we would want to look at the prior history of an optimum coupling function for dayside reconnection voltage but would need to use a different coupling function to best predict, for example, the geomagnetic disturbance.
A number of other physical mechanisms have been proposed as ways of further preconditioning the magnetosphere.They include: mass loading of the near-Earth tail with ionospheric O + ions from the cleft ion fountain (Yu & Ridley, 2013); the formation of thin tail current sheets (Pulkkinen & Wiltberger, 2000); the development of a cold dense plasma sheet (Lavraud et al., 2006).Another proposed preconditioning effect is the effect on the reconnection rate in the cross-tail current sheet of enhanced ring current, as has been proposed by Milan et al. (2008Milan et al. ( , 2009) ) and Milan (2009).The magnetosphere sometimes responds to continued solar wind forcing (over a period of tens of minutes) by generating a substorm, or a string of substorms and sometimes with a steady convection event (e.g., Kissinger et al., 2012).Studies (e.g., Gleisner & Lundstedt, 1999) have demonstrated that the response of the auroral electrojet indices depends on the current Dst value.O'Brien et al. (2002) studied two intervals in which the solar wind coupling function was similar, one of which resulted in an isolated substorm and the other in a steady convection event.They noted the main difference was the pre-existing state of the magnetosphere in that prior to the substorm, the magnetosphere was quiet but whereas before the steady convection event the magnetosphere was already undergoing enhanced activity.McPherron et al. (2005) estimate that about 80% of steady convection events are associated with a substorm onset but thereafter the magnetospheric behavior diverges.The work of Juusola et al. (2013) strongly suggests that enhanced ring current is the reason that a steady convection event forms as opposed to a substorm, quite possibly through the mechanism proposed by Milan and co-workers.
Hence, preconditioning of the magnetosphere undoubtedly occurs through at least one mechanism, and this will be an inherent noise factor in the derivation of a simple correlative coupling function and hence a major limitation on the performance of that coupling function.The problem is that not only are the effects of the various mechanisms on the response different, the time constants of the prior activity that is influencing the response will be different in each case.This means that the time profiling of any preconditioning quantification factor in a coupling function using the prior history of the interplanetary parameters will depend on the mechanism.
To underline this point about the importance of the mechanism that is causing pre-conditioning, note that some mechanisms, such as the cold dense plasma sheet, would emphasize prior periods of quiet, northward IMF conditions as giving higher activity for a given input (Borovsky & Denton, 2006, 2010;Lavraud et al., 2006), whereas others, such as the ring current enhancement mechanism would emphasize prior periods of enhanced solar wind magnetosphere coupling.The time constants for forcing in the build-up to ring current enhancements (Lockwood et al., 2016) are different to those for the development of a cold, dense plasma sheet (Fuselier et al., 2015).Yet another proposed preconditioning mechanism involves the effect of solar wind dynamic pressure and thus would introduce yet another different precursor time profile (Xie et al., 2008).Some of these preconditioning effects have been predicted by numerical modeling (e.g., Lyon et al., 1998;Wiltberger et al., 2000) and it is quite possible that we may need numerical simulations to isolate the preconditioning effects and determine how best to allow for them.
The first column gives the basis of the formulation in each case, which is given in the second column.Columns 3-6 give the exponents a, b, c, and d and column 7 the F(θ) function used (h.w.r. Figure1ademonstrates that it is far from a valid assumption to take <C f > τ and [C f *] τ to be the same, using the example of theVasyliunas et al. (1982) energy transfer coupling function P α for a coupling exponent α = 1/3 (hence this P α is an example of C f with a = 2/3, b = 1/3, c = 5/3, and we here have used d = 4).The specific exponents do not change the general principles demonstrated by Figure1.The raw data in Figure1are the 9,930,183 valid Figure 1.
Figure 1b compares the combine-then-average values and the average-the-combine values for G (for the same example as shown in Figure 1a and in the same format), <G> τ , with a corresponding average-then-combine value [G] τ = <B ⊥ > a <ρ SW > b <V SW > c : again, all values have been normalized by dividing by the overall mean, G o .Note that we here use <B ⊥ > a and not [B ⊥ ] τ a (where [B ⊥ ] τ is defined by Equation Figure 1.Comparison of combine-then-average, average-then-combine and our compromise hybrid procedure for averaging 1-min data into 1-hr data (τ = 1 hr).In all panels, the horizontal axis gives the result of the combine-then-average approach which is what we ideally would wish to use to mimic solar wind forcing of the magnetosphere.The vertical axes in (a-e) give the result of an average-then-combine procedure.In each case the fraction of samples n/Σn is color-coded on a logarithmic scale, where n is the number of samples in small bins.The raw data used are 9,930,183 valid 1-min integrations of estimated power input to the magnetosphere, P α , and 11,646,678 valid 1-min values of the IMF clock angle θ and tangential component B ⊥ observed between 1995 and 2020 (inclusive).(a) The coupling function P α for α = 1/3 and d = 4 (the normalizing factor P o is the arithmetic mean of P α for all datapoints) in bins of P α /P o of size 0.08.The x-axis shows the means of 1-min values of P α , <P α > 1hr and the y axis the values [P α *] 1hr computed from 1-hr averages (including computation of the clock angle[θ]  1hr and the transverse magnetic field [B ⊥ ] 1hr from hourly means of the IMF components <B Z > 1hr and <B Y > 1hr ).(b) The corresponding plot for G, which is P α without the IMF orientation factor.(c) The IMF clock angle (in the GSM frame of reference) θ in bins that are 2° × 2°.(d) The tangential IMF component B ⊥ = (B y 2 + B x 2 ) 1/2 in bins of 0.5 nT × 0.5 nT.(e) For sin d (θ/2) in bins 0.01 × 0.01.(f) Compares <B ⊥ > a with <B ⊥ a > (where a = 2α for the P α coupling function).(g) Compares <sin(θ/2)> d with <sin d (θ/2)>.(h) The y-axis is the result of our hybrid averaging procedure for P α , [P′ α ] 1hr , defined by Equation 15.

Figure 2 .
Figure 2. Data density plots of normalized geomagnetic indices as a function of normalized transpolar voltage, Φ PC /<Φ PC > (a) the am index and (b) the SML index.The fraction of samples (on a logarithmic scale) in bins that are 0.03 wide in the x dimension and 0.06 in the y dimension.The black points are means in bins of Φ PC /<Φ PC > that are 0.1 wide and the black error bars are between the 1σ points of the distribution of normalized geomagnetic index in the bin.The mauve line is a third-order polynomial fit to the data.

Figure 3 .
Figure 3. (Top) Lag correlograms (linear correlation coefficient, r, as a function lag, δt) of predicted variations using 61-point boxcar (running) means of the coupling function C f from 1-min interplanetary parameters with hourly observations of the transpolar voltage Φ PC (in mauve), the interpolated am geomagnetic index (in blue) and hourly means of the SML index (in green).Note.That unless otherwise stated, C f in this and later figures refers to hourly means [C′ f ] 1hr , derived from our hybrid formulation, Equation15.The Φ PC , am, and SML data are all for the full 25-yr dataset, but only for hours when the number of SuperDARN radar echoes n e exceeds the threshold n min .This yields N = 65,133 data points.The hourly am data are derived from the observed three-hourly am values using PCHIP interpolation to the mid-points of the hourly integration periods for the radar data.The lag δt = 0 means that the radar data and the Omni interplanetary data are averaged over the same one-hour interval and positive δt corresponds to the interplanetary data leading the terrestrial data.The exponent d is assumed to be 3 but tests of values between 1 and 6 made negligible differences to the optimum values of δt, δt p , derived.The dark gray, lighter gray, and lightest gray areas define, respectively, the 1σ, 2σ, and 3σ uncertainty bands in the lag δt p and are defined using the Meng-Z test (see text for details).The vertical dashed lines give the lag δt p that yields the peak r, r p , which is 0.862 at δt p = 18.5 ± 1.3 min for Φ PC , 0.818 at δt p = 31.5 ± 4.0 min for am, and 0.803 at 45.3 ± 7.0 min for SML, the quoted uncertainties being at the 2σ level.(Bottom) The best-fit exponents a, b, and c as a function of δt (lines marked by squares, triangles, and circles, respectively), derived using the Nelder-Mead search algorithm to maximize r.
Figure4shows the distributions of various parameters relevant to the IMF orientation factor, all panels being for 1-min integrations of data and in the GSM frame of reference.Figure4is for 11,646,678 1-min Omni data samples from 1995 to 2020, inclusive.The vertical axis is the fraction of samples n/Σn in 100 bins of width that are 1% of the range of the horizontal axis.The sequence of Figures 4a-4e are fromLockwood, Bentley, Owens, Barnard, Scott, Watt, et al. (2019) and explain how strange, highly asymmetric distributions of 1-min samples of

Figure 5 .
Figure 5. Distributions of the IMF orientation factor F(θ) = sin d (θ/2) for d = 4, where θ is the IMF clock angle in GSM coordinates, for data averaging timescales τ of: (a) 1 min; (b) 15 min; (c) 1 hr (used in this paper); (d) 2 hr; (e) 6 hr; (f) 1 day;(g) a solar rotation period of 27 days, and (h) 1 yr.The numbers of samples, n, as a fraction of the total number Σn, in bins 0.01 wide are shown in each case and the dataset used is the same as in Figure4.The vertical mauve dashed lines are for the overall average of all samples.The vertical green line is at θ = 90° for which the IMF lies the GSM equatorial plane.Note.That the lowest bin in sin 4 (θ/2), which is 0-0.01,corresponds to a range in θ of 0°-36.9°whereas the highest bin (0.99-1) corresponds to 171.9°-180°.

Figure 7 .
Figure 7. Analysis of the effect of the exponent of the d of the F(θ) = sin d (θ/2) IMF orientation factor for all N = 65,133 samples which meet the criterion of the hourly mean number of radar echoes n e > n min = 255.For each value of d, the value of the other three exponents a, b, and c are derived by the Nelder-Mead simplex search method to maximize the correlation coefficient r between the hourly lagged coupling function C f .The results for observed Φ PC are in mauve, interpolated hourly values of am are in blue and hourly means of SML in green.The vertical dashed lines mark the peak correlation in each case, the vertical solid lines the optimum d (that gives linearity and determined from Figures 9-11) and the gray areas the 1σ, 2σ, and 3σ uncertainty bands of the optimum d.(a) The correlation coefficients, r, as a function of d.(b) The best fit values of the exponents a (identified by squares), b (triangles), and c (circles) as a function of d.
17)which we can test for.Equation 17 stresses the point that d is not an independent fit variable from the other exponents because for a given a, b, and c and set of interplanetary data, G is proscribed which means there is a unique exponent d in F(θ) = sin d (θ/2) that ensures the linearity of C f = G.F(θ) with T. The supplementary material toLockwood, Bentley, Owens, Barnard, Scott, Watt, et al. (2019) showed that this test yields F(θ) = sin 4 (θ/2) for

Figure 8 .
Figure 8. Data density plots of normalized coupling function C f /<C f > as a function of normalized transpolar voltage Φ PC /<Φ PC > in a similar format as Figure 2 (except mean values are shown in red and without standard deviations for clarity and the color scale is linear in fraction of samples, rather than logarithmic).The orange dashed line in each panel is the best linear regression to the individual data pairs and the green dashed line is the best second-order polynomial fit.The panels are for (a) d = 1.1.(b) d = 2.2.(c) d = 6.5.(d) d = 2.5.In each panel, the best-fit exponents a, b, and c are given for the d used (as in Figure 7), as is the correlation coefficient, r and the rms deviation of the normalized C f and Φ PC value pairs, Δ rms .

Figure 9 .
Figure 9. Tests of the IMF orientation term, F(θ) = sin d (θ/2) for the transpolar voltage Φ PC .(a-c) show plots of the means of R Φ = (Φ PC −i Φ )/G as a function of mean F(θ), both averaged for 25 bins of F(θ) that are 0.04 wide.G is given by Equation 14, where C f is the optimum coupling function for the optimum exponents a, b, and c for the d in question, as shown in Figure 7. (a) For d = 1.5.(b) For the derived optimum d of 2.50.(c) For d = 5.The green and red lines are linear and quadratic fits, respectively, to the mean values.The values of the linear regression coefficients s Φ and i Φ (see Equations 16 and 17) are given in (b), where the s Φ values are for B ⊥ in nT, N SW in 10 6 m −3 , V SW in km s −1 , and m SW in kg.The mauve line in (d) is the coefficient of the quadratic term of the second-order polynomial fit to the means, a Φ , shown as a function of d: the optimum d gives a proportional relationship between <R Φ > and <F(θ)>, that is, when a Φ = 0, marked by the vertical dashed line.Under the mauve line in three shades of gray are the 1σ, 2σ, and 3σ uncertainty bands in a Φ , the limits to which define the corresponding uncertainty bands in the optimum d, giving a 2σ uncertainty in the optimum d of ±0.08.Note.That in this case for Φ PC the differences between the uncertainty bands are often so small that they cannot be discerned; they are more clearly seen in Figure10for am and Figure11for SML.(b) Confirms this proportional relation at this optimum d = 2.50 for which the exponents are given in Table2.
Figure9shows that the required d is 2.50 ± 0.08 (at the 2σ uncertainty level) for Φ PC , Figure10shows that it is 3.97 ± 0.18 for am and Figure11shows that it is 5.20 ± 0.41 for SML.Hence the optimum IMF orientation factors for Φ PC , am, and SML are not the same within 2σ uncertainties and in all three cases are larger than the value derived by correlation.Essentially SML requires a function, that is, most like a half-wave rectified function and Φ PC requires a function that is least like one.The optimum d and their uncertainty bands for Φ PC , am, and SML are also shown in Figure7which reveals that the uncertainties do not overlap even at the 3σ uncertainty level.Note that the commonly used value of d = 4 is too large for Φ PC and am but too small for SML.Some agreement between the behavior of am and SML is to be expected because both are dominated, at high activity at least, by the effect of the substorm current wedge and so do show considerable agreement(Adebesin, 2016; supplementary information toLockwood, Bentley, Owens, Barnard, Scott, Watt, & Allanson, 2019).However, they are different indices and, as indicated by Figure2, they have a different relationship to the transpolar voltage.The values of s T and i T for the optimum d are given in Figures9b, 10b, and 11b.

Figure 10 .
Figure10.The same as Figure9for the am index.The blue line in (d) is the best-fit coefficient of the quadratic term of the polynomial fits illustrated in (a), (b), and (c), a am , under which the three gray areas define the 1σ, 2σ, and 3σ uncertainty bands in a am , the limits to which define the vertical uncertainty bands in the optimum d as shown.The optimum d giving the proportional relationship is d = 2.87 ± 0.18 for which the exponents a, b, and c are given in Table2.

Figure 11 .
Figure 11.The same as Figures 9 and 10 for the SML index.The green line in (d) is the best-fit a SML under which the three gray areas define the 1σ, 2σ, and 3σ uncertainty bands in a SML , the limits to which define the vertical uncertainty bands in the optimum d shown.The optimum d giving the proportional relationship is d = 5.20 ± 0.41 for which the exponents a, b, and c are given inTable 2.

Figure 12 .
Figure 12.Data-point density plots of predicted against observed values of (a) the transpolar voltage Φ PC , (b) the am geomagnetic index, and (c) the −SML index, each for their optimum d value defined in Section 3.These data are for the fit dataset which is for 2012-2020.In both cases, the optimum fit of C f has been scaled to the data by ordinary least squares linear regression.The numbers samples n (as a faction of the total number Σn) in bins, which are 1 kV × 1 kV wide in (a), 1 nT × 1 nT wide in (b), and 5 nT × 5 nT wide in (c), are color-coded on the logarithmic scales given.The diagonal mauve lines mark perfect agreement of observed and predicted values.The correlation coefficient r and the rms deviation Δ of observed and predicted values are given in each panel, along with the total number of valid data-point pairs, Σn.

Figure 13 .
Figure13.Same as Figure12but for the independent test dataset from 1995 to 2011, computed using the best-fit exponents, regression coefficients, and optimum lags derived as used for the fit dataset (2012-2020).The correlation coefficients r and the rms deviations Δ are very similar to the corresponding values for the fit dataset shown in Figure12.For these plots, the data had no role at all in deriving the fit exponents and coefficients.

Figure 14 .
Figure 14.Distributions of fitted values of exponents a (left panel), b (middle panel), and c (right panel) for fits to the transpolar voltage, Φ PC , drawn from the entire 25-yr dataset of 65,133 values with n e > n min = 255.The fraction of samples n/Σn in bins of width (1/30) of the maximum range of each exponent are plotted.In each case, three histograms are shown: (a) The light gray histogram bounded by the mauve line is for (1/25) of the whole dataset (N = 2,606 samples, on average corresponding to 1 yr of data).(b) The darker gray bounded by the blue line is for (1/10) of the whole dataset (N = 6,513 samples, on average corresponding to 2.5 yr of data).(c) The darkest gray bounded by the black line is for (1/2.5) of the whole dataset (N = 26,503 samples, on average corresponding to 10 yr of data).The standard deviation of the distribution is given in each case with the generic name σ xi where x is the exponent in question and i is the number of the dataset number.The distributions are generated by taking 1,000 random selections of   samples from the total of 65,133 samples with n e > n min = 255 available.The vertical dashed lines give the values for the full set of 65,133 samples.

Figure 15 .
Figure 15.Distributions of observed normalized transpolar voltage, Φ PC /<Φ PC >: the gray area bounded by the black line is for all 65,133 Φ PC hourly estimates for 1995-2019 for which the mean number of SuperDARN radar echoes exceeds 255; the cyan line is for those for which the lagged hourly IMF was northward (B Z ≥ 0 in GSM) and the yellow line for those for which it was southward (B Z < 0 in GSM).The green line is the distribution predicted by the Boyle et al. (1997) formula, C BEA and the mauve line the prediction for the best fit empirical coupling function, C f , given by Equation 9 with d = 2.50 and the other exponents given in Table 2.The blue line is the distribution of the best-fit C f to hourly means of the SML index (d = 5.20 with the other exponents also given in Table2).
and the Borovsky and Birn (2014) formulae for interplanetary Mach number M A < 6 and M A > 6.For am we use the best-fit version of the Vasyluinas et al. (1982) energy input formulation, P α (with d = 2 and coupling exponent α = 0.34) and for SML we shown the P α formulation with best fit values of d = 4 and α = 0.26.

Figure 17
Figure17is aimed at understanding the difference between the dependences of am and SML on the solar wind dynamic pressure P sw .It shows the (normalized) ratios of the geomagnetic indices per normalized transpolar voltage for (top panels) am and (bottom panels) SML, as a function of the normalized dynamic pressure P sw .divides the data up into subsets for Φ PC ≤ 20 kV and Φ PC > 20 kV which roughly corresponds to northward and southward IMF, but more importantly is above and below the change of gradient in Figure2b.For am there is an addition dependence of am, compared to Φ PC that varies as (P sw /<P sw >) e where e = 1 for Φ PC ≤ 20 kV and e = 0.61 for Φ PC > 20 kV (as shown by the dashed mauve lines).This is consistent with Figure11ofLockwood  and McWilliams (2021).On the other hand, for SML there is no additional dependence beyond that of Φ PC (e ≈ 0-in fact the best fit is e = −0.1)for Φ PC ≤ 20 kV and e = 0.25 for Φ PC > 20 kV case.The difference for the Φ PC > 20 kV data is Δe = 0.61-0.25 = 0.36 which is the additional P sw dependence for am (compared to that for SML) derived above from the exponents in Table2.

Figure 17 .
Figure 17.Data density plots for (top) the normalized am index per unit transpolar voltage (am/<am>)/(Φ PC /<Φ PC >) and (bottom) the normalized SML index per unit normalized transpolar voltage, (SML/<SML>)/(Φ PC /<Φ PC >) both as a function of normalized solar wind dynamic pressure (P SW /<P SW >) and in the same format as Figure 2. The data are divided into two subsets by transpolar voltage with Φ PC ≤ 20 kV in the left-hand panels and Φ PC > 20 kV in the right-hand panels.The mauve lines are the variations of P SWe /<P SWe >) for best-fit exponents e of 1, 0.61, 0.01, and 0.25 in (a-d).

Table 1 A
List of Proposed Coupling Functions That Share the General Functional Form B a Uncertainties in a, b, and c allow for both the fit uncertainties at a given d and the uncertainty caused by the uncertainty in d.The correlation coefficients are for all available data for 1995-2020: for Φ PC this means the hourly 65,133 samples with the mean number of radar echoes exceeding 255; for am this means the 69,028 three-hourly means with simultaneous interplanetary data yielding a valid hourly coupling function; and for SML this means the 241,848 hourly means with simultaneous interplanetary data yielding a valid hourly coupling function.The best-fit exponents are derived always from the 65,133 samples (using the optimum lag), using interpolated values in the case of am and simultaneous means for SML.