On multifractality of high-latitude geomagnetic fluctuations

In order to contribute to the understanding of solar wind/magnetosphere interactions the multifractal scaling properties of high-latitude geomagnetic fluctuations observed at the Thule observatory have been studied. Using the local observatory data and the present experimental knowledge only it seems hard to characterize directly the, presumably intermittent, mesoscale energy accumulation and dissipation processes taking place at the magnetotail, auroral region, etc. Instead a positive probability measure, describing the accumulated local geomagnetic signal energy content at the given time scales has been introduced and its scaling properties have been studied. There is evidence for the multifractal nature of the so defined intermittent field ε, a result obtained by using the recently introduced technique of large deviation multifractal spectra. This technique allows us to describe the geomagnetic fluctuations locally in time by means of singularity exponents α, which represent a generalization of the local degree of differentiability and characterize the power-law scaling dependence of the introduced measure on resolution. A global description of the geomagnetic fluctuations is insured by the spectrum of exponents f(α) which represents a rate function quantifying the deviations of the observed singularities α from the expected value. The results show that there exists a multifractal counterpart of the previously reported spectral break and different types of f(α) spectra describe the fluctuations in direct dissipation or loading/unloading regimes of the solar wind-magnetosphere interaction. On the time scale of substorms and storms the multifractal structure of the loading-unloading mode fluctuations seems to be analogous to the simple multiplicative P-model, while the f(α) spectra in direct dissipation regime are close but not equal to the features of a uniform distribution. Larger deviations from the multiplicative model are observed when the influence of the solar wind fluctuations is examined. On this basis it is expected that an extended multifractal analysis of the singularity structure of near-Earth plasma system fluctuations would lead to improved geomagnetic diagnosis of the magnetospheric dynamics.

Abstract. In order to contribute to the understanding of solar wind±magnetosphere interactions the multifractal scaling properties of high-latitude geomagnetic uctuations observed at the Thule observatory have been studied. Using the local observatory data and the present experimental knowledge only it seems hard to characterize directly the, presumably intermittent, mesoscale energy accumulation and dissipation processes taking place at the magnetotail, auroral region, etc. Instead a positive probability measure, describing the accumulated local geomagnetic signal energy content at the given time scales has been introduced and its scaling properties have been studied. There is evidence for the multifractal nature of the so de®ned intermittent ®eld , a result obtained by using the recently introduced technique of large deviation multifractal spectra. This technique allows us to describe the geomagnetic¯uctuations locally in time by means of singularity exponents a, which represent a generalization of the local degree of dierentiability and characterize the power-law scaling dependence of the introduced measure on resolution. A global description of the geomagnetic¯uctuations is insured by the spectrum of exponents f a which represents a rate function quantifying the deviations of the observed singularities a from the expected value. The results show that there exists a multifractal counterpart of the previously reported spectral break and dierent types of f a spectra describe thē uctuations in direct dissipation or loading-unloading regimes of the solar wind±magnetosphere interaction. On the time scale of substorms and storms the multifractal structure of the loading±unloading mode¯uctuations seems to be analogous to the simple multiplicative P-model, while the f a spectra in direct dissipation regime are close but not equal to the features of a uniform distribution. Larger deviations from the multiplicative model are observed when the in¯uence of the solar wind¯uctuations is examined. On this basis it is expected that an extended multifractal analysis of the singularity structure of near-Earth plasma system¯uctuations would lead to

Introduction
After decades of eort the understanding of fully developed turbulence occurring in geophysical, astrophysical or laboratory¯ow systems is still far from complete. One of the promising trends is represented by hydrodynamic or MHD models describing turbulent intermittency in terms of inertial range canonical cascades of various nonlinearly conserved¯uxes. An adequate description of the deviations from Gaussian statistics in intermittent ®elds, however, requires higher order moments to be investigated. The departure from Gaussianity at small scales and the scaling of higher order moments is usually interpreted in terms of multifractal processes (Frisch, 1995;Biskamp, 1993). A widely used approach in solar wind studies is based on q-th order structure functions computing ensemble averaged q-th powers of absolute increments from velocity, magnetic, temperature, etc. ®elds (Burlaga, 1991;Carbone, 1994;Tu et al., 1996), which typically exhibit power-law scalings in the inertial range of scales. Another approach is based on the examination of thē uctuations of energy dissipation ®elds deriving nonnegative measures from the data ®rst . The relation between the two approaches is discussed in Frisch (1995).
As far as the intermittency of¯uctuations of the Earth's magnetosphere is concerned the structure function method was used to describe the statistical properties of collisionless plasmas in the vicinity of the Earth's bow shock (Dudok de Wit and Krasnoselskikh, 1996) and a generalization of the method known as extended self-similarity (ESS) was applied to seek for scaling laws from geomagnetic time series (VoÈ roÈ s et al., 1998). In these cases, however, due to the nonlinear dispersion eects, presence of coherent wave®elds, short time recordings or triviality of ESS in a great variety of stochastic ®elds, it was not possible to interpret the results in terms of nonlinear cascade models. The experimental structure function approach also has its own limitations (Muzy et al., 1993).
In parallel, a multifractal measure derived from AE index time series was characterized by the scaling features of its coarse-grained weight and the occurence of intermittence in``magnetospheric turbulence'' was proven on the basis of anomalous scaling of the corresponding partition function (Consolini et al., 1996). The underlying multiplicative nature of the AE index scaling seems to be well ®tted by the p-model of Meneveau and Sreenivasan (1987).
In the paper of Consolini et al. (1996) less than two months of AE index data with 1-min time resolution was analyzed. Here we analyze four years of 1-min data from Thule observatory (77:47 N, 69:23 W) using the method of continuous large deviation multifractal spectrum (LDMS). This method was tested on both synthetic and real data sets and has been proven to be more precise than the Legendre spectrum and more regular than the Hausdor spectrum (Canus et al., 1998;VeÂ hel and Vojak, 1998). We note that the LDMS method provides robust estimations of singularity spectra even in a case of non pure multiplicative processes.
Let us explain now in more traditional terms what are we looking for. Our goal is to proceed with the quantitative description of the irregular or``erratic'' nature of high-latitude geomagnetic¯uctuations. We are going to show that a positive measure, constructed from Thule observatory time series, fails to have a density and this singularity appears due to the irregularity of¯uctuations present over a range of scales. Focusing on external ®elds, the irregularity itself appears in consequence of the non-equilibrium state of the magnetosphere±ionosphere (M±I) system caused by bursty changes in solar wind forcing or induced by internal M±I instabilities. It is appropriate to mention geomagnetic substorms and storms in this respect.
Let us consider at least several examples when irregular¯uctuations observed in a high-latitude observatory may be generated in the remote regions of the magnetosphere. One of the causes of the near-Earth signatures of substorm onset is the Earthward high-speed¯ow in the neutral sheet. Shiokawa et al. (1998) reported that these¯ows undergo abrupt changes in speed at the boundary between tail-like and dipolar ®elds. It was pointed out by Daglis et al. (1999), however, that some bursty fast¯ows, often with short duration, may remain undetected on the ground because of their localized nature. Of course, there are other known sources which may cause erratic substorm appearance in a high-latitude observatory such as M±I coupling reinforced by ionospheric feedback (Kan, 1991), nonlinearity of the loadingunloading cycle (Klimas et al., 1996), etc., all intensi-®ed during bursty space weather conditions. Truly global substorms can be associated with a number of nonstationary and nonlinear plasma processes in magnetospheric tail, inner magnetosphere and ionosphere (Daglis et al., 1999). Hence, all these processes, even those in remote regions, may contribute to the spectrum of geomagnetic¯uctuations which can be observed at a high-latitude observatory. Considering that the local observatory signatures of the contributing sources may dier from each other in amplitude, frequency and other characteristics, accordingly, the resulting local signal exhibits irregular and intermittent uctuations and we are searching for generic laws which could be related to the scaling properties of these¯uctuations.
Without doubt, any kind of averaging somehow distorts the scaling properties of a signal, laying stress only upon dominant contributions from the most signi®cant sources and discarding the less representative components. It was emphasized by Coles (1990) that standard practices in geomagnetic observatories include the formation of hourly, daily, etc. averages which lead to considerable distortions of the signals. In this respect the de®nition of the widely used geomagnetic indices seems to be unsatisfactory, too. Let us consider, for example, the auroral electrojet (AE) index which was introduced to describe the global activity of the auroral zone electric currents (Davis and Sugiura, 1966). The AE-index is derived, after the substraction of base line values, from evaluation of the variations measured at 12 stations located near the northern auroral zone. From the recordings of all these stations the greatest (upper envelope) and smallest (lower envelope) values are taken at intervals of one minute and their dierence de®nes the AE index. The¯uctuations with values between the upper and lower envelopes are not taken into account at all. Again, our expectation is that the singularity spectra which characterizes the¯uctuations over a range of scales would provide a more complete description of the processes involved. The fact that we analyze here data only from one observatory can also lead to some loss of important physical information. However, we do it using the technique of large deviation multifractal spectra which is more powerful than the methods used for derivation of geomagnetic indices. As we are going to analyze the distribution of local HoÈ lder exponents (see later) over a longer period of time the description provided by multifractal spectra is statistical which focuses on long term generic characteristics of high-latitude¯uctuations.

Large deviation multifractal spectrum (LDMS)
In order to study the local irregularities of measures derived from geomagnetic data we apply the method of large deviations which provides a statistical description of a singular measure. This method can describe interwoven fractal subsets with singularities a and dimensions f g a and it also allows us to characterize irregular but otherwise arbitrary (non-fractal) signals (VeÂ hel and Vojak, 1998). A geometrical description of a multifractal measure is based on the estimation of the Hausdor spectrum (f h ) de®ned for the subsets having the same HoÈ lder exponents. Last but not least, the scaling of the q-th order moments of a measure allows us to characterize the multifractality (f l ) through Legendre transforming of the exponent functions called ReÂ nyi exponents. The latter approach was also used by Consolini et al. (1996). Generally f h f g f l . In a number of compound processes with non-concave singularity spectra f l overestimates the true dimensions while f g still provides reliable estimates (VeÂ hel and Vojak, 1998). Other methods based mainly on codimension formalism are outlined in (Schertzer and Lovejoy, 1993).
In the following we use the algorithm proposed by VeÂ hel and Vojak (1998) and by Canus et al. (1998). It is not commonly used in geophysical literature, therefore, we outline the basic assumptions of this approach. One can ®nd further details in their works. In the next paragraph we will introduce a probability measure l using a time series from the Thule geomagnetic observatory. Our intention is to provide robust estimations of singularity spectra describing the scaling properties of l. In order to introduce the continuous large deviation multifractal spectrum (LDMS) (VeÂ hel and Vojak, 1998;Canus et al., 1998) we consider P P n n!1 ; the sequence of partitions P n of the unit interval 0; 1. P n is constructed as P n I k n and 0 k 2 n À 1 1 with I k n k Á 2 Àn ; k 1 Á 2 Àn : 2 As a following step we use the concepts of HoÈ lder exponents and singularity spectra. The HoÈ lder exponent (or singularity exponent, also crowding index) is de®ned at a point x 0 Suppl as a limit (Muzy et al., 1993): where B x 0 g is the ball centered at x 0 , the size being g, and Suppl denotes the support of the measure l. If l is de®ned e.g. as a probability measure, Eq. (3) expresses the power law dependence of probability measure on resolution g with a set of HoÈ lder exponents ax 0 (TeÂ l, 1988). In our notation n 3 I as g 3 0. In a number of cases a measure can be characterized by its density An erratic behaviour appears in the absence of a density for a singular l. In such a case l decays e.g. roughly exponentially fast when g 3 0 (the size of the ball shrinks down to a point) and the exponent function ax 0 could be thought of as a generalization of the local degree of dierentiability at x 0 (Riedi and Ribeiro, 1999). For a monofractal ax const: for all x, while in a case of multifractal measure (non-uniform distribution) a changes from a point to point. For instance, the fractional Brownian motion (f Bm) or continuous Itô processes represent self-ane¯uctuations governed by a single HoÈ lder exponent (Riedi and Ribeiro, 1999). The f a singularity spectrum of a measure l associates to any given a, the Hausdor dimension of the set of all the points x 0 , which are such that ax 0 a: with the number of subsets N a with the same a (TeÂ l, 1988): Now, for every interval I of size jIj g let a g I : log lI= log g be the coarse grain HoÈ lder exponent of I. Then a Lebesgue measure is computed being the reunion of all intervals of the same size for which the coarse grain HoÈ lder exponent is equal to a HoÈ lder exponent a. The Lebesgue measure is de®ned as (Canus et al., 1998): where the set E g is It is possible to introduce an e precision for coarse grain HoÈ lder exponent by e ! ja g I À aj 8 and compute the Lebesgue measure as p c;e g Z ae aÀe p c g bdb ; 9 for which the continuous (superscript c) LDMS is (Canus et al., 1998) f c g a : 1 À lim e30 lim g30 log p c;e g a log g : 10 We have already introduced in Eq. (4) the f a singularity spectrum as a set of non-integer dimensions describing the`sizes' of interwoven fractal subsets with singularity a. Now we recall the interpretation of the same thing in terms of the large deviation principle (Riedi and Ribeiro, 1999). Here f g a (Eq. 10) represents a rate function which measures the deviation of the observed a from the``expected value'' a for which f g a 1. The smaller f g a, the fewer points display singularity a. It is obvious from the de®nition of a (Eq. 3) that the singular measure at x 0 is of Dirac d tipe for ag 3 0 0. Also, the measure is more singular for ax 0 < 1. In other words, smaller values of ax 0 correspond to the more singular measures, that is, burst of events around x 0 , while regions where events occur sparsely are less singular with higher value of a around x 0 (Riedi, 1999).

Data analysis
We study here the 1-min mean Thule geomagnetic Hcomponent time series from 1994±1997 consisting of 2 102 400 data points. There were only short data gaps ($0.04%) which were interpolated. The raw data were transformed in several steps ®rst to create the basic time series for our investigation. In the ®rst step we dierenced the data according to the relation This equation represents a kind of high pass ®ltering process necessary to set o the¯uctuations in the higher frequency range. For a certain value of s a break in scaling characteristics of geomagnetic data can be observed which is related to the spectral break at about several hours VoÈ roÈ s et al., 1998).
In turbulence studies usually the velocity dierence inside eddies is a positive and additive quantity de®ning a measure. In these studies the time structure between instants separated by s of the velocity ®eld measured in a point of the¯ow used to be related to the spatial structure (eddy size) assuming Taylor's frozen turbulence hypothesis to be valid (Frisch, 1996).
In our case the energy content of the signal is estimated by the squared value of DH: E DH 2 , then the data are normalized through Finally, the integrated quantity de®ned by j X jiT i EiDt 13 is considered to be the probability measure of the accumulated/dissipated signal energy during a time interval T ; (Dt 1). It is important to note that Eq. (12) is not substituted directly into Eq. (13) but the data normalization is performed in the ®rst step. Then the probability measure j is computed which is one data value at the end of every interval T . This is why Eq. (13) represents the aggregated signal energy accumulated during T . We can justify the introduction of the parameter T and Eq. (13) in the following manner. The majority of the bursty events in the near-Earth magnetospheric tail with clear signatures in high-latitude observatory time series appear due to the energy and mass accumulation-release processes (loading±unloading cycle) on time scales larger than 20 (min) (Bargatze et al., 1985) while the sampling interval of Thule data is 1 (min). On the other hand, no bursty energy release is expected on the time scales of <20 (min). Therefore, we expect to observe dierent multifractal scaling laws, if any, as the integration parameter T changes. Naturally, it would be possible to introduce dierent kinds of measures being more or less``sensitive'' to the various near-Earth plasma processes of interest. For example, taking a positive cubic quantity E jDH j 3 , used by a number of authors for the turbulent energy transfer rate in hydrodynamic or MHD¯ow systems, another family of measures could be de®ned. The solar wind±magnetosphere±ionosphere system, however, is not a``simple'' MHD¯ow system and we have no precise a priori knowledge about the energy transfer rates or cascade processes interconnecting multi-scale phenomena in the near-Earth magnetosphere. Therefore, by way of precaution, we use here the multifractal analysis ®rst of all as a powerful signal processing tool which could help, in analogy with turbulent¯ows, to understand the scaling/singularity features of¯uctuations allowing a better geomagnetic diagnosis of the multi-scale magnetospheric processes.
The introduction of the free parameters by the above transformations (Eqs. 11±13) needs further physical argumentation. Therefore, let us examine the eect of these transformations on power spectral density (psd). The upper curve on Fig. 1 represents a robust estimation of the averaged psd computed from 1994±1997 Thule raw data with N 2 21 points, sampling frequency f s 1=60 Hz using a blackman window of 2 17 points. The spectral break present on Fig. 1 was previously reported from AE and AL index time series (Tsurutani et al., 1990) as well as medium latitude observatory data studies (VoÈ roÈ s et al., 1998). The spectral break (or transition region) is thought to be within the range of 1/(1.5 h) and 1/(5 h). It divides the whole spectrum into the higher frequency part which is more intrinsic to the magnetosphere (smaller values of s) and the lower frequency part being in¯uenced mainly by turbulent Fig. 1. Power spectral density (psd). The upper curve represents an averaged psd for 1994±1997 Thule raw data (N = 2 21 points total, blackman window of 2 17 points, f s 1=60 Hz). The denoted spectral indices (<3) indicate that the ®eld has stationary increments. The lower curve, which has been shifted to ease comparison, represents an averaged psd of dierenced (s 20) and downsampled ( f s 1=600 Hz) data driving of the solar wind (larger values of s, . Within the transition region the external and the internal in¯uences can be mixed (VoÈ roÈ s et al., 1988) and towards higher frequencies (substorm time scale) the best predictor of the magnetospheric response is represented by nonlinear (local-linear) ®lters (Klimas et al., 1996).
For our purposes is important to recognize that if the analogy with the energy transfer in turbulent¯ows works then from the values of the power law spectral exponent b the following information can be drawn: (a) if b < 1 the data are stationary; (b) if 1 < b < 3 then a local spectrum corresponds to the stochastic intermittent processes in multiplicative energy cascade models, which means that most of the contributions to the dynamically signi®cant quantities are due to structures with neighbouring wave numbers in Fourier space (Schertzer and Lovejoy, 1993). If b < 3 the ®eld has stationary increments and the small-scale gradient ®eld will be stationary (Davis et al., 1994). Indeed, the spectral exponent for the dierenced data (lower curve on Fig. 1)  Let us examine now the eect of the transformations (Eqs. 11±13) on empirical probability distribution functions (EPDF). The whole raw data set can be characterized by scale (s) dependent EPDF (lack of global self-similarity) and the results of Hinich hypothesis testing (Hinich, 1982) we applied to the raw data reveal that Gaussianity and linearity assumptions should be rejected. Without going into the details here, we note that, if the bispectrum of a process is not zero then the process is non-Gaussian and if the process is linear as well as non-Gaussian the computed bicoherence is a non-zero constant. The EPDFs of the raw data exhibit non-Gaussianity with long-tail wings (see also VoÈ roÈ s et al., 1998), but the``mean'' or the expected value of it still characterizes well the geomagnetic data, especially during quiet conditions. The EPDFs of the integrated data exhibit similar characteristics as the EPDFs of the raw data but all values are positive with a clear absence of an intrinsic scale or a mean value. Figure 2 illustrates the power-law scaling behaviour of the EPDF as calculated from 1994± 1997 Thule data for s 10 and T 20 (min). The slope derived from a log±log scale plot of Àq is q $ 2:7 which slowly changes with s. A similar exponent (q 2:5 AE 0:2) for a large enough threshold discriminating the Gaussian part from the long-tailed one has been found previously from the EPDF plots of some geomagnetic indices and other physical processes (VoÈ roÈ s 1998). A possible interpretation of these scalings was oered by Chang (1992Chang ( , 1999 who conjectured that the dynamic nonlinear behaviour of the magnetotail being the consequence of its self-organization into criticality is responsible for the observed power laws. Self-organized criticality (SOC) models explaining the scale invariance property of global magnetospheric processes were also oered by a number of authors (Consolini, 1997;VoÈ roÈ s, 1998;Chapman et al., 1998). We emphasize that multifractal processes reach self-organized criticality through generic ®rst order phase transitions (Schertzer and Lovejoy, 1994) which was demonstrated for geo-magnetic¯uctuations in VoÈ roÈ s (1998). Nevertheless, much more eort is needed to provide unambiguous evidences for SOC in magnetospheric dynamics. In any case, the SOC model remains very promising because of its robust scaling properties that provide a possible explanation of the bursty nature of the transport in magnetosphere in terms of plasma physics (Chang, 1999;Watkins, personal communication, 1999).
The f a singularity spectrum can be calculated theoretically for a number of multifractal models. We use the analytical expression for the P-model (Halsey et al., 1986) which was introduced by Meneveau and Sreenivasan (1987) a À log 2 p w À 1 log 2 1 À p w f a À w À 1 log 2 w À 1 À w log 2 w w where w is a free parameter.
It was widely used in solar wind studies  and it was also used to explain the AE-index intermittence (Consolini et al., 1996). The P-model does not allow multifractal phase transitions. This model describes energy cascade processes in turbulent¯ows. We outline only the basic assumptions of this model. The largest turbulent eddy is assumed to be built up by a speci®c energy¯ux per unit length. Then a scaleindependent space-averaged cascade rate is considered and the¯ux density is transferred to the two smaller eddies with the same length but dierent¯ux probabilities p 1 and p 2 (p 1 p 2 1). This process with randomly distributed p 1 and p 2 is repeated again and again. Then the parameter p p 1 1 À p 2 governs the asymmetric breakdown in the fragmentation process. The value of p p 1 p 2 0:5 corresponds to the homogeneous energy transfer rate with no intermittency eects while p > 0:5 corresponds to an intermittent ow.
We do not expect that the P-model will exactly describe the energy cascading processes within the magnetosphere. The P-model is suitable, however, for a rough comparison of our results with the previously obtained values in solar wind and magnetosphere turbulence studies. Moreover, the corresponding f a spectrum is a smooth concave function with a parabolic shape like the symbol , describing measures with multiplicative rescaling structure. Any deviation from a simple shape would contain important physical information about the system under study.
To characterize the multifractal behaviour, Eq. (10) was employed to obtain the f a spectrum for dierent values of s and T . The f a spectrum was examined for ®xed values s, T 6, 60, 600 (min) respectively. When s was ®xed T was changed and vice versa. Figure 3a shows the dependence of singularity spectrum on s, while T is ®xed to 6 (min). For the sake of perspicuity, the parameter s has been changed from 10 to 100 (min) by 10 (min) steps, and from 100 to 1000 (min) by 100 (min) steps. All the curves are also projected to the (a; f a; s 3000 min) plane. The point symbols refer to the time lags s 300 (min), while the other curves correspond to the lags 1000 ! s ! 300 (min). It is discernible that these two subsets of curves dier, namely for the ®rst set (s 300 min) there is a clear increasing part (as a increases from right to left), having practically no decreasing part. Vice versa, the second set of the curves (s > 300 min) has a larger decreasing part. Similar``left-sided'' and``right-sided'' spectra were also useful in computer network studies for distinguishing incoming and outgoing tracs (Riedi and VeÂ hel, 1997) on the basis of the multifractal technique. On Fig. 3a, there is a dierence between the expected values of a, too. as 300 1:218 AE 0:001 and as > 300 1:241 AE 0:007, respectively. The probability of the singularity exponents within the interval a P 0:8; 1:6 is close to 1, and it is not possible to ®t a P-model to this distribution. For a uniform distribution ax 1 for all x, therefore the observed curves on Fig. 3a are close but not equal to a uniform l. For s > 1000 (min) and s 300 (min) the spectra are similar (not shown). From these facts the conclusion can be drawn that, as s changes, the left-sided rate function turns into the right-sided one at about 100Ä300 (min), that is there, where the spectral exponent changes (see the break point on Fig. 1). There were two spectral exponents needed for describing the power-laws in psd. In contrast with it a whole spectrum of singularity exponents is required for more complete description of the geomagnetic¯uctuations. This``additional'' information, however, is hidden for second order statistics (psd, autocorrelation).
The ®rst shaped curve depicted with points on Fig. 3a refers to the parameters T 60 and s 20 (min) [actually s is shifted to the (a; f a; s 0 min)] plane and projected to the (a; f a; s 3000 min) plane, too. This curve is added to Fig. 3a in order to show how the spectrum changes for a larger value of T . To this curve the multiplicative P-model has been ®tted using the Levenberg±Marquardt nonlinear algorithm (Press et al., 1996) resulting in p 0:700 AE 0:005. In this case, the rate function decreases more rapidly, and the deviations of the observed a-s from a $ 1:19 are smaller than in previous case for smaller T .
The simple multiplicative P-model, however, fails to describe precisely the cases when T fixed is 60 (min), but s ) 20 (min). Figure 3b shows these cases with a step of s as earlier. The point symbols refer to the time lags s 100 (min) and the other curves correspond to the lags 1000 > s > 100 (min). The projection of the curves to the (a; f a; s 3000 min) plane reveals two dierent subsets of curves, well visible especially for larger values of a. The projections of these two subsets of curves are Fig. 3a±c. Large deviation multifractal or singularity spectra as derived for a range of the parameters s and T . The individual curves are projected to the (a; f a; s 3000) plane. a T fixed 6 min; the shaped curves (point symbols, shifted to the plane s 0; 3000 min) correspond to the values T 60; s 20 min and the ®t (line) is the P-model with p 0:700 AE 0:005. b T fixed 60 min; the subset A comprises the curves for s P 30; 100 and B for s P 300; 1000 (min); c T fixed 600 min marked by squares as well as capital letters A and B. The expected value is a 1:110 AE 0:005 for the A and a 1:13 AE 0:01 for the B family of curves. The curves exhibit the same qualitative behaviour as the spectra on Fig. 3a, namely, there is a transition from A to B type of curves at about 100Ä300 (min), which is the multifractal counterpart of the spectral break. As described earlier, the second curve from the front, for s 20 (min) on Fig. 3b, can be modelled by the P-model taking p $ 0:7. It is clearly visible that for the restricted range of a P 0:7; 1:4 the f a spectra exhibit the parabolic, shape and the P-model represents a good approximation of the observed curves with eective parameters p $ 0:69 for the A and p $ 0:71 for the B group of curves. Outside this range the f a spectra are not exactly shaped and larger deviations from a simple multiplicative model occur. The easiest way of modelling this property is to consider sums of multiplicative measures of disjoint supports in which case the large deviation multifractal spectrum is simply the maximum of the individual spectra (Riedi and VeÂ hel, 1997;Radons and Stoop, 1996).
The non-parabolic shape of f a in our case gives some evidence for the phenomenon of phase transition. At the a values where the f a spectrum is out of parabolic shape the major contributor to the observed singularities changes from one measure to another. As dierent physical processes may generate dierent measures (distributions), possible models with similar characteristics as the observed spectra may contain physical information on the contributing remote magnetospheric sources. For example, the non-parabolic shape of the less singular wing (a > 0:8 of the B-family of curves in Fig. 3b could be described by a P-model with p $ 0:735, but it already overestimates the spectra at the more singular wing (a < 0:8), that is, the observed spectra are too complex and sums of simple multiplicative measures usually fail to describe fully the observed phenomena. Nevertheless, exactly the rich structure of the large deviation spectra for changing parameters T and s may reveal new types of relationships between the geomagnetic data and remote magnetospheric or solar wind processes.
For larger values of T the dependence of large deviation spectra on s does not change substantially. Figure 3c shows the curves for T fixed 600 (min). As Eq. (13) downsamples the signal, the dierences between the Fig. 3b and c can be attributed to the decreased sample size for large T .
Figures 4a±c shows how the spectra change as T increases for some ®xed values of s. The parameter T changes from 6 to 18 (min) by 2 (min) steps (depicted by point symbols), from 20 to 100 (min) by 10 (min) steps (lines) and from 100 to 1000 (min) by 100 (min) steps (lines). First, the transition region at about T 100 Ä 300 (min) is absent. However, a sudden change in the shape of the curves is present at about T 15 (min). This change is also visible in Fig. 3a, where the spectra for T fixed 6 (min) are compared with the singularity spectrum computed for T 60 and s 20 (min). The sudden change at about T 15 (min) can be explained in the following way. Equation (13) represents a system with a kind of memory because at the end of every interval T the aggregated value over this interval is computed. Moreover, T is changed, that is larger and larger intervals are considered. In this respect Eq. (13) resembles linear ®lters which are de®ned in discrete case as where It, Ot represent the input and output time series, H j is a linear ®lter (also weight function or impulse response). Equation (14) brings into connection the weighted inputs preceding the output by time lag j Á dt. Bargatze et al. (1985) have analyzed the magnetospheric impulse response for many levels of geomagnetic activity using the solar wind V Á B s (V ± solar wind speed, B s the southward component of the magnetic ®eld) as the input data, and the magnetospheric AL index as the output data. Two peaks in impulse response were identi®ed. One at about 20 (min) and another about an hour. The ®rst peak at 20 (min) arises in consequence of the direct dissipation of the solar wind energy within the M-I system. The second peak which corresponds to the loading-unloading (or storage-release) mechanisms arises or disappears depending on the geomagnetic activity level.
In Fig. 4a, b very dierent singularity spectra belong to the parameter values T < 15 (min) and T > 15 (min). The separation of these two regimes for s fixed 600 (min) in Fig. 4c is smaller. On the basis of the comparison with linear ®lter studies, we believe that the f a spectra for T < 15 (min) describe the singularity properties of the direct dissipation processes, while the curves for T > 15 (min) correspond to the loadingunloading regime.
Again, using the Levenberg±Marquardt algorithm the P-model has been ®tted to the singularity spectra corresponding to the loading-unloading regime (T > 15 (min) ± thick curves at T 0 and T 3000 (min) on Fig. 4a, b). The intermittency parameters are p 0:685 AE 0:005 for s fixed 6 min and p 0:710 AE 0:005 for s fixed 60 (min). The deviations from the parabolic shape are small and could be interpreted as earlier.
Lastly, let us return to the manner in which the measure used in this study is de®ned. We re-counted all the spectra as earlier, but this time the positive cubic quantity E jDH j 3 was considered instead of E DH 2 . No qualitative changes were observed. As an example, Fig. 5a shows the``cubic'' counterpart of Fig. 3b, that is, the dependence of the singularity spectra on s and T fixed 60 (min). The only dierence is that the spread of a around a is larger on Fig. 5a and the corresponding eective intermittency parameters are p $ 0:79 for A and p $ 0:81 for B family of curves. As previously, these values due to the deformed shapes of the spectra represent rough estimates valid for the most probable singularities. Figure 5b shows the cubic counterpart of Fig. 4b. Again the same qualitative behaviour, the sudden transition from direct dissipation to the loading unloading regime is clearly visible at about T 15 (min). The P-model intermittency parameter computed for the curves within the interval 1000 min ! T ! 20 min is p 0:790 AE 0:005 (thick lines at T 0, 3000 min).
On this basis we conclude that geomagnetic¯uctuations seem to be more intermittent for the measure derived from the positive cubic quantity E jDH j 3 , nevertheless, the singularity spectra exhibit the same qualitative behaviour as for the measures derived from the quadratic quantity E DH 2 . Further physical investigations and comparative studies including the analysis of the solar wind magnetosphere input±output system are needed for making decisions in favour of a proper measure.

Discussion
We presented a multifractal analysis of singularity spectra of the measures derived from high-latitude observatory time series. Quiet and disturbed periods of geomagnetic activity were not separated. There exist several reasons for such an approach. First of all, our goal was to understand the long term generic dynamical behaviour of``geomagnetic turbulence''. A separation of quiet and disturbed periods could destroy the global character of intermittence. In fact, as the results obtained by Consolini and De Michelis (1998) suggest, the Earth's magnetosphere regardless of its closed or open con®guration, always exists in a non-equilibrium dynamical con®guration characterized by non-uniform and spotty dissipation in time.
Our intention was to prove the eectiveness of the multifractal technique as a powerful signal processing tool for describing the singularity spectra of highlatitude geomagnetic¯uctuations. To this end a measure derived from the raw geomagnetic data was introduced and its scaling properties and bursty structure were examined using the LDMS. In analogy with spectral and linear ®lter studies, the multifractal counterparts of the spectral break and of the transition from direct to the loading±unloading dissipation regimes were identi®ed, namely for the measures derived from quadratic (Figs. 3a±c; 4a±c) and for the measures derived from positive cubic quantity (Fig. 5a, b). In this respect we note, that the AE index of geomagnetic activity is also known to be a compound index which mixes driven and Fig. 5a, b. Singularity spectra for the measure derived from E jDH j 3 . a T fixed 60 min; b s fixed 60 min. The thick curve on b (shifted to T 0; 3000 min represents the P-model ®t with p 0:790 AE 0:005 loading-unloading eects (Kamide and Baumjohann, 1991).
In addition, it was shown that within the considered interval of parameters, a whole spectrum of the singularity exponents a with probabilities of occurrence (rate function) f a (or dimensions of sets with equal a) is required for a more complete characterization of the geomagnetic¯uctuations. In a number of models the local HoÈ lder exponent can only take a ®nite number of values. For example, a simple fractal model, the fractional Brownian motion (f Bm) contains a simple local HoÈ lder exponent which rules the long range dependence of the¯uctuations. In fact,  have introduced, mainly on the basis of second order statistics, an f Bm type model, against which the hypothetical low dimensional nonlinear (chaotic) behaviour of the magnetosphere should be tested. They have found a clear resemblance between the geomagnetic AE index data and bicoloured noise with two exponents. The multifractal analysis, however, provided extended characteristics of the geomagnetic data which are hidden for second order statistics. The results here show that the local HoÈ lder exponents are time-dependent and the bicoloured noise model is not relevant.
Multifractal models are used with success in turbulence studies. The simple P-model describes how energȳ ux can be transferred between scales with multiplicative rescaling structure. In analogy with turbulence, but with the mentioned restrictions in the mind, the P-model has been ®tted to the observed f a spectra.
First of all, the P-model is not suitable for the description of the singularity spectra of¯uctuations in direct dissipation regime when T < 15 min.
The P-model is close to the observed curves for higher frequency (smaller s)¯uctuations in loadingunloading regime 1000 min ! T > 15 min. Thesē uctuations are presumably of magnetospheric origin corresponding to the physical processes on the time scale of substorms and storms. The estimated intermittency parameters depend on s and on the way how the measure is derived. For example p 0:685 AE 0:005 for s 6 min, T > 15 min and E DH 2 (Fig. 4a), while p 0:790 AE 0:005 for s 60 min, T > 15 min and E jDH 3 j (Fig. 5b).
As s increases (lower frequency range ± mainly of solar wind origin) the observed curves start to exhibit non-parabolic shapes and these deviations from a multiplicative process could be partially, but presumably not fully, explained by sums of multiplicative measures. We conjecture that this notion could be helpful in identi®cation of the sources of the geomag-netic¯uctuations. The estimated eective P-model parameters are roughly between 0:7 Ä 0:8, e.g. p $ 0:735 for a > 0:8, B family of curves in Fig. 3b; p $ 0:81 for a > 0:8, B family of curves in Fig. 5b. The former is close to the value obtained from the multifractal analysis of AE index data pAE 0:746 AE 0:002 (Consolini et al., 1996). In addition, in high-speed Helios-2 solar wind streams the analysis of f a spectra for turbulent energy transfer rate led to the value p $ 0:87 for the time shift of s 81 (s) and p $ 0:74 for s 684 (s) . Structure function analyses of low-speed Helios-2 bulk velocity data recorded in the inner heliosphere (0.33±0.94 AU) provide a best ®t P-model parameter which is 0:714 < p < 0:743 with no clear radial evolution with heliocentric distance (Carbone et al., 1996). According to Tu et al. (1996) the intermittency parameter computed from the solar wind velocity components concentrates in the range 0:7 < p < 0:8, while for the AlfveÁ n velocity components the p-values are distributed widely from p 0:64 to 0:89. Near the current sheet at 1 (AU) for s P 81; 1000 (s) p $ 0:88 .
Naturally, it would be thoughtless to over-estimate the coincidences between the corresponding parameters describing the intermittency in the solar wind¯ow and the singularity structure of high-latitude geomagnetic ®eld¯uctuations. In the mentioned solar wind studies the P-model parameter is usually obtained for the range of s P 100; 1000 (s) while here the time shift s changes from 10 to 1500 (min). The nature of spatial¯uctuations in the solar wind is usually analyzed invoking the Taylor hypothesis supposing that turbulence is convected past the sensors at nearly constant speed of several hundred km/s. The application of the Taylor hypothesis within the magnetosphere is limited (Dudok de Witt and Krasnoselkhikh, 1996), besides, the appearance of geomagnetic ®eld¯uctuations is connected with a formation of structures, magnetosphere-ionosphere coupling, etc., with characteristic speeds which are very dierent from those observed in the solar wind. For example, during a substorm expansion phase the plasma convection speed toward the equatorial plasma sheet is $13 (km/s) (Kan, 1991) which is much less than the average velocity of the solar wind. Therefore, we expect that the interaction of the solar wind plasma with the geomagnetic ®eld would change the singularity structure of the penetrating solar wind¯uctuations substantially. This question seems to be analogous with the problem of the penetration of upstream waves into the inner magnetosphere leading, after several transformations, to the generation of quasi-periodic Pc pulsations. The understanding of the transformations, of wave-like structures has led to improved diagnostics of the MHD processes and structures in near-Earth plasma system, though it remains an open ®eld of research for the future. Our expectation is that a similar analysis of the singularity spectra describing solar wind and magnetospheric intermittent¯uctuations and their evolution would contribute to the understanding of the near-Earth plasma processes as well.