Identiﬁcation of Fractal Properties in Geomagnetic Data of Southeast Asian Region during Various Solar Activity Levels

: The fractal properties of geomagnetic northward component data (H-component) in the equatorial region during various phases of solar activity over Southeast Asia were investigated and then quantiﬁed using the parameter of the Hurst exponent ( H ). This study began with the identiﬁcation of existence of spectral peaks and scaling properties in international quiet day H-component data which were measured during three levels of solar activity: low, intermediate, and high. Then, various cases of quiet and disturbed days during different solar activity levels were analyzed using the method that performed the best in the preceding part. In all the years analyzed, multifractal scaling and spectral peaks exist, signifying that the data have fractal properties and that there are external factors driving the ﬂuctuations of geomagnetic activity other than solar activity. The analysis of various cases of quiet and disturbed days generally showed that quiet days had anti-persistence tendencies ( H < 0.5) while disturbed days had persistence tendencies ( H > 0.5)—generally a higher level of Hurst exponent compared to quiet days. As for long-term quiet day H-component data, it had a Hurst exponent value that was near H (cid:39) 0.50, while the long-term disturbed day H-component data showed higher values than that of the quiet day. Author Contributions: Conceptualization, F.N.R. N.S.A.H.; Methodology, F.N.R. N.S.A.H.; Validation, N.S.A.H. and A.B.R.; Investigation, F.N.R.; Resources, A.Y.; Data curation, A.Y.; Writing— draft preparation, F.N.R.; Writing—review editing, N.S.A.H. and A.B.R.; Visualization, F.N.R.; Supervision, N.S.A.H.; Project administration, N.S.A.H.; Funding acquisition, N.S.A.H.


Introduction
The geomagnetic field of the Earth is an astonishing feat of nature as it protects us from the incoming forces from the Sun, such as charged particles in the solar wind [1,2] that may be the harbinger of catastrophe on earth, if the geomagnetic field were to never exist in the first place. This so-called catastrophe is the geomagnetic storm; its occurrence of which can be attributed to the compression of the magnetosphere due to the solar wind [3]. When such a storm occurs, it can disrupt the telecommunication systems and electrical grids on earth, among many other things. There are two classes of geomagnetic activity, namely 'quiet day' and 'disturbed day'. A quiet day is the period where the fluctuations of geomagnetic activity are only affected by solar and lunar daily variations, while a disturbed day is the period where the fluctuations are affected predominantly by the occurrence of a geomagnetic storm.
In understanding the characteristics of geomagnetic activity, one of the attributes analyzed is its fractal properties [4][5][6], which are the characteristics of self-similarity that a particular object has when subjected to varying degrees of magnification [7,8]. A fractal parameter, the Hurst exponent, is used to quantify these properties by defining the longterm correlation of the time series of a geomagnetic activity [9]. The value of Hurst exponent

Data Procurement
The geomagnetic northward component, H, or H-component, was chosen as the basis for this study as it is very susceptible to the geomagnetic volatility of the Earth at the geomagnetic equator-the scope of the area investigated in this study-and virtually considered to show the same characteristics as the total field component [14,24]. The data used were of 1-s resolution for the year 2009, 2013, and 2015 which covered three different solar activity levels of low, intermediate, and high solar activity, respectively.
The H-component data used in this study was taken from the geomagnetic stations of Davao, Philippines (7.00 • N, 125.40 • E) and Langkawi, Malaysia, (6.30 • N, 99.78 • E) -henceforth abbreviated as DAV and LKW, respectively-all of which are a part of the Magnetic Data Acquisition Systems (MAGDAS) authorized by the International Center for Space Weather Science and Education (ICSWSE) of Kyushu University, Japan [25]. Figure 1 illustrates a sample of the H-component data analyzed in this study.
of Davao, Philippines (7.00° N, 125.40° E) and Langkawi, Malaysia, (6.30° N, 99.78° E)henceforth abbreviated as DAV and LKW, respectively-all of which are a part of the Magnetic Data Acquisition Systems (MAGDAS) authorized by the International Center for Space Weather Science and Education (ICSWSE) of Kyushu University, Japan [25]. Figure  1 illustrates a sample of the H-component data analyzed in this study. Dst index which is an aggregate index of H-component, A-index which is an index of Kp index average over the period of 24 h, and the International Quiet and Disturbed Day Database were utilized to identify quiet and disturbed days of each year. Datasets of the Dst index and International Quiet and Disturbed days were obtained from WDC Kyoto web services, while datasets of A-index were obtained from NASA's OMNIWeb service.
This study was divided into three parts. The first part involved the identification of the fractal properties of International Quiet Day H-component data which required the selection of H-component data of 60 quietest days according to the International Quiet Day data for the year 2009, 2013, and 2015 from DAV station. The selected data were later analyzed using the PSA method. It was found that the data from LKW station were scattered with noises if it were to be analyzed in the long-term fashion. Therefore, it was decided that the data from LKW station would be omitted for the first part of the study and to be used exclusively for the short-term event analyses.
The second part involved the testing of various fractal methods to find the best method in determining the Hurst exponent in the geomagnetic data. This required all the aforementioned methods-namely the PSA, RRA, DFA, and r-DFA-to be tested multiple times against various synthetic signals that simulated the characteristics of the geomagnetic data. Each method was tested against five different types of signals, measured by Hurst exponent. For each type of signal, three distinctive datasets were synthesized to be Dst index which is an aggregate index of H-component, A-index which is an index of Kp index average over the period of 24 h, and the International Quiet and Disturbed Day Database were utilized to identify quiet and disturbed days of each year. Datasets of the Dst index and International Quiet and Disturbed days were obtained from WDC Kyoto web services, while datasets of A-index were obtained from NASA's OMNIWeb service.
This study was divided into three parts. The first part involved the identification of the fractal properties of International Quiet Day H-component data which required the selection of H-component data of 60 quietest days according to the International Quiet Day data for the year 2009, 2013, and 2015 from DAV station. The selected data were later analyzed using the PSA method. It was found that the data from LKW station were scattered with noises if it were to be analyzed in the long-term fashion. Therefore, it was decided that the data from LKW station would be omitted for the first part of the study and to be used exclusively for the short-term event analyses.
The second part involved the testing of various fractal methods to find the best method in determining the Hurst exponent in the geomagnetic data. This required all the aforementioned methods-namely the PSA, RRA, DFA, and r-DFA-to be tested multiple times against various synthetic signals that simulated the characteristics of the geomagnetic data. Each method was tested against five different types of signals, measured by Hurst exponent. For each type of signal, three distinctive datasets were synthesized to be tested against each method. The results considered for comparison were the averaged results.
The third and final part involved the characterization of both short-and long-term cases involving quiet and disturbed days. For the short-term cases studied, all the major solar event days with Dst < −200nTin the year 2013 and 2015 were required to be concatenated into a single dataset for each year and then analyzed by using data from DAV and LKW stations. Next, two major solar events-namely, the moderate solar storm and intense solar storm which happened in 2013 and 2015 respectively-were analyzed using data from DAV and LKW stations. The results of the preceding analyses were then compared with a couple of quiet days in their respective period, following the period of study of the second analysis which happened to be June solstice for 2015 and March equinox for 2013. All of the dates for the short-term cases analyzed in this study are listed in Table 1. Next, for the long-term cases studied, the 60 quietest days according to the International Quiet Day data for the year 2009 and the 60 most disturbed days according to the International Disturbed Day data for the year 2013 and 2015 were analyzed using data from DAV station. The datasets were then refined by eliminating days with A-index below 25 (A-index < 25) for each year to be analyzed.

Methods
Four methods were utilized for identifying fractal properties in this study-namely, the Power Spectrum Analysis (PSA), Rescaled Range Analysis (RRA), Detrended Fluctuation Analysis (DFA), and Robust Detrended Fluctuation Analysis (r-DFA). Alongside all of the abovementioned methods, the concept of fractional Brownian motion was also heavily utilized throughout this study, which is explained extensively in this section.

Power Spectrum Analysis (PSA)
The first method utilized was the power spectrum analysis. It is a method that evaluates a time series in a way that enables us to point out random inflections within it [26], thus allowing for a thorough examination of the scaling features [14]. The analysis starts with the Fourier transformation of the geomagnetic time series where y n is the geomagnetic time series and m = 0, 1, 2, N − 1. The time series may indeed be a self-affine fractal if the correlation of the power spectrum, P(f ) and the frequency, (f ) shows a linear relation where the power spectrum decreases alongside the frequency, P( f ) ∼ f −β where β indicates the degree of correlation [26] or in our case, the spectral exponent [14]. From the spectral exponent, β, the Hurst exponent can be determined. The Hurst exponent for a stationary fractional Gaussian noise can be calculated with H = (β + 1)/2, for −1 ≤ β < 1. Meanwhile, for a non-stationary fractional Brownian motion, it can be calculated with H = (β − 1)/2 for 1 < β ≤ 3. On another note, if β happens to be equal to 1, the noise shall be in the form of 1/f [26,27]. The second method used was the rescaled range analysis (also known as R/S analysis). It is a method devised by Mandlebrot and Wallis [28], based on the pioneering hydrological study of Hurst [29], which measures the parameter of Hurst exponent, H. This parameter allows the measurement of the intensity of long-range dependence in a time series [10].
The analysis starts with the division of the geomagnetic time series into n subintervals, y n , with size τ. The range of each subinterval is then calculated Next, the standard deviation of each subinterval is calculated The Hurst exponent is finally obtained using the correlation between R/S and τ; R/S n ∼ τ H where H is the slope of the double logarithmic plot of <R/S> n against τ.
As of recently, the method of RRA is still being utilized in analyzing medical data [30], geological data [31][32][33], internet traffic data [34,35], and solar activity data [36]. Additionally, it is also utilized in the field of engineering [37].

Detrended Fluctuation Analysis (DFA)
The third method utilized was the detrended fluctuation analysis. It is a method devised by Peng et al. [38] to detect scaling features that a particular time series has. It is exceptionally created such that the spurious detection of scaling and long-range correlation can be avoided; thus, a generally more accurate result can be produced.
The detrended fluctuation analysis starts with the integration of the geomagnetic time series, y(t), of size N. Next, the time series is divided into n subintervals, y n (t), of size τ. The subsequent time series is then put through a detrending process in which the local trend is subtracted in each subinterval, which resulted in the fluctuation of the time series, F(τ). This calculation is repeated for each of the subintervals; the average fluctuation can be represented with the correlation of F(τ) n ∼ τ α where α is the scaling exponent, representing the slope of double logarithm plot of F(τ) against τ. From the scaling exponent, α, the Hurst exponent can be determined. The Hurst exponent for a stationary fractional Gaussian noise is H = α. Meanwhile, for a non-stationary fractional Brownian motion, the Hurst exponent is H = α − 1.

Robust Detrended Fluctuation Analysis (r-DFA)
The last method used was the robust detrended fluctuation analysis. It is a method devised by Habib et al. [47] as an improvement to the already established detrended fluctuation analysis (DFA) method. The improvement involves additional steps of analysis to the established method of DFA.
After a particular geomagnetic time series is computed using the regular DFA method, surrogate datasets are made on which the regular DFA will be performed upon. The results of the DFA analysis on the surrogate datasets would be used to modify the original DFA result, thus diminishing the bias inherent when using DFA of higher order [15]. Next, piecewise linear regression is used to optimize the location of the data where there is a change in scaling properties, which in this case can be referred to specifically as Universe 2021, 7, 248 6 of 14 'crossover' [48,49]-the existence of which is predefined-and also to minimize the least square error between data that is bound to happen.
If crossovers are indeed identified to exist in the data, a further test is performed in the form of ANCOVA (analysis of covariance) and multiple comparison procedure which would detect the exact number of the crossovers that the data has until the point of saturation and produces the local scaling exponent, α L for each section according to the division made by the crossovers. Instead, if crossovers are not identified, a robust regression is used to determine the global scaling exponent, α. As for the determination of the Hurst exponent, the method remains the same as the regular DFA method and can be implemented seamlessly with this method.

Fractional Brownian Motion (fBm)
Fractional Brownian Motion is devised by Mandlebrot and Van Ness [50] as a form of stochastic process which is continuous in nature and consists of Gaussian increments. It is defined by where B H (t) is fractional Brownian motion at time t, Γ is gamma function, H is Hurst exponent and B(s) is standard Brownian motion at time s [51]. At H = 0.5, it generalizes into a standard Brownian motion. The concept of fractional Brownian motion (fBm) was used extensively in the second part of this study; it involved the testing of the various methods using various synthetic signals. The various synthetic signals were in fact made using fBm as fBm can provide simulations that are quite close to real geomagnetic data.   One common attribute that can be observed from the periodograms is the scaling properties. The change in data scaling rule indicates that the data has fractal properties [8] which were present in all of the periodograms analyzed in this study. The markings at the bottom of Figure 1-namely M1 and M2-indicate that all of the three periodograms have virtually the same scaling properties, with two distinct areas having different scaling rule.  One common attribute that can be observed from the periodograms is the scaling properties. The change in data scaling rule indicates that the data has fractal properties [8] which were present in all of the periodograms analyzed in this study. The markings at the bottom of Figure 1-namely M1 and M2-indicate that all of the three periodograms have virtually the same scaling properties, with two distinct areas having different scaling rule.

Fractal Properties of Quiet Day Geomagnetic Data
While the identification of the point at which the scaling properties change is purely a qualitative process [52,53] and the use of power spectrum analysis as a method is not optimized to indicate precisely the scaling properties of a data [54], the fact remains that all the periodograms showed a change in scaling rule, a multifractal scaling [55,56]. Hence, all of the long-term quiet day data analyzed in this study exhibited fractal properties.
Aside from the scaling properties of the data, another attribute that can be observed from the periodograms is the spectral peaks. From Table 2  While the identification of the point at which the scaling properties change is purely a qualitative process [52,53] and the use of power spectrum analysis as a method is not optimized to indicate precisely the scaling properties of a data [54], the fact remains that all the periodograms showed a change in scaling rule, a multifractal scaling [55,56]. Hence, all of the long-term quiet day data analyzed in this study exhibited fractal properties.
Aside from the scaling properties of the data, another attribute that can be observed from the periodograms is the spectral peaks. From Table 2, this study found that the year 2009, 2013, and 2015 showed the same number of peaks at the same hour marks; 6 h, 8 h, 12 h, and 24 h. As for the study of Rabiu et al. [23], all the years studied-1996, 2000, and 2002-also showed the same number of peaks at the same hour marks (albeit different from our findings of peaks) of 8 h, 12 h, and 24 h, with an exception for the peak at the 6 h mark which was only present in the KOU station H-component data for the year 2000. The existence of peaks itself is often attributed to several factors such as solar daily variation [14,23,57], lunar daily variation [58], and atmospheric tides [57,59]. Solar daily variation is caused by the rotation of the Earth, causing the Earth's changes from night to day, which primarily causes the existence of the spectral peak at the 24 h mark. Lunar daily variation is caused by the motion of the moon revolving around the Earth, which causes the spectral peaks at the 12 h mark. The other factor that drives the existence of other peaks is the atmospheric tides caused by the rotation of the Earth, creating a global scale of oscillations which naturally have harmonic periods [59]. These harmonic periods correspond directly to the period in which the spectral peaks exist, which in turn, causes the remaining spectral peaks of 6 h and 8 h to exist.
It is concluded that all the long-term geomagnetic data over Southeast Asia during low, intermediate, and high solar activity showed fractal properties, and the existence of multiple spectral peaks showed that solar activity was not the sole driver of geomagnetic fluctuations. The findings in this study generally agreed with Rabiu et al. [23] study which used data from a different region. We came to this decision by considering the stations from which the data were used and the period of the study. 2000 and 1996 were set during the solar maximum while the years 2002 and 2015 were set during the solar minimum. It is worth noting that this study did not subtract the baseline of H-component as implemented by Rabiu et al. [23]. In addition, their study analyzed the 1 h resolution data of H-component for the whole year while this study analyzed 1 s resolution data of H-component for selected quietest days of the year according to the international quiet day data. Nevertheless, the spectral peaks present in both H-component data are virtually the same. Table 3 shows the result of the accuracy test for each method using various synthetic fBm signals. The Hurst exponent values were obtained by using linear regression between parameters specified to be utilized in each method (Section 2.2), while the errors were based on the standard error of the regression which indicates how far the observed values are from the regression line. The color gradient of orange-yellow-green-dark green indicates the degree of accuracy of each method, with orange being the least accurate and dark green being the most accurate. The PSA method came out to be the least accurate out of all the methods, while the DFA came out to be the most accurate. The PSA method is perfectly capable of identifying in detail the inflections and the scaling properties within a time series as demonstrated by the preceding part; a task that the other methods tested here were not that capable of. However, in measuring the overall Hurst exponent of a time series, it turns out that the PSA method was not that adequate at it. While it offered a great precise measurement by offering a small amount of standard error, the measurements just fell far off the mark. When subjected to fBm signals of H = 0.5 and above, it failed to produce a precise value for each noise and instead returned H 0.4-a point of saturation-for the noise H = 0.5, H = 0.7, and H = 0.9.

Fractal Methods to Determine the Hurst Exponent of Geomagnetic Data
The RRA method may be the oldest method tested in this study, but according to our test, it is still adequate at measuring the overall Hurst exponent of a time series. The results produced by this method did not go that far off the mark, with a standard error level that left a lot to be desired-at worst, the amount of error was ±0.07-but was just acceptable given the age of this method. Moreover, the developments and improvements of various fractal methods over the period of half a century leaves us with plenty of methods to choose from, among them the DFA method.
Despite being the newest method in our test, r-DFA did not turn out to be the best method in this test. The r-DFA method was designed from the ground up, intended to be used for hydrogeological data. The incapability of the method in measuring the exact Hurst exponent of a synthetic geomagnetic time series does not come as a surprise since hydrogeologists undoubtedly have a different set of goals in using r-DFA as a method. In addition, there were several adjustments that had to be made in this test including the number of surrogate sets, the r-DFA order number, and the number of crossovers predetermined.
The number of surrogate sets utilized in this test was 5 sets instead of the supposed 100 sets the originators of the method had performed earlier. As for the r-DFA method, the results of the sixth order were compared with the other methods as the results yielded the nearest reading of the simulated signals. Additionally, there were no crossovers predetermined to exist in our testing. All these adjustments were made for the sake of comparison with the other methods and to reduce biases as the other methods had fewer variables to consider. Still, this study does not discredit the fact that it is an established method in the field of hydrology; however, it is just not adequate for our use in this study. The DFA method, over the years, has been stated as one of the most accurate fractal methods [10,12,14], and apparently it still is. It was found that the DFA method offered the most accurate measurement of Hurst exponent of time series compared to the other methods in this test. Not only that, but it also boasts the least amount of standard error compared to the others, at ±0.01 across all measurements. Therefore, going forward, we utilized the DFA method to characterize the various cases of quiet and disturbed days, in the form of short-term datasets, as well as long-term.  Table 4 shows all the Hurst exponent measurement results of various short-term data of disturbed and quiet days. There are three classifications of events of which results we would compare against one another, namely Disturbed Day (all major events in one particular year), Disturbed Day, and Quiet Day. All H-component data were sourced from DAV and LKW stations and analyzed using DFA imposed with spectral limit of 2.7788< log τ < 4.334 (τ = 6 h until τ = 10 min) following the approach made by previous studies [14,24].   The first class, the 'Disturbed Day (all major events in one year)', is the analysis of days with major geomagnetic events, namely days with Dst < −200 nT-in the year 2013 and 2015, utilizing the H-component data from DAV and LKW stations. All of the disturbed days' data were concatenated into one single dataset according to its year and station, subject to each station's availability of data, causing the variety of the datasets' length to be analyzed using DFA. All of the measurements of the Hurst exponent for all datasets in this class showed persistence tendencies, with data from LKW in particular yielding far higher levels of Hurst exponent compared to DAV station.

Characterization of Geomagnetic Data during Various Cases of Quiet and Disturbed Days
The second class, the 'Disturbed Day', is the analysis of two days in which the lowest point happened to be during a particular major geomagnetic event in the year 2013 and  The first class, the 'Disturbed Day (all major events in one year)', is the analysis of days with major geomagnetic events, namely days with Dst < −200 nT-in the year 2013 and 2015, utilizing the H-component data from DAV and LKW stations. All of the disturbed days' data were concatenated into one single dataset according to its year and station, subject to each station's availability of data, causing the variety of the datasets' length to be analyzed using DFA. All of the measurements of the Hurst exponent for all datasets in this class showed persistence tendencies, with data from LKW in particular yielding far higher levels of Hurst exponent compared to DAV station.
The Overall, the results clearly established that disturbed days showed persistence tendencies while quiet days showed anti-persistence tendencies. This finding reaffirms the long standing theory that disturbed days generally has a higher level of scaling exponent compared to quiet days [12,14,60,61], resulting in a higher level of Hurst exponent. The higher level of Hurst exponent also indicates that the entropy of geomagnetic data strongly decreases during disturbed period [62,63]. It was also found that the results from LKW station were slightly higher than that of DAV station, with different levels of deviation for each class. The results of LKW station tend to deviate further from the results of DAV station in multiple day analyses compared to single day analyses, of which deviations still exist but negligible compared to the former. This deviation may be caused by the difference in the geographic state of each station [64,65], in that DAV and LKW stations are located in different parts of the Southeast Asia region, with LKW station located in the west while DAV station is located in the east. Further study utilizing geomagnetic and local wind data from different longitude sector, such as South America and India, is needed to draw a conclusion to confirm this hypothesis. Table 5 shows all the Hurst exponent measurement results of various long-term data of quiet and disturbed period. For this part of the long-term data analysis, only DAV station data were utilized as this study only focused on comparing two conditions: the long-term quiet day data during low solar activity, the quiet period; and the long-term disturbed day data during intermediate and high solar activity, the disturbed period. The DFA result of 60 quiet days in the year 2009 yielded the Hurst exponent measurement of 0.51 ± 0.03, showing persistence tendencies that bordered on random tendencies (H = 0.50). This finding is in line with the findings of Wanliss [12] and a couple of studies preceding it [66,67] in which these authors argued that the value of the Hurst exponent for the quiet time shall be around H 0.50. It is surmised that the random tendencies of the H-component data are caused by the nature of the geomagnetic data itself which is completely random and barely changes unless there is an external influence that acts upon it. The year 2009 itself was a year of low solar activity, suggesting that external influences were scarce and inconsequential compared to that of high solar activity. Therefore, geomagnetic activities during this time were generally random.
The DFA of 59 disturbed days in the year 2013 and 42 disturbed days that had A-index level above 25 in the year 2013 yielded the same level of Hurst exponent at 0.55 ± 0.02. For the DFA of 60 disturbed days in the year 2015 and 58 disturbed days that had Aindex level above 25 in the year 2015, both analyses also yielded the same level of Hurst exponent at 0.55 ± 0.01. While both years showed persistence tendencies, the periods of solar activity in which both years were situated-intermediate solar activity for the year 2013 and high solar activity for the year 2015-caused little to no difference in the Hurst exponent. These findings reaffirm that disturbed days generally has a higher level of Hurst exponent compared to that of quiet days.
One last thing worth noting in our observations is that the value of Hurst exponent for analyses with shorter length of H-component datasets tend to skew extremely to either having persistence or anti-persistence tendencies, while longer H-component datasets tend to have value nearing random tendencies. Hence, it is surmised that as the H-component dataset gets longer, it tends to have more random tendencies, as evidenced by our findings. This particular characteristic may have a potential to be utilized in geomagnetic storm monitoring, particularly in the Southeast Asia region, using fractal analysis of real-time geomagnetic data as fractal analysis on short-term real-time geomagnetic data may indicate whether geomagnetic storm is still ongoing or otherwise, or showing potential to subside in the near future by referring to the Hurst exponent of the geomagnetic data.

Conclusions
In this work, we have analyzed the long-term quiet day H-component data of the Southeast Asia geomagnetic equatorial region during low, intermediate, and high levels of solar activity. We found that multifractal scaling exists in all of the data analyzed, signifying that the data exhibit fractal properties. We also found the existence of spectral peaks that were present at 6, 8, 12, and 24 h, which were marked throughout all of the years analyzed. This shows that there are external factors that drive the fluctuations of geomagnetic activity that are not limited to solar activity. The Detrended Fluctuation Analysis (DFA) method was found to be the best method to characterize fractal behavior of geomagnetic data. It was found that disturbed days generally showed a higher level of Hurst exponent compared to quiet days. As for the long-term data, disturbed days generally showed a higher level of Hurst exponent compared to quiet days, which had a value near H 0.50. This particular characteristic indicates that the Hurst exponent may have a potential to be utilized with geomagnetic data in monitoring geomagnetic storms in the Southeast Asia region. Additionally, we suggest future studies to identify whether the value of Hurst exponent will vary during different events that perturb daily variations of geomagnetic data such as earthquakes (lithospheric events) and solar flares (space events).