Analysis of fractality and complexity of the planetary K-index

The objective of this research is to explore the inherent complexities and multifractal properties of the underlying distributions in the daily Planetary K-index time series collected from NOAA Space Weather Prediction Center. In this article, non-stationary and nonlinear characteristics of the signal have been explored using Smoothed Pseudo Wigner–Ville Distribution and Delay Vector Variance algorithms, respectively, while Recurrence Plot, 0–1 test, Recurrence Quantification Analysis and correlation dimension analysis have been applied to confirm and measure the chaos in the signal under consideration. Multifractal detrending moving average has been used to evaluate the multifractality and also recognise the singularities of the signal. The result of these analyses validates the nonstationary and nonlinear characteristics of the Planetary K-index signal, while a significant presence of deterministic chaos in it has also been noticed. It has also been confirmed that the Planetary K-index exhibits multifractal nature with positive persistence. The long-range temporal association and also the large pdf are discovered to be the primary factors that contribute to the multifractal behaviour of the Kp-index.


Introduction
Immense disturbances in the geomagnetic field occur mainly due to the impact of solar storms (like solar flares, solar wind, coronal mass ejection (CME)) on the Earth's magnetosphere. This geomagnetic field perturbation, commonly termed as geomagnetic storms, has a significant impact on the backbone of our modern day civilisation, power grids, electronic devices, navigation systems, spacecraft operations and global communications [1,2]. A geomagnetic storm (GMS) is the most evident manifestation of magnetospheric phenomena, and it can be conceived of as a series of nonlinear electromagnetic processes propagating from the Earth into the entire magnetospheric zone. The magnitude of the geomagnetic storms can be characterised and quantified by the value of Planetary K-index [3].
K p -index is a measure to quantify the fluctuation of the magnetic field of our planet, more specifically it is horizontal component. K index is an integer number lying between 0 and 9. For a value of less than 5, our planet's magnetic field is considered calm; otherwise, it is a geomagnetic storm. The device used to evaluate the variations of the horizontal component in geomagnetic field is a magnetometer. The official K p -index is estimated by capturing a weighted average of K-indices from a specified system of geomagnetic observatories, where the K-index is a code that refers to the maximum horizontal component variations found on a magnetometer over a three-hour period. The planetary K-index can be considered as a global indicator of geomagnetic activity around the globe.
Due to the severe influence of the geomagnetic storm on human civilisation, there is a need to learn about the dynamical behaviour of it. The study of the different characteristics and complexity of Planetary K-index signal with an investigative approach may help to infer about the actual nature of the geomagnetic storm's dynamics. NOAA Space Weather Prediction Center derives the estimated 3-h Planetary K-index using the data from the ground-based magnetometers located in different countries. The daily magnitude is calculated from the averages of the eight 3-h K p -indices per day. Here, regular K p -index data over the period of February 1999 to December 2007 has been taken for analysis. 1 There are about 88 geomagnetic storm events that have occurred in solar cycle 23. Only a few of these 88 storms have a significant impact on our Earth's environment. On the 4-7 October 2000, there are records of extreme GMSs that had a major effect on the ionosphere. Though the ion temperatures of the terrestrial magneto-tail generally used to remain consistent, it has been noticed that in October 2000 the estimated temperature of ion was 2-3 times more compared to the average value [4]. Another geomagnetic storm event was noticed on 21 October 2001, and its effect has been explained in Jordanova et al. research paper [5]. They have observed massive loss of electrons in the radiation belt into the atmosphere and electron flux dropout due to the outward radial diffusion. The massive ionospheric disturbance has been registered for geomagnetic Sudden Storm Commencement (SSC) occurred on [29][30] October 2003 which has yielded a 'swirling' effect on these days. On 20-22 November 2003, a super magnetic storm hit the Earth's atmosphere that generated ionospheric disruptions in the southeast region of Asia, causing GPSnavigation to stop for several hours. Also, these GMSs on 29 October 2003, 30 October 2003and 20 November 2003 produce very high value geomagnetically induced currents (GIC) of 57.05 Amp, 48.57 Amp, 23.86 Amp, respectively [6]. It is well known that high-value GIC is always a threat to our power grids. One such well-known GIC event happened in March 1989 caused a province-wide blackout in Quebec, Canada. An intense GMS triggered on 24 August 2005 with a maximum K p -index value of ≈ 9 while another strong geomagnetic storm triggered on 11 September 2005 with a peak K p -index value of ≈ 8 . These disturbances modulated the galactic cosmic rays (GCR) which produced a severe Forbush effect (FE) on 24 August and 11 September 2005. Papaioannou et al. [7] discuss the major effects in the Earth's magnetosphere from 22 August to 17 September 2005 and also calculate cosmic ray gradients. On 14-15 December 2006 the geomagnetic storm events were triggered when Earth's magnetosphere was affected by the CME-associated interplanetary shock. Ionospheric disturbances have been recorded due to these geomagnetic storm events which produce fluctuations in electron densities and also produce GIC [8]. In space research, the Earth's magnetosphere is always an exciting domain for researchers to explore. Various researchers have conducted various studies to determine the type of interrelationship between GMSs and several solar activities [9][10][11]. However, in this article, an attempt has been made to discover the nature of GMS fluctuation by the analysis of the K p -index signal. The prime objective of this work is to surface the behavioural phenomena like stationarity, nonlinearity, chaos, complexity, self-similarity and multifractality of the daily K p -index as represented in Fig. 1 and hence the geomagnetic storm. For this purpose, appropriate statistical tools have been used.
Non-stationary components are inherently present in real-world signals when they change over time. From the spectral approach, the stationary signal is one where the frequency content does not change with time, while it changes in the case of a non-stationary signal. Classical spectral analysis can only detect the frequencies within the signal, but the specific time of presence of these frequencies cannot be located and can therefore not be used to characterise the non-stationary signal. Therefore, a technique based on time-frequency representation (TFR) can be used for natural signals that are usually non-stationary in nature. Various well-established TRF based method like short-time Fourier transform (STFT) [12], wavelet transform (WT) [13], and Wigner-Ville distribution (WVD) [14], etc., can be used to confirm the non-stationary characteristics of the signal. The STFT method generates time-frequency spectrum by taking FT, while this method restricts temporal and spectral resolution by a pre-specified window length [15]. SPWVD [16], a TFR-based technique that is basically a modified WVD algorithm, has been applied in this analysis to validate the non-stationary nature of the signal under investigation. WVD has variable resolution along the time-frequency plane by having decent temporal resolution and frequency resolution at high and low frequencies, respectively, while the WVD method generates cross-term for a multi-component signal. These unwanted components can be suppressed by applying the SPWVD method where comparatively high resolution is achieved with independent time and frequency function as kernel function. Therefore, the SPWVD method has been chosen here to validate the nonstationary nature of the signal under investigation by time-frequency energy mapping. The features like nonlinearity and chaosity give an indication about the behavioural complexity of the signal. Nonlinear characteristics of a time series are traditionally being judged by signal processing tools like deterministic versus stochastic method, correlation exponents, etc. In this work, the nonlinear characteristics of the signal have been revealed by applying DVV algorithm [17]. It is necessary to detect the presence of any chaotic behaviour in the signal in order to obtain a full understanding of the K p -index's nonlinear dynamics. '0-1 test' which is a better alternative compared to the traditional maximal Lyapunov exponent method, has been chosen to investigate the trace of the chaotic phenomenon in the signal [18,19]. Whether the chaos, which may be present in the K p -index dynamics, is deterministic or stochastic is being determined by the correlation dimension analysis [20,21]. The structural complexity (the scaling properties and self-similarity) of the dynamics of the geomagnetic storm dynamics has been explored by fractal analysis of the K p -index signal [22]. The widely used methodologies for fractal analysis are rescaled adjusted range analysis (R/S) method by Hurst [23], the detrended fluctuation analysis (DFA) method introduced by Peng et al. [24]. However, R/S and DFA techniques are not useful for the multifractal signal. The drawbacks of these techniques have been overcome by methods, like wavelet transform modulus maxima (WTMM) [25] and multi-fractal detrended fluctuation analysis (MFDFA) [26] where MFDFA is more robust compared to the former [27]. Recently, a modified detrending moving average (DMA) method [28] has been developed by Gu et al. [29], known as multifractal detrending moving average (MFDMA). In this investigation, the backward MFDMA method which is advantageous than MFDFA [29][30][31] has been used for the fractal analysis of the signal.
The possible cause of the multifractality has been identified by performing fractal analysis of the shuffled and surrogated form of the signal. The major contributory factors for which a dynamical system is nonstationary, chaotic and multifractal is the hidden recurrent patterns and structural changes in it. Such recurrences are the fundamental characteristic of many dynamical systems. Recurrences provide all pertinent details on the behaviour of the system. To visualise recurrences of the K p -index, recurrence plot (RP), proposed by Eckmann et al. [32], is applied in this paper. A recurrence plot is the graphical representation of the recurrence matrix of the trajectory of the dynamical system in its phase space. RP gives the first notion about the patterns of recurrences which will let to study the underlying dynamics of the process and its trajectory. With the intention to go beyond the visual realisation by RPs, some other measures of complexity that quantify the small scale structures in RPs have been proposed in [33], and their analysis is recognised as recurrence quantification analysis (RQA). The quantification procedure is based mostly on concentration of recurrence points as well as the RP's diagonal/vertical line formations. The RQA parameters of the K p -index have been estimated for the understanding of its scale complexity. The SPWVD method, DVV, correlation dimension analysis, 0-1 Test method, MFDMA algorithm and RP analysis are discussed in Sect. 2. Section 3 provides the relevant findings and observations derived from the Planetary K-index signal applying the signal processing tools described above, while this paper concludes with inferential comments made in Sect. 4.

SPWVD analysis
In this work, Smoothed Pseudo Wigner-Ville Distribution (SPWVD) has been chosen over various other TFR-based methods to evaluate the dependency of the frequency of the signal on time, i.e. stationary or non-stationary. The widely used TFR-based methods like short-time Fourier transform or wavelet transform do not offer sufficient crisp resolution like SPWVD analysis [34]. In 1932, the WVD method had been proposed by Wigner and later modified by Ville in 1948 which is expressed as where s(.) represents the signal achieved by using the Hilbert transform to the input signal, f represents frequency, τ is the lag variable.
Interferences in the time and frequency domain have been identified in the WVD because of its bilinear structure [35]. Time domain cross terms are attenuated by using the Pseudo-WVD (PWVD) method which is the convolution of WVD and a regular window function h(t) [36]: whereas the interferences in the frequency domain are attenuated by smoothing the PWVD where it is made to pass through a low-pass filter having window function (t) . This yields smoothed-PWVD (SPWVD) [37,38].
For stationary signal, the frequencies content of the signal will be found to exist continuously along the entire time axis of the SPWVD time-frequency spectrum, whereas for nonstationary, the frequencies of the signal will exist sporadically with time.
The frequency components of a stationary signal would be observed to appear continually over the entire time axis of the SPWVD spectrum, while the frequency components of a nonstationary signal would occur sporadically over time.

DVV analysis
For a given optimal embedding dimension m and embedded lag τ, a delay vector set (DVs) The DVs which are clustered within the pair-wise Euclidian distance of d from DV x(t) , produce the sets t d . d lies between the ranges of min 0, d − n d d to d + n d d where d denotes standard deviation, d is the mean value and n d represents controlling parameter. For every set t d , the variances of the corresponding targets 2 t d are computed. The target variance * 2 d is determined by taking the mean of these variances and normalising them to the variance of signal x(k) . The plot of the target variance as a function of the standardised distance is known as a "DVV plot". In order to standardize the distance axis, d is replaced by ( d − d )∕ d . d and d are the mean and the standard deviation of each DV, respectively. Randomness can be measured by determining the least value of target variance * 2 . A surrogate data series has been computed using iterative amplitude-adjusted Fourier transform (IAAFT) technique [10] as well as "DVV plots" for both surrogate and natural signal are compared. In case, "DVV plots" shows no difference between surrogates and original signal and then, original signal is considered linear. "DVV scatter diagrams" is the graphical mapping of the target variances of the surrogate signal to the target variances of the original signal for corresponding standardised distances. If the dynamical system is linear, the mapping locus bisects the plot diagonally; otherwise, it will be significantly different from the diagonal bisector line then. The deviation of the locus from the bisector line in terms of root mean square error (RMSE) value defines the degree of nonlinearity [39]. [40][41][42][43][44] technique has been applied to confirm the existence of chaotic nature in the signal under investigation. For measured data, this tool gives better result compared to the long-established maximal Lyapunov  [18]. This binary test method is straightforward and effective to judge the chaosity of a system. If the output of the test is towards "1", then the chaotic behaviour of the signal has been confirmed, whereas if the output has an inclination towards "0", the nature is not chaotic.

0-1 Test
Let, x(k) is a time series of length L and X l is its Fourier transform. The inherent characteristics of the x(k) can be judged as regular or chaotic depending upon the relation of smoothed mean square displacement d c (l) with l. In case signal has chaotic dynamics then d c (l) will continuously increase with l and X l going to produce a Brownian motion in complex plain. The smoothed mean square displacement d c (l) can be expressed as x(k) and l ≤ m = L∕10 < L . 100 values of c ranging from π/5 to 4π/5 are selected [43].
To have an optimum result for noise effected signals, The term indicates the test sensitivity. The correlation between l and d * c (l) gives the asymptotic growth rate K c = corr l, d * c (l) . The median of K c is the final output for 0-1 test analysis.
Based on the value of the final output K, the time series may be judged as regular or chaotic. The signal can be categorised as normal or chaotic depending on the value of the final output K.

Correlation dimension analysis
Correlation dimension analysis [20] is an important indicator which not only informs the type of chaos (deterministic or stochastic) existing in a dynamical system but also enlightens about the dimension of the system. The correlation dimension D(r, m) quantifies the ability to occupy the phase space of a time series (with embedding dimension of m and time lag τ) by the cells of size r where correlation integral C(r, m) is the probability of any two randomly chosen points within a given distance r.
where r is the distance among coordinates i and j for ith and jth points, respectively, in phase space for a set of N number of points and H is the Heaviside function.
The mean of D(r, m) for all r values is calculated to yield D(m) for every m = 1, 2, … , L while L ≤ 10 . For a stochastic time series D(m) exhibits a monotonically increasing trend towards infinity with m. However, if the process is deterministic, this continuous increment of D(m) ceases at some small m where it saturates in a plateau [45]. This saturated value of D(m) gives the correlation dimension. So, stochastic systems are 'infinite-dimensional' and deterministic systems are 'finite-dimensional' .

MFDMA
A signal can be divided into several sub-signals. Now, if each sub-signal is a statistical copy of the original signal, then these sub-signals can be referred as fractals, i.e. fractal systems can be typified by self-similarity. MFDMA can be considered as an effective method for evaluating the fractal characteristics for any nonstationary time series signal [29,46].
Let U(k) be the cumulative sums for the signal x(k) of length N. The moving average functions Ũ (k) , with the moving window length of l, can be expressed as [47] where ⌊ ⌋ are the highest integer and ⌈ ⌉ smallest integer above . Position parameter denoted as θ lies between 0 and 1. For the values of = 0, 0.5 and 1 , MFDMA method can be categorised into three specific cases like backward moving average method, centred moving average method, forward moving average method, respectively. It is established in a paper [48] that out of these methods, backward moving average ( = 0) based MFDMA outperforms the other two. In this work, the backward moving average-based MFDMA has been used for fractal analysis of the Planetary K-index signal.
The difference between Ũ (k) and U(k) gives the residual series (k) as The root mean square function F v (l) can be obtained as where v (k) = (n + k) for 1 ≤ k ≤ l and n = (v − 1)l.
The q th -order overall fluctuation function F q (l) is If the h(q) is the Hurst index, then the power-law relation between the function F q (l) and the size scale l can be established as: F q (l) ≈ l h(q) . Also, the multifractal behaviour of a signal can be characterised by determining the scaling exponent (q) which can be defined as: (q) = qh(q) − 1 . Using the Legendre transform, the singularity strength function and the multifractal spectrum f ( ) are expressed as Eqs. (14) and (15), respectively [29].

RP and RQA
The recurrence plot (RP) and recurrence quantification analysis (RQA) are summarised as follows:

RP
All RP method is the graphical representation of a twodimensional square matrix where a matrix element is located at (i, j) if a state at time i recurs at a different time j. If a signal trajectory ⃗ i N i=1 is in its phase space, then the mathematical equation of RP can be defined as [39] where H(•) represents Heaviside step function and threshold distance is denoted as . For i = j , RP i,j = 1 which generate main diagonal line of the RP called line of identity (LOI).
The system dynamics can be explored by visualising the structural pattern of the RP. The typical pattern and their interpretation are given in Table 1.

RQA
After getting an idea about the underlying system dynamics by visualising the RP plot, there is a need to determine different recurrence variables which are revealed by the plot to measure the complexity [33,49]. Among many other recurrence variables, four variables are essential for the quantification of recurrences, like %REC, LMAX, %DET and ENTR.
The density of recurrence points in the RP can be measured by calculating the recurrence variable %REC, given as Signal frequency changes with time, i.e. the signal generated is from a non-stationary system. The non-stationarity is due to existence of a trend or a drift in the system 3 Disruptions The system is non-stationary. The non-stationarity is due to the existence of some rare states in the system. Transition may have occurred 4 Periodic/quasi-periodic patterns The system has cyclicities. For quasi-periodic system, the distances between long diagonal lines is varying 5 Single solitary points There is strong fluctuation in the signal, i.e. random process 6 Diagonal lines (parallel to LOI) Signal generating system is deterministic. If there are isolated single points along with the diagonal lines, the process can be judged as chaotic 7 Diagonal lines (orthogonal to LOI) Here, also the states evolve equally at various times but in reverse order. It may be to insufficient embedding 8 Vertical and horizontal lines/clusters Time evolution of the states is either constant or varies slowly indicating the presence of laminar states 9 Long bowed line structures Here also the states evolve equally at various times but with differing velocities Another RQA measure is the LMAX which can be defined as the longest diagonal line present in the RP plot. If ith diagonal line length is l i, then LMAX can be expressed as where the number of diagonal lines is denoted by N l . Shorter LMAX signifies the rapid divergence of the trajectory segments, whereas longer LMAX denotes that the trajectory segments diverge slowly [50].
%DET is the measure of determinism of the system quantified by the density of lines, with length l l ≥ l min , which are aligned diagonally [51]. If there are no diagonals or if the diagonals are very short nearly resembling isolated recurrence points, the system can be judged uncorrelated or poorly correlated, stochastic and %DET is close to 0%, while deterministic system yields relatively longer diagonals and %DET is quite away from 0%.
where hist(l) is the histogram of the diagonal lines with lengths l.
ENTR gives the Shanon entropy of the probability distribution of the trace of diagonal lines with exact length l. Mathematically ENTR can be defined as where p(l) is the probability to have diagonals of exact length l.
The ENTR is a measurement of how much data are required to recover the process. A low ENTR means that less information is required to define the process and it has less complexity, whereas a high ENTR suggests that more information is necessary and that the system is more complex. In Fig. 2, it is clearly visible that the frequency components which are present in the signal do not exist continuously along the entire time axis which confirms that K p -index fluctuation is nonstationary in nature. Here, it can be noticed that the frequency of K p -index signal is continuously varying for every time instant, while the colour bar represents the energy strength of every frequency components.

Results and discussion
The detection of nonlinear behaviour of the Planetary K-index using DVV analysis has been illustrated in Fig. 3. A total of 99 IAAFT-based surrogates have been generated for this purpose. Figure 3a represents that the DVV plots of surrogated and the original signal do not precisely match but deviate from each other which confirms the nonlinearity of the Planetary K-index signal. Figure 3a also reveals the least value of the target variance, * 2 min = 0.3612 which is substantially less than that of the surrogate. This signifies that K p variation is not much affected by noise and hence has a deterministic behaviour.
The DVV scatter diagram as shown in Fig. 3b clearly depicts the variation of the target variances mapping locus of the surrogate and original signals from the bisector line. This deviation validates the nonlinear property of the K pindex profile.
A nonlinear system may exhibit chaotic behaviour. To trace the existence of chaos within the time series 0-1 Test has been performed whose outcomes are demonstrated in Fig. 4. The computed D c (n) and D * c (n) for the Planetary K-index time series rise continuously as n increases in Fig. 4a, c, which indicates that the source of the signal is affected by the chaos. From Fig. 4b, Brownian motion is observed in the complex plane of the Fourier spectrum X n which indicates that the signal's possible underlying geometrical structure is random. The variation of the asymptotic growth rate K c is displayed in Fig. 4d. The final binary output for the Planetary K-index signal is computed to be 0.9745. This nearly unit value of the 0-1 test output corroborates about the chaotic behaviour of the signal.
To identify whether chaos is deterministic or stochastic in nature, correlation dimension analysis has been done and is illustrated in Fig. 5. The correlation dimension D(m) saturates at a plateau consequently as seen in the graph. It can be inferred that the chaos of the geomagnetic storm is a deterministic phenomenon, i.e. less disturbed by external noises. Besides the nature of the curve of the D(m) , the magnitude of the D(m) is also an important parameter for evaluating the nature of the system under scrutiny. If the value of the correlation dimension tends towards infinity, more specifically if it is 5 or above, generally the signal generating process is considered as a stochastic one, i.e. heavily administered by noises [52]. As the D(m) attains the plateau at the value of 0.7348, the correlation dimension of the geomagnetic storm is 0.7348 which is a very low value compared to 5. This low value re-establishes the claim of deterministic character of the geomagnetic storm. Moreover, if the correlation dimension is found to have finite integer value, the system dynamics exhibit non-chaotic and strongly periodic deterministic quality. But its fractional value with smaller magnitude implies that the process can be considered as a low dimensional deterministic chaotic one [52][53][54]. So, the fractional correlation dimension of 0.7348 gives the notion that variation of the geomagnetic storm is absolutely a chaotic process having low-dimensional determinism. If the process is regular and deterministic, by Takens embedding theorem [55], the correlation dimension quantifies minimum the number of equations required to express the process's dynamics which is known as degrees of freedom (dof), whereas for chaotic and deterministic process, the integer value just above the fractional correlation dimension gives Fig. 6 a F(q) versus q, b h(q) versus q, c (q) versus q, d f ( ) versus But it must be kept in mind that the existence of a plateau and a finite fractional correlation dimension is not enough to give a conclusive inference about the process, but it gives a probable inclination towards the inferences made about the process. Figures 6 and 7 represent the output of MFDMA analysis based on the backward moving average. In Fig. 6, the reliance of F(q) , h(q) , (q) on q as well as nature of the singularity spectrum f ( ) is noticed. The scale parameter n is set to a range of 10 to L/10 [34] and the exponent q is selected with steps of 0.5 from − 20 to + 20. In Fig. 6a, b, the significant variation of the F(q) and nonlinear dependence of h(q) function with q unveil multifractal behaviour of the signal under scrutiny. The computed value of h(2) = 0.93891 ± 0.0059 suggests the existence of longrange positive correlation within data series. The existence of several structures of different scales is also suggested by the nonlinear shape of the (q) curve. The multi-fractal characteristic of Planetary K-index signal can be explored more thoroughly from the singularity spectrum f ( ) in Fig. 6d. f ( ) gives an estimate about the degree to which the signal is occupied by singularities of varying strengths. The parameters min and max are lowest and highest values of the Hölder exponent for which f ( ) = 0 . The width of the spectrum Δ max − min is often used to measure the strength of multifractality which indicates the 'affluence' of the signal structures. The computed values of max , min , Δ and 0 are tabulated in Table 2.
The multifractal behaviour generally observed in time series is due to (1) flat probability density function (pdf ) of the data series and/or (2) various temporal structures (nonlinearity and long-range correlations) for different fluctuations [56]. Fractal analysis of surrogated and shuffled data produced from the original data is used to evaluate the true cause of multifractality. If upon shuffling the series all its long-range correlations are destroyed making h(q) = 0.5 for all q, but pdf remain unchanged, the multifractality is due to reason (2), i.e. temporal structure. The signal is not influenced by the shuffling process due to the multifractality existence of the signal due to reason (1), i.e. fatness of the pdf. In this case, surrogation (amplitudeadjusted Fourier transform) of the signal will change its pdf to Gaussian distribution, while correlation is not disturbed. The shuffled series' h(q) on q dependence would be almost equivalent with that of the natural signal, while h(q) of surrogate may not be relying on q. While both factors contribute to the signal's multifractality, h for surrogated and shuffled signals would be dependent on q, and the shuffled sequence will have weaker multifractality than the initial one. Figure 7 represents the comparative fractal analysis for original, shuffled and surrogated of the Planetary K-index data series achieved using backward MFDMA. The dependence of h(q) on q for the surrogated and shuffled data, as seen in Fig. 7b, claims that the multifractality present in the Planetary K-index signal is attributed to wideness of the pdf as well as temporal structure. The magnitudes of h(2) for the three signals as shown in Table 2 indicate the existence of long-range association which further computed applying MFDMA algorithm for all three signals nearly are equal to or even greater than 1 reveals that the signal under investigation is nonstationary and that this property is preserved even after randomisation. The skewness ( r = max − 0 0 − min ) is also determined for K p -index signal which is 1.1359. There is an inverse relationship between the values of and the strength of the singularity whereas for a rough signal having high magnitude singularities are characterised by low value [46]. For right skewed profile:r > 1 and it signifies the singularities of lower strength are predominant in the signal. For left skewed: r < 1 and it advocates the presence of higher strength singularities. If the singularities with high and low strengths are equally distributed in the signal, its singularity spectrum will be symmetric for which r = 1 . However, r = 1.1359 for K p -index it can be assumed that the fluctuation of the geomagnetic storm is not much interrupted by singularities of higher strength but exhibit as a smooth signal having low singularity strength [57]. It is visualised from the recurrence plot (RP) in Fig. 8 that the patterns which are noticed match with that mention in point number 2, 3 and 6 in Table 1. So, it can be inferred that the Planetary K-index signal is absolutely non-stationary which may be attributed to the existence of some trends or drift that are either unusual or irregular. The presence of diagonal lines parallel to LOI and isolated single points indicate that the signal has deterministic chaotic feature.
The low value of the %REC for the time series as tabulated in Table 3 suggests that recurring points are rare and that denotes the absence of cyclicities in the process. The considerable magnitude of %DET validates the claim that the Planetary K-index signal is deterministic. The low value of ENTR is obtained for the signal under investigation which reveal the chaotic feature of the signal. The shorter LMAX signifies that the rapid divergence of the trajectory segments of the signal. As a whole, the RQA states that Planetary K-index signal has deterministic chaotic  behaviour with the fast diverging trajectory in the phase space. Also, the temporal evolution graph of four parameters (REC, LMAX, %DET and ENTR) for the quantification of recurrences are represented in Fig. 9

Conclusion
Non-stationary behaviour of the K p -index signifies that the frequency (or periodicity) of occurrence of the geomagnetic storm of nearly the same strength is not uniform all through time. The time intervals between the geomagnetic storms of nearly equal strengths are varying with time. The geomagnetic storm may be recurrent or non-recurrent type. The first one is caused due to the high-speed solar wind from co-rotating interaction region (CIR) of the Sun, and it has a uniform periodicity of nearly 27 days. This periodicity corresponds to the rotational period of the Sun. The non-recurrent type storms, on the other hand, are generally caused due to the high speed Coronal Mass Ejection (CME) whose periodicity of occurrence is not uniform. As the K p -index is found to be nonstationary, it can be concluded that the non-recurrent phenomenon like CME is more dominating over the CIR solar wind to cause the fluctuation of the geomagnetic storm.
Since the K p variation exhibits deterministic characteristic, it can be inferred that the geomagnetic storm is not much affected by the causes which may introduce noise or randomness in its fluctuation. The intensity of the storm due to the CIR or CME is so strong that the effect of the other causes on the storm fluctuation seems noisy and negligible. This non-randomness feature of the K p -index fluctuation indicates a very well-defined "cause-effect" relationship. So, if the initial state of the geomagnetic system is known, the K p -index (effect) in its future state can be precisely predicted for a certain change of the inputs like CIR or CME (cause) and hence helps in geomagnetic forecasting.
The nonlinearity of the K p -index establishes that this 'cause-effect' association is not governed by some simple, proportional polynomial equations of degree one but is ruled by some set of polynomials with degrees more than one. From the Takens' theorem (1980), we can determine the actual dimension (d) of the any process using the relation d ≤ (m − 1)∕2 . Since the computed value of m for our signal is 13, the actual dimension of signal generation process is d ≤ 6 . It suggests that a maximum of six nonlinear polynomial differential equations will be used to model process dynamics. Since a nonlinear system may have multiple attractors, there remains the chance that the dynamics of a nonlinear system may depend on its initial state. The direction of evolution of process parameters in state space is represented by the trajectories. The way the trajectories change with the variation of the system parameters determines the robustness of the system against small perturbation in the system parameters. If the trajectory is found to be irregular but remain confined in region of state space where there is no stable limit cycle or fixed points, the system is said to lie on a strange attractor. If any two arbitrary initial points which are close to each other diverge from each randomly after any number of iterations, and again come closer to each other after any number of iterations making the trajectory confined within a certain region of the state space generates the strange attractor. This repetition of rapid divergence and subsequent closeness of any two points within a confined space can be interpreted that a system with strange attractor is locally unstable but globally stable. The phenomenon of having strange attractor | https://doi.org/10.1007/s42452-021-04622-4 by a system is better known as chaos. The proof of the presence of chaos in the K p -index, as established by 0-1 test, indicates that there is strange attractor in its phase space. It means that the variation of the geomagnetic storm is not much robust against any small perturbation at locally but as a whole it remains stable.
The chaos which is there in the variation of the geomagnetic storm is found to be deterministic with low dimension. The presence of deterministic chaos reconfirms that the system is "sensitive dependence on initial conditions". The aperiodic variation of the trajectory that is generated from an "autonomous" low-dimensional nonlinear dynamical system is generally referred as deterministic chaos. The term "autonomous" signifies the system which has no inputs neither deterministic nor noisy.
The presence of this chaos is due to the fluctuations in the primary states which can cause significant deviation in future. The system is neither a complete chaos nor complete deterministic but it is in between these two states, i.e. deterministic chaos, the state where the system acquires fractal structure. As the system is not completely deterministic, the system must be energy dissipative and hence the system generating geomagnetic storm is a nonlinear dissipative dynamical system (NDDS). Again the absence of complete chaos advocates that there are nonlinear interactions between different temporal scales of the geomagnetic storm time series. These interactions enable the geomagnetic storm generating system to build and rebuild medium to large substructures from short-scale chaotic fragments. This particular state between chaos and determinacy is characterised by the presence of self-similar substructures called fractals. The presence of multifractal indicates that the variation in the strength of the geomagnetic storm characterised by many clusters of self-similar structures. The physical reason of varying slope of (q) is the presence nonlinear interaction between these structures or fractals. The attributes like multifractality, long-range positive correlation present in the Planetary K-index signal is because of both the wideness of its probability density function and the temporal structure, whereas positive persistence implies that the nonlinear interaction between the fractals is strong and their evolution moves in tandem. The superiority of scaling properties of small fluctuations over large fluctuations of K p -index signifies that the small structured fractals are playing more important role to determine the evolution pattern of the geomagnetic storm dynamics.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.