A comparative study on chaoticity of equatorial/low latitude ionosphere over Indian subcontinent during geomagnetically quiet and disturbed periods

In the present study, the latitudinal aspect of chaotic behaviour of ionosphere during quiet and storm periods are analyzed and compared by using GPS TEC time series measured at equatorial trough, crest and outside crest stations over Indian subcontinent, by employing the chaotic quantifiers like Lyapunov exponent (LE), correlation dimension (CD), entropy and nonlinear prediction error (NPE). It is observed that the values of LE are low for storm periods compared to those of quiet periods for all the stations considered here. The lowest value of LE is observed at the trough station, Agatti (2.38 ◦ N, Geomagnetically), and highest at crest station, Mumbai (10.09 N, Geomagnetically) for both quiet and storm periods. The values of correlation dimension computed for TEC time series are in the range 2.23–2.74 for quiet period, which indicate that equatorial ionosphere may be described with three variables during quiet period. But the crest station Mumbai shows a higher value of CD (3.373) during storm time, which asserts that four variables are necessary to describe the system during storm period. The values of non linear prediction error (NPE) are lower for Agatti (2.38 N, Geomagnetically) and Jodhpur (18.3 ◦ N, Geomagnetically), during storm period, compared to those of quiet period, mainly because of the predominance of non linear aspects during storm periods The surrogate data test is carried out and on the basis of the significance of difference of the original data and surrogates for various aspects, the surrogate data test rejects the null hypothesis that the time series of TEC during storm and quiet times represent a linear stochastic process. It is also observed that using state space model, detrended TEC can be predicted, Correspondence to: K. Unnikrishnan (kaleekkalunni@gmail.com) which reasonably reproduces the observed data. Based on the values of the above quantifiers, the features of chaotic behaviour of equatorial trough crest and outside the crest r gions of ionosphere during geomagnetically quiet and disturbed periods are briefly discussed.


Introduction
In the real world system, pure determinism is rather unlikely to be realised, since all systems somehow interact with their surroundings and deterministic picture may be regarded only as a limiting case of a more general framework involving fluctuations in the environment and in the system itself (Hegger et al., 1999). In recent years, nonlinear time series methods were employed to study the magnetospheric chaos, and these studies strongly support the evidence of magnetospheric chaos (Vassiliadis et al., 1990;Shan et al., 1991;Chen and Sharma, 2006).
During the last decade, several interesting studies and modeling efforts on magnetospheric dynamics based on theories of nonlinearity, chaos and stochasticity were successfully carried out. Bhattacharyya (1990) had studied the chaotic behaviour of ionospheric density fluctuations, using amplitude and phase scintillation data, and found the existence of low-dimensional chaos. Also, Wernik and Yeh (1994) have studied chaotic behaviour of ionospheric turbulence using scintillation data and numerical modeling of scintillation at high latitude. Also, Kumar et al. (2004) reported the evidence of low dimensional chaos in a set of TEC data, obtained by Faraday rotation technique, measured at a high-latitude station, Goose Bay (47 • N, 286 • E). Recently, Unnikrishnan et al. (2006) have analyzed the deterministic chaotic behaviour of GPS TEC fluctuations at mid-latitude. However, comprehensive studies on non linear 766 K. Unnikrishnan: Chaoticity of equatorial/low latitude ionosphere aspects and chaotic behaviour of equatorial/low latitude ionosphere during quiet and disturbed periods, especially using GPS TEC time series are not conducted so far.
Most of the Indian region encompasses the equatorial and low-latitude ionospheres. The morphology of the equatorial ionosphere is quite different from that of other latitudes because the magnetic field (B) at the equatorial region is nearly parallel to the Earth's surface. During daytime, the E-region dynamo electric field (E) is eastward. This field in the E-region and at off-equatorial latitudes maps along the magnetic field to F-region altitudes above the magnetic equator, resulting in E × B drift, which transports F-region plasma upward over the magnetic equator. The uplifted plasma over the equator then moves along magnetic field lines in response to gravity, diffusion, and pressure-gradient forces. The physical mechanism, which is traditionally involved in order to explain the diurnal variation of the equatorial anomaly, is the so-called fountain effect and this fountain is controlled by the E × B drift (Abdu, 1997). As a result, the equatorial ionization anomaly (EIA) is formed with reduced F-region ionization density at the magnetic equator and increased ionization at the two anomaly crests around ±15 • in magnetic latitude to the north and south of the magnetic equator (Rama Rao et al., 2006).
In the present study, the latitudinal aspect of chaotic behaviour of ionosphere during quiet and storm periods are analyzed and compared by using GPS TEC time series measured at equatorial trough, crest and outside crest stations over Indian subcontinent. For this the chaotic quantifiers like Lyapunov exponent, Correlation dimension, entropy and nonlinear prediction error are estimated for three stations.

Data and methodology
In the Indian region the augmenting of GPS is planned through a regional Satellite Based Augmentation System (SBAS), called GPS Aided Geo Augmented Navigation (GAGAN) from the Indian Space Research Organisation (ISRO) and Airport Authority of India (AAI). The Slant Total Electron Content (STEC) is the measure of the total number of free electrons in a column of the unit cross section along the path of the electromagnetic wave between the satellite and the receiver. The total number of free electrons is proportional to the ionospheric differential delay between L1 (1575.42 MHz) and L2 (1227.60 MHz) signals. The Vertical TEC is obtained by taking the projection from the slant to vertical using the thin shell model assuming a height of 350 km, following the technique given by (Klobuchar, 1986): where Re = 6378 km, h max = 350 km, = elevation angle at the ground station. Rama Rao et al. (2006) observed that the IPP (Ionospheric Pierce Point) altitude of 350 km is valid In the present study three Indian stations, namely Agatti, Mumbai and Jodhpur are considered, which represent the equatorial trough, crest and outside the crest regions of EIA feature (Table 1). GAGAN-TEC data for a period, 11-14 June 2005, with 1-min resolution, measured at 3 stations given in the Table 1 are used as the time series representing storm period for this study. During this period of study, a Storm Sudden Commencement (SSC) was observed on 12 June 2010 at 12:42 IST (07:42 UT), with Dst max = -106 nT, which is a major geomagnetic storm event. To compare with the storm time chaotic behaviour, GAGAN-TEC data for a period 7-10 June 2005 is chosen as the quiet period.

Time series and phase space
A time series, S n is the sequence of scalar measurements of some quantity, which depends on the current state of a system, taken at multiples of a fixed sampling time ( t), and Fig. 1a represents the typical time series of GAGAN-TEC measured at Agatti. Figure 1a shows that the dominant dynamics could be the diurnal variations, in addition to the fluctuations due to the internal dynamics. Since, we are interested in the nonlinear TEC variabilities, it is essential to minimize the influence of diurnal variations observed in the data. To reduce the influence of the diurnal variations of the data, we have carried out "epoch analysis", which is as follows.
Here, the number of data points in one day is 1440. Let a i , i = 1,2,3...1440 denotes the average of the values t i , t i+1440 , t i+2880 , etc., where t i , i equals 1, 2, 3, etc., denotes the observed time series of TEC. Then the diurnal variation reduced time series is given by where i = 1,2,3..., j = mod(i, 1440), if mod(j , 1440) = 0, and j = 1440 if mod(j , 1440) = 0. By using this method on the measured TEC, we get the detrended time series of TEC ( Fig. 1b), and the dynamical results of the detrended TEC are explained further (Kumar et al., 2004;Unnikrishnan et al., 2006). The detrending filter will transfer the nonlinear TEC variabilities/fluctuations, by minimizing the influence of diurnal variations, as much as possible, and Fig. 1c depicts the variations of a i . The detrended time series of TEC measured at various stations of equatorial/low latitude region, considered in the present study, are subjected for analysis to obtain mutual information, false nearest neighbours, delay coordinates, correlation dimension, Lyapunov exponent, entropy, and non linear prediction error.
To reveal the multidimensional aspects of the system, phase space reconstruction of the time series is required. For this, the proper choice of embedding dimension (m) and delay time (τ ) are essential (Fraser and Swinney, 1986;Kennel et al., 1992). If the time delayed mutual information shows a marked minimum, that value can be considered as a reasonable time delay. Likewise, the minimal embedding dimension, which corresponds to the minimum number of false nearest neighbours can be found out, and is treated as the optimum value of m (Unnikrishnan et al., 2006). With the help of time delayed mutual information ( Fig. 1d) and the false nearest neighbour method (Fig. 1e), it is observed that any delay ≥ 18 and any dimension ≥ 6 are the suitable choices for τ and m, respectively. By repeating the similar analysis for all the time series considered in the present study, these choices for τ and m are found to be true. Hence, the choice of τ as 20 and m as 8 for further analysis of this study are reasonable. Based on embedding theorem, a multidimensional state space can be reconstructed as follows: where Y n are the vectors. Figure 2a represents the reconstructed phase space for the time series of TEC observed at Agatti (Fig. 2a), using time delay embedding for τ = 20 and m = 8.
To check the stationarity, probability density of the time series of TEC was calculated, by dividing the range of values of TEC into short intervals, and then by counting the values of the series fall in each interval (Kumar et al., 2004;Unnikrishnan et al., 2006;Unnikrishnan and Sudha Ravindran, 2010). The probability density for the original TEC time series (continuous line) and those of the first half of the series (dashed line) were calculated and shown in Fig. 2b. From this Figure, it is evident that the trends of curves for the probability density of the original TEC time series, and that of the first half of the TEC time series are very similar, which supports stationarity of the data.
It is known that Lyapunov exponent is a measure of the rate of attraction to or repulsion from a fixed point in the state space. One of the most prominent evidences of chaotic behaviour of a dissipative deterministic system is the existence of positive Lyapunov exponent. A positive Lyapunov exponent indicates divergence of trajectories in one direction, or alternatively, expansion of an initial volume in this direction, and a negative Lyapunov exponent indicates convergence of trajectories or contraction of volume along another direction.
The average rate of divergence can be estimated by the first Lyapunov exponent λ 1 (Wolf et al., 1985).
The Lyapunov exponent for the TEC time series measured at Agatti, was computed and its evolution, as the state space trajectory is scanned with τ = 20, m = 8, is shown in Fig. 3a.
The value of Lyapunov exponent was calculated for various time delays by keeping embedding dimension a constant. A plot of Lyapunov exponent versus time delay, by keeping embedding dimension a constant (Fig. 3b) for TEC time series measured at Agatti shows that, its values initially increase up to a time delay 15, and thereafter remain almost constant.
The spectrum of the Lyapunov exponents ( Fig. 3c) can be estimated from a time series by flowing the evolution of small perturbations of the reconstructed orbit, making use of a linearized approximation. The evolution of the displacement vector between the neighboring points X(i) and X(i) + w(i) in the reconstructed phase space is given by the equation where DF denotes the derivative matrix of F. A local approximation of the matrix DF can be found by minimizing the following expression where k is the number of the neighbours to X(i) regarding k different perturbations w j , j = 1,...,k, which are used to estimate A i ≡ DF at the point X(i). The Lyapunov spectrum is found by repeating this process for all N reconstructed       points X(i), i = 1,...,N, that is where { e i j } new set of orthogonal vectors produced by orthonormalisation of the vectors at time i in order to retain the local orthogonal spanning of the state space (Athanasiu et al., 2003). Another significant quantifier used in chaotic analysis is the correlation dimension, D which is defined as: lnr (8) where C(r) is the correlation sum for radius r. Also, it can be seen that C(r) ∼ r d for r → 0. The dependence of correlation sum of the time series of length N, on the embedding dimension m of the reconstructed phase space is given by where is the Heaviside step function, (a) = 0 if a ≤ 0 and (a) = 1 for a > 0.
The parameter, mτ , (i.e. the product of m and τ ) is the time span represented by an embedding vector, which is nearly equal to w, called Theiler window. We have calculated the Theiler window as w = 160, which is approximately equal to mτ . Figure 4a shows the curves of correlation sum C(r) verses r in log scale, where the scaling region is prominent. We have estimated the correlation dimension by varying Theiler window, w from 170 to 300 in small steps and there is no appreciable difference. Hence we have selected ω as 200 for estimating the correlation dimension. Generally, the best approximation to the actual dimension of the attractor is the value of correlation dimension (Grassberger and Procaccia, 1983), and is given by y intercept value of the horizontal line, drawn along the plateau of the local slopes of correlation sum for various embedding dimensions. In experimental data, getting a clear plateau is very difficult in a deterministic case. At the same time, the local slopes converge to some value instead of getting a plateau for all embedding dimensions, m greater than the actual dimension. Local slopes of the logarithms of correlation sum for normalised and detrended TEC (Fig. 4b) exhibits a clear plateau, indicating a fractal dimension of the attractor, whose value is obtained as, D = 2.74 ± 0.02.
Recurrence Plot (RP) first described by Eckmann, et al. (1987), and further modified by Zbilut and Webber (1992) is a powerful analytical tool for the study of nonlinear dynamical systems. With RP, one can graphically detect hidden patterns and structural changes in data or see similarities in patterns across the time series under study. If the time series is truly random and has no structure, the distribution of colors over the RP will be uniform, and has no identifiable patterns. On the other hand, if there is some determinism in the signal generator, it can be detected by some characteristic, distinct distribution of colors. The Shannon entropy of the probability distribution of the diagonal line lengths p(l) is given as: where l is the lowest number of upward diagonal recurrent points required to define a deterministic line. We have performed similar nonlinear time series analysis for all the TEC time series measured at different stations considered in the present study, and it is seen that the choice of τ as 20 and m as 8 are appropriate for all of them. Figure 4c and d presents the RPs of the typical TEC time series for quiet and storm periods, respectively, which exhibit different characteristic patterns. By preparing the Recurrence Plots (RPs) for all the time series studied here, their entropies were estimated using Eq. (10).

Implementation of surrogate data test
The method of surrogate data (Theiler et al., 1992) obtained by using Monte Carlo techniques enable us to assess the significance of the results, we have derived in this analysis. It has been accepted that, the method of surrogate data test could be a successful tool for the identification for nonlinear deterministic structure in an experimental data (Pavlos et al., 1999). According to this method, typical characteristics of the data under study are compared with those of stochastic signals (surrogates), which have the same autocorrelation function and the power spectrum of the original time series. If the behaviour of the original data and the surrogate data are significantly different, in such characteristics, it may be safely concluded that, a stationary linear Gaussian stochastic model cannot describe the process under study. In the present study, we used the improved algorithm of Schreiber and Schmitz (1996) to generate 40 surrogate data. In this case, the generated "nonlinear" surrogate data mimic the original time series in terms of the autocorrelation function and the probability density function. Then, we compared the geometrical and dynamical characteristics of the measured data and those of surrogates. The significance of difference of the characteristics can be defined as where α and σ are the mean and standard deviation values of the characteristic of surrogates and α original is the value of the same characteristic for the original data.
In the present study, the significance of difference of the original data and surrogates for mutual information, fraction of false nearest neighbors, and Lyapunov exponent (Fig. 5ah), for quiet period are greater that 2, which enables to reject the null hypothesis, that the time series of TEC represent a linear stochastic process, with confidence greater than 95% (Pavlos et al., 1999). The nonlinear prediction error for surrogates of quiet time TEC for Agatti is, 0.3823 ± 0.008052, while that for detrended TEC is, 0.2485. In this way, the geometrical and dynamical characteristics of all the time series considered in the present study and those of surrogates are compared and the significance of difference are found to be greater than 2. This implies that determinism is associated with the GPS-TEC time series along with its chaoticity. Values of correlation dimension (CD) computed for TEC time series during storm period show appreciable variabilites over the latitude sector of the present study (Fig. 6b). Storm time values of this parameter are lowest at Agatti (2.51) and Jodhpur (2.58), which are equatorial trough and outside the crest stations, respectively. However, correlation dimension value of the crest station, Mumbai is highest compared to other two stations, during storm period. It is interesting to note that, the value of CD observed at Mumbai during quiet time is appreciably less than that during storm period, where as the values of this parameter for quiet and storm periods are almost same for Agatti and Jodhpur.

Latitudinal variations of chaotic quantifiers
Latitudinal dependence of entropy is not so prominent and its value observed at Mumbai during storm period is lesser than that of quiet period (Fig. 6c). However, Agatti and Jodhpur show no such trends, and having almost same values for quiet and storm periods.
Non linear prediction error (NPE) has strong latitudinal dependence during quiet period, which is absent for storm period (Fig. 6d). These values during quiet period observed at Agatti (trough station) and Jodhpur (out side the crest station) are comparatively higher than those of storm period. The quiet time value of NPE observed at Mumbai (crest station) is lowest, which is almost identical with the corresponding storm time value.

Prediction using phase space model
In a low dimensional dynamical system, the trajectories are well defined, although over long periods, they may be chaotic. This makes short term predictions of the dynamics in such systems, reasonable and practical. In the case of a system specified by the observational data of one of its variables, the reconstructed phase space can be used to make predictions (Farmer and Sidorowich, 1987;Sharma, 1997).
The basic principle is illustrated in Fig. 7a. The filled circles within the neighbourhood (represented by the bigger circle of radius r) denote the locations of the trajectories and those outside denote the known locations of these points at the next time step. Knowing how the neighbouring trajectories evolve, the location of the current state X(t) at time X(t + τ ) can be predicted (Farmer and Sidorowich, 1987;Sharma, 1997).
Using this technique, the Detrended TEC(t) value at the next time interval Detrended TEC(t + 1) is expressed as where the functional F of the state vector X(t) is obtained from the dynamical features of the reconstructed phase space. The local value of F is obtained by a Taylor expansion around X(t) and the coefficients computed by fitting procedure using the evolution of the nearest neighbours of X(t). When the Taylor expansion is limited to the linear term, this yields a local linear technique. Using this technique, detrended TEC is predicted, which reasonably reproduces the observed data (Fig. 7b).

Discussion
It is known that E × B drift brings plasma from low to high altitudes during daytime and from high to low altitudes (downward) during night. In the morning and evening, electron density and temperature peaks are observed (Fukao et al., 1991). The local time variations of the neutral wind, electric fields, production-recombination rates produce inherent irregularities leading to chaotic variations, even at 772 K. Unnikrishnan: Chaoticity of equatorial/low latitude ionosphere quiet time ionosphere. This shows that, the response of ionosphere to the above variabilities is irregular, showing a disorder with determinism. However, during storms, this behaviour is modulated and influenced by storm time energy drivers like solar wind and subsequently disturbed thermospheric wind electric field systems. Anomalous night time enhancement in TEC producing TEC fluctuations during local night times (Balan and Rao, 1987), and collisional Rayleigh-Taylor instability (Huba et al., 1996), playing a crucial role in the development of equatorial spread F, are another points to account for ionospheric irregularities. The evolution of the magnetospheric system will be the result of the combined effects of local couplings of the magnetic and plasma structures, through the nonlinearities 774 K. Unnikrishnan: Chaoticity of equatorial/low latitude ionosphere   The magnetosphere thus exhibits the signatures of both self-organization and self organized criticality. This point of view also supports the recent results of  that the substorm activity resembles the non equilibrium (first and/or second order) phase transitions. The phase transition analogy provides a conceptual frame work for combining the global and multi scale aspects of substorm dynamics in a unified model. The global behaviour of magnetosphere has clear dynamical behaviour and can be predicted using dynamical models. The multiscale features on the other hand has a power law behaviour and can not be modeled in the dynamical sense and can be described/modeled in terms of statistical properties (Sharma et al., 2005). Like magnetosphere, ionosphere may also exhibit global coherent dynamics and multi-scale features (first and/or second order) at times of geomagnetic disturbances.
Both Borovsky et al. (1993) and Huang et al. (2003) explained that, the combination of the storm and substorm caused some unique and well-correlated phenomena in the auroral-subauroral ionosphere. Infact, magnetosphere has an intrinsic oscillation period of approximately 2-3 h during disturbed times and can maintain such periodic oscillations even after substorms have stopped. Particle precipitation from ring current or plasmasphere during the magnetospheric oscillations causes simultaneous increases in electron density and temperature in the equatoriallow/mid latitude ionosphere (Huang et al., 2003). Hence, the superposition of a large number of active degrees of freedom can produce extremely complicated signals during a magnetic storm period, which in turn modify the stability/instability conditions of magnetosphere. Hence during geomagnetic disturbances, the above effects via substorms may increase/modulate the periodicity in TEC fluctuations and thereby decrease/modulate the chaoticity, which is evident by the decrease/modification in values of Lyapunov exponent during storms, compared to quiet times.
The station, Agatti is well within the equatorial trough, and the TEC fluctuations at this region are expected to be lower, compared to the stations at/near the boundary of the anomaly crest. Infact, chaos arises from the exponential growth of infinitesimal perturbations in a system, and Lyapunov exponent is the estimate of its rate. This may be a reason for lower values of Lyapunov exponent observed at Agatti, which indicates small exponential divergence of the initially closed trajectories in the phase space of TEC time series. But station like Mumbai, lies near/at the edge of anomaly crest, which may suffer more fluctuations, because of the spatial movement of EIA, based on local time and the influence of EEJ/counter electrojet. This could be a reason for higher value of LE observed at Mumbai. Unnikrishnan et al. (2006) analyzed the deterministic chaotic behaviour of GPS TEC fluctuations at mid-latitude, and their study emphasis that the influence of an external stochastic driver (solar wind) could alter the inherent dynamics of a system (ionosphere), if the coupling is powerful. This could be a possible reason for lower values of Lyapunov exponent during storms, with respect to the corresponding quiet time values (Unnikrishnan, 2006) observed at mid-latitude. In the present study also, it is observed that the values of LE are low for storm periods compared to those of quiet periods for all equatorial/lowlatitude stations considered here.
The values of correlation dimension computed for TEC time series are in the range 2.23-2.74 for quiet period, which indicate that equatorial ionosphere may be described with three variables during quiet period. But the crest station Mumbai shows a higher value of CD (3.373) during storm time, which asserts that four variables are necessary to describe the system during storm period. It implies that storm time crest station Mumbai needs one more variable to describe/model TEC time series compared to trough and outside the crest stations during quiet and storm periods. The values of non linear prediction error (NPE) are lower for Agatti and Jodhpur during storm period, compared to those of quiet period, mainly because of the predominance of non linear aspects during storm periods. It is understood from recent studies that the magnetospheric dynamics is neither clearly low dimensional nor completely random, but exhibit combinations of these two aspects, which yields an improved and effective concept for forecasting space weather Ukhorskiy et al., 2004). Like magnetosphere, ionosphere is also a complex system, where the basic fluctuations are of internal origin, though the bursts can be triggered by an external perturbation (solar wind).

Summary
In the present study, the chaotic behaviour of ionosphere during quiet and storm periods over equatorial trough, crest and outside the crest stations of Indian sub continent are analysed and compared, by using GPS-TEC time series.
The station, Agatti is well within the equatorial trough, and the TEC fluctuations at this region are expected to be lower, compared to the stations at/near the boundary of the anomaly crest. This may be a reason for lower values of Lyapunov exponent observed at Agatti, which indicates small exponential divergence of the initially closed trajectories in the phase space of TEC time series. But station like Mumbai, lies near/at the edge of anomaly crest, which may suffer more fluctuations, because of the spatial movement of EIA, based on local time and the influence of EEJ/counter electrojet. This could be a reason for higher value of LE observed at Mumbai. It is to be noted that, during geomagnetic disturbances, the effects of substorms may increase/modulate the periodicity in TEC fluctuations and thereby decrease/modulate the chaoticity, which is evident by the decrease in values of Lyapunov exponent during storms, compared to quiet times.
The values of correlation dimension computed for TEC time series are in the range 2.23-2.74 for quiet period, which indicate that equatorial ionosphere may be described with three variables during quiet period. But the crest station Mumbai shows a higher value of CD (3.373) during storm time, which implies that four variables may be required to describe the system during storm period. The values of non linear prediction error (NPE) are lower for Agatti and Jodhpur during storm period, compared to those of quiet period, mainly because of the predominance of non linear behaviour of ionosphere over these stations, during storm periods.
To further confirm the deterministic nature of the time series, the surrogate data test was carried out on time series considered in the present study. It is observed that the significance of differences of the original data and surrogates for mutual information, fraction of false nearest neighbors, Lyapunov exponent and local slopes of the correlation sums are greater than 2, which enable to reject the null hypothesis, that the time series of TEC represent a linear stochastic process, with confidence greater than 95%. But still then it is not possible to confirm that equatorial/low latitude ionosphere is purely low dimensional under all conditions. A more reasonable and logical concept is that, like magnetosphere, ionosphere is neither clearly low dimensional nor completely random, but exhibit combinations of these two aspects, which yields an improved and effective concept for forecasting space weather.