COMPARISON OF CHAOTIC AND FRACTAL PROPERTIES OF POLAR FACULAE WITH SUNSPOT ACTIVITY

To investigate the nonlinear behavior underlying a scalar time series $x=\{{x}_{1},{x}_{2},\cdots ,{x}_{{\rm{n}}}\},$ the first step is to reconstruct an m-dimensional phase space (Packard et al. 1980). Two different coordinate sets can be used, i.e., the delay and the derivative coordinates (Takens et al. 1981). According to Takens' theory, the reconstructed phase space y, constructed mathematically from the time series x, can be defined as:

$$\begin{eqnarray}&&y=\\ (\begin{array}{cccc}x(1) & x(2) & \cdots & x(n-(m-1)\cdot \tau )\\ x(1+\tau ) & x(2+\tau ) & \cdots & x(n-(m-2)\cdot \tau )\\ \cdots & \cdots & \cdots & \cdots \\ x(1+(m-1)\cdot \tau ) & x(2+(m-1)\cdot \tau ) & \cdots & x(n)\end{array})\end{eqnarray} \tag{ 1 }$$

where τ is the time delay and m is the embedding dimension. It is necessary to calculate the appropriate values of these two parameters before reconstructing the phase space. As pointed out by Takens et al. (1981) and Sauer et al. (1991), if the sequence x consists of the scalar measurements of the state of a dynamical system, then the time delay provides a one-to-one image of the original time series x. It should be emphasized that this technique deals with the infinite, noise-free trajectories of a dynamical system, but is not guaranteed to hold for the finite series of the measured data (Schreiber 1999). Therefore, we have to assume that the results of the analysis in this paper are valid for the limited, and not the noise-free, time series.

In this paper, we apply the mutual information method (Fraser & Swinney 1986) to calculate the time delay of each time series. To estimate the embedding dimension, we use an algorithm written by Cao (1997) based on the false nearest neighbor method.

The mutual information metod, which was first proposed by Fraser & Swinney (1986), is used to estimate the time delay of a chaotic time series. According to the information theory, the mutual information of two variables X and Y can be defined as

$$\begin{eqnarray}I(X,Y) & = & {\displaystyle \int }_{X}{\displaystyle \int }_{Y}{f}_{{XY}}(x,y)\\ & & *{\mathrm{log}}_{2}\displaystyle \frac{{f}_{{XY}}(x,y)}{{f}_{X}(x)*{f}_{Y}(y)}{dxdy},\end{eqnarray} \tag{ 2 }$$

where f_XY (x, y) is the joint probability density function of X and Y , and f_X (x), f_Y (y) are the probability density functions of X and Y, respectively. Assuming a partition of the domain of X and Y, the double integral becomes a sum over the cells (Kantz & Schreiber 1997). For a scalar time series $x=\{{x}_{1},{x}_{2},\cdots ,{x}_{{\rm{n}}}\}$ sampled at equal distant interval τ_s, the mutual information between x(i) and x(i + τ_s) is defined as a function of the time delay τ. It can be described as

$$\begin{eqnarray}I(\tau ) & = & \displaystyle \sum _{x(i),x(i+{\tau }_{s})}P[x(i),x(i+{\tau }_{s})]\\ & & *{\mathrm{log}}_{2}\displaystyle \frac{P[x(i),x(i+{\tau }_{s})]}{P[x(i)]*P[x(i+{\tau }_{s})]},\end{eqnarray} \tag{ 3 }$$

where P[x(i)] and P[x(i + τ_s] are the probability density functions of the time series x(i) and x(i + τ_s), respectively. An appropriate choice of the time delay for the embedding is one such that the mutual information is at a minimum. That is to say, if the mutual information exhibits a marked minimum at a certain value of τ, then it is a good candidate for a reasonable time delay. More information can be found in the papers of Fraser & Swinney (1986), Kantz & Schreiber (1997), and Hegger et al. (1999).

Figure 2 displays the results of the mutual information method as a function of time delay. The graphs are given for polar faculae (left panel), sunspot numbers (middle panel), and sunspot areas (right panel). From this figure, one can see that the profiles of the mutual information method for three time series slightly differ from each other. The reasonable time delays are found to be 31 months for polar faculae, 29 months for sunspot number, and 31 months for sunspot areas.

Zoom In Zoom Out Reset image size

As the original attractor of the considered system is unknown, one should attempt to reconstruct an attractor that preserves the invariant characteristics of the original attractor. Based on Taken's theorem, an m-dimensional phase space can completely capture the nonlinear behavior of a given time series under the condition m ≥ d, where d is the embedding dimension of the original attractor (Karakasidis & Charakopoulos 2009). For an infinite and noise-free time series, the embedding dimension m can be determined by Takens' theorem. As proved by Takens' theorem, it is sufficiently to reconstruct the attractor under the condition m ≥ 2d + 1. However, it should be emphasized that the condition m ≥ 2d + 1 is sufficient but not necessary, because an original attractor may sometimes be reconstructed triumphantly with an embedding dimension as low as d + 1. If the embedding dimension m is too low compared to the dimension d of the attractor, the attractor displays self intersections, spatially nearby points on the attractor are either real neighbors due to the system's dynamics or false neighbors due to the self intersections. In a higher dimension, where the self intersections are resolved, the false neighbors are revealed as they become more remote. One tries to find a threshold value for m where no false neighbors are identified as one moves to increased dimensions. Here, we used the false nearest neighbor method, which was proposed by Kennel et al. (1992), to determine the minimal sufficient embedding dimension m. Using Takens' theorem one can estimate the dimension of a possible attractor.

The fundamental principle of the false nearest neighbor method is to increase the dimension of the phase space up to the point where there are no longer any self intersections of the trajectory. Based on Takens' theorem, Cao (1997) defined a parameter E₁(d)—the relative change in the average distance between the two contiguous points—when the dimension is increased from d to d + 1. If the relative changes in the average distance saturate around 1, then the embedding dimension is equal to the minimal dimension. Moreover, another parameter E₂(d) was proposed to discriminate the deterministic system from the stochastic process. For a random process, the future values are independent of the past values, so the values of E₂(d) are equal to one for any dimensions. For a chaotic system, the values of E₂(d) are not constants for all dimensions and they are related to the embedding dimension of the considered system. It should be pointed out that different time delays may bring out different embedding dimensions, especially for the time series from continuous temporal–spatial systems (Cao et al. 1998). An appropriate choice of the time delay will decrease the embedding dimension, which is necessary for phase-space reconstruction, as suggested by Letellier et al. (2006). That is to say, choosing a range of the time delays for which the embedding dimension does not change constitutes a suitable parameter that actual properties of the dynamical behavior are identified (Rhodes & Morari 1997). For more details of the false nearest neighbor method, the interested reader is referred to Kennel et al. (1992), Cao (1997), and Cao et al. (1998).

Figure 3 displays the variations of the parameters E₁(d) and E₂(d) as a function of the dimensions for polar faculae (left panel),  sunspot numbers (middle panel), and sunspot areas (right panel). From this figure one can see that the profiles of E₂(d) for three time series are not straight lines, providing strong evidence that the three solar time series analyzed in this paper exhibit a chaotic behavior. Moreover, their embedding dimensions are all equal to three, implying that both the high- and low-latitude solar activity are governed by a low-dimensional chaotic attractor, i.e., a three-dimensional phase space should be sufficient to properly unfold the underlying nonlinear behavior of the Sun. Previous studies have shown that the attractor dimension (D₂) of the low-latitude solar activity (e.g., sunspot numbers and sunspot areas) is around 3, such as ${D}_{2}\approx 2.3$ (Mundt et al. 1991), D₂ = 2.4  ±  0.2 (Kremliovsky 1994) and D₂≈3 (Letellier et al. 2006). Therefore, the embedding dimension of sunspot activity is consistent with that reported by previous authors. Furthermore, our result provides convincing evidence that the chaotic behavior of the high-latitude solar activity should also be low-dimensional.

Zoom In Zoom Out Reset image size

After estimating the time delay and embedding dimension of each time series, in Figure 4 we plot the reconstructed phase spaces for  polar faculae (left panel), sunspot numbers (middle panel), and sunspot areas (right panel). From this figure, one can see that the phase portraits for the three time series are not like "balls of wool" (i.e., totally unfolded), and there seem to be some intrinsic structures underlying the considered systems. Here, we must keep in mind that the smoothing procedure can remove part of the original dynamics (the ideal situation is to remove only the stochastic component), but it cannot inject the nonlinear dynamics into the time series (Letellier et al. 2006). The stranger attractors characterized by the infinite self-similarity structures are clearly revealed, which are the one of the typical properties of the chaotic system. By comparing their complexities shown in the reconstructed phase spaces, we find that the chaotic dynamics of polar faculae are the most complicated, followed by that of the sunspot areas, and then the sunspot numbers. This is the first clue to support the hypothesis that the chaotic dynamics of high-latitude solar activity is more complicated than low-latitude solar activity.

Zoom In Zoom Out Reset image size

Historically, the first approach to quantify the long-range persistence of a time series was proposed by Hurst (1951), who spent his life studying the hydrology of the Nile River, in particular, the record of floods and droughts. To better reveal the long-term correlation of his empirical data, he introduced R/S analysis. This approach is based on the random-walk theory. Its predominant concept was developed at a time (1) when computers were in their early stages, so calculations had to be done manually, and (2) before fractional noises or motions were introduced. Currently, the use of R/S analysis is still popular and often applied. In this section we apply R/S analysis to investigate the long-range correlation and the fractal property of the above three time series.

R/S analysis is a simple but a strong nonparametric method for fast fractal analysis (Das et al. 2009). It is performed on the discrete time series X(i) of dimension N by calculating three factors: first the mean $\langle X\rangle ,$ second, R, which is the range of commutative difference of X(t, s) at discrete time t over a time span s, and third, the standard deviation S(i), estimated from the observed values of X(t, s), where

$$\begin{eqnarray}&&\langle X\rangle =\displaystyle \frac{1}{N}\displaystyle \sum _{t=1}^{N}X(i)\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&S(i)={\{\displaystyle \frac{1}{s}\displaystyle \sum _{i=1}^{s}[X(i)-\langle X\rangle \}}^{2}{}^{1/2},s=1,\ldots ,N.\end{eqnarray} \tag{ 5 }$$

The range of cumulative difference of the time series is given by

$$\begin{eqnarray}&&R(i)=\mathrm{Max}[Y(i)]-\mathrm{Min}[Y(i)]\end{eqnarray} \tag{ 6 }$$

where Y(i) is defined as

$$\begin{eqnarray}&&Y(i)=\displaystyle \sum _{k=1}^{i}[X(k)-\langle X\rangle ].\end{eqnarray} \tag{ 7 }$$

Hurst found that the ratio (R/S) is well described for the long-term time series of natural phenomena by the following experimental relationship (Hurst 1951):

$$\begin{eqnarray}&&\langle {\rm{R}}/{\rm{S}}\rangle ={s}^{H},\end{eqnarray} \tag{ 8 }$$

where s and H are the time span and the Hurst exponent, respectively. The slope of the plot of R/S versus the time span s on log–log plot gives rise to H. The value of H indicates whether a time series is random or successive increments in time series are not independent. For H equals to 0.5, the time series represents a random walk, or is uncorrelated, and each observation is fully independent of all prior observations. The value of 0 < H < 0.5 indicates that it exhibits an anti-persistence behavior, i.e., its trend will likely reverse in the future. The value of H lying between 0.5 and 1 implies persistent time series characterized by long-range effects. The successive increments are positively correlated with the preceding observations (Kilic & Golbasi 2011).

To characterize the type of the self-affinity and the fractal property of solar time series, the Hurst exponents of the polar faculae, the sunspot numbers, and the sunspot areas are calculated by R/S analysis. Figure 5 shows the log–log plots of R/S as a function of the time span s, which estimates the values of the Hurst exponent H = 0.7692  ±  0.0525 for polar faculae, H = 0.7400  ±  0.0684 for sunspot numbers, and H = 0.7729  ±  0.0672 for sunspot areas, respectively. Since the H values of the three time series are larger than 0.5, this indicates that both the high- and low-latitude solar activity show persistence on large scales. The Hurst exponents of the three time series analyzed in this paper are agreement with those reported by previous authors, whose studies focused on the sunspot numbers (H = 0.86 ± 0.05; Mandelbrot & Wallis 1969; H = 0.8033; Zhou et al. 2014), the sunspot areas (H = 0.7834; Zhou et al. 2014), the solar radius (H = 0.72; Kilic & Golbasi 2011), Hα flare activity (H = 0.73 ± 0.01; Lepreti et al. 2000), and the cosmogenic radiocarbon data (H = 0.84; Ruzmaikin et al. 1994).

Zoom In Zoom Out Reset image size

More importantly, the Hurst exponent H of a time series is directly related to its fractal dimension D, and the mathematical relationship between these two parameters is D = 2 − H (Das et al. 2009). The fractal dimension D is a statistical quantity that measures the dynamical complexity of the chaotic time series. Using this equation, the fractal dimensions are found to be 1.2308 ± 0.0525 for polar faculae, 1.2600 ± 0.0684 for sunspot numbers, and 1.2271 ± 0.0672 for sunspot areas. Our results suggest that both the high- and the low-latitude solar activity have distinctly fractal natures. Furthermore, the values of D show that the fractal nature of the solar activity is due to its long-range correlation.

After calculating the fractal dimensions D of the three time series, their corresponding predictability indicators P can be determined from the statistical relationship $P=2| D-1.5| .$ If the value of P of a time series is close to zero, the considered system will be unpredictable. However, if the value of P becomes close to one, the dynamical system can be taken as a highly predictable process (Ghosh et al. 2014). Using this equation, the predictability indicators can be evaluated as 0.5384 ± 0.1050 for polar faculae, 0.4800 ± 0.1368 for sunspot numbers, and 0.5458 ± 0.1344 for sunspot areas, respectively. The P values of the three time series are around 0.5, indicating that the solar activity cycle is predictable in nature. It should be emphasized that the analysis results do not provide any concept on the predictable timescale. Subsequently, the method named the largest Lyapunov exponent is utilized to describe the dynamical complexity and to calculate the predictability timescale of the three time series.

A chaotic system is a deterministic system that is usually characterized by the fractal structure and sensitivity to small changes in initial conditions. The chaotic property of an attractor can be quantitatively described by many chaotic variables, such as the largest Lyapunov exponent, the correlation dimension, the Kolmogorov entropy, and so on (Hegger et al. 1999). In the Lyapunov exponent analysis (Eckmann et al. 1986), the quantities investigated are assumed to follow a set of differential equations in a high-dimensional phase space. These unknown equations are usually nonlinear and the apparent randomness of the generated series may reflect chaotic solutions of this set of nonlinear equations. The distance between two initially close orbits evolves exponentially, with a characteristic exponent λ, known as the Lyapunov exponent. The Lyapunov exponent is a measure of the rate at which nearby trajectories in phase space diverge. Moreover, it measures how predictable or unpredictable the system is. For chaotic orbits, there is at least one positive Lyapunov exponent. The positive Lyapunov exponent is a necessary condition for deterministic chaos but not sufficient. For periodic orbits, all Lyapunov exponents are negative. The Lyapunov exponent is zero near a bifurcation, which reflects the fact that the orbits evolve slower than exponentially.

The largest Lyapunov exponent, defined as the divergence of nearby trajectories in the phase space, is a convincing indicator to characterize the chaotic behavior and the dynamical complexity of a time series. Rosenstein et al. (1993) proposed an algorithm to calculate the largest Lyapunov exponent: after reconstructing the phase space, choose a point x_no and all neighbor points x_n, then calculate their average distance from that point, and finally estimate the average quantity S (stretching factor) by repeating the above procedures for N numbers of those points along the orbits:

$$\begin{eqnarray}&&S=\displaystyle \frac{1}{N}\displaystyle \sum _{\mathrm{no}=1}^{N}(\mathrm{ln}\displaystyle \frac{1}{| {u}_{{x}_{\mathrm{no}}}| }\sum | {x}_{\mathrm{no}}-{x}_{{\rm{n}}}| ),\end{eqnarray} \tag{ 9 }$$

where $| {u}_{{x}_{\mathrm{no}}}|$ is the number of neighbor points around x_no. The plot of the stretching factors S versus the number of points N (or time interval t = NΔ t) will produce a curve with a rapid increase at the beginning and followed by a relatively flat region. Hegger et al. (1999) described the implementation of methods of nonlinear time series analysis, which are based on the paradigm of deterministic chaos and fractal behavior. A range of these algorithms has been made available in the form of computer programs by the TISEAN package⁸ (for details, see Hegger et al. 1999).

In Figure 6, we plot the variations of the divergence as a function of the time interval for polar faculae (left panel), sunspot numbers (middle panel), and sunspot areas (right panel), respectively. From this figure, one can see that the profiles of the three curves display expected linear increases and relatively flat regions with some fluctuations superimposed on the part of the curves. The slopes (corresponding to the largest Lyapunov exponent for the three time series) are estimated through the least-squares line fitting method. We find that reasonable values of the largest Lyapunov exponent are 0.1429 ± 0.0304 month⁻¹ for polar faculae, 0.0654 ± 0.0151 month⁻¹ for sunspot numbers, and 0.1198 ± 0.0054 month⁻¹ for sunspot areas. These positive values indicate the possible existence of chaotic behavior in the solar activity. This result is consistent with the estimation of a low-dimensional chaotic attractor from the phase-space reconstruction and the false nearest neighbor method, as already discussed above. According to the information theory method, the larger the largest Lyapunov exponent, the more complex the dynamical behavior of the system. Our results suggest that the dynamic behavior or the chaotic degree of the polar faculae and the sunspot areas is more complex than that of the sunspot numbers. The predictability timescale, defined as the inverse of the largest Lyapunov exponent, is taken as the length of the strange attractor. In this sense, the predictability timescales are calculated to be 0.6108 ± 0.1299 year for polar faculae, 1.3443 ± 0.3091 year for sunspot numbers, and 0.6996  ±  0.0.0552 years for sunspot areas. Therefore, the effective timescale used for forecasting the solar activity cycle is smaller than three years, i.e., the high-accuracy prediction of the solar magnetic activity should only be done for short- to mid-term due to its intrinsically dynamical complexity.

Zoom In Zoom Out Reset image size

Understanding the nonlinear behavior of a time series plays an integral role in the detection of deterministic chaos in physical systems (Schreiber & Schmitz 2000). One approach is to formulate a null hypothesis for a specific class of processes and compare the system output to this hypothesis. The surrogate data method is a classical way to establish such a null hypothesis and can be undertaken in two distinct ways. Typical realizations are Monte Carlo generated surrogate time series from a model that provides a suitable fit to the data; constrained realizations are surrogate time series generated from the data values to conform the certain properties of the data (Suzuki et al. 2007). The latter approach is more suitable for hypothesis testing as it does not require the definition of a pivotal test statistic (Theiler & Prichard 1996). Theiler et al. (1992) suggested some algorithms for generating surrogate data, such as the Fourier transform (FT) algorithm, the windowed Fourier transform algorithm, the amplitude-adjusted Fourier transform (AAFT) algorithm, and so on. In this paper, we generate 30 surrogate time series, using the AAFT algorithm, that have the same mean, variance, and autocorrelation function as the original time series. The AAFT algorithm can be described by the following three steps:

(1). We produce Gaussian random time series η(t)  ∼  N(0,1) and shuffle the temporal order to obtain the same rank order as with the original time series X(i). Here, the rank order is defined as an order of state values of a time series (Suzuki et al. 2007). In this process, the shuffled time series ω(i) corresponds to a linear stochastic process and X(i) is considered to be a realization of a transformation of ω(i) through a static monotonic nonlinear function h( · ). That is, X(i) = h(ω(i)).

(2). Since the shuffled time series ω(i) is a realization that can be characterized by its power spectrum, we apply the FT algorithm to ω(i). The FT algorithm is as follows: the FT is applied to a target time series, in this case ω(i), and its power spectrum is obtained. Next, preserving a symmetry of phase components of the FT, the phase components are randomized. Then, the inverse FT generates an FT surrogate, which preserves the power spectrum of the target time series ω(i). Namely, in this process, we make an FT surrogate $\hat{\omega }(i)$ of ω(i), which has exactly the same power spectrum as ω(i).

(3). We shuffle the original time series X(i) to obtain the same rank order as the FT surrogate time series $\hat{\omega }(i)$ of ω(i) given in the previous step. Here, we again consider a static monotonic nonlinear function h( · ). Namely, $\hat{X}(i)$=h($\hat{\omega }(i)$), and $\hat{X}(i)$ is called the AAFT surrogate data version of X(i).

Based on the AAFT algorithm, we generate 30 surrogate time series each of polar faculae, sunspot numbers, and sunspot areas.  We then calculate the Hurst exponents and the largest Lyapunov exponents of each time series, as shown in Tables 1 and 2. From Table 1, the average values of the 30 Hurst exponents are H = 0.7765  ±  0.0444 for polar faculae, H = 0.7786  ±  0.0496 for sunspot numbers, and H = 0.7945  ±  0.0491 for sunspot areas. Using the mathematical relationship among the Hurst exponent, the fractal dimension, and the predictability indicator, the fractal dimensions and the predictability indicators are calculated to be D = 1.2235  ±  0.0444 and P = 0.5530  ±  0.0888 for polar faculae, D = 1.2214  ±  0.0496 and P = 0.5572  ±  0.0992 for sunspot numbers, and D = 1.2054  ±  0.0491 and P = 0.5891  ±  0.0982 for sunspot areas. From Table 2 aj521731t2, the average values of the 30 largest Lyapunov xponents are 0.1148 ± 0.0292 month⁻¹ for polar faculae, 0.0675 ± 0.0236 month⁻¹ for sunspot numbers, and 0.0950 ± 0.0225 month⁻¹ for sunspot areas. The predictability timescales are calculated to be 0.7761 ± 0.1974 year for polar faculae, 1.4065 ± 0.4918 year for sunspot numbers, and 0.9293 ± 0.2201 year for sunspot areas.

All analysis results are summarized in Table 3. From this table, one can see that the analysis results obtained by the surrogate data method are very similar to those given by the above sections. Therefore, it is safe to say that the quantitative analysis results through the nonlinear time series analysis techniques used in this paper have a high degree of credibility. Based on the analysis results shown in Table 3, we make the following conclusions: (1) both the high- and the low-latitude solar activity show persistence on large scales and have distinctly fractal natures. Moreover, the fractal nature of the solar activity is due to its long-range correlation; (2) the solar magnetic activity cycle is predictable in nature, but the high-accuracy prediction should only be done for short- to mid-term due to its intrinsically dynamical complexity; and (3) the dynamical behavior or the chaotic degree of polar faculae is the most complicated, followed by sunspot areas, and then sunspot numbers because the largest Lyapunov exponent of polar faculae is the largest, which is in agreement with the results shown in the reconstructed phase spaces (Figure 4).