Complexity in the scaling of velocity ﬂuctuations in the high-latitude F-region ionosphere

. The temporal scaling properties of F-region velocity ﬂuctuations, δ v los , were characterised over 17 octaves of temporal scale from τ =1 s to < 1 day using a new data base of 1-s time resolution SuperDARN radar measurements. After quality control, 2.9 (1.9) million ﬂuctuations were recorded during 31.5 (40.4) days of discretionary mode soundings using the Tasmanian (New Zealand) radars. If the ﬂuctuations were statistically self-similar, the probability density functions (PDFs) of δ v los would collapse onto the same PDF using the scaling P s ( δ v s , τ )= τ α P ( δ v los , τ ) and δ v s = δ v los τ − α where α is the scaling exponent. The variations in scaling exponents α and multi-fractal behaviour were estimated using peak scaling and generalised structure function (GSF) analyses, and a new method based upon minimising the differences between re-scaled probability density functions (PDFs). The efﬁciency of this method enabled calculation of “ α spectra”, the temporal spectra of scaling exponents from τ =1 s to ∼ 2048 s. The large number of samples enabled calculation of α spectra for data separated into 2-h bins of MLT as well as two main physical regimes: Population A echoes with Doppler spectral width < 75 m s − 1 concentrated on closed ﬁeld lines, and Population B echoes with spectral width > 150 m s − 1 concentrated on open ﬁeld lines. For all data there was a scaling break at τ ∼ 10 s and the similarity of the ﬂuctuations beneath this scale may be related to the large spatial averaging ( ∼ 100 km × 45 km) employed by SuperDARN radars. For Tasmania Population B, the velocity ﬂuctuations exhibited approximately mono fractal power law scaling between τ ∼ 8 s and 2048 s (34 min), and probably up to several hours. The scaling exponents were generally less than that expected for basic MHD turbulence ( α =0.25), except close to magnetic dusk where they peaked towards the basic MHD value. For Population A, the scaling exponents were larger than for Population B, having values generally in the range expected for basic MHD and Kolmogorov turbulence ( α =0.25–0.33). The α spectra exhibited complicated variations with MLT and τ which must be related to different physical processes exerting more or less inﬂuence.


Introduction
In accordance with scientific method, the space physics community has adopted simple homogenous models of the ionosphere, yet accepted more complicated dynamics and morphology when compelled to do so by new observations. For example, consider an early, basic model of high-latitude ionospheric convection (Volland, 1975), the greater detail contained in more recent statistical models (e.g. Papitashvili and Rich, 2002), and the real time variability actually observed by DMSP spacecraft (http://cindispace.utdallas.edu/ DMSP/) and the Super Dual Auroral Radar Network (Super-DARN) Chisham et al., 2007). Diverse behaviour is superimposed on the underlying statistical average models.
The space environment is as complex as the systems studied in biology, economics, geology, meteorology, etc. By complexity, here we refer to the collective multi-scale dynamics which emerge from, but are more than, the sum of their parts. Note that many systems are complicated, but not necessarily complex. It is accepted that a complex systems approach, as well as the traditional reductionist approach, will help us to fully understand the Geospace environment and its response to the solar wind (Wanliss and Dobias, 2007). Chang et al. (2006) provided a definition of "dynamical complexity" as "a phenomenon exhibited by a nonlinearly interacting dynamical system within which multitudes of different sizes of large scale "coherent structures" are formed, resulting in a global nonlinear stochastic behaviour for the dynamical system, which is vastly different from that could be surmised from the original dynamical equations." The present study will make use of higher order statistical moments which are more sensitive to the large-scale fluctuations associated with the large-scale coherent structures forming within a system.
The present paper will quantify the complexity of ionospheric electric fields using scaling techniques which yield fractal dimensions suggesting complexity greater than found in the solar wind (Parkinson et al., 2007a). Although the scaling of electric field fluctuations in the ionosphere and solar wind may be related, the unusual complexity of the ionosphere must arise because it is an electrified, magnetised, multi-component fluid (ion, electron, and neutral gases) supporting a menagerie of complicated plasma instabilities. The behaviour of the ionosphere is also regulated by feedback between the coupled thermosphere and magnetosphere, and possibly the lower atmosphere (e.g. Tinsley et al., 2007).
Spatial and temporal series of measurements of natural systems fluctuate in characteristic ways which provide clues about the physical processes controlling the system. The characteristics of the fluctuations can change in space and time along with the dynamics of the system. Techniques of modern statistical analysis have been developed to quantify the self-similarity (Benzi et al., 1993;Bershadskii and Sreenivasan, 2004) or lack there-of in fluctuating measures. Herein we make use of peak scaling and generalised structure function (GSF) analyses (e.g. Chapman et al., 2005;Kiyani et al., 2006) to characterise the temporal scaling of electric field fluctuations in the high-latitude ionosphere inferred using SuperDARN observations of F-region Doppler velocity made at 1-s time resolution.
Peak scaling and GSF analyses have been used to identify power law scaling exponents of fluctuations in solar wind parameters (Hnat et al., 2002a(Hnat et al., , b, 2003. Solar cycle variations in these scaling exponents were subsequently studied by Parkinson et al. (2007a). For the range of temporal scales τ =64 s to ∼2 h, the scaling exponent for solar wind speed vary between the value expected for MHD turbulence, α=0.25 (Irosnikov, 1964;Kraichnan, 1965) (IK hereafter) to the value expected for hydrodynamic turbulence, α=0.33 (Kolmogorov, 1941(Kolmogorov, , 1962) (K41 hereafter). The scaling exponents for fluctuations in magnetic energy density B 2 vary between the value expected for K41 turbulence (0.33) to less than the value expected for Gaussian fluctuations (0.5).
SuperDARN radars often observe a boundary separating regions of small spectral width at lower latitude from regions of large spectral width at higher latitude. The spectral width boundary (SWB) is a proxy for the open-closed magnetic field line boundary Lester et al., 2001;Parkinson et al., 2002;Chisham and Freeman, 2004). Parkinson et al. (2003) separated measurements into Popu-lation A echoes with spectral width <75 m s −1 and Population B echoes with spectral width >150 m s −1 . Although both kinds of echoes can be found at any latitude, Population A (B) echoes are concentrated on closed (open) field lines. Parkinson (2006) used SuperDARN radar measurements to estimate temporal scaling exponents of zonal electric field fluctuations in the high-latitude ionosphere. For the range of temporal scales τ =1 to 256 min, the scaling exponents for Population A echoes were α=0.29±0.06 (0.21±0.06) in the dayside (nightside) ionosphere, values similar to those expected for IK turbulence. Population B echoes are more directly linked to the solar wind, yet the scaling exponents were only α=0.17±0.06 (0.12±0.06), suggesting complexity beyond that of IK and solar wind turbulence. As with the solar wind, the origin and meaning of this complexity requires a physical explanation. Uritsky et al. (2001) presented results relevant to the issue of solar wind versus stochastic components of variability in magnetospheric current activity as measured by the AE index. They found that AE index burst size was independent of the integrated solar wind energy input for time scales <3.5 to 5.0 h. For major solar wind energy input, the AE index burst size showed considerable scatter, but a convergence toward an upper bound imposed by the available solar wind input. Their results provide evidence for self-organised activity at low solar wind input, becoming forced and constrained at high solar wind energy input. Abel et al. (2006) used SuperDARN radar measurements to estimate spatial scaling exponents of electric field fluctuations. They found evidence for scale-free velocity fluctuations for separations from 45 km to ∼1000 km. The power law scaling exponent for fluctuations equatorward (poleward) of the open-closed boundary (OCB) was estimated to be α=0.39 (0.31). They suggested the scaling exponent poleward of the OCB was related to spatial structure in the solar wind, whereas the scaling exponent equatorward of the OCB was related to internal dynamics of the magnetosphereionosphere system. Their spatial scaling exponents need to be reconciled with the smaller temporal scaling exponents found by Parkinson (2006).
The preceding studies used SuperDARN common time observations made at 1 or 2 min resolution. This paper extends the results of these earlier studies by analysing a data base of experimental radar measurements made at 1-s resolution. This represents a 60-fold decrease in the minimum temporal scale used to characterise the fluctuations. The scaling exponents were also estimated using a new efficient method based upon minimising the differences between re-scaled probability density functions (PDFs). This enabled rapid calculation of the temporal spectra of scaling exponents ("α spectra") between 1 s and 1 day. The α spectra resolve the minimum temporal scale in the quasi-linear scaling regimes identifiable using SuperDARN, as well as substantial changes in the scaling exponents with temporal scale due to multi-fractal behaviour.

SuperDARN measurements
SuperDARN radars are over-the-horizon HF backscatter radars designed to map middle to large-scale high-latitude convection Ruohoniemi and Baker, 1998;Chisham et al., 2007). They do so by measuring the line-of-sight (LOS) Doppler velocity of 10-m scale F-region ionospheric irregularities known to drift with the E×B/B 2 velocity (Villain et al., 1985); hence the velocity measurements represent electric field measurements. During common time operations , SuperDARN radars scan 52 • of azimuth by electronically stepping the radar beam through 16 beam directions separated by 3.24 • . The transmitter pulse width usually corresponds to 45 km of group range, and the first range gate is 180 km, with 70 or more ranges separated by 45 km. The integration time is usually 3 or 6 s per beam so that one full scan is completed per minute or two.
The Tasman (Dyson and Devlin, 2000;Dyson et al., 2003Dyson et al., , 2005. Figure 1 illustrates their default azimuthal scans. For this study, the Tasmanian (New Zealand) radar was run in a 1-s camping beam mode for 31.5 (40.4) days of discretionary time evenly spread throughout the ∼1 year sunspot minimum interval 4 December 2005 to 21 November 2006 (Table 1). Beam 4 of the Tasmanian radar was used because it points down the magnetic meridian, thus measuring zonal electric fields, whereas beam 12 of the New Zealand radar was used because it passes above Macquarie Island, the location of an ionosonde and magnetometer. Approximately 2.8 million 1-s integrations were performed per radar throughout the campaigns.
The LOS Doppler velocity (m s −1 ) of echoes were estimated in real time using a standard pulse set scheme and analysis algorithm known as "FITACF" Ponomarenko and Waters, 2006). FITACF version 1.22 was used in this study and different results might be obtained when new pulse sets and analysis algorithms are developed. The quality of the 1-s radar velocities was partly controlled by rejecting echoes with signal-to-noise ratio less than 3 dB, velocity errors greater than 200 m s −1 , and all echoes with the ground scatter flag set. Echoes from range gates less than 630 km were rejected because they tend to be from meteors and the E region. Echoes from ranges greater than 3330 km were also rejected because of bad lags in the autocorrelation functions (ACFs). A filter was also used to reject velocities clustering near the ion acoustic speed.
An integration time of 1 s is shorter than the 3 s and 6 s employed during normal SuperDARN operations. The in- tegrity of the 1-s radar data was tested using 4-day discretionary mode data which interleaved 1-s and 6-s integration times. There was no systematic decrease in lag zero power (signal-to-noise ratio) when integrating for 1 s as opposed to the standard 6 s because the pulse width and pulse power remained unchanged. Spatial and temporal variations in the concurrent velocities were very similar, as were the PDFs of the velocities. Although the velocity "errors" were statistically larger for the 1-s radar observations, the 60-fold increase in the number of samples obtained when continuously dwelling for 1 s compensated for any loss of accuracy in the individual samples. Note that more reliable measurements of irregularities and their velocities are obtained using 1-s data when the ACFs are non-stationary on time scales of ∼3 s. Abel et al. (2006) calculated spatial structure functions using radar data recorded simultaneously at many echo ranges; hence temporal variations were completely separated from spatial variations. Here we use FITACF data to calculate temporal structure functions. Velocity changes observed at the smallest temporal scales represent pure temporal variations to a very good approximation. However, at temporal scales of 2 h and more the radar beam rotates a significant distance in MLT and spatial variations in the underlying twocell convection pattern begin to influence the observed velocity changes. This problem could be solved by a set of meridional radar beams fixed in MLT for the duration of the experiment (impossible).
SuperDARN radars do not make point-like measurements of the LOS velocity flow. The sampling volume of Super-DARN measurements are determined by the 300-µs (45-km) pulse width and the ∼3-4 • half power beam width which corresponds to ∼100 km at typical F-region group delays. The vertical dimension is controlled by where the refracted radio waves become orthogonal to the magnetic field lines. SuperDARN radars do not resolve velocity fluctuations at scale sizes less than the ∼100 km×45 km size of individual samples. The LOS velocities are also highly conditioned because the fundamental lag length effectively aliases velocities greater than 2-3 km s −1 (Parkinson, 2006). Figure 2 illustrates the spatial coverage of the 1-s velocities estimated using the Tasmanian radar. These magnetic clock dial plots show the number of valid echoes found in 2 h bins of magnetic local time (MLT) centred on each of the 24 h of MLT. The 45-km bins in group delay were mapped to magnetic latitude using the AACGM coordinate system (Baker and Wing, 1989). The two panels separate the counts for Population A and B echoes with spectral width <75 m s −1 and >150 m s −1 , respectively. The SWB often has an ephemeral quality and its location can be identified using various algorithms (e.g. Chisham and Freeman, 2003). It is not known whether regions of large and small spectral width are more dynamically fundamental than regions found above and below the SWB (Parkinson et al., 2008).

Spatial coverage of velocity samples
The overall spatial distribution of velocity counts shown in Fig. 2 are consistent with previous statistical studies (Villain et al., 2002;Parkinson et al., 2003). The peak count rates for Population B echoes were concentrated at higher latitude than for Population A echoes, consistent with the ten-dency for Population A (B) echoes to occur on closed (open) magnetic field lines. For the Tasmanian radar, 2 558 234 (2 996 816) Population A (B) echoes were selected for further analysis. The individual counts for the bins used in Fig. 2 ranged from ∼1000 to 20 000. Assuming shot noise, these sample sizes suggest estimates of statistical averages will have errors less than 3%.

High-latitude velocity fluctuations
The signatures of magnetic reconnection in the ionosphere have been observed to reoccur on time scales of several minutes (Pinnock et al., 1995). Transient reconnection may be a fractal process occurring at many spatial and temporal scales (Coleman and Freeman, 2005). Time series of LOS Doppler velocity, v los , measured by high-latitude SuperDARN radars are often exceptionally bursty. Digital ionosonde measurements of polar cap drift share a similar quality, especially when measurements are made in proximity to the ionospheric cusp (Parkinson et al., 1999). The presence of barely resolved flow bursts can make the time series seem amorphous. However, measurement noise must also contribute to the velocity scatter.
A velocity fluctuation at time t on scale τ is defined as δv los (t,τ )=v los (t)−v los (t−τ ). Here the smallest possible time scale is τ 0 =1 s. For convenience, we will consider the fluctuations occurring at the octaves of temporal scale, τ =2 n τ 0 , found between 1 s and 1 day. For example, Fig. 3 shows the velocity fluctuations, δv los , calculated using Tasmania beam 4 observations made during 10:00 to 18:00 UT on 17 December 2005. From top to bottom, the fluctuations were sorted according to representative temporal scales τ =1, 8, 64, 2048, 4096, and 8192 s. Each of the 6 panels actually shows a mass plot of the fluctuations at 8 representative, F-region group delays.
As will be demonstrated, Fig. 3 reveals a gradual change in the character of the velocity fluctuations with temporal scale, but there are three distinct regimes. There were numerous velocity fluctuations of amplitude ∼200 m s −1 between τ =1 and 8 s, giving the mass plots an amorphous quality. The velocity fluctuations between τ =64 and 2048 s were similar in character, but they were more structured than those for τ =1 to 8 s. By τ =8192 s, the character of the fluctuations was different again, only partly because of the influence of spatial variations.
Similar results were obtained for New Zealand beam 12, except the amplitude of the fluctuations were smaller at the shortest τ , suggesting the fluctuations were anisotropic.   Figure 4 show the Probability Density Functions (PDFs) of the velocity fluctuations calculated for each octave of temporal scale between τ =1 s and 1 day. These PDFs were calculated using the entire set of filtered Tasmania beam 4 fluctuations (2 915 381 samples at τ =1 s). As for hydrodynamic turbulence, the PDFs had more small fluctuations at small τ , and more large fluctuations at large τ . The PDFs were nearly indistinguishable for τ =1 to 8 s. However, the main feature of the PDFs was their strong leptokurtic form with weak convergence towards a Gaussian shape at the largest scales. The smooth black curves represent Gaussian functions with the same standard deviation as the measured PDFs at τ =1, 64 and 86 400 s. A Gaussian function is fatter than the equivalent leptokurtic function at small δv los ; the standard deviation of the latter is strongly influenced by its heavy tail. The leptokurtic PDFs are consistent with the exceptionally bursty time series data, Fig. 3.  Table 1. The fluctuating coloured curves correspond to temporal scales τ =1 s (black), 2 (purple), 4 (blue), . . . , 86 400 s (1 day) (back to purple). The smooth black curves represent Gaussian functions with the same standard deviation as the measured PDFs at τ =1, 64 and 86 400 s. The number of δv los samples ranged from a maximum of 2 915 381 at τ =1 s to a minimum of 189 483 at τ =65 536 s. A bin size of 1 m s −1 was used throughout this study.

Statistical scaling of the velocity fluctuations
The behaviour of fluctuations in a data record encode information about the underling dynamics of the system. A decomposition of the fluctuations according to scale size can help reveal this information. A familiar example is the use of power spectral analysis to identify a power law scaling regime with a scaling exponent of β=5/3, thus suggesting the presence of K41 hydrodynamic turbulence. The power spectrum represents a frequency domain decomposition of the 2nd order statistical moment (Sornette, 2000). Higher order moments provide further clues about the underlying dynamics. Note the interpretation of scaling is not limited to turbulence models; the statistical behaviour might implicate another dynamical process, for instance, some variant of self-organised criticality (Bak et al., 1987).

Data conditioning
Unless observations are made on very long scales, the estimation of higher order statistical moments becomes unreliable because of sensitivity to rare large fluctuations -the few large fluctuations observed may not be representative of the underlying population of large fluctuations. Conditioning in GSF analysis usually refers to, paradoxically, the rejection of the most extreme events to help stabilise the estimation of the statistical moments (Kiyani et al., 2006). SuperDARN radar fluctuations are highly conditioned to the 3 σ level by the ef-fective Nyquist velocity of ±2-3 km s −1 (Parkinson, 2006;Abel et al., 2007).
Another form of conditioning affecting the reproducibility of our results was used here. Time series of v los measured using 3 and 6 s integration times revealed small but significant numbers of outliers concentrated in side bands located ±(300-400 m s −1 ) about the main concentration of velocities. The speed of the outliers was similar to the ion-acoustic speed of type I instabilities (Kelley, 1989). They became more prevalent when using 1-s integration time and exaggerated the number of large fluctuations observed at small τ . The outliers do not belong to the dominant population of v los associated with irregularities drifting at the E×B/B 2 velocity and they were rejected as follows.
For a sample to be included in the data base there had to be valid samples located within the 10 s intervals immediately before and after. This simple rule had the effect of rejecting the majority of outliers which tended to occur intermittently during intervals of patchy scatter with weak power. Next, a sample was rejected if the absolute difference between the sample found immediately before and after exceeded a threshold of 250 m s −1 . This meant the sample would still be included if the velocity change was >250 m s −1 , yet the change persisted for more than one sample. This method was more effective at linearising the structure functions at small τ than application of a simple, running median filter.
Finally, Hnat et al. (2005) discuss the use of overlapping and non-overlapping intervals for generating differenced time series. Here we used overlapping intervals to greatly increase the number of available δv los samples at large τ . This had little or no affect on the estimation of statistical moments at short τ . However, the statistical independence of samples decreased with larger τ and the associated errors were larger than the number of samples might otherwise indicate. Note that for 32 days of observation, there are ≤768 non-overlapping intervals for τ ∼1 h, and ≤32 nonoverlapping intervals for τ ∼1 day. Thus we expect our results to become unreliable for τ >2 h.

Peak scaling and GSF analysis
Peak scaling and Generalised Structure Function (GSF) analysis are standard techniques which have been discussed extensively in the context of space physics (e.g. Hnat et al., 2002aHnat et al., , b, 2003Hnat et al., , 2005Chapman et al., 2005;Kiyani et al., 2006;Parkinson et al., 2006Parkinson et al., , 2007a. The PDFs, Fig. 4, revealed a temporal regime where the peak probability density decreased with increasing τ . If the peak probability density obeys a power law then P (0,τ )∝τ −α Peak and a plot of log 2 (P (0,τ )) versus log 2 τ will be linear with the gradient giving the peak scaling exponent α Peak . Peak probability densities represent the probability of the time series remaining the same, or of δv los →0. Unlike GSF analysis, estimates of α Peak are very sensitive to the statistical moments of negative order (1/δv los →∞).  5. Scaling of the PDF peaks, log 2 P (0,τ ) versus log 2 τ for the δv los curves shown in Fig. 4. The linear best fit to the data yields a slope of α Peak =0.234±0.007 using a PDF bin size of 1 m s −1 . The error bars assume Gaussian noise in the peaks, and are drawn at the 3 σ level. The vertical dashed line indicates the nominal outer temporal scale for the power law regime. Figure 5 shows log 2 (P (0,τ )) versus log 2 τ for the PDFs of δv los shown in Fig. 4. The quasi-linear regime between τ ∼1 s and 8 s has a near zero scaling exponent and was a stronger feature prior to data conditioning. This "flat" regime might arise because the radars do not fully resolve the fluctuations occurring at scale sizes less than the sampling volume, but measurement noise must also contribute. Another quasi-linear regime can be identified between τ ∼8 s and ∼2048 s (34 min). The straight line represents the results of a weighted linear regression with the error bars drawn at the 3 σ level. The estimated scaling exponent is α Peak =0.234±0.007 (1 σ error).
In Fig. 5, the point P (0, 86 400 s) is purely temporal in character because the radar beams samples the same spatial location every day. P (0, 86 400 s) is not co-linear with the points at shorter τ which are approximately temporal. This suggests a fundamental transition in the temporal scaling of fluctuations located somewhere between τ =2048 and  Table 1). The results show the variation of the moments S m (τ )= |δv sw (t,τ ) | m with temporal scale τ for orders m=−1 to 6. The results are quasi-linear on a log-log plot between τ ∼8 s and ∼2048 s (34 min). The variation of the standard deviation σ (τ )=[S 2 (τ )] 1/2 with τ (squares) is also shown. Error bars were drawn at the 3σ level but are negligible. Weighted least squares fits have been applied to all the moments. 86 400 s. However, the location of the outer scaling break is uncertain because the error bars and scatter in log 2 (P (0,τ )) increase with τ because of the decreasing number of truly independent samples (Fig. 4 shows significant asymmetry in the PDFs at large τ for the same reason). Corresponding results obtained using New Zealand beam 12 measurements were very similar (α Peak =0.220±0.009), but they are not shown for brevity. All the main scaling exponents are summarised in Table 2. GSF analysis makes use of higher order moments and provides a more comprehensive decomposition of the data than peak scaling or power spectra. GSF analysis involves calculating the time average of the statistical moments, S m (τ )= |δv los (t,τ ) | m where m is any real number, not necessarily positive. A variogram is a plot of S m versus τ for the different order m (Fig. 6). For single exponent scaling of δv los (t,τ ), S m (τ )∝τ ζ (m) , and a log-log plot of S m versus τ will reveal a straight line for each m with gradients ζ (m). If the time series of δv los is mono-fractal, then ζ (m)=α GSF m with a single scaling exponent α GSF . A plot of ζ (m) versus order m is known as a "ζ plot" (Fig. 7). Figure 6 shows the variogram for all the velocity fluctuations measured on Tasmania beam 4. The results are shown for τ =1 s to 86 400 s and for orders m=−1 to 6. The higher order moments are the least accurate because they are increasingly influenced by the few extreme events actually observed. Note that the moments for different m have been shifted in the y-direction to faithfully reveal their full dynamic range. Similar to peak scaling (Fig. 5), the variogram reveals three basic scaling regimes. From τ =1 s to ∼8 s there is a distinct variation of S m with τ ; this was especially the case prior to application of the ion acoustic filter. Thereafter there is a linear increase in the growth of S m between τ ∼8 s and ∼2048 s. Lastly, there is a divergence of the gradients for τ >2048 s. Figure 6 also shows the variation of the standard deviation σ (τ )=[S 2 (τ )] 1/2 with τ (squares) which provides an estimate of the Hurst exponent (Sornette, 2000). between PDFs and re-scaled PDFs at the next two largest τ versus temporal scale τ and a plausible range of scaling exponents α PDF . The best estimates of α P DF were given by the minima in the sums of the differences across the entire PDF (purple). These results are for all the δv los samples observed on Tasmania beam 4 (Table 1). Figure 7 shows the ζ plot where each ζ (m) is the gradient of the m-th best fit curve between τ =8 s and 2048 s (Fig. 6). The weighted least squares fit to the ζ plot was to the points for m=0 to 6, with the solution constrained to pass through the origin because ζ (0)=0 by definition. The use of weighted least squares fitting also weights the fit to the lower order moments which are more statistically significant. The result is α GSF =0.126±0.004 (1 σ error), but α GSF may have been ∼0.16 when considering only ζ (m) between m=0 and 3. Departure from linearity across the full range of m might be due to under-sampling of the largest fluctuations, or multi-fractal behaviour .
If the ζ plot curves were purely linear, then the least squares fit errors would be interpreted as statistical errors associated with measurement noise. In this study, the ζ plot curves are non-linear to a greater or lesser extent depending on the particular data set. Hence we know larger 1σ errors tend to indicate stronger non-linearity and multi-fractal behaviour, and the errors have physical meaning in the context of our results. Departures from linearity in the ζ plots were one of the main motivations for developing the "α spectrum" analysis presented in the following section.
3.3 PDF scaling collapse: "α spectra" It is well known the presence of power law scaling can be tested by collapsing the PDFs at separate τ onto a single curve using the transformations P s (δv s ,τ )=τ α P (δv los ,τ ) and δv s =δv los τ −α (e.g. Hnat et al., 2002a, b). Scaling collapse was used here to provide alternative estimates of the scaling exponents. The basic idea was to estimate α PDF by minimising the sum of the differences between the re-scaled PDF elements. A simple algorithm involved storing every PDF(τ , δv los ) in computer memory and then re-scaling them for a range of plausible α PDF values. The sum of the differences were weighted according to the peaks of the PDFs or the wings of the PDFs to obtain results more closely related to those of either peak scaling or GSF analysis, respectively. Computational efficiency was achieved because no statistical moments were calculated and fitted by linear regression.
Quasi-linear scaling regimes were identified in the variograms (e.g. Fig. 6), but the gradients of the GSFs gradually change with every octave of τ . The preceding minimisation algorithm was applied to a short running window of τ to construct an "α spectrum" revealing the variation of α PDF with τ . Figure 8 shows the relative variation of the sums of the differences ( Diffs) between every PDF at temporal scale τ and the re-scaled PDFs at the next two largest τ . The best estimates of α PDF were inferred from the minimum sums of differences across the entire PDF (purple). Similar results (±0.01) were obtained when using various methods of weighting the sums of differences, including weights inversely proportional to errors assuming Gaussian noise.
In effect, Fig. 8 constitutes an "α spectrum" for all the velocity fluctuations observed at 1-s resolution on Tasmania beam 4. Reasonable scaling collapse was achieved for all τ <8192 s. Consistent with the indistinguishable PDFs for small τ (Fig. 4), α PDF was <0.1 for τ =1 to 4 s. Thereafter, α PDF increased, reaching a maximum value of 0.27 for τ =64 s to 128 s, and then decreasing to a local minimum of α PDF =0.17 at τ =2048 s. For τ >4096 s, no meaningful scaling collapse was achieved, possibly due to inadequate sample sizes. The scaling exponents peaked at the values expected for intermittent IK turbulence, but clearly the behaviour was more complicated than a single mono-fractal from τ =1 to 4096 s.

MLT dependence of scaling exponents
Population A consists of echoes with spectral width ≤75 m s −1 concentrated on closed field lines, whereas Population B consists of echoes with spectral width ≥150 m s −1 concentrated on open field lines. For both populations and radars, approximately 2-3 million velocities survived rejection. Figure 2 showed that after further separation into 2-h bins of MLT, thousands of velocities per bin were available for further analysis. Peak scaling, GSF analysis, and PDF rescaling were applied to each 2-h bin of MLT. Figure 9 shows the MLT dependence of the different scaling exponents for Population A (a) and Population B (b) velocities measured on Tasmania beam 4, and Fig. 10 shows the same for New Zealand beam 12.
In Fig. 9a the scaling exponents α PDF (τ =64 to 256 s) and α Peak (τ =8 to 2048 s) are similar but different to α GSF and H (τ =8 to 2048 s) which are also similar. When estimating α PDF by re-scaling the PDFs the peaks of the PDFs exerted a stronger influence on the sums of differences because no Tasmania Beam 4 39 Fig. 9. Tasmania Beam 4 Fig. 9. (a) The variation of scaling exponents α PDF (black), α Peak (purple) α GSF (red), and H (orange) with MLT for Population A echoes measured using Tasmania beam 4. The α PDF results are shown for τ =64 to 256 s. All other scaling exponents were calculated using the quasi-linear regimes identified between τ ∼8 s and ∼2048 s (34 min) in log-log plots. Error bars are at the 1 σ level and scaling exponents estimated using less than 23 284 samples were rejected. (b) The same except for Population B echoes. In both parts, the continuous curves show the number of δv los samples identified at τ =1 s, with corresponding labels at right (e.g. 50 K=50 000 samples).
attempt was made to weight the calculations in favour of the wings of the PDFs. Also, both techniques must be sensitive to the chosen bin size of the PDFs (Parkinson, 2006). Hence α PDF and α Peak values are sometimes more alike than α GSF and H . The latter should also be similar because H is estimated using the m=2 structure function.
The continuous curves shown in Figs. 9 and 10 indicate the number of δv los samples identified at τ =1 s in each 2h bin of MLT. Obviously, the results are more significant where the sample sizes are larger: for smaller sample sizes, the 1 σ errors were generally larger and there was a greater divergence in the scaling exponents estimated using different methods. This is especially the case for New Zealand Population A (Fig. 10a). Although ∼1000 samples is adequate for estimating the low order moments at small τ , the δv los samples were identified at many octaves of τ , and few truly independent samples were identified at the largest τ . As will be shown, the scaling exponents estimated using less than approximately 23 284 samples were unreliable and thus rejected in Figs. 9 and 10.
Ironically, our understanding of the analysis methods was improved because limited sample sizes in some MLT sectors stressed our estimates beyond the limits of statistical significance. Nevertheless, Figs. 9 and 10 point towards some interesting results in need of confirmation with larger sample sizes. Figure 9 shows that the 1σ errors for Population B scaling exponents were much smaller than for Population A. Also, the different kinds of scaling exponents were similar for Population B, but not Population A. The corresponding ζ plots for individual 2-h bins of MLT (not shown) indicate Population B results exhibited greater linearity up to higher order. The scaling collapse was also generally better for Population B velocities. Combined, this suggests Popula- tion B (A) velocities were more mono-fractal (multi-fractal) (Parkinson, 2006).
Scaling exponents estimated for Tasmania Population B had the best consistency with MLT and estimation technique (Fig. 9b). Hence this data (at all MLT) was used to probe the effects of decreasing sample size on the reliability of the scaling exponents. Figure 11 plots estimates of the scaling exponents using all techniques and their 1 σ errors versus log 2 N where N is the number of δv los samples at τ =1 s. Only the first N samples were used as N was decreased by the same amount at all τ ; hence statistical stationarity throughout all the Table 1 intervals was assumed. The results show that α GSF and H were the most robust exponents for small sample sizes, whereas α PDF and α Peak rapidly became meaningless for small N. Examination of the PDFs, spectra, and structure functions confirmed the scatter and errors became too large for N <23 284 (log 2 N <14.5). Ideally, many millions of independent samples should be used for a statistical analysis of this kind.
Superimposed in Figs. 9 and 10 are horizontal dashed lines indicating the scaling exponents expected for K41 hydrodynamic (α=0.33) and IK MHD turbulence (α=0.25). For Tasmania Population A, the H and α GSF exponents are in the range expected for MHD turbulence (especially ∼00:00-01:00 MLT), and MHD to K41 turbulence (especially ∼16:00-23:00 MLT). There were strong minima in all the scaling exponents near 12:00 MLT, but the exponents were rejected because of small sample sizes. The PDFs for large τ will have too much probability density concentrated As in Fig. 9, scaling exponents were rejected when there were fewer than 23 284 samples at τ =1 s (dotted regions).
at small δv los because of the small number of large fluctuations observed during ∼31.5 (40.4) days. Figure 9b shows that all the Population B scaling exponents were in the range ∼0.1-0.25, but mostly smaller than the value expected for basic IK turbulence (0.25). There was a weak local maxima in the scaling exponents at ∼03:00 MLT, but more clearly defined maxima occurred at 12:00 MLT and 18:00-19:00 MLT, possibly suggesting different levels of intermittency in the velocity fluctuations near magnetic noon and dusk. Depending on the physical model which ultimately explains these results, the intermittency was stronger on open field lines (Population B) than on closed field lines (Population A) for both radars.
The results obtained for New Zealand beam 12 (Fig. 10) were generally less significant and many more scaling exponents were rejected. Like the results for Tasmania beam 4, they suggest the scaling exponents were larger for Population A. However, the different results obtained for the two radars point toward possible anisotropies in the scaling exponents in some MLT sectors. Table 2 lists the scaling exponents estimated using data integrated across all MLTs; these values are far more reliable.

MLT dependence of α spectra
The scaling exponents α PDF shown in Figs. 9 and 10 were "localised" values estimated by re-scaling the PDFs across τ =64 s to 256 s, but the α spectrum (Fig. 8) revealed a complicated variation of α PDF with temporal scale τ . Figure 12 shows the variation of α PDF with τ and MLT for Tasmania Population A (a) and Population B (b). As with Fig. 9, results have been excluded where there were too few samples (dotted regions). The scaling exponents were very low, mostly <0.12 (black and purple), for τ <8 s. They were also close to or less than zero (black), for τ >1024 s. As mentioned, there was a different scaling regime for τ <8 s and the velocities do not re-scale for τ >2048 s. Thus we focus on the MLT variation of α PDF for τ ∼8 s to 1024 s. Consistent with Figs. 9 and 10, Fig. 12 shows the scaling exponents were larger for Population A (a) than Population B (b). The scaling exponents were basically MHD-like (0.25) (green) for Population A (a), yet predominantly less than MHD (blue) for Population B (b). However, more features are revealed by this analysis. For Population A, the scaling exponents were most Kolmogorov-like at ∼17:00-18:00 MLT for τ ∼64-128 s. They even became like those expected for Gaussian fluctuations in the pre-midnight sector ∼22:00-23:00 MLT for τ ∼256-512 s. For Population B (Fig. 12b), there were two strong MHD maxima in the scaling exponents, one at ∼11:00-13:00 MLT for τ ∼32-64 s, and a stronger feature at ∼18:00-19:00 MLT for τ ∼64-128 s. The scaling exponents do not supercede 0.15 in other MLT sectors.
Results for the New Zealand radar are not shown because they are reliable over a smaller range of MLT.

Discussion and summary
The present results have been interpreted in the framework of the scaling exponents expected for turbulence, but other dynamical theories may ultimately explain our results. If we think of larger fractal dimension (D) indicating greater complexity (D=2−α), then smaller scaling exponents indicate greater dynamical complexity. Multi-fractal behaviour indicates different physical processes operating at different spatial and temporal scales, and for different event sizes (Chang and Wu, 2008), adding another aspect to the complexity. For velocity fluctuations in the high-latitude ionosphere, we found scaling exponents less than the value expected for basic IK turbulence (0.25). This suggests complexity greater than found for fluctuations in solar wind (Parkinson et al., 2007a), perhaps related to complicated coupling between the thermosphere, ionosphere, and magnetosphere.
The solar wind ε parameter (Perreault and Akasofu, 1978) incorporates the solar wind speed, magnetic pressure, and the IMF clock angle, and is a measure of energy coupling to the magnetosphere via magnetic reconnection. The ε parameter has been widely used but it is not the optimum parameter (Newell et al., 2007). Fluctuations in ε are usually in the range expected for K41 turbulence (α=0.33), except during fast streams when they become MHD-like (Parkinson et al., 2007a). It is thought that comparing measures of solar wind and magnetospheric activity might reveal the extent to which the magnetosphere is self-organised or directly driven by the solar wind (Freeman et al., 2000;Uritsky et al. 2001).
It is well known that the behaviour of the noon sector ionospheric cusp is strongly driven by fluctuations in the solar wind via magnetic reconnection (e.g. Pinnock et al., 1993). Intense ionospheric flows associated with individual flux transfer events can be related to local-scale features in the ionosphere (Pinnock et al., 1993). The scaling of velocity fluctuations in the noon sector high-latitude ionosphere might be related to fluctuations in the driving solar wind. The present radar observations were made throughout solar minimum, an interval when ε fluctuations were expected to be Kolmogorov-like. Population B velocities tend to occur on open magnetic field lines. Figures 9b and 12b suggest the electric field fluctuations had scaling exponents beyond (less than) the values expected for basic MHD turbulence in the polar cap ionosphere. However, Figs. 9b and 12b revealed local maxima in the scaling exponents at ∼12:00 MLT, with the values peaking towards those expected for MHD turbulence. Perhaps these noon sector enhancements are a signature of the magnetospheric cusp and solar wind. Figure 12b showed another stronger peak in α PDF at ∼18:00-19:00 MLT. Perhaps this feature is related in some way to the Kelvin-Helmholtz instability acting in the dusk sector ionosphere.
Population A velocities tend to occur on closed magnetic field lines. Figures 9a and 12a showed the scaling exponents were generally in the range expected for MHD and K41 turbulence, with values larger than found for Population B. Population A fluctuations also exhibit stronger multi-fractal behaviour, as indicated by the larger error bars and discordant values estimated using different techniques (Fig. 9a). Population A velocities are concentrated deeper inside the magnetosphere than Population B velocities, and it would be surprising if their fluctuations were more directly controlled by the solar wind.
Here we presented "α spectra" showing the variation of the scaling exponents with τ . The spectra were calculated using data separated into 24 bins of 2-h in MLT. Variations in the α spectra are related to data conditioning and the available number of samples. However, when the calculations use a large number of samples, they provide information about MLT variations in the relative contribution of different physical processes. For example, there are well known MLT variations in auroral activity, flux transfer events, plasma instabilities, turbulence, ULF waves, etc. The electric field fluctuations approximate to a single quasi-linear power law scaling regime over τ ∼8 s to 2048 s (34 min) (Figs. 5 to 7). This might indicate the different physical processes are related via the same multi-scale energy cascade. Figures 5 and 6 revealed a separate scaling regime for velocity fluctuations beneath τ ∼8 s, yet the author is unaware of a fundamental transition in the plasma physics at τ ∼8 s. For example, the gyro period of an oxygen ion is only ∼18 ms in the high-latitude ionosphere. In general, fluctuations in fluid flow at small temporal scales tend to occur at small spatial scales. For example, the speed of the Atlantic Gulf Stream does not fluctuate over its entire length at τ ∼1 s! In the case of the ionosphere, many velocity fluctuations at τ ∼1 s will be concentrated at small spatial scales. Hence a plausible explanation for the transition in velocity fluctuations beneath τ ∼8 s is spatial averaging of small-scale velocity fluctuations across the ∼100 km×45 km radar pixels.
The effects of spatial averaging might suggest there is no advantage in making SuperDARN measurements using an integration time less than ∼8 s. However, SuperDARN radars routinely use 3-s integration times during 16-beam full scan observations to achieve the Nyquist frequency necessary to prevent aliasing of convection cells and ULF wave signatures changing on time scales less than 2 min. Parkinson (2006) identified power law scaling regimes from τ =1 to 256 min using SuperDARN velocities measured using 3 and 6 s integrations once every minute or two during year 2000. Allowing for experimental error, the results listed in Table 1 of Parkinson (2006) are similar to the results reported here (Table 2) using 1-s resolution data. However, there may be significant solar cycle changes in the characteristics of electric field fluctuations and these will be the subject of future studies.
However, the results of the present study extend beyond those of Parkinson (2006) because the scaling behaviour is quantified at much smaller temporal scales, as well as the scaling exponents and multi-fractal behaviour versus MLT for echoes found predominantly on open or closed field lines. Furthermore, the α-spectra quantified the variation of the scaling properties with temporal scale, whereas the complementary "rank-order multifractal spectrum" recently proposed by Chang and Wu (2008) evaluates the variation of scaling properties with event size. Abel et al. (2007) reported the results of a structure function analysis of spatial fluctuations in ionospheric plasma velocities measured using the Halley SuperDARN radar. They presented a thorough discussion on the effects of data conditioning. The results obtained without data conditioning were difficult to reconcile with existing turbulence models. After conditioning the data by removing all velocity fluctuations with magnitudes >3 standard deviations, they found the structure function scaling was most consistent with IK versions of the P and log-normal models of turbulence.
In this study, the TIGER radar data were conditioned in the same way as Abel et al. (2007). The results obtained, combined with those of Parkinson (2006), imply that especially for TIGER Population B, there was reasonable conformity to power law scaling from τ =8 s (and possibly less) to ∼256 min (and greater?). However, the scaling exponents for Population B were consistently less than the basic IK value of 0.25. The discrepancies between the results obtained in the spatial (Abel et al.) and temporal domains is not understood. Nor is the author aware of a physical theory which explains the small values of α and their complicated variation with MLT and temporal scale.
Here we summarise the design and main results of this study:  (Table 1). This high time resolution (HTR) mode generated 60 times more valid samples on the chosen beam than during HTR common time. The 1-s resolution also provided a 60-fold decrease in the minimum temporal scale used to characterise the electric field fluctuations.
2. The velocities were analysed separately for Population A echoes with low spectral width (<75 m s −1 ), and Population B echoes with high spectral width (>150 m s −1 ). Population A echoes are concentrated on closed field lines and Population B echoes on open fields, but there is no consensus about the physical meaning of the SuperDARN spectral widths (Parkinson et al., 2008). The velocities were also analysed using separate 2-h bins in MLT, but too few samples were available for accurate calculations in some MLT sectors.
3. The velocity fluctuations were analysed using peak scaling, GSFs, and a new PDF re-scaling technique which involved minimising the sums of differences between re-scaled PDF elements. When applied to a running window of τ , PDF re-scaling permitted calculation of "α spectra" revealing detailed variations in the scaling properties with τ . The different analyses were repeated for Populations A and B, and 2-h bins of MLT.
4. Tasmania beam 4 is a magnetic meridional beam and New Zealand beam 12 points poleward and westward. The two beams intercept at −78 • magnetic latitude. Fewer echoes were recorded with the New Zealand radar and the results are less reliable. However, after rejecting the least significant results, the scaling exponents estimated using the Tasmanian and New Zealand radars were similar but not identical at all MLTs (Figs. 9 and 10). The possibility of anisotropic scaling remains an open question.
5. The variation of the statistical moments with τ were divided into three main regimes: (i) a quasi-linear regime from τ ∼1 to 8 s, probably related to spatial averaging of velocity components across the ∼100 km×45 km radar pixels, (ii) a quasi-linear power law scaling regime from τ ∼8 to at least 2048 s (34 min), probably indicative of an inertial energy cascade, and (iii) a different scaling regime for τ >2048 s associated with large-scale convection. The latter results had large errors due to a decreasing number of truly independent samples, and thus our choice of 2048 s for the outer scale is probably conservative.
6. For TIGER Population B echoes, the scaling exponents estimated using the various techniques were similar, and the error bars were relatively small (Fig. 9b). There was reasonable conformity to a single exponent power law scaling regime across τ ∼8 to 2048 s (34 min). Most of the scaling exponents were less than the value expected for IK MHD turbulence (0.25), except near magnetic dusk where the α PDF exponents peaked toward IK values.
7. Most of the scaling exponents for Population A echoes were larger than the values found for Population B, and generally in the range expected for IK and K41 turbulence (∼0.25-0.33). However, there were also larger discrepancies between the scaling exponents estimated using the different techniques, and consistent with the non-linearity of the ζ plots, the error bars were relatively large (Figs. 9a and 10a). This suggests stronger multi-fractal behaviour in regions of low spectral width and on closed field lines.
8. The α spectra revealed some interesting variations in the scaling exponents with MLT: the scaling exponents for Population A were most Kolmogorov-like at ∼17:00-18:00 MLT for τ ∼64-128 s, and there were localised enhancements to Gaussian-like values in the premidnight sector ∼22:00-23:00 MLT for τ ∼256-512 s. For Population B, there were two strong MHD maxima in the scaling exponents, one at ∼11:00-13:00 MLT for τ ∼32-64 s, and a stronger feature at ∼18:00-19:00 MLT for τ ∼64-128 s. These variations must be related to different physical processes exerting more or less influence at different MLT and τ .
9. Electric field fluctuations in the polar cap ionosphere (Population B) can be approximated by single, beyond MHD scaling exponent, but there was significant nonlinearity present in the ζ plots. A more detailed analysis reveals variations in the behaviour of electric field fluctuations with τ and MLT. The variations are more strongly multifractal for Population A echoes on closed field lines. Future studies will need to quantify the variation of the scaling exponents with solar activity, season, and geomagnetic activity, and understand their numerical values in terms of the best available physical theory.