Double rank-ordering technique of ROMA ( Rank-Ordered Multifractal Analysis ) for multifractal fluctuations featuring multiple regimes of scales

Rank-Ordered Multifractal Analysis (ROMA), a technique capable of deciphering the multifractal characteristics of intermittent fluctuations, was originally applied to the results of a magnetohydrodynamic (MHD) simulation. Application of ROMA to measured fluctuations in the auroral zone, due to the dominant physical effects changing from kinetic to MHD as the scale increases, requires an additional level of rank-ordering in order to divide the domain of scales into regimes. An algorithm for the additional step in this double rank-ordering technique is discussed, and is demonstrated in the application to the electric field fluctuations in the auroral zone as an example. As a result of the double rank-ordering, ROMA is able to take into account the nonlinear crossover behavior characterized by the multiple regimes of time scales by providing a scaling variable and a scaling function that are global to all the time scales.


Introduction
Electromagnetic field fluctuations in space plasmas are often observed to be intermittent.Such has been the case, for instance, for the measured fields in the solar wind, magnetospheric cusp, plasma sheet and the auroral zone (Burlaga, 1991;Bruno et al., 2001;Weygand et al., 2005;Tam et al., 2005;Echim et al., 2007).The origin of intermittent fluctuations in magnetized plasmas was interpreted as the result Correspondence to: S. W. Y. Tam (sunwytam@pssc.ncku.edu.tw) of the sporadic mixing and/or nonlinear interactions of localized coherent structures (Chang et al., 2004; and references contained therein).An example of such structures, in the solar wind for instance, was suggested to be in the form of fieldaligned magnetic flux tubes of various sizes, convected by the flow (Bruno et al., 2001).In the case of the auroral zone, coherent structures probably include nearly two-dimensional oblique potential structures, such as those generated in simulations based on the reduced MHD formulation of the inertial Alfvén fluid equations (Seyler, 1990), along with other small-scale kinetic coherent structures.In general, coherent structures in space may also be in other forms such as double layers, ion holes, electron holes, current filaments, etc. Intermittent fluctuations of electromagnetic fields are generated when those structures interact stochastically.Such signatures of the interactions are Doppler shifted when detected in frames moving relative to the plasma.Spatial fluctuations and spatial scales of intermittent turbulence appear as temporal when recorded in time-series measurements by spacecraft.Hence, time series of intermittent fluctuations observed in space plasmas contain information about the interactions of the structures.The multifractal behavior associated with the intermittent nature of the fluctuations, in particular, is related to the different fractal dimensions of those interactions at various scales.Thus, the ability to analyze multifractal behavior is essential for the understanding of multiscale interactions in the turbulent medium.
S. W. Y. Tam and T. Chang: Double rank-ordering technique of ROMA first proposed by Chang and Wu (2008), who applied it to the spatial fluctuations extracted from the results of a twodimensional magnetohydrodynamic (MHD) simulation.The method was later applied to time series measured in the solar wind environment (Chang et al., 2008) as well as other solar-terrestrial environments (Chang et al., 2010).At the spatial and temporal scales considered in those studies, the fluctuations were primarily governed by the physics of MHD.Recently, an extension of the ROMA technique has been applied to electric field fluctuations observed in the auroral zone (Tam et al., 2010).Unlike the cases of the previous ROMA studies, the observed fluctuations in the auroral zone carried the signatures of both MHD and kinetic effects, with the former effect dominating at large scales, the latter effect at small scales, and a crossover behavior at scales in between.In other words, the different governing physics divides the scales into different regimes.The extension of the ROMA method, using a double rank-ordering technique, was designed to obtain global scaling properties over all the regimes of temporal or spatial scales.
In Sect.2, we review the ROMA technique applicable to an individual regime of temporal or spatial scales.Our discussion will focus on the similarities and differences between this technique and the traditional structure function analysis along with the one-parameter scaling of monofractals, and explain why ROMA has certain advantages over the traditional method when analyzing multifractal fluctuations.In Sect.3, we shall discuss in detail how ROMA can be systematically extended to analyze fluctuations featuring crossover behavior in the domain of temporal or spatial scales.We shall demonstrate the extension of the ROMA technique in Sect. 4 by applying it to measured electric field fluctuations in the auroral zone.

ROMA and structure function analysis: similarities and differences
The traditional structure function analysis combined with the idea of single-parameter scaling is useful for describing monofractal behavior.However, as we shall discuss below, the analysis becomes inadequate in the case of multifractals.
The ROMA technique remedies such inadequacy.ROMA and the structure function analysis share the same initial steps in their applications to a time (or spatial) series.Given a time series of a quantity X, both methods begin with the calculation of the increment of the quantity over a prescribed time scale τ : δX ≡ X(t +τ )−X(t).One can then obtain the PDFs of the absolute value of δX at scale τ , P (|δX|,τ ), with the normalization condition: For the traditional analysis, one calculates the structure functions of time scale τ and moment order q > 0 as: where ... denotes averaging over t.One then finds the fractal dimension ζ q of the moment order q by looking for the scaling behavior S q (τ ) ∼ τ ζ q . (3) For the special case where ζ q is proportional to q, i.e.
all the fractal properties of the fluctuations can be characterized by a single number ζ 1 , meaning that the fluctuations are monofractal.One may then scale the PDFs for different τ 's with one scaling function P s and one power-law scaling parameter s = constant as follows (Chang et al., 1973;Hnat et al., 2002): where τ 0 is a reference time scale.Equation ( 5) implies a scale-invariant functional relation between P (|δX|,τ ) τ τ 0 s and By substituting Eqs. ( 5) and (6) into Eq.( 2), one can show that which, from Eq. (3), gives With the monofractal condition, Eq. ( 4), s takes on the constant value ζ 1 .For fluctuations that exhibit multifractal behavior, Eq. ( 4) is no longer true and s cannot be a constant, thus invalidating the scaling relation Eq. ( 5).
For multifractal fluctuations, however, there may still be portions of the PDFs for different τ 's where Eq. ( 5) is valid with a constant scaling exponent s, in which case such portions of the PDFs would share the same fractal behavior.The basic idea of ROMA is to rank order the domain of the PDFs such that each rank can be scaled according to Eq. ( 5).In the original application of the ROMA technique (Chang and Wu, 2008), the rank-ordering parameter was therefore chosen as A scaling exponent s was then determined for each rank.The function s(Y ) can thus be considered as a spectrum of the scaling exponent varying over the different ranks.As a result of this rank-ordering, the PDFs for different τ 's can be scaled with the scaling function P s (Y ) and the power-law scaling exponent s(Y ) as follows: The methods of solving for s(Y ) and P s (Y ) were detailed in the literature (Chang and Wu, 2008;Tam et al., 2010;Wu and Chang, 2011) and will not be repeated here.
3 ROMA with double rank-ordering

Motivation
For the application of ROMA to the electric field fluctuations observed by the SIERRA sounding rocket in the auroral zone (Fig. 1 of Tam et al., 2010), however, there are reasons to believe that the rank-ordering scheme based on Eq. ( 9) alone is not adequate enough to cover all the time scales τ .As discussed in Section 1, the time series of the electric field fluctuations measured by the rocket could be interpreted as spatial fluctuations, given that the rocket was moving relative to the turbulent plasma.With the assumptions that a significant fraction of the fluctuations were electrostatic and transverse (Chang, 2001), and that the horizontal speed U of the rocket was much larger than that of the movements of the turbulent fluctuations and the geomagnetic field was essentially vertical, the time scales τ in the study and in the discussion below may be viewed approximately as spatial scales ≈ U τ , where U ≈1.5 km s −1 .With structures ranging in size from kinetic to MHD scales in the auroral zone, the measured fluctuations included a mixture of the effects due to the different governing physics.From the physical point of view, there is no particular reason why fluctuations due to the MHD and kinetic effects, and a mixture of the two should be characterizable by the same scaling exponent.To support the above argument, let us approach the problem from another perspective by examining Fig. 1, plots of S q (τ ) vs. τ in logarithmic scales (referred to as the log-log plots hereafter) generated from the electric field fluctuations in the auroral zone for a number of moment order q at time scales ranging from 5 to 1280 ms.At the smaller time scales, the plot for any given q appears as a straight line, implying a power-law relationship between S q (τ ) and τ .When τ increases to about 160 ms, some of the results start to deviate from the straight lines, particularly those for moment order q ≤ 2, indicating that the above power-law relationship starts to break down.As τ further increases to 320 ms and beyond, it is clear from all the plots that the relation between S q (τ ) and τ is totally different from the power-law relationship exhibited at the small scales.Based on Eq. ( 3), we can interpret the fractal dimension ζ q as the slope of the plot logS q (τ ) vs. logτ .The local slope at τ thus corresponds to the fractal dimension at that time scale.As the local slopes of the plots in Fig. 1 change with τ and become drastically different at the S q (τ) q = 0.5 q = 1 q = 1.5 q = 2 q = 2.5 q = 3 q = 3.5 q = 4 q = 4.5 q = 5 Fig. 1.Plots of S q (τ ) vs. τ generated from measured electric field fluctuations in the auroral zone at selected moment order q.
two ends of the time-scale range, we expect the fractal behavior of the fluctuations to change significantly across the range as well.
The changing fractal behavior with the time scale motivated us to introduce another level of rank-ordering for ROMA, which involves τ (Tam et al., 2010).The domain of time scale is to be divided into regimes, each of which to be further examined with the original rank-ordering scheme of ROMA (see Sect. 2) to find a spectrum of the scaling exponent s(Y ).Hence, τ is the first parameter to be rankordered in this double rank-ordering ROMA scheme; each rank at this level corresponds to a regime of time scales.The rank-ordering is based on how well the data points in the loglog plots can be fitted with straight lines.Time scales that are fitted well together belong to the same regime.In other words, the time scales within each regime is characterized by a power-law relationship between S q (τ ) and τ .

Rank-ordering algorithm on time scales
To discuss the rank-ordering procedure involving τ in more detail, let us use the notation We assume that each regime of time scales is continuous and covers at least two of the τ i 's.The regimes are also assumed to be contiguous, with the non-power law crossover ranges between adjacent regimes being very narrow such that the changes across regimes are essentially characterized by abruptly changing but piecewise continuous power-laws.The boundary between two adjacent regimes is taken to be at one of the τ i 's, which is a common time scale that belongs to both regimes.The rank-ordering procedure is based on applications of an algorithm that determines the upper limit of a regime, provided that the lower limit of the regime is given.The algorithm utilizes the idea that linear fitting of three of more data points in the log-log plots is more accurate for the case where the τ i 's at those points are all from the same regime than in the scenario where they belong to different regimes.This idea is illustrated in Fig. 3, whose two panels show the log-log plot for q = 2.5 along with a few fitted straight lines for subsets of the data points.One can see, in particular, that the data points from τ 1 to τ 5 , corresponding to time scales from 5 to 80 ms, line up close to a straight line, while the data points at τ 7 (τ = 320 ms) or larger clearly deviate from that line.Now if we try to fit a straight line through a subset of the first five data points, as we do for three of the subsets to result in the three straight lines in the top panel, we can see that the discrepancy between the fitted values and the actual values of logS q at the τ i 's that belong to the subset is small.In contrast, the corresponding discrepancy becomes larger, as shown in the bottom panel, when the fitting includes also the data point at τ 7 , a point that clearly does not line up with those at smaller τ i 's.The contrast between the two scenarios suggests that the discrepancy discussed above can be utilized to determine whether a specific τ i shares the same regime with other time scales.To explain how this idea is applied to our algorithm to determine the domain of a regime, let us use the example of Regime 1, which is defined to be the regime of the smallest time scales in our study.The lower limit of Regime 1 is set at τ 1 .The regime must also include τ 2 , but may or may not include more time scales.Hence, the upper limit of Regime 1 is at τ n , for some n ≥ 2. Because the regime is continuous by assumption, for any i that satisfies n ≥ i ≥ 1, the time scale τ i must belong to the regime.To find what the value of n is, we check consecutive τ i 's one at a time, starting with τ 3 in the case of Regime 1, examining the average accuracy of the fitted straight lines in the log-log plots involving that specific τ i , together with all the subsets of at least two smaller time scales that are in the regime.If the average accuracy is too low (i.e. the average discrepancy of the fitting, to be defined below, is larger than an acceptable level), then the upper limit of the regime would be τ i−1 .On the contrary, if the average discrepancy of the fitting is within an acceptable level, then Nonlin.Processes Geophys., 18,[405][406][407][408][409][410][411][412][413][414]2011 www.nonlin-processes-geophys.net/18/405/2011/ the τ i being examined would be added to the regime.We would then move on to use the same procedure to examine τ i+1 and so on, until we find the upper limit of the regime.
Because the upper limit of Regime 1 corresponds to the lower limit of Regime 2, and so on, the algorithm can be readily applied to the next adjacent regime, eventually distinguishing the domain of all the regimes in terms of the τ i 's.
To complete the description of the rank-ordering algorithm on the time scales, we shall discuss how the average discrepancy involving a particular τ i is defined for the purpose of checking whether the time scale belongs to a certain regime.First of all, in any given log-log plot of S q (τ ) vs. τ , there are nine data points corresponding to the nine τ i 's in our study.If one chooses a subset of three or more data points to fit a straight line, the fitted values of logS q may not be exactly the same as the actual values at the corresponding subset of τ i 's.Thus, there is a discrepancy associated with every such fitted line.For our purpose to examine whether a τ i belongs to a certain regime, we take into consideration all the possible combinations of data points that include τ i along with at least two of the smaller time scales that are in the regime, and calculate an average discrepancy over the fitted lines of all such combinations at all the moment order q.In order to keep track of the data point combinations for presentation purpose, we find it convenient to use the binary form of nature numbers.The binary digits, reading from the right, would represent the order of the τ i 's; a "1" in the digit would indicate that the corresponding τ i is included in the combination and a "0" means that it is not.The solid red line in the top panel of Fig. 3, for instance, is determined by fitting data points of the subset of time scales (τ 3 ,τ 4 ,τ 5 ).The line thus corresponds to the binary number 11100, which is the natural number N = 28.Having quantified data point combinations with the label N, we now proceed with the mathematical formulations for finding the average discrepancy.Let the equation of the line fitted for a subset of data points, represented by N, at moment order q be y q = a N,q x +b N,q , where x = logτ and y q = logS q . (11) Suppose the number of data points in the subset is n N .Provided n N ≥ 3, we take the discrepancy associated with the linear fit line for a data point to be: where the summation is over all of the n N time scales in the subset, and y j,q = logS q (τ j ) is the actual value of y q at the time scale τ j .We then obtain the root-mean-square discrepancy over all the moment order for the subset of data points: S S N vs. log 2 N for every natural number N < 2 9 that represents a subset of three or more τ i 's.(Bottom) Same as above but for N < 2 6 only.
We would like to relate y to the discrepancy from the structure function due to the fitting.The differential form associated with Eq. ( 11) suggests the following approximation for the average fractional difference between the values of the structure functions and the results of the fitting based on the data point subset labeled N: S S N = (ln( 10))( y) N . ( The quantity in Eq. ( 14) constitutes one of the two requirements for adding a specific time scale τ i to a specific regime: S S N ≤ 0.1 for all the data point subsets N that are relevant to that specific step of examining τ i in the algorithm.With this criterion, the structure functions for each subset of τ i 's within a regime deviate from the respective fitted lines in the log-log plots by less than 10 % on average.If multiple subsets of data points are relevant, then there is an addition requirement of S S τ i ≤ 0.05, where www.nonlin-processes-geophys.net/18/405/2011/ Nonlin.Processes Geophys., 18, 405-414, 2011

S S τ i
≡ is the weighted average of the relevant S S N with the weighing factor being the number of points used in the linear fitting of the log-log plots.With the additional requirement, our study essentially aims for an average uncertainty of 5 % or less in the fitting of the log-log plots when we examine the possibilities of four or more τ i 's in the same regime.Such a stricter requirement in uncertainty is imposed because Eq. ( 15) involves multiple subsets of data points, which makes the concept of average uncertainty statistically more meaningful.

Results of ROMA application to the auroral zone electric field
We have applied the rank-ordering algorithm described in Sect. 3 to the electric field fluctuations in the auroral zone, as shown in Fig. 1 of Tam et al. (2010).Below, we discuss how to determine the domain of Regime 1 as an example.The top panel of Fig. 4 shows the results of S S N for all subsets consisting of three or more data points, i.e. n N ≥ 3. The same results, but limited to N ≤ 64, is shown in the bottom panel.Note that a given τ m is the largest time scale in a subset of data points if and only if the corresponding label N is in the range [2 m−1 ,2 m − 1].Thus, to determine whether τ m belongs to Regime 1, the results in Fig. 4 with m − 1 ≤ log 2 N < m are all that is relevant.We notice from the top panel of the figure that as N increases, S S N in general becomes significantly larger in each consecutive range of N = [2 m−1 ,2 m − 1] for m ≥ 7 (log 2 N ≥ 6), meaning that the approximation of a linear relationship in the loglog plots becomes increasingly inaccurate when each of the time scales of τ 7 or larger is added into consideration.In particular, it is clear that τ 7 does not belong to Regime 1, because when it is added as the largest time scale for the fitting (6 ≤ log 2 N < 7), the corresponding S S N may reach as high as 0.2.On the other hand, all the values of S S N at log 2 N < 6 appear to be smaller than 0.1 in the plot.This is confirmed by the plot in the lower panel of Fig. 4, which shows S S N only at such a range of log 2 N. When it comes to determining the upper limit of Regime 1, as discussed earlier, the first time scale to examine is τ 3 .But for τ 3 to be the largest time scale for the fitting of three or more data points, there is only one possibility, namely the subset (τ 1 ,τ 2 ,τ 3 ) corresponding to N = 7.Thus, for τ 3 to be included in Regime 1, there is only one criterion: S S N=7 ≤ 0.1, which is satisfied as shown in Fig. 4. The subsequent steps in the algorithm would be to examine τ 4 and, if necessary, τ 5 and τ 6 .There are multiple subsets of data points that would be relevant to the procedure of checking these time scales.In fact, within a given regime of time scales, there are 2 m−1 − m subsets of three or more data points where the m-th smallest τ i is the largest time scale of the subset.We have shown in Fig. 4 that the criterion based on S S N does not exclude τ 4 , τ 5 and τ 6 from Regime 1.In Fig. 5, we show the results of S S τ i for these three time scales (40, 80 and 160 ms), obtained by applications of the averaging scheme in Eq. ( 15) to the 2 m−1 − m relevant subsets in the range of N = [2 m−1 ,2 m −1] with m = 4, 5 and 6, respectively.As S S τ i is smaller than 0.05 for only τ 4 and τ 5 , but not τ 6 , the upper limit of Regime 1 would be at τ 5 , which corresponds to a time scale of 80 ms.For Regime 2, it would then consist of at least τ 5 and τ 6 .To check whether τ 7 belongs to Regime 2, there is only one relevant subset of data points, namely (τ 5 ,τ 6 ,τ 7 ), corresponding to N = 112.We find that S S N=112 = 0.125, which is too large for τ 7 , which corresponds to a time scale of 320 ms, to be in the regime.Regime 2, therefore, ranges from 80 to 160 ms.The results for Regime 2 provide the lower limit of Regime 3 at τ 6 = 160 ms.Applying the same rank-ordering procedure to Regime 3, we find that S S N=224 = 0.274, where N = 224 corresponds to the subset (τ 6 ,τ 7 ,τ 8 ).The range of Regime 3 is thus determined to be from τ 6 to τ 7 , that is, from 160 to 320 ms.Thus, τ 7 is the lower limit of Regime 4, which also includes τ 8 .At this point, only one τ i in our study remains unaccounted for, namely τ 9 = 1280 ms.To determine whether τ 9 would belong to Regime 4, we only need to calculate S S N=448 for the subset (τ 7 ,τ 8 ,τ 9 ).It turns out that S S N=448 = 0.089, small enough for τ 9 to be incorporated into Regime 4. Hence, the auroral zone electric field fluctuations are rank-ordered into four different regimes in terms of the time scales from 5 to 1280 ms: Regime 1 from τ ≈ 5 to 80 ms; Regime 2 from τ ≈ 80 to 160 ms; Regime 3 from τ ≈ 160 to 320 ms; and Regime 4 for τ 320 ms.
For each of the regimes, we have applied the original ROMA technique to rank order the domain of |δE|, the absolute value of the increment of the electric field that plays the role of |δX| in the discussion of Sect.2, to study the multifractal characteristics of each regime separately (Tam et al., 2010).The reference time scale τ 0 is taken to be the smallest τ i of the regime.Figure 6 shows the results of the scaling exponent s i (Y i ), where the subscript i is added to "s" and "Y " (see Eq. 9, for instance) to denote the regime number.The variation of s i as a function of Y i can reveal certain information about the fluctuations at the time scales of the regime.For example, Tam et al. (2010) has pointed out the similarity between s i and the Hurst exponent; thus, the persistency (anti-persistency) of the fluctuations is indicated by s i of values larger (smaller) than 0.5, following the classical demarcation for the Hurst exponent.In addition, how fast s i changes with Y i suggests how developed the turbulence is at the time scales of the regime.Based on these general indications by s i , we can see from the top left panel of Fig. 6 that the fluctuations are persistent at the time scales of Regime 1.The persistency may be due to kinetic effects, which are probably important in this regime of small scales.At small values of Y 1 , which correspond to small sizes of |δE|, the scaling exponent s 1 increases rapidly, an indication of possible developing instabilities and turbulence.The fluctuations seem to settle down to a more stable and developed turbulent state at larger Y 1 , as the values of s 1 seem to become more and more slowly varying.Such interpretations of the turbulent state based on the behavior of s 1 (Y 1 ), in particular, are consistent with the effects of rapidly growing linear or nonlinear instabilities when the amplitudes of the turbulent fluctuations are small (corresponding to small Y 1 ) but has smaller influence on the nearly stationary statistical processes (including such effects as nonlinear saturation) for large turbulent fluctuations (corresponding to large Y 1 ).For Regime 2, s 2 exhibits fluctuations in the range around 0.5 at small values of Y 2 , as shown in the top right panel of Fig. 6.Such fluctuating behavior is rather similar to that of s 1 at small values of Y 1 , except that the values for s 2 are considerably lower.Thus, the developing turbulence seems to be of a mixture of persistent and anti-persistent nature, probably as a result of effects beyond the kinetic range starting to play a non-negligible role at the scales of this regime.As Y 2 becomes larger, the values of s 2 become more stable, indicative of the turbulence settling down to a more stable and developed state, similar to the case for Regime 1.The apparent persistent nature of the fluctuations suggested by the values of s 2 at large Y 2 is perhaps due to kinetic effects still being more dominant than those of larger scales.Regime 3 features large fluctuations in s 3 at small Y 3 (lower left panel of Fig. 6).The significant decrease in the values of s 3 from 0.677 to 0.285, which covers the range of Y 3 from 5 to 19, indicates that the turbulence is highly unstable at the time scales of this regime at those values of Y 3 .At larger Y 3 , s 3 seems to settle at around 0.5 when the turbulence becomes more stable, suggesting a mixture of persistent and anti-persistent fluctuations similar to the case of Regime 2. However, the lower settling value of the scaling exponent in Regime 3 compared with Regime 2 is probably due to the kinetic effects becoming even less dominant with the increase in time scale.For Regime 4, as shown in the bottom right panel of Fig. 6, there is an increase in s 4 at small values of Y 4 .After reaching a peak value close to 0.5, s 4 then decreases monotonically as Y 4 increases further.The monotonically decreasing trend of the rank-ordered spectrum and the anti-persistent nature of the fluctuations are qualitatively similar to those obtained for the solar wind (Chang et al., 2008) and MHD numerical simulations (Chang and Wu, 2008).Thus, it is tempting to conjecture that in Regime 4, where the time scales are larger compared with the other regimes, the fluctuations bear the signatures of MHD turbulence.With the solution s i (Y i ), we have mapped the PDFs of the time scales in each regime according to Eq. ( 10).As shown in Fig. 7, the PDFs collapse quite well into the corresponding scaling function P si (Y i ) of the regime, particularly at values of Y i that are not too large so that the samples are sufficient for the statistics to be meaningful.We should note that the scaling exponent of one regime is generally not applicable to the time scales of other regimes.For instance, if we map the PDF of a time scale in Regime 2 using the scaling exponent s 1 (Y 1 ), as shown in Fig. 8 for τ = 160 ms, the mapped PDF does not collapse into the scaling function P s1 (Y 1 ) of Regime 1.This discrepancy between the mapped PDF for τ = 160 ms and those for the time scales in Regime 1 justifies the necessity of our approach to rank order the time scales by regime.We have shown that application of ROMA to the four regimes allows us to find the power-law scaling exponent s i (Y i ) and the scaling function P si (Y i ) for each regime i.The scaling relation for the PDFs in the i-th regime (i goes from 1 to 4 in our case) is: where τi is the smallest time scale of the i-th regime and Y i is the "scaled" parametric scaling variable implicitly provided by the equation: The four scaling variables as well as the scaled PDFs for i = 1,2,3,4 are related due to the assumed piecewise continuous property across the contiguous regimes.To derive the relationship among the scaling variables, let us consider fluctuations of an arbitrary magnitude |δE| = X 0 at an arbitrary time scale τ = τ * in the (i+1)-th regime.The fluctuations correspond to a value of Y i+1 , which, according to Eq. ( 17), is implicitly given by: But based on Eq. ( 17), the fluctuations would share the same value of Y i+1 with fluctuations of magnitude |δE| = X 0 τ * τi+1 −s i+1 (Y i+1 ) at the time scale τi+1 , which also belongs to the i-th regime and thus corresponds to: Equations ( 18) and ( 19) imply a recursion relation for the scaling variables: To derive the relationship among the scaled PDFs, let us consider the time scale τ = τi+1 , which belongs to both the i-th and the (i+1)-th regimes.The PDF at this time scale can be scaled by the power law of either regime according to Eq. ( 16): Equations ( 21) and ( 22) together lead to the following recursion relation for the scaled PDFs:

Summary
We have discussed how the original ROMA technique was able to decipher the multifractal properties of temporal and spatial fluctuations by building on the ideas of traditional structure function analysis and one-parameter scaling for monofractals.Extending the idea of ROMA a step further by introducing an additional level of rank-ordering, we have shown that the technique is applicable to fluctuations that feature crossover behavior due to different governing physics over the domain of time scales.Using the electric field fluctuations in the auroral zone as an example, we have discussed an algorithm that enables us to separate the time scales into regimes, each of which shown to follow its own power-law scaling.Those scaling behaviors are essential for the description of the crossover behavior over all the regimes.And a global scaling variable and a global scaling function that characterize the transition over the regimes can be obtained based on the scaling exponents of the individual regimes.The value of the global scaling variable Y global changes within and across the time scales of different regimes, and would serve as useful guidelines for theoretical studies of physical processes or effects that span multiple regimes.

Fig. 3 .
Fig.3.Linear fitting of the log-log plot of S q (τ ) vs. τ at q = 2.5 based on subsets of data points that correspond to the following time scales: (Top) solid red line: (20, 40 and 80 ms); dashed blue line: (10, 20 and 80 ms); dot-dashed magenta line: (5, 10 and 80 ms).(Bottom) The data point at 320 ms is added to the subsets of the top panel, represented by the corresponding line styles.
Fig. 4. (Top)S S N vs. log 2 N for every natural number N < 2 9 that represents a subset of three or more τ i 's.(Bottom) Same as above but for N < 2 6 only.

Fig. 6 .
Fig. 6.Rank-ordered spectra for the scaling exponents s i (Y i ) of the four regimes, with i =1, 2, 3 and 4 corresponding to the regime number.The extent of the horizontal lines indicates the ranges in Y i over which the scaling exponent s i is obtained.

Fig. 7 .
Fig. 7. Scaling function P si (Y i ) obtained from the PDF at each time scale of Regime i, where i =1, 2, 3 and 4.

Fig. 8 .
Fig. 8. of the PDF at the time scale τ = 160 ms using the scaling exponent of Regime 1.Comparison with the same mapping for the PDFs at the time scales of Regime 1 shows that the PDF at 160 ms does not collapse to the scaling function of Regime 1.

Fig. 9 .
Fig. 9. (Top) The rank-ordered spectra s 1 (solid black), s 2 (dotted red), s 3 (dot-dashed blue), and s 4 (dashed green) as functions of Y global .(Bottom) Global scaling function P s1 (Y global ) obtained from the PDFs at all the time scales of the four regimes.
P s(i+1) (|δE|) = τi+1 τi −s i (Y i ) P si |δE| τi+1 τi −s i (Y i ) .(23) Applying Eq. (20) to the four time regimes, we find a global scaling variable Y global across the four regimes: with Y global , P s1 is the scaling function that is global to all four regimes: allows us to express s i for i = 1,2,3,4 as a function of Y global , as shown in the top panel of Fig. 9. From the figure, we can clearly see that except for highly unstable turbulence, the fluctuations exhibit a generally decreasing trend of s i at given Y global as i increases from 1 to 4. The fluctuations become increasingly anti-persistent as the regimes cross over from kinetic to MHD.Based on the fractal exponents s i , P (|δE|,τ ) of all the time scales of the four regimes can be mapped to collapse into one profile, namely the global scaling function P s1 Y global as shown in the bottom panel of Fig. 9.