Annales Geophysicae (2001) 20: 175–183 c ○ European Geophysical Society 2001

Abstract. A detailed nonlinear time series analysis of the hourly data of the geomagnetic horizontal intensity H measured at Kodaikanal (10.2° N; 77.5° E; mag: dip 3.5° N) has been carried out to investigate the dynamical behaviour of the fluctuations of H . The recurrence plots, spatiotemporal entropy and the result of the surrogate data test show the deterministic nature of the fluctuations, rejecting the hypothesis that H belong to the family of linear stochastic signals. The low dimensional character of the dynamics is evident from the estimated value of the correlation dimension and the fraction of false neighbours calculated for various embedding dimensions. The exponential decay of the power spectrum and the positive Lyapunov exponent indicate chaotic behaviour of the underlying dynamics of H . This is also supported by the results of the comparison of the chaotic characteristics of the time series of H with the pseudo-chaotic characteristics of coloured noise time series. We have also shown that the error involved in the short-term prediction of successive values of H , using a simple but robust, zero-order nonlinear prediction method, increases exponentially. It has also been suggested that there exists the possibility of characterizing the geomagnetic fluctuations in terms of the invariants in chaos theory, such as Lyapunov exponents and correlation dimension. The results of the analysis could also have implications in the development of a suitable model for the daily fluctuations of geomagnetic horizontal intensity. Key words. Geomagnetism and paleomagnetism (time variations, diurnal to secular) – History of geophysics (solar-planetary relationships) Magnetospheric physics (storms and substorms)


Introduction
The geomagnetic field pervades the region around the Earth, extending to several times of the radius of the Earth.Solar output, in terms of solar plasma and magnetic field, ejected Correspondence to: G. Renuka (dlcampus@vsnl.com)out into the interplanetary medium contributes greatly to the perturbation in the geomagnetic field.Episodes of extra ordinary fluctuations in Earth's magnetic field were detected as storms in the mid 1800's.The analysis of the fluctuations in the geomagnetic field has many practical applications in magnetic navigation, orientation control, geophysical exploration, etc. (Newitt, 1993;Kerridge, 1993;Gonzalez et al., 1994;Sutcliffe, 2000).The analysis of storm morphology has been undertaken by several authors.Different definitions of geomagnetic storms have been given by Gonzalez et al. (1994).According to the classical definition, a geomagnetic storm occurs when the daily A p index exceeds 29, a minor storm occurs when 30 ≤ A p < 50; a major storm occurs when 50 ≤ A p < 100 and a severe storm occurs when A p ≥ 100 (Lundstedt, 1996).
A continuous recording of any of the components of the geomagnetic field typically exhibits two types of variations: a smooth, regular variation, known as Sq, and the solar quiet day variation, which arises as the magnetic signature of the E-region ionosphere current driven by a dynamo action (Campbell, 1989) and a rapid irregular fluctuation, referred to as a geomagnetic disturbance or storm, the magnitude of which may be such that the regular Sq variation is swamped and thus, not easily discernible.Although the Sq variations are the most regular of all the geomagnetic field variations, tending to repeat itself with a periodicity of 24 h, significant day-to-day differences do occur (Hibberd, 1981;Sutcliffe, 2000).At low and middle latitude stations H , is known to change drastically during a geomagnetic storm.Wright (1962) found that the H field at Ibadan (dip lat. 3 • S) for any hour of the day was lower on international disturbed (ID) days than on international quiet (IQ) days.Bharghava and Subramanyan (1964) found that at Kodaikanal (10.2 • N; 77.5 • E; mag: dip 3.5 • N), there is little variation in the daily range of H and Z on disturbed days when compared to that under quiet conditions.Vestine et al. (1947) and Wright (1962) also found that the daily range of H at low-latitudes remains unchanged on magnetically active days.Bhargava and Yacob (1969) found a systematic decrease in H at low-latitudes, with an increasing A p index during extremely quiet periods and suggested the effect is associated with the interaction of the solar wind on the magnetosphere.One of the objectives of investigation of the dynamical behaviour of the fluctuations of the geomagnetic field is to predict storms.Wang (1996) has applied fractal theory in a quantitative analysis of geomagnetic storms and has estimated the correlation dimension of the attractor for storm data from the Beijing observatory (40.0 • N, 116.2 • E).
Nonlinear dynamic methods have been applied to magnetospheric data in order to study the underlying dynamics (Sharma, 1995;Klimas et al., 1996).Studies using these methods have given results supporting the concept of magnetospheric chaos (Vassiliadis et al., 1990(Vassiliadis et al., , 1992;;Robert et al., 1991;Shan et al., 1991;Pavlos et al., 1992Pavlos et al., , 1994Pavlos et al., , 1999a, b, c), b, c).However, several studies have given evidence against the hypothesis of magnetospheric chaos and indicate the significant role of the stochastic solar wind driver (Prichard and Price, 1993;Price et al., 1994;Takalo andTimonen, 1994, Prichard, 1994).Klimas et al. (1996) have given an excellent review of the studies on these aspects.The criticism about the magnetospheric chaos has been recently addressed in a series of papers by Pavlos et al. (1999a, b, c).They have tested the null hypothesis that the observed AE index signal is generated by a linear stochastic signal, possibly perturbed by a static nonlinear distortion.In the first paper (Pavlos et al., 1999a) of the series, they have used four distinct geometric parameters derived from the slope of the correlation integral as discriminating statistical procedures in order to test the null hypothesis of the nonlinear stochastic surrogate data, which have the same power spectrum and amplitude distribution as the original data.In the second paper (Pavlos et al., 1999b), dynamical characteristics, such as Lyapunov exponents, nonlinear dynamic models and mutual information, were used to test the null hypothesis.The results of these tests suggest the rejection of the null hypothesis that the AE index signal belongs to the family of stochastic signals undergoing a static nonlinear distortion, i.e. the results of these studies strongly support the hypothesis of nonlinearity and chaotic behaviour of the underlying dynamics of the magnetospheric system.In continuation of these studies, they have introduced significant theoretical concepts about the magnetospheric system and its dynamical interaction with the solar wind (Pavlos, 1999c).Based on the comparison of the observational behaviour of the magnetospheric system with the results of the analysis of the different types of stochastic and deterministic input-output systems, they have observed that the hypothesis of low-dimensional chaotic behaviour of the magnetospheric dynamics remains a possible and fruitful concept which must be developed further.Hence, we feel that the tools of nonlinear time series analysis can be used with confidence in order to obtain useful information about the internal deterministic component of a magnetospheric time series.
Our main objective in this work is to carry out a detailed nonlinear time-series analysis of the time series of the measurements of the geomagnetic horizontal field H .The data we used in this analysis represent the geomagnetic horizontal intensity H , measured during the year 1991 at a one hour interval at the Kodaikanal observatory and published by the Indian Institute of Geomagnetism, Bombay.The importance of this data is that H varies drastically during the year and in addition, we noted 33 storms, out of which five were severe.The time series of H is plotted in Fig. 1a and its delay representation in Fig. 1b.The origin of the values of H has been shifted to 39 000 nT.

Nonlinear time series analysis
In a purely deterministic system, the states of all future times are determined once its present state is fixed.Thus, we can study the dynamics of the system by studying the dynamics of the corresponding state space for the study of systems with deterministic properties.However, at times, the only information of the system available is a series of univariate measurements equidistant in time, i.e. a time series.Under these circumstances, one has to construct a new state space in which the mapping from one point of the trajectory to the successive one is unique.This is accomplished using time delay coordinates.Embedding theorems (Sauer et al., 1991) guarantee that for an appropriate delay depending on the data, at most 2d + 1 delay coordinates are enough, when d is the fractal dimension of the attractor.A time series is a sequence of scalar measurements of some quantity which depends on the current state of the system, taken at multiples of a fixed sampling time: where η n is the measurement noise.A delay reconstruction in m dimensions is then found by the vectors y n , given by The time difference in number of samples v between adjacent components of the delay is referred to as the lag or delay time.One of the problems of the nonlinear time-series analysis is to find an optimal embedding dimension m.However, for many practical purposes, the most important embedding parameter is the product mτ of the embedding dimension m and the delay time τ , since mτ is the time span represented by an embedding vector.A precise knowledge of m is only required when we want to exploit determinism with minimal computational effort (Kantz and Schreiber, 1997).Nevertheless, there are several indicators of an optimal embedding dimension m.One such indicator is the correlation dimension D, defined by where C(r) is the correlation sum for radius r, which reveals a scaling profile as C(r) ∼ r d for r → 0. The correlation sum depends on the embedding dimension m of the reconstructed phase space and the length of the time series N as where is the Heaviside step function, (a) = 0 if a ≤ 0 and (a) = 1 for a > 0. The scaling exponent d increases with m and saturates to a final value of D for sufficiently large embedding dimension m o .In most cases, m o may be the smallest integer larger than D (Ding et al., 1993).When the slopes d of the correlation integral for various embedding dimensions reveal a plateau at low values of r and the plateau converges for increasing m, then there is strong evidence for a low-dimensionality of the underlying dynamics of the system.Calculation of the correlation dimension from the time series is easy, compared to the calculation of other dimensions and, in most cases, gives a good approximation to the actual dimension of the attractor.Another technique to estimate an optimal value of m is to look for the closed false neighbors (FNN) in the phase space at a given value of m (Kennel et al., 1992).As the embedding dimension increases, the number of false neighbours decreases.Thus, one can detect the minimal embedding dimension which corresponds to the minimum number of false neighbours.These techniques give the optimal number of independent variables required to reconstruct the underlying attractor of the dynamical system from the time series.Most of the significant characteristics of the original phase space are carried over to the reconstructed phase space for the sufficient large embedding dimension obtained in the above discussed methods.Mathematically speaking, the geometrical characteristics of the original phase space are topologically equivalent to its mirror dynamical flow in the reconstructed phase space.Consequently, it is possible to make both qualitative and quantitative statements about the system based on the time series by using the tools of a nonlinear time-series analysis.Thus, it is relevant to carry out a detailed, nonlinear time-series analysis of the fluctuations of the geomagnetic horizontal intensity H.The results of the analysis are discussed in the following section.

Results and discussion
We started the analysis by estimating an optimal value of the embedding dimension m, using the method of the closest false neighbours.The fraction of the closest false neighbours as a function of the embedding dimension is plotted in Fig. 2. It is evident from the figure that the fraction of false neighbours decreases drastically after the first few values of m.In order to establish the low dimensionality of the attractor, we have also calculated the correlation dimension.
Attractors of dissipative chaotic systems generally have very complex geometry and hence, are called the strange attractors.One of the parameters that characterizes an attractor is its dimension, which can be regarded as a measure of the amount of information necessary to specify the position of a point on the attractor within a given accuracy.The correlation dimension is one such parameter which depends on the spatial correlation of points on the attractor.
The temporal correlation present in a time series can lead to a spurious estimation of the correlation dimension of the Initially, the curves increase quickly and later, the stable diurnal oscillations repeat.
attractor (Theiler, 1991).In order to remove this spurious effect, the temporally correlated points are excluded from the pair counting in the correlation sum, i.e. the sum in Eq. 4 is restricted to pairs of points whose time indices differ by at least ω, called the Theiler window (Theiler, 1991).Hegger et al. (1999) have suggested that the value of ω should be chosen generously.The space-time separation plot introduced by Provenzale et al. (1992) provides a good means of determining a sufficient value of ω (Hegger et al., 1999).In the presence of temporal correlations, the probability that a given pair of points has a distance smaller than r depends on both r and also on the time t between the two points.In a space-time separation plot, the number of pairs is plotted as a function of two variables, the time separation t and the spatial distance r (Hegger et al., 1999).The space-time separation plot of the geomagnetic horizontal field H is given in Fig. 3.The diurnal variability is dominant in this system and is reflected in Fig. 3.The temporal correlation is evident within the first 12 time steps where the lines increase consistently.However, the oscillations saturate after a few cycles.To be safe, we have chosen the minimum correlation time to be 100.However, we have varied the value of the Theiler parameter ω from 100 to 1000 for the calculation of the correlation dimension, but there was no significant difference.The value of the correlation dimension calculated is 3.70 ± 0.04 which corresponds to a good plateau, as shown in Fig. 4.This fractional correlation dimension indicates the existence of a strange, low dimensional attractor.
We then analyzed the deterministic nature of the time series.This is important since a finite time series with a broad band power spectrum may be a realization of a stochastic process governed by an autoregressive moving average model or of a low dimensional deterministic chaotic process (Eckman and Ruelle, 1985).It has also been noted that some geometrical and dynamical characteristics (low correlation dimension and positive Lyapunov exponent, etc.) of the low dimensional chaotic dynamics can also be observed from a particular linear stochastic process.We employed the method of surrogate data and also the Recurrence Plots (Eckman et al., 1987) to distinguish between linear stochastic and deterministic dynamics.
The method of surrogate data has been reported to be a successful tool of choice for the identification of a nonlinear deterministic structure in an experimental data (Mitschke and Dämmig, 1993).It is specifically devised to contrast a data set under study with data that are similar with respect to linear correlations, but otherwise purely random.According to this method, the geometrical and dynamical characteristics of the data under study is compared with those of stochastic signals which have the same Fourier amplitudes and the same distribution of values.If the behaviour of the original data and the surrogate data are significantly different in such characteristics, then it may be safely concluded that the process under study cannot be described by a stationary, linear Gaussian stochastic model.The significance of the difference of the local slopes of the correlation sums can be defined as where µ is the mean value of the local slopes of the correlation sums of different surrogates, and σ is their standard deviation (Mitschke and Dämmig, 1993;Pavlos et al., 1999a).
For our analysis, a group of 40 surrogates of the time series H were generated based on the algorithm of Schreiber and Schmitz (1996).The plot of the local slopes of the correlation sums of the surrogates, along with that of H as a function of ln(r), is given Fig. 5a. Figure 5b gives its mean value and the standard deviation (S.D.).For the region where S > 2, we can reject with a confidence above 95% the null hypothesis that the time series of H is a realization of a linear stochastic process (Pavlos et al., 1999a).As shown in Fig. 5c, the significance of the difference S reaches above 7 for small values of r, and remains above 2 for a larger interval of small r values.
We have also used the recurrence plot of the time series of H to identify the structure in the data.Recurrence plots are useful to graphically detect hidden patterns and structural changes in data or to see similarities in patterns across the time series under study (Eckman et al., 1987).Once the dynamical system is reconstructed by means of delay coordinates, the distances between all pairs of vectors y i and y j are computed and various colour-codes are assigned to different distances.In a two-dimensional recurrence plot, a colour-code at position (i, j ) specifies the distances between the vectors y i and y j .For random signals, a uniform distribution of colours over the entire plane is obtained, whereas for deterministic signals, the recurrence plot contains beautiful structures.The recurrence plot in the grey scale of the time series in this study is given in Fig. 6, which provides evidence of the determinism in the data.
In addition to these techniques, we have also computed the spatiotemporal entropy.This quantity compares the distribution of distances between all pairs of vectors in the reconstructed phase space with distances between different orbits evolving in time.For random signals, the value of the spatiotemporal entropy will be 100%, whereas for deterministic signals the value lies in between 0 and 100.In our case the value computed was close to zero, showing perfect structure in the data.From the above observations, we can safely conclude that the geomagnetic horizontal intensity H does not belong to the family of linear stochastic signals.Thus, the fluctuations observed in the time series of H are due to a low dimensional deterministic process.
It has been observed that a time series of coloured noise can also exhibit a power law power spectrum ω −α which is typical of a chaotic time series.The correlation dimension D is related to α by D = 2/(α − 1) and hence, a coloured noise time series can also lead to a low correlation dimension (Osborne and Provenzale, 1989;Theiler, 1991).A time series of coloured noise can be generated by letting where the coefficients C k , related to the power law power spectrum P (k) ∝ k −α by and the phases φ k , are randomly distributed in [0, 2π].Pavlos et al. (1992) have discussed the methods to compare the chaotic characteristics of an experimental time series with the pseudo-chaotic characteristics of a coloured noise time series.
The autocorrelation coefficients of H and that of the coloured noise time series with same fractal dimension is plotted in Fig. 7.For small lag times, the autocorrelation coefficient of H drops quickly, whereas that of the coloured noise decreases very slowly.There is also a significant difference in the decorrelation time.The autocorrelation coefficient of H drops abruptly by days 14 and 15 and reaches zero by day 40.However, the coloured noise takes about 72 days to become decorrelated.
The randomization of the phases of a deterministic chaotic signal can destroy its low dimensional chaotic profile.However, this profile of a coloured noise time series remains invariant under the phase randomization (Pavlos et al., 1992).We have compared the low dimensional profile of the time series of H to that of the time series obtained by the randomization of phases of the time series of H.The phase randomized time series of H was obtained from the Fourier series representation of H where the phases φ k are distributed randomly on the interval [0, 2π]. Figure 8a shows the local slopes of the logarithm of the correlation sums of the phase randomized time series of H and Fig. 8b shows that of the original time series of H .It is evident from these figures that the low dimensional chaotic profile is destroyed by the phase randomization.
We then compared the stationarity of these two times series.The fractional Brownian motions with coloured noise profiles can have a low dimensional character, but they are not stationary signals unlike motions on a strange attractor (Pavlos et al., 1992).For a stationary process, the probability density in phase space must be invariant with time.The probability density of the time series of H is calculated by dividing the range of the values of H into short intervals and counting the values of H which fall in those intervals.Figure 9a shows the probability density calculated for the original times series H.The solid line shows the probability density calculated for the entire time series and the filled circles corresponds to the first half of the time series.The coincidence is very clear in this figure.Figure 9b shows the probability calculated for the phase randomized time series of H.In this figure, the solid line corresponds to the entire time series and the dashed line with filled circles corresponds to the first half of the series.The coincidence in this figure is very low and this reveals the nonstationary character of the phase randomized signal.These results also show the deterministic character of the time series of H.
The power spectrum of the time series of H shows exponential decay which excludes the possibility of random behaviour and thus, indicates the chaotic behaviour of the time series (Fig. 10).It is also evident from Fig. 10 that there are five distinct peaks and they corresponds to A p = 300 nt, 300 nt, 179 nt, 236 nt, 300 nt, respectively, and hence, represent severe storms.
One important part of the nonlinear time series analysis is the calculation of Lyapunov characteristic exponents which, if positive, by definition, are the most striking evidence of chaos (Kantz, 1994).The first algorithm to a compute Lyapunov exponent for a time series was introduced by Wolf et al. (1985).One of the drawbacks of the algorithm of Wolf et al. (1985) is the strong dependence of the embedding dimension (Kantz, 1994).In this work, we used the algorithm developed by Kantz (1994) to calculate the maximum Lyapunov exponent.According to this method, one has to compute for a point y n o of the time series in the embedding space, where U(y n o ) is the neighborhood of y n o with diameter r.This is repeated for many values of n 0 so that the fluctuations of the effective expansion rates are averaged out.For an intermediate range of values of n, S( n) increases linearly with slope λ, which is our estimate of the maximal Lyapunov exponent.In this work, we have repeated the computation of S( n) for different values of the embedding dimension m and the diameter of the neighbourhood r.The maximum Lyapunov exponent was estimated to be λ = 0.25 ± 0.006 for the time series of H (Fig. 11).Sensitive dependence on initial conditions is an important feature of a chaotic system.The average exponential rate at which initial trajectories diverge is described by the Lyapunov exponent.Thus, the average error made when forecasting the future measurements of a chaotic system increases exponentially with time, even though it is governed by a deterministic law.The length of the period, over which accurate, short-term predictions of successive fluctuations of the signal, is determined by the accuracy of the initial condi- tions and the Lyapunov exponent.The prediction of subsequent values of the times series of H , using a simple but robust zero-order nonlinear prediction method (Kantz and Schreiber, 1997), is given in Fig. 12, along with the normalized root mean-square errors.It may be noted that the root mean-square errors increases exponentially.This could be seen as further evidence of the deterministic nature of the fluctuations and also of the underlying chaotic behaviour.

Conclusion
In this work, an investigation of the fluctuations of the time series of the geomagnetic horizontal intensity H , using the tools of nonlinear time series analysis, has been carried out.The recurrence plots and the results of surrogate data method, in addition to the estimate of spatiotemporal entropy, shows the deterministic nature of the data.We have also compared the chaotic characteristics of the time series of H with the pseudo-chaotic characteristics of coloured noise time series generated by the phase randomization of the original time series.The results of the comparison also reveal the determinism in the data.In addition to the estimated value of the correlation dimension, the analysis of the data, according to the method of the closest false neighbours, also shows the low dimensional character of the underlying dynamics.The positive value of the maximum Lyapunov exponent and the exponential decay of the power spectrum shows the chaotic behaviour of the dynamics.Thus, the physical process underlying the fluctuations of H is deterministic, low dimensional and chaotic.We have also shown that the error involved in the short-term prediction of the successive values of H , using a simple but robust, zero-order nonlinear prediction method (Kantz and Schreiber, 1997), increases exponentially and this indicates the sensitive dependence on the initial conditions.Wang (1996) has proposed that a combination of the correlation dimension of the attractor of the storm data and the magnetic index k could perhaps better describe the degree of solar disturbance than the single parameter k.We have estimated both the geometric and dynamical invariants of the attractor for the geomagnetic horizontal intensity, such as the correlation dimension and Lyapunov exponents, and we feel that it would be more advantageous to include the dynamical invariants, such as Lyapunov exponents, into this combination to describe the solar disturbance.In general, the quantities involved in chaos theory could be used to characterize the fluctuations of the geomagnetic intensity and hence, to describe the associated phenomena more accurately.The results of the analysis could also have implications in the development of a suitable model for the daily fluctuations of geomagnetic horizontal intensity.
Fig. 1.(a) Time series of the geomagnetic field H.(b) Delay representation of the time series of H.

Fig. 2 .
Fig. 2. The fraction of the closest false neighbours as a function of the embedding dimension m for the time series.

Fig. 3 .
Fig.3.The space-time separation plot for the time series of H. Initially, the curves increase quickly and later, the stable diurnal oscillations repeat.
Fig. 5. (a) The local slopes of the correlation sums of the surrogates as function of ln(r) for m = 12, τ = 10 and ω = 100.(b) The mean values of the local slopes of the surrogates with standard deviation for m = 12, τ = 10 and ω = 100.(c) Plot of the significance of the difference S versus ln(r).

Fig. 6 .
Fig. 6.The recurrence plot of the time series of H.

Fig. 7 .
Fig. 7.The autocorrelation coefficient of the time series of H and that of the coloured noise with the same correlation dimension D = 3.7.
Fig. 9. (a) The probability density function of the original time series of H based on the entire time series (solid line) and on the first half of the series (filled circles).(b) The probability density function of the phase randomized time series of H based on the entire time series (solid line) and on the first half of the series (dashed line with filled circles).

Fig. 10 .Fig. 11 .
Fig.10.The power spectrum of the times series with w = 0.001.The first part of the spectrum corresponds to an abrupt decay and the second part corresponds to a slow decay, as shown by the discontinuous lines.

Fig. 12 .
Fig.12.The predicted (using the zero-order nonlinear prediction method) values of the times series of H is plotted along with the actual values and the normalized error.