ROMA (Rank-Ordered Multifractal Analyses) of intermittency in space plasmas – a brief tutorial review

Intermittent fluctuations are the consequence of the dynamic interactions of multiple coherent or pseudocoherent structures of varied sizes in the stochastic media (Chang, 1999). We briefly review here a recently developed technique, the Rank-Ordered Multifractal Analysis (ROMA), which is both physically explicable and quantitatively accurate in deciphering the multifractal characteristics of such intermittent structures (Chang and Wu, 2008). The utility of the method is demonstrated using results obtained from large-scale 2-D MHD simulations as well as insitu observations of magnetic field fluctuations from the interplanetary and magnetospheric cusp regions, and the broadband electric field oscillations from the auroral zone.


Introduction
As plasmas interact nonlinearly in the electromagnetic environment of space, coherent structures are generally created. For example, in the solar wind, field-aligned magnetic flux tubes of varied sizes have been observed (Bruno et al., 2001). Other examples of coherent structures observed in space are double layers, ion holes, electron holes, drift Alfvén vortices, current filaments, stationary and convecting structures. When these (sometimes partially formed) structures interact stochastically, intermittent fluctuations are generated.
Correspondence to: T. Chang (tsc@space.mit.edu) For intermittent turbulence, one may visualize the fluctuations to be composed of many types, each being characterized by a particular fractal dimension. Two questions arise: (i) What are the different types of fractal dimensions? and (ii) How are they distributed in the turbulent medium? Recently, a new technique of analyzing intermittent fluctuations has been developed to specifically address these questions . The technique, known as Rank-Ordered Multifractal Analysis (ROMA), retains the spirit of the traditional structure function analysis and combines it with the idea of one-parameter scaling of monofractals.
We briefly review this method below in terms of largescale two-dimensional (2-D) magnetohydrodynamic (MHD) simulations and then apply it to the study of several observational examples in the space environment.

Fractals
We begin with an elemental introduction of the concept of fractals. In ordinary Euclidean geometry, the area of a two-dimensional square may be expressed as the square of the length of one of its sides, and the volume of a threedimensional cube is the cube of the length of one of its sides. Generalizing, the volume V of a hyper-dimensional cube of dimension d is V = L d , where L is the length of one of the sides of the hypercube. Consider now the fluctuations δX of some physical quantity X of a turbulent plasma. We may construct some convenient scale-dependent measure µ of the fluctuations -equivalent to the volume of the hypercube -and study how the measure varies with the scale size δ -equivalent to the length of one side of the hypercube.   If the variation between µ and δ follows that of a power law, µ ∼ δ d f , then the value d f may be visualized as the dimension of the fluctuations for the particular chosen measure µ.
If the obtained dimension of the measure of the fluctuations may be understood in terms of straight forward dimensional analysis, one visualizes the fluctuations fully occupying the space in which they reside. On the other hand, if the obtained dimension of the fluctuations differs from that surmised from dimensional analysis, the fluctuations do not fully occupy the space in which they exist. Such a dimension has been termed a "fractal" or fractal dimension by Mandelbrot (1977).
We generalize this idea to the study of the probability distribution functions (PDFs) of fluctuations. To be more specific, we refer to the results of a 2-D MHD simulation for homogeneous turbulence (Chang et al., 2004). In the simulation, ideal compressible MHD equations expressed in conservative forms were solved numerically in a doubly periodic (x,z) domain of length 2π in both directions with 512 by 512 grid points ( Fig. 1) or 1024 by 1024 grid points (Fig. 2) using the WENO code (Jiang and Wu, 1999) so that the total mass, energy, magnetic fluxes and momenta were conserved. The initial conditions consisted of random magnetic and velocity fields with a constant total pressure for a high beta plasma. (We note that the ideal MHD equations are scale independent. They can be made dimensionless with the specification of the units of length, mass and time. For reference, the initial parameters in the simulation for the maximum perturbation velocity v and magnetic field B were chosen such that the sound wave and Alfvén wave traveling times through a distance 2π are about 4.4 and 60, respectively.) After sufficient elapsed time, the system evolved into a set of randomly interacting multiscale coherent structures exhibiting classical aspects of intermittent fluctuations. Spatial values of the magnetic field energy density B 2 were collected over the entire (x,z)-plane at a given time and PDFs P (δB 2 ,δ) were constructed for the spatial fluctuations (left panel). The PDFs are seen to be non-Gaussian and become more and more heavy-tailed at smaller and smaller scales.
We now attempt to understand the scaling behavior among the family of PDFs of the turbulent fluctuations in analogy to the idea of fractals. We assume that the δB 2 's and the PDFs vary with δ as power laws: δB 2 /δ a = I and P /δ b = J , where (a,b) are the exponents and (I , J ) are constants -i.e. invariants with respect to the scale δ. If, in addition, the form of P (δB 2 ,δ) is invariant with respect to the changes of the scale δ, then there must be a functional relation between the two invariants, (I , J ), (Chang et al., 1973). Therefore, along with the normalization condition for the PDFs, we obtain the following scaling form for the PDFs: where s = a = −b is the lone scaling exponent, and P s is a scaling function of the invariant δB 2 s = δB 2 /δ s ≡ Y .
The scaling exponent s may be interpreted as the fractal (monofractal) measure for (1). If the PDFs are self-similar Gaussian distributions for all scales representing random diffusion, the scaling exponent s is known to be equal to 0.5. For other monofractal distributions, the scaling exponent may take on any positive real value. An attempt to collapse the unscaled PDFs of Fig. 1 (left panel) according to expression (1) indicated approximate scaling with an estimated scaling exponent s = 0.34, Fig. 1 (right panel). The scaling was quite robust for small values of δB 2 s , but became less and less satisfactory as the value of δB 2 s increased. When the scaling law (1) is fully satisfied for the full range of the scale invariant, the scaling of the PDFs is called self-similar or monofractal. Otherwise, the fluctuating phenomenon is multifractal.

Multifractals and ROMA
One established procedure for measuring the multifractal nature of fluctuations is by studying the moments of the PDFs, S m -commonly called structure functions. If the scaling exponents ζ m for S m vary linearly with the moment order m, one visualizes that the various portions of the PDFs representing the different sizes of the fluctuations have the same fractal characteristics and therefore all the PDFs collapse onto one scaled PDF according to (1). In fact, one may easily demonstrate that PDFs satisfying (1) will have the monofractal property ζ m = sm, and s has been called the Hurst exponent.
For multifractal fluctuations, the fluctuations within a small differential range dY about the value Y may be assumed to exhibit approximate monofractal behavior. Thus, the differential structure function dS m will vary with the scale as δ sm within the differential range according to (1), i.e., dS m δB 2 m P δB 2 ,δ d δB 2 = δ sm Y m P s (Y )dY (2) The solution s(Y ) for the above functional differential expression (if the ansatz is valid) may be obtained approximately by integrating (2) over small contiguous ranges of Y with the assumption that within each incremental range the scaling exponent s is approximately constant. Figure 2 shows the results of such a rank-ordered spectrum using the numerical data from a MHD simulation similar to that discussed in the previous section. Because of the approximate symmetry of the PDFs with respect to δB 2 , we have defined Y ≡ δB 2 /δ s . We note from Fig. 2 that the scaling exponent s appears to approach the self-similar value of 0.5 representing Gaussian random diffusion as Y → 0. As the size of the scaled fluctuation Y increases the scaling exponent decreases indicating that the fluctuations are becoming less and less space-filling. This method of analyzing intermittent fluctuations, ROMA (rank-ordered multifractal analysis), was first suggested by .

Application to space plasma observations
Recently, ROMA has been applied to the analyses of turbulence data observed within several regimes of the space environment. We briefly review these results below.

Solar wind
Most of the turbulence data in space came from measurements obtained by one single spacecraft. For solar wind, the collected data are generally analyzed based on the Taylor hypothesis (1938) which assumes that the transit time of eddies (bundled fluctuations) is much less than the characteristic evolution time. For fully developed MHD turbulence, the characteristic eddy turnover (evolution) time may be estimated as where λ is the typical size of the eddy, v A is the Alfvén speed based on the mean magnetic field B 0 and δB is the average fluctuating magnetic field of the eddy in question. The transit time of the eddy is λ/v s where v s is the relative streaming speed between the plasma medium and the spacecraft. Thus, as it has been demonstrated by Matthaeus and Goldstein (1982) that the requirement of t(evolution) be much larger than t(transit) becomes ν S /ν A 2π (δB/B 0 ). This condition is generally satisfied for solar wind turbulence.
Because the spacecraft speed is much smaller than the solar wind speed, the data gathered by instruments aboard the spacecraft essentially sample the spatial fluctuations in the direction of the solar wind, particularly at times when the solar wind turbulence is homogeneous at large scales and fully developed. One may then construct, for example, the probability distribution function P (δB 2 ,τ ) with τ being the time scale of the fluctuations in analogy to the analysis of the simulation results discussed above.
Based on a particular set of WIND solar wind magnetic field data, Hnat et al. (2002) demonstated that the PDFs of the magnetic field energy density fluctuations P (δB 2 ,τ ) nearly collapsed onto a single scaling curve according to (1), indicating that the observed solar wind turbulence was essentially monofractal.
Using a different set of WIND data,  obtained a ROMA spectrum for the observed solar wind turbulence, Fig. 3. The spectrum had the signature of that of Fig. 2. However, it was much flatter, i.e., the range of values of the fractal index s (roughly between 0.44 and 0.38) for the full observed range of the scale invariant Y ≡ δB 2 /τ s was quite small. In fact, it was found that the observed PDFs nearly collapse unto one single curve according to (1) by setting the fractal index s as a constant (0.44), Fig. 4, quite similar to the result obtained by Hnat et al. (2002). We now proceed to consider turbulent fluctuations within other space environments where multifractal structures are more evident.

Magnetospheric cusp turbulence
One of its regions of the magnetosphere most exposed to the solar wind is the magnetospheric cusp. It has been demonstrated by Echim et al. (2007)  There was strong indication that the statistical properties of the fluctuations in the cusp are inhomogeneous and did not represent fully developed turbulence. How inhomogeneity may be included explicitly and more rigorously in a ROMA analysis is discussed in the concluding section.
A demonstration of the validity of the Taylor hypothesis is not straightforward for the cusp region. In Echim et al. (2007), we provided an argument for the justification of the hypothesis to the data set that was analyzed. Firstly, Lundin et al. (2003) found that the plasma bulk velocity is of the order of 100-200 km/s for the same Cluster data as those discussed in this paper. Thus the plasma speed took on values larger than the Cluster velocity indicating that temporal variations could be transformed into spatial variations. Secondly, a correlation analysis was performed. Both the dynamic linear autocorrelation factors for the B-field fluctuations on-board each Cluster satellite and the dynamic cross correlation factors for the same B-field component (or amplitude or B 2 ) between pairs of Cluster spacecraft were computed. It was found that the dynamic autocorrelation coefficients are similar to the cross correlation factors, which could be interpreted as a verification that temporal fluctuations correspond to spatial fluctuations (Weygand et al., 2005).
With this assumption, a ROMA analyses had been performed for the chosen data set in the cusp region, Fig. 6 (Lamy et al., 2008). When the turbulence is fully developed, the ROMA spectrum generally has a smooth monotonic behavior. On the other hand, if the turbulence is unstable and not yet fully developed, its ROMA spectrum may vary rapidly over small ranges of Y . Under such circumstances, the increment Y in the ROMA analysis need to be adjusted sufficiently small to avoid the spurious occurrence of double roots. Since such circumstances occur generally for small Y values, there usually will be sufficient number of data points to determine the spectrum values even for relatively small Y increments (Tam et al., 2010). We have employed this procedure of variable Y in determining the ROMA spectra in this and the following section.
We note that for small values of Y , the fluctuations were persistent (s > 0.5), indicating the turbulence was unstable and probably not yet completely fully developed. For larger values of Y , the fluctuations became anti-persistent (s < 0.5) and the turbulence was probably well-developed and became sparser and sparser as the value of Y increased similar to the MHD simulation result. The ROMA spectra for all 4 spacecraft were very similar; indicating the turbulence there was essentially isotropic.  Tam et al. (2010) have applied the ROMA technique to analyze the multifractal characteristics of the auroral zone electric field fluctuations observed by the SIERRA (Sounding of the Ion Energization Region: Resolving Ambiguities) rocket. Since the dominant fraction of the electric field variations in the observational region were known to be nearly non-propagating intermittent turbulence composed mainly of multiple interacting inertial Alfvénic potential structures, the observed time series were interpreted in terms of spatial fluctuations with ≈ U τ , where U ≈ 1.5 km/s is the horizontal speed of the rocket, τ the time scale and the corresponding spatial scale.

Auroral electric field fluctuations
In analyzing the data, it was discovered that the observed fluctuations (commonly called the broad band extremely low frequency or BBELF fluctuations) span across four contiguous regimes of scales (from kinetic to inertial Alfvénic, to intermediate, and eventually to the MHD scales) with different multifractal characteristics. A ROMA spectrum was found for each regime, s i (Y i ) for i = 1,2,3,4. The dividing times between the adjacent regimes were determined to be roughly τ ≈ 80 ms, τ ≈ 160 ms, and τ ≈ 320 ms.  Thus, the observed intermittent turbulence of the electric field fluctuations were found to be characterized by two rankordered parameters: the index i and the power-law scaling variable Y i . Assuming the nonlinear crossover regions between the adjacent regimes were narrow, the ROMA spectra were mapped onto one global scaling variable with a single global scaling function, Fig. 7. Such a mapping served the purpose of demonstrating the ROMA spectrum as multivalued and provided a convenient display of the individual multifractal characteristics for each rank-ordered index i. It also demonstrated the generalization of the idea of monofractal mapping (1) to the concept of multi-parameter mapping involving crossover effects among multifractals. 5 Summary, conclusion and speculative remarks 1. We briefly reviewed the ROMA technique of analyzing multifractals. The method was developed using the results of large scale 2-D MHD simulations. Applications of the method to the interplanetary medium, the magnetospheric cusp, and the auroral region were briefly described.
2. Such an implicit multifractal spectrum has several advantages over the results obtainable using the conventional structure functions. Firstly, the utility of the spectrum is to fully collapse the unscaled PDFs. Secondly, the physical interpretation is clear. It indicates how intermittent (in terms of the value of s) the scaled fluctuations are once the value of Y is given. Thirdly, the determination of the values of the fractal nature of the grouped fluctuations is not affected by the statistics of other fluctuations that do not exhibit the same fractal characteristics. Fourthly, it provides a natural connection between the one-parameter monofractal scaling idea (1) and the multifractal phenomenon of intermittency.
3. Most of the existing analyses on multifractals are in terms of either spatial or temporal fluctuations. On the other hand, fluctuations that exist in the real world are generally spatiotemporal in nature. In an interesting paper, Pulkkinen et al. (2006) discussed the scaling of the second order structure functions that are intrinsically spatiotemporal and monofractal or self-affine.
In principle, ROMA may be generalized to consider the scaling of multifractals in a similar manner. For monofractal processes, the underlying stochastic equation may generally be expressed in terms of a Langevin equation with Gaussian noise, or equivalently a Fokker-Planck equation. For multifractal fluctuations, the underlying stochastic equation would be a spatiotemporal Langevin-like equation driven by a generalized noise. And the solution of such an equation would require the meticulous application of the methods of path integrals and the nonlinear dynamic renormalization-group (Chang et al., 1978(Chang et al., , 1992(Chang et al., , 2004. The results of a multiparameter ROMA study with nonlinear crossover effects based on the concept of spatiotemporal fluctuations can shed light to the description of such a theoretical analysis. These studies are currently being pursued by the authors.
4. Traditional multifractal analyses are based on the statistics of incremental changes. An alternative popular method of complexity analysis of fluctuations, the concept of self-organized criticality (SOC) (Bak et al., 1987), considers the statistical description of avalanches structured around the ideas of threshold cutoffs and finite-size scaling effects. Recently, Uritsky et al. (2007Uritsky et al. ( , 2009a analyzed an extended data set of extreme ultraviolet images of the solar corona using both methods. The data exhibited SOC power-law avalanche statistics as well as multifractal properties of structure functions for spatial activity. It will be interesting to investigate if this simultaneous existence of SOC and multifractal phenomenon is universal for dynamical complexity processes or if there are multifractal processes that do not possess the hallmark of SOC behavior. Threshold cutoffs and finite-size effects as well as power-law scaling are related to the effects of different invariant sizes. Since ROMA, by construction, can differentiate monofractal and multifractal behavior in terms of the invariant scale sizes, it is a method tailor-designed for the determination of this intriguing question related to complexity processes. The authors and their collaborators have been actively pursuing the answer to this important intellectual question using ROMA along with the method of the dynamical renormalization-group.
5. Space plasmas are inhomogeneous at all scales. Examples are found in the auroral zone (Uritsky et al., 2008(Uritsky et al., , 2009b, the magnetospheric cusp as indicated above and elsewhere. As demonstrated by the calculations described for the auroral electric field fluctuations, the rank-ordering idea may be generalized to the phase space of multiple parameters representing the physical characteristics of the real space environment such as inhomogeneity. The authors are currently applying such an idea to the cusp data observed by Cluster.