Multifractal structure of turbulence in the magnetospheric cusp

Magnetospheric cusps are regions which are characterized by highly turbulent plasma. We have used Polar magnetic field data to study the structure of turbulence in the cusp region. The wavelet transform modulus maxima method (WTMM) has been applied to estimate the scaling exponent of the partition function and singularity spectra. Their features are similar to those found in the nonlinear multifractal systems. We have found that the scaling exponent does not allow one to conclude which intermittency model fits the experiment better. However, the singularity spectra reveal that different models can be ascribed to turbulence observed under various IMF conditions. For northward IMF conditions the turbulence is consistent with the multifractal p-model of fully developed fluid turbulence. For southward IMF experimental data agree with the model of non-fully developed Kolmogorov-like fluid turbulence.

have found that fluctuations at different frequencies form wave trains suggesting multiscale, intermittent processes operating in the TBL.They have shown that the three-wave, nonlinear interaction is responsible for strongest wave trains.However, magnetic turbulence in TBL is dominated by, albeit weak, random-like fluctuations with smooth, continuous spectrum, which cannot be adequately studied with the bi-spectral, or even tri-spectral analysis.In this case we have to resort to more sophisticated methods of data analysis, which make use of the higher order statistics.
Experiments provide ample and unquestionable evidence that turbulence is not adequately and fully specified by spectral analysis alone (e.g.Frisch, 1995;Paladin and Vulpiani, 1987).The power spectrum, being related to a second moment of the probability distribution function (PDF), fully describes fluctuations if they have a Gaussian PDF.In this case, the turbulence is scale-invariant and self-similar.In the intermittent turbulent media, the PDF is increasingly non-Gaussian at smaller scales, turbulence is no longer scaleinvariant, and higher order moments are needed to characterize properties of the fluctuating field.
Intermittency is observed in turbulent fluid flows (e.g.Frisch, 1995), as well as in the magnetohydrodynamic turbulent media, such as a solar wind (e.g.Marsch and Tu, 1997;Horbury et al., 1997;Sorriso-Valvo et al., 2001).Several intermittency models have been proposed.The multifractal pmodel by Meneveau andSreenivasan (1987, 1991) seems to reproduce fluid experiments better than other models.Carbone (1993) adopted the p-model to the case of developed MHD turbulence.Motivated by the MHD turbulence observed in the solar wind, this model has been extended to the case of still developing, evolving turbulence (e.g.Tu et al., 1996;Marsch and Tu, 1997).The model successfully describes the power spectrum evolution in high-speed streams.Unlike in the p-model, which depends on a single intermittency parameter, the so-called extended structure-function model introduces a second parameter characterizing the scaling properties of the space-averaged cascade rate.In this paper we analyze the scaling properties of turbulence in the polar cusp of the magnetosphere, as revealed by magnetic field fluctuations measured on board the Polar satellite.In Sect. 2 the method of analysis is described.We used the wavelet transform modulus maxima method, which was reviewed, for instance, by Muzy et al. (1994) and Arneodo et al. (1999).The method allows one to determine the scaling properties of measured fluctuations for negative moments.Results of analysis are discussed and compared with the models of energy cascade for turbulence in Sect.3. In Sect. 4 we present the physical interpretation of results and conclusions.

Method of analysis
The multiscaling properties of fully developed turbulence are conventionally investigated by calculating the q-th order structure functions of a measured fluctuating parameter g(x) where L l is the length of the signal.It should be noted that the structure function exhibits exponential scaling only if the turbulence is "locally" self-similar at a scale near l.
By Legendre transforming the scaling exponents ζ (q) of the structure functions one can obtain the singularity spectrum D(h), defined as the fractal dimension of the set of points with the Hölder exponent h.According to Parisi and Frisch (1985) D(h)=min q (qh−ζ (q)+1).To quantify the fractal properties of a function one should find a set of locations of the singularities and estimate the values of the Hölder exponent h.
As has been shown by Muzy et al. (1991)(for a comprehensive review, see Muzy et al., 1994), the wavelet transform is a very suitable tool to detect singular behavior of self-similar functions.In the wavelet transform, one approximates a function g as a sum of properly weighted basis functions: where a is the scale (or inverse frequency), b is the dilation or translation parameter.For our purpose we use L 1 norm wavelet transform and a real valued transforming function ψ called mother wavelet.The mother wavelet is chosen to be well localized in both space and frequency.It is also required that ψ have a certain number of vanishing moments.For instance, the N -th order derivative of the Gaussian function has N vanishing moments, while the Haar wavelet, which is the equivalent wavelet for the structure function, has only one vanishing moment.The wavelet transform can be considered as a decomposition of the function g into space-scale contributions.Mallat and Hwang (1992) have shown that singularities in g produce a maximum in the modulus of the wavelet transform coefficients and that at small scales this maximum gives a location of the singularity in the signal.Muzy et al. (1994) call modulus maximum of the wavelet transform, "any point (x 0 , a 0 ) of the space-scale half-plane which corresponds to the local maximum of the modulus of T ψ [g](x, a 0 ) considered as a function of x", and the curve which connects the modulus maxima is the maxima line.The coefficients of these maxima (which are a small fraction of the total number of coefficients) are enough to encode the information contained in the signal.This allows for the calculation of the singularity exponents by a power law fit of the wavelet coefficients along the maxima line.This approach, introduced by Mallat and Zhong (1992), is called a wavelet transform modulus maxima (WTMM) method.
The main disadvantage of using the structure function method in characterizing the singular structure of a function is that it often diverges for q<0.In the WTMM method, in order to avoid this effect, at a given scale, one calculates the partition function as a sum over local maxima of the modulus of the wavelet transform (Muzy et al., 1991).The waveletbased partition function is defined as: where L(a) is a set of all the maxima lines l existing at a scale a, and b l (a) is the position, at a, of the maximum belonging to the line l.Each line l={b l (a), a} is pointing (when a goes to 0) towards a point b l (0) which corresponds to a singularity of g.Because one does not sum over places where the wavelet modulus is zero, the partition function is also defined for q<0.If g(x) is self-similar, then along the maxima line the partition function behaves like: Z (q, a) ∼ a τ (q) .(4) The singularity spectrum D(h) is obtained by the Legendre transform of the function τ (q): Thus the relation between the scaling exponent of the structure function ζ (q) and the WTMM exponent τ (q) is: The spectrum of the scaling exponent is an important statistical characteristic of the turbulent field.When derived experimentally it can be compared to that given by models of the turbulence.In the case of local, fully developed, isotropic fluid turbulence (Kolmogorov, 1941) the structure function scaling exponent ζ (p) behaves like p/3.In the MHD analog of the turbulence (Kraichnan-Iroshnikov theory) (Kraichnan, 1965;Iroshnikov, 1963) the scaling exponent is p/4.But the experiments on fluid turbulence and observations show that the spectrum of the scaling exponent is nonlinear.This nonlinear behavior is interpreted as an intermittency phenomenon and as a direct consequence of the existence of spatial fluctuations in the local regularity of the velocity field.
The intermittency can be simply visualized by plotting the probability functions of fluctuations for various scales.If for small scales the PDF is spiky, with stretched wings, and as scale increase, it becomes closer to Gaussian, then we say we are dealing with the intermittency.

Results of data analysis and comparison with models
Magnetic field data from NASA Polar satellite (Russell et al., 1995) are used in this study.The Magnetic Field Experiment provides three components of the magnetic field sampled with the frequency 8.333 Hz.The total magnetic field B for 9 October 1996 is plotted as a function of universal time in the top panel of Fig. 1.The measurements were taken in the northern cusp region, at MLT from 12:17 to 12:50, magnetic latitudes between 55.46 and 70.58 • , and distance to the reference magnetopause 3-4 R E .The bottom panels show the differences B(t)=B(t+ t)−B(t) for the time delays t=7, 29, 117, and 612 s and probability distribution function of B normalized to its standard deviation.
The equivalent Gaussian distribution is plotted for comparison.Time delays are chosen arbitrarily, just to show the differences in the behavior of the PDF for different delays.It is seen that the larger the time delay, the closer to the Gaussian the PDFs are.At small delays, distributions are spiky and have extended wings.The degree to which PDFs depart from the Gaussian can be quantified with kurtosis and skewness.One can see that for small delays the kurtosis is considerably different from zero, meaning that higher order statistics is needed for a full characterization of the PDFs of fluctuations and that the cusp turbulence structure appears intermittent, singular.Note, however, that the skewness differs from zero only slightly.
To find the singularity spectrum we used the WTMM method described in the previous section.The data set has been divided into seven segments, each 8192 samples long.By using the short data segments we tried to avoid the effect of non-stationarity.At this point we wish to note that we made use of the ergodic and Taylor hypotheses.The first one is necessary to replace the ensemble averages by spatial averages.Taylor's "frozen turbulence" hypothesis allows one to replace the spatial statistics with temporal statistics.According to this hypothesis, the entire spatial pattern of turbulence is transported past the probes with a constant speed.To verify experimentally the validity of Taylor hypothesis one should have access to the spatiotemporal data, for which spaced-probe measurements are necessary.For the case under consideration, the mean drift speed V 0 , calculated from the electric and magnetic field Polar data, is of the order of 100 km/s.We assume that turbulent eddies of all sizes considered here have an intrinsic speed much smaller than V 0 , and are convected with this high speed past the probe moving with the velocity ∼2 km/s.We should note that coherent structures present in the turbulent flow and associated with the velocity bursts might invalidate both ergodicity and Taylor hypotheses.In spite of that, most experimental results assume ergodicity and frozen flow, and we follow this practice being aware, however, of the problem.
The top panel of Fig. 2 shows a segment no. 5 of magnetic field fluctuations.In the middle panel of Fig. 2 the spectrogram (modulus of the wavelet transform) is given.The mother wavelet used here is the 4th order derivative of the Gaussian function.This wavelet makes possible the accurate determination of the location of singularities (Muzy et al., 1994).Long "fingers" extending from low to high frequencies are signatures of singularities.The bottom panel shows the lines of modulus maxima of the wavelet transform.
Figure 3 compares the Fourier (solid line) and wavelet (dashed line) power spectra.Both spectra agree quite well.Note that the spectrum is feature-less and continuous.We should mention, however, that segments nos. 1 and 4 exhibit a small but distinct maximum around 0.3 Hz.The spectrum seems to follow a power law dependence f −α on the frequency with α=2.36±0.04 over the frequency range 0.06-0.78Hz.The spectral indices α for other segments fall into the range from 1.87 to 2.62, with the mean 2.27.Smallest α index is observed for a segment no. Figure 4 shows the behavior of τ (q) exponents (dots) derived from WTMM for power q in the range −4 to +4.Given a relatively small number of samples in our data segments, we did not attempt to perform calculations for larger q.We note that the partition function (3) was calculated only for those lines of modulus maxima that extended in frequency over more than 2.5 octaves.To avoid the effect of noise, we considered only frequencies lower than one-fifth of the Nyquist frequency, i.e. ≈0.8 Hz.The rms error of τ is small, of the order of 10 −2 , and error bars are not discernible on the τ (q) plot.However, the test computer runs show that the error increases if shorter lines are included in the partition Fig. 4. Comparison between the experimental τ (q) (dots) and bestfit model scaling exponent.The bottom plot shows the results of fitting for τ (q) limited to q>0 and, for comparison, the scaling exponent of the structure function (squares).The error bars give the standard deviations of the least-squares fits.In spite of different parameters, curves for Kraichnan-like and Kolmogorov-like models are practically indistinguishable on this plot.function calculations and the range of q is increased.The points in Fig. 4 appear to form a curved line rather than a straight line, which means that the turbulence is indeed intermittent.
We attempted to compare experimental results with several intermittency models.The simplest multifractal cascade model is the p-model (Meneveau andSreenivasan, 1987, 1991), in which the energy cascade flux transfers from the larger eddy to two smaller eddies with the same scale size, but with different flux portions defined by randomly distributed probabilities P 1 and P 2 =1−P 1 .The intermittency parameter P 1 describes the spatial inhomogeneity of the cascade rate.For the case without intermittency P 1 =0.5, while for fully intermittent turbulence P 1 =1.For the p-model the partition function scaling exponent is given by: τ (q) = − log 2 P q/3 1 + (1 − P 1 ) q/3 .
(7) Since α=τ (2)+2, the maximum spectral index α=2 allowed by the model is achieved for P 1 =1.However, as pointed out by Tu et al. (1996), the spectral index directly calculated from the power spectrum does not need to be the same as that derived from the structure function.We wish to add that the same is true for τ (2).This is due to the fact that two methods give different weighting across time scales.For instance, in our example τ (2)=−0.02±0.02,therefore α=1.98±0.02,which differs from the power spectrum index 2.36±0.04.The dashed line in Fig. 4 represents the best-fit τ (q) for the p-model.We can see that the p-model departs from observation for the positive power q.As a quantitative measure for the goodness-of-fit we used the chi-square test (Press et al., 1986).Table 1 shows the values of P 1 with their rms error, and the probability Q that the computed fit would have a value χ 2 (the sum of squared differences between the fit and data) or greater.If Q is small, then the differences between observations and model are unlikely to be random fluctuations and the model can be rejected.Except for the data set 3, the theoretical p-model cannot be fitted satisfactorily to the observational results.Tu et al. (1996)(see also Marsch and Tu, 1997) introduced an intermittency model that applies to the turbulence not fully developed.They derived the following scaling functions for the Kolmogorov-like cascade: and for the Kraichnan-like cascade: These extended intermittency models depend on two parameters, the intrinsic spectral index α and intermittency parameter P 1 .For the case without intermittency P 1 =0.5, α=α and τ (q)= − 1+(α−1)q/2, for both kinds of cascade.For P 1 =1, we have α=α +1/3 and α=α +1/2, for the Kolmogorov-and Kraichnan-like cascade, respectively.The scaling exponent is τ (q)=(−1+α/2)q in both cases.
The least-squares fitting of experimental τ (q) to ( 8) and ( 9) gives P 1 =0.80, α =2.18, and P 1 =0.87, α =2.23, respectively.The rms errors of parameters are of the order of 10 −3 .To calculate α we used the experimentally derived τ (2), instead of the power spectral index.Both models fit data very well, except at positive q.In spite of different parameters, curves for Kraichnan-like and Kolmogorov-like models are practically indistinguishable on this plot.Best-fit results for all analyzed Polar data segments are given in Table 1.In general, we may say that a good agreement between data and intermittency models is achieved, which indicates that the basic concept of the multifractal energy cascade applies to the turbulence in the cusp region.However, using just the scaling function we cannot distinguish which of the two models, Kraichnan-or Kolmogorov-like, better describes the measurements.Tu et al. (1996) reached the same conclusion when investigating the structure function of the solar wind velocity fluctuations.We also see that the intermittency parameters P 1 obtained from the Kraichnan-like model are systematically higher than that from the Kolmogorov-like model, again confirming the Tu et al. result.The bottom plot of Fig. 4 shows the τ (q) dependence restricted to q>0.Except for the p-model, for which the chisquare probability Q=0, the fit to the theoretical curves is very good (Q=1).In Table 1, for the segment no. 5, entries in italics are the model parameters derived from the truncated τ (q).One can see that for the extended models they are very close, albeit systematically smaller, to those calculated when both positive and negative powers are taken into account.However, when attempting to calculate the singularity spectra one should include in the partition function (3) the negative values of q.If one uses only positive qs then the analysis is restricted to the strongest singularities characterized by the Hölder exponent h smaller than the most "frequent" one.
In Fig. 4 we also show, for comparison, the results derived from the structure function (squares).The τ (q) has been computed by linear least-squares fit to the doublelogarithmic plot log S q -log t shown in Fig. 5.The lower limit of the time delay t≈1 s is chosen to correspond roughly to the highest frequency used in the WTMM.The upper limit of t≈31 s assured a reasonable fit.At this point we note that since the fractal functions may have, at any scale, increments close to zero, the structure function will diverge for q<0.Thus, the structure function method does not provide a reliable generalization of the multifractal formalism to fractal functions (Muzy et al., 1994).In addition, the structure function is very sensitive to any outliers present in the data.Jaffard (1994) proved mathematically that a slightly modified WTMM method is superior to the structure function method in giving the correct singularity spectrum for all q.
From τ (q), through the Legendre transformation (5), we have derived the singularity spectrum D(h), which is a measure of the local scale-invariance.Local scale invariance means that for each h there is a fractal set with the dimension D(h) near which scaling with the exponent h holds.
The relation ( 5) can be rewritten as: The Hölder exponent h has been calculated by numerical differentiation of the function τ (q) and used to derive the singularity spectrum D(h) from the second Eq. ( 10).An alternative method in which h and D(h) are calculated directly from the wavelet transform modulus maxima (Arneodo et al., 1992), without explicitly Legendre transforming, has also been used.Both approaches have been tested on artificially generated data sets representing p-model.In Fig. 6 the derived singularity spectra are compared with the theoretical spectrum for the intermittency parameter P 1 =0.8.Errors were calculated and are shown only for the Legendre transform method.To quantify the disparity between the theoretical and numerical singularity spectra we have used the mean squared deviation 2 h of the theoretical and numerical scaling exponent h.We have found that for the singularity spectrum computed directly the disparity measure h ≈0.037, while for the Legendre transform (10) h ≈0.023.We see that both methods give quite consistent results with slightly better agreement for the Legendre transform.Therefore, in our estimation of the experimental singularity spectra we used (10), which is more straightforward and less computationally involved.
In Fig. 7 we depicted the function D(h) for our sample segment.The experimental singularity spectrum has a characteristic parabolic shape, typical of other nonlinear systems, and its support extends from 0.40 to 0.89 over the range of qs from −4 to 4 with a maximum at h max ≈0.62 and D(h max )≈1.00.
In Fig. 7, we also plotted the D(h) singularity spectra for three models considered here with parameters derived from the τ (q) dependence (see Table 1).It is clear that, in spite of a good τ (q) fit, the experimental singularity spectrum considerably departs from the models.One may observe, however, that the disagreement depends on the model.This fact was used to discriminate between extended models, which was not achievable from the best fit to τ (q).In the case of our sample data segment no. 5 the disparity measure h is smallest for the Kolmogorov-like model (cf.Table 2).Indeed, Fig. 7 shows that the Kolmogorov-like spectrum has h max and support closest to those found in the experiment.
Singularity spectra parameters for all data sets are given in Table 2.The h max varies between 0.47 and 0.62.It is interesting to note that its mean value ≈0.55 is not very much different from that found for the solar wind (Marsch et al., 1996).When varying q from −4 to 4, the widest range of h is observed for the segments nos. 2 and 3.If one uses as a criterion of the agreement between experiment and model the values, then several types of intermittency can be distinguished within the analyzed data set.
Data segments nos.2, 3, and 4 seem to conform to the p-model, and segments nos.5, 6, and 7 resemble the Kolmogorov-like model.For the data segment no. 1 it is not possible to judge which of the two extended models of intermittency fit the experimental singularity spectrum better.

Discussion and conclusions
We have investigated the scaling properties of magnetic field fluctuations as measured in the turbulent boundary layer.The wavelet transform modulus maxima (WTMM) method has been used to estimate the scaling behavior of the partition function and the multifractal spectrum of turbulence.We have found that their features are similar to those found in the nonlinear multifractal systems.The experimental scaling exponent τ (q) and singularity spectra D(h) have been compared with models of the intermittent turbulence: a) pmodel, which was introduced to describe the intermittency in the fully developed neutral fluid turbulence (Meneveau and Sreenivasan, 1987), b) extended model, which takes into account the average energy cascade rate changes with scale and simulates non-fully developed turbulence.Two versions of the extended model have been considered: Kolmogorov-like cascade in the neutral fluid turbulence, and Kraichnan-like cascade in the MHD turbulence (Tu et al., 1996).We have found that the scaling exponent does not allow one to conclude which of the two extended models fits the experiment better.However, comparison of the experimental and model singularity spectra reveals that different data segments can be described by different models.The fact that the singularity spectra better differentiate between models is apparently due to the fact that the singularity spectra, effectively dependent on the derivative (gradient) of the scaling function, are more sensitive to the model and its parameters.
The physical situation, which corresponds to different types of turbulence, is difficult to describe in detail.Figure 8, in which the solar wind Geotail magnetic field vector is plotted, will help further the discussion.One can see that the time interval corresponding to the first four data segments is characterized by variable y and z components of interplanetary magnetic field (IMF).Yet for set no. 3 B z is definitely positive (northward).For this set we have found that the singularity spectrum extends over a wide range from 0.05 to 0.68 and agrees with the p-model spectrum describing a fully developed turbulence.The corresponding intermittency parameter P 1 =0.81 is relatively high.Taylor and Cargill (2002) discussed recently plasma flows when the magnetosheath interacts with the magnetopause indentation at cusp under northward IMF conditions.They have shown that when the plasma velocity is in excess of the fast mode magnetosonic wave speed, a highly turbulent, albeit thin, boundary layer forms which enters the cusp indentation.
At times coinciding with the last three segments, the IMF is stable with large positive B y and negative B z , which suggests that Polar spacecraft senses plasma on open field lines flowing toward the magnetotail.For such configuration of IMF the reconnection in the vicinity of the sub-solar point affects the cusp structure.In this case the scaling behavior of the partition function and the singularity spectra reveal that the magnetic field has a multifractal structure compatible with the non-fully developed Kolmogorov-like (fluid) turbulent cascade.This leads to the conclusion that the turbulence is dominated by flow eddies.
The conclusions drawn here on the basis of limited data need to be confronted with results of a more complete study.

E
. Yordanova et al.: Multifractal structure of turbulence in the magnetospheric cusp

Fig. 1 .
Fig. 1.Total magnetic field B (upper panel) and differences B(t)=B(t+ t)−B(t).Right plots show the probability distribution function of B normalized to its standard deviation.The dotted curve represents the equivalent Gaussian distribution.Kurtosis and skewness are also given.The data segments used in the analysis are marked with vertical lines in the upper panel.

Fig. 2 .
Fig. 2. Segment of data shown in Fig. 1 (top panel), its spectrogram (middle panel), and the map showing locations of the modulus maxima (lower panel).

Fig. 5 .
Fig. 5.The structure function of the magnetic field shown in Fig. 2. Curves are labeled with the power q.

Fig. 6 .
Fig. 6.Singularity spectra of the p-model calculated directly from the wavelet transform modulus maxima (upper) and by the Legendre transform (lower).Errors were calculated and are shown only for the Legendre transform method.Dashed curves represent the theoretical singularity spectrum.

Fig. 8 .
Fig. 8.The Geotail solar wind magnetic field.Time intervals corresponding to the analyzed Polar data segments are marked with vertical lines.

Table 1 .
The best fit intermittency model parameters.