Dynamical characteristics of magnetospheric energetic ion time series: evidence for low dimensional chaos

. This paper is a companion to the ﬁrst work (Pav-los et al., 2003), which contains signiﬁcant results concerning the dynamical characteristics of the magnetospheric energetic ions’ time series. The low dimensional and nonlinear deterministic characteristics of the same time series were described in Pavlos et al. (2003). In this second work we present signiﬁcant results concerning the Lyapunov spectrum, the mutual information and prediction models. The dynamical characteristics of the magnetospheric ions’ signals are compared with corresponding characteristics obtained for the stochastic Lorenz system when a coloured noise perturbation is present. In addition, the null hypothesis is tested for the dynamical characteristics of the magnetospheric ions’ signal by using nonlinear surrogate data. The results of the above comparisons provide signiﬁcant evidence for the existence of low dimensional chaotic dynamics underlying the energetic ions’ time series.


Introduction
The hypothesis of the chaotic nature of the magnetospheric behaviour has been supported in the last decade by a large number of theoretical studies (Pavlos, 1988;Baker et al., 1990;Klimas et al., 1991Klimas et al., , 1996Pavlos et al., 1994). These results were further supported by a series of studies including chaotic analysis of ground measured magnetospheric signals with AE index data (Vassiliadis et al., 1990(Vassiliadis et al., , 1992Roberts et al., 1991;Shan et al., 1991;Pavlos et al., 1992aPavlos et al., , b, 1994Klimas et al., 1996). A series of significant studies (Price and Prichard, 1993;Price et al., 1994;Prichard, 1995) showed the weakness of the nonlinear analysis of the magnetospheric AE index time series, especially in relation to the strong null hypothesis of Theiler (Theiler et al., 1992a, Correspondence to: M. A. Athanasiu (mathanas@ee.duth.gr) b). In a recent series of papers by Pavlos et al. (1999a, b, c), Athanasiu and Pavlos (2001), the hypothesis of magnetospheric chaos, supported by the analysis of the magnetospheric AE index, was reestablished. Moreover, the application of chaotic analysis to signals observed by more than one spacecraft was examined for the first time by Pavlos et al. (1999c). In this paper a more thorough study of the magnetospheric ion signals is presented. In the first work (Pavlos et al., 2003) the correlation dimension and other geometrical quantities were estimated by using the magnetospheric ions' time series, observed during days 7-8 December 1994 at the dawn magnetosheath of the Earth's magnetosphere, and they were used as discriminating statistics between nonlinear dynamics and linear stochastic signals. The null hypothesis which was tested, concerns the observed time series that arises by a static nonlinear distortion of a Gaussian signal x(t) = h(s(t)), where h is a monotonic nonlinear function. (Theiler 1992a, b;Schreiber and Schmitz, 1996;Schreiber, 1998). Statistics of geometrical characteristics showed that the magnetospheric ions' time series was clearly distinguishable from a Gaussian, linear stochastic signal that had the same power spectrum and amplitude distribution. We must point out here that the geometric characteristics are measures of the spatial distribution of the sample points along a system orbit in the reconstructed phase space of the system. In this case there is no information about the dynamic evolution of the system in the phase space. Dynamic characteristics that connect current and future states of the system are the Lyapunov exponents, the average mutual information, the local linear prediction and nonlinear modeling. The results presented in Pavlos et al. (2003) indicate clearly the existence of a nonlinear, low-dimensional dynamical process underlying the magnetospheric energetic ions' time series. In order to test this hypothesis here we use dynamical characteristics as discriminating statistics between the magnetospheric ions' time series and stochastic (surrogate) data. The results obtained indicate the existence of low-dimensional chaotic dynamics underlying the energetic ions' time series.
It is important to note here that the prediction models are used only as indicators of the nonlinear and chaotic behaviour of the magnetospheric dynamics. Moreover, in a series of previous papers, prediction methods have been applied for understanding the physical process of solar windmagnetosphere coupling (Bargatze et al., 1985;Goertz et al., 1993;Klimas et al., 1996Klimas et al., , 1997Vassiliadis et al., 2000). In Sect. 2, the theoretical part of our study is presented in relation to the dynamical characteristics of the energetic ions' time series. In Sect. 3 a comparison between the energetic ions' time series and the surrogate data is carried out based on the dynamical characteristics, i.e. spectrum of Lyapunov exponents, mutual information, local linear prediction and global linear and nonlinear polynomial fitting. Finally, in Sect. 4 a summary is given and the results of this work are discussed.

Theoretical framework
In this section we present some theoretical concepts concerning the spectrum of the Lyapunov exponents, the mutual information, the modeling and prediction, which constitute the main tools of our analysis.

Spectrum of the Lyapunov exponents
The spectrum of Lyapunov exponents measures the rate of convergence or divergence of close trajectories in all d directions of the phase space. 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. For flows, there is always a zero Lyapunov exponent corresponding to the direction of the flow. The Lyapunov exponents can thus be ordered as λ 1 ≥, ..., ≥ λ d , and a positive λ 1 indicates the existence of chaos for a dissipative deterministic system.
The spectrum of the Lyapunov exponents can be estimated from a time series by following 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 } is a new set of orthogonal vectors produced by orthonormalization of the vectors at time i in order to retain the local orthogonal spanning of the state space (Sano and Sawada, 1985;Eckman et al., 1986;Holzfus and Lauterborn, 1988;Karadonis and Pagitsas, 1995). For the estimation of the maximum Lyapunov exponent (L max ) we use the equation where d(t) = |x 2 (t) − x 1 (t)| measures the distance between neighbouring points in the reconstructed phase space (Wolf et al., 1985). It follows, for finite data, that the initial d(0) is limited by the distance of the closest neighbors and the time t is limited by the time period of the observation.

Dynamics and mutual information
Chaotic or stochastic dynamical systems can be described by using the concept of information. For this scope we suppose that the random behaviour of the system is a realization of Shannon's concept of an ergodic information source (Shaw 1981(Shaw , 1984Abarbanel et al., 1993). If S is some property of the dynamical system and s i , i = 1, 2, . . . possible values of S, then the average amount of information gained from a measurement that specifies S is given by the entropy H (S) where P (s i ) is the probability that S equals s i and the logarithm is taken with respect to base 2. An estimate of P (s i ) is given by n(s i )/n T , where n(s i ) is the number of times that the value s i is observed, and n T is the total number of measurements. The same concept can be used to identify how much information we obtain about a measurement of an observable S quantity from measurement of another observable Q quantity. This concept is the basis for the definition of mutual information. For a time series, we consider a general coupled system (S, Q) with Q = {x(i)} and S = {x(i + τ )}, where x(i), x(i+τ ) corresponds to scalar samples from a dynamical system at discrete times t i and t i+τ . The conditional uncertainty of S given that Q = q i is defined as where P (s j /q i ) is the conditional probability of S = s J given that Q = q i . Thus, we define the conditional uncertainty of S given Q as a weighted average of the uncertainties H (S/Q = q i ), that is Using the fact that P (q i , s J ) = P (q i )P (s j /Q i ) we have The amount by which a measurement of Q reduces the uncertainty of S (average mutual information) is given by the relation For more details about the definitions of the above quantities refer to Ash (1990) and Papoulis (1991). If this relation is applied to time series leads to The mutual information between the two samples independent, then the mutual information will vanish for this value of τ , i.e. knowledge for the second sample can not be gained by knowing the first. On the other hand, if the first sample uniquely determines the second sample, then I (τ ) = I max , which is most likely to be true when τ = 0. In this paper, we follow the work of Fraser and Swinney (1986) for the estimation of the mutual information (according to Eq. 7) of an experimental time series, which is used as discriminating statistics between the surrogate data and the ions' time series.

Modeling and prediction
The observable information x(ti) ≡ x(i) on the temporal evolution of the orbit in the reconstructed phase space can be used for prediction or modeling purposes, building the map F (x(i)). The estimation of the map F , for experimental data and for T time steps ahead, F T is reduced to determine a group of parameters a, given a class of functional forms for F T (x(i), a). When the functional form of F T is known, then the parameters a are chosen by using some form of a cost function which measures the matching of the observed future sample x(i + T ) with the predictedx(i + T ) = F T (x(i), a) (see Abarbanel et al., 1993). Our purpose is to construct both linear and nonlinear maps, as well as local or global maps by looking for the best approximation of F T . Modeling and prediction make use of phase space reconstruction, implicitly (like the autoregressive (AR) models) or explicitly (like the local linear maps) (see Weigend and Gershenfeld, 1993). Moreover, the map F T (x(i), a) may be approximated with different functional forms of the global, local or semi-local type (Lillekjendlie et al., 1994). For any functional form, the parameters involved are estimated directly from the data. In the case of modeling the whole set of available data is used in fitting, while in the case of prediction a subset of the data is used initially for fitting an appropriate model (training set) and the rest is used for testing the predictability of the proposed model (test set).
To verify the performance of the model, two measures of the modeling or prediction error are often computed. The first is the normalized root mean square error (NRMSE), which is defined as the root of the mean square differences of the modeled (or predicted) values from the actual values, normalized by dividing by the standard deviation of the data. If NRMSE ∼ = 0, perfect model performance is achieved and if NRMSE > 1, the modeling or prediction is worse than this obtained using the mean value as a model. The other measure is the correlation coefficient (CC), which gives the correlation between the modeled (or predicted) data and the actual data. An estimate of CC is obtained from the ratio of the covariance over the root of the product of the variance of the two data sets. It is known that CC takes values in [-1, 1]. When CC ∼ = 1 best correlations are obtained, i.e. the performance of the model is excellent, while for values of CC close to zero or negative, the performance is very poor. In this work, F T is approximated with global polynomials (linear and nonlinear) and local linear models.

Local models
In the case of local models which describe deterministic systems, nearby trajectories evolve similarly, at least for a short time if the system is chaotic. Thus, on the reconstructed attractor, for any point x(i) we can locally approximate F T which leads to the estimation of x(i + T ), taking into account the k nearest neighbours of x(i), {x(i(1)), ..., x(i(k))}. It should be noted that for each target point x(i) and time step T a different model is computed.
The local approximation of F T may be done with a linear map of the formx(i + T ) = a 0 + a T x(i). Assuming that this model is good enough for the neighborhood of x(i), i.e. the above equation also applies to the mapping of the k neighbours, where k > m, we can solve a system of k equations with m + 1 unknown parameters {a 0 , a} using the ordinary least squares (OLS). The values of the estimated parameters are then used to find the mappingx(i + T ) (Farmer and Sidorowich, 1987;Casdagli et al., 1992).

Global polynomial models
A simple approximation of F T is a polynomial, which may involve only linear terms (this is actually the autoregressive (AR) model) or linear terms plus nonlinear terms of a degree q as well. Certainly, a polynomial of a small degree q in m delay variables, call it p T q , cannot model complicated dynamics and for pure chaotic systems is more likely to be insufficient. However, when we are dealing with real data,  this may turn out to be an advantage, because the evident dynamics (linear or nonlinear) are only modeled and the hidden are not modeled, as they may be covered by noise as well (Barahona and Chi-Sang Poon, 1996).
The general form of p T q , wherex i+T = p T q (x i ), x ∈ R m , is given by the Volterra-Wiener series of degree q and memory m . In our case we use q = 2 because we are only interested in investigating the existence of nonlinearity in the data. Moreover, we construct all M polynomials, starting with the simplest which contains the first linear term and go on by adding one term of the Volterra series at a time.

Data analysis and results
In the following we present the results about the dynamical characteristics of magnetospheric ions' time series and the corresponding discriminating statistics which involve surrogate data. The energetic ions' time series (35-46.8 keV) was observed by the experiment EPIC/ICS during the days 7-8 December (days 341-342) 1994 at the dawn magnetosheath of the Earth's magnetosphere. As was explained in the first paper of this study (Pavlos et al., 2003), it is reasonable to suppose that these particles were accelerated in the inner magnetosphere during periods with strong coupling of the magnetospheric system and the solar wind, simultaneously with strong bursts of electrons and O + and a clear enhancement of the AE index. Therefore, it can be supposed that the dynamics of the energetic ions mirror the internal magnetospheric dynamics, similar to the AE index during periods with a strong coupling of the magnetosphere and the solar wind (Pavlos et al., 1999c). The energetic particle differential fluxes are provided via the Energetic Particle and Ion Composition (EPIC) instrument of the GEOTAIL spacecraft and essentially remained close to the ecliptic plane (Williams et al., 1994). The sampling time for the energetic ions analyzed here was 6 s.
Moreover, in order to understand the underlying dynamical process of the energetic ions' time series, we study the first component (V 1 -component) and the reconstructed time series (V 2−10 component), according to the theoretical concepts of singular value decomposition (SVD) analysis (Broomhead and King, 1986;Pavlos et al., 2003).
As it is shown in Athanasiu and Pavlos (2001) the original time series x(t) can be reconstructed by using n new time series V i (t) according to the relation where n is the embedding dimension and V i (t), i = 1, ..., n are given by the first column of the matrix (Xc i )c T i , X is the trajectory matrix constructed by the observed time series and c i are the eigenvectors of the covariance matrix X T X (Broomhead and King, 1986). The reconstructed time series corresponds to the component V 2−10 = V i (t), i = 2 − 10. It was shown in Pavlos et al. (2003) that the V 1 -component of the ion time series corresponds to a coloured noise external component of the magnetospheric dynamics, while the V 2−10 component corresponds to the internal magnetospheric.
3.1 Lyapunov exponents Figure 1 shows the first 6 Lyapunov exponents estimated according to Eqs. (1)-(3) of Sect. 2.1, for the energetic ions, and the reconstructed and V 1 -component time series as a function of embedding dimension. Figure 1a clearly shows that the first Lyapunov exponent of the energetic ions' time series is positive (∼0.26 bit/s), while the second is marginally positive, the third is marginally negative and the others are negative. The spectrum of Lyapunov exponents for the reconstructed time series (see Fig. 1b) also presents one positive Lyapunov exponent (∼0.27 bit/s), one zero exponent, while the others are negative. This observed profile of the spectrum of the Lyapunov exponents corresponding to the reconstructed time series is in good agreement with the spectrum of a strange attractor, where positive exponents signify mechanisms of instability, negative exponents correspond to mechanisms of convergence, while at least one zero exponent must exist corresponding to the expansion along the trajectory. Figure 1c shows the spectrum of the Lyapunov exponents corresponding to the V 1 -component. In this case the first Lyapunov exponent is much smaller than the first Lyapunov exponent obtained from the original ions' time series and the V 2−10 reconstructed time series. The observed profile of the V 1 -component is similar to the profile which is expected for a coloured noise signal.
In addition, the maximum Lyapunov exponent L max has been estimated independently according to Eq. (4). Figure 2 shows the L max for the energetic ions, the reconstructed and V 1 -component time series as the state space trajectory evolves. In this case the L max was found to be approximately the same (∼0.075 bits/s) for both the energetic ions and the reconstructed time series, and many times higher than the L max of V 1 -component, which tends to a zero value. These results provide further evidence for the existence of chaotic behaviour in the energetic ions and the reconstructed time series, and suggest that the V 1 -component corresponds to different, possibly linear, contamination with noise dynamics.
For a purely deterministic system the existence of positive L-exponents implies chaotic dynamics, but for a signal contaminated by noise it is possible for some L-exponents to be positive due to the stochastic perturbation. It is known that stochastic data often give positive L-exponents without the underlying dynamics necessarily being chaotic (Osborne et al., 1986;Provenzale et al., 1991;Argyris et al., 1998). When the embedding dimension is larger than the degrees of freedom of the underlying system, spurious exponents occur as an artifact of the embedding, but the existence of more than one positive L-exponents for the energetic ions' time series constitutes significant evidence that at least one of them  could correspond to the underlying deterministic dynamics.
In order to verify the validity of this statement, we estimate the L max for the surrogate data according to the null hypothesis that the energetic ions' time series belongs to a family of linear stochastic signals transformed by a nonlinear static distortion. As in Pavlos et al. (1999a), we use in the following the surrogate data, as described in Schreiber et al. (1996), which mimic the time series in relation to the amplitude distribution and the autocorrelation function. In order to obtain convincing results, we created a rich sample including 40 surrogate data. The results of this comparison between the energetic ions' time series and the surrogate data are shown in Fig. 3a. The discrimination between the original time series and the surrogate data for L max is possible, because the significance of the statistical test remains higher than two sigmas, fluctuating at the value of ∼ = 4 sigmas (see Fig. 3d). These results for the L max as a discriminating statistic clearly permit the rejection of the null hypothesis with a confidence greater than 95%. Figure 3b presents the maximum Lyapunov exponent for the reconstructed time series and the surrogate data, while Fig. 3d shows the significance of the statistical test which approaches the value ∼20 sigmas, much greater than the original one of ∼4 sigmas. This result rejects the null hypothesis with high confidence. Figure 3c is similar to Fig. 3b but it shows the L max for the V 1 -component and the surrogate data. In this case the discrimination between the V 1 -component and the surrogate data is impossible because the significance of the statistical test remains lower than 2 sigmas, as seen in Fig. 3d. The above results suggest that the reconstructed time series is nonlinear in a more evident way than the original time series of the energetic ions, including the main magnetospheric dynamics, in contrast to V 1 -component, which corresponds to external dynamics of a coloured noise. The low value of L max for the V 1 -component, as well as the previous resultsrelated to the high correlation dimension, the long  decorrelation time and the nonstationarity, set more clearly the question of whether the V 1 -component could be a form of coloured noise (Athanasiu and Pavlos, 2001). In order to support further this argument about the V 1 -component, we estimated the L max for the x(t) variable of the Lorenz chaotic system, as well as for its V 2−10 reconstructed component and its V 1 -component. In this case the Lorenz system was perturbed by external coloured noise of low dimensionality (Provenzale et al., 1992). Figure 4 shows the maximum Lyapunov exponent for the x(t) variable, its reconstructed and V 1 -component time series. It is observed that the L max of the reconstructed and V 1 -component approximately equals the L max of the original time series, which stabilizes at the value of ∼2.2, as it is known. This result is expected because the reconstructed and V 1 -component time series are a linear transformation of the original time series according to SVD analysis (Athanasiu and Pavlos, 2001). Figure 5a presents the L max for the variable x(t), contaminated with 110% observable colored noise, as well as for the surrogate data. The value of L max for the noisy time series was estimated to be ∼4.5, larger than the corresponding value of ∼2.2 of the original Lorenz times series, due to the coloured noise influence (Argyris et al., 1998). Nevertheless, as shown in the same figure, the value of L max ∼ = 10 for the surrogate data is about two times that of the original time series. Figure 5b-c are similar to Fig. 5a but correspond to reconstructed and V 1 -component time series of the noisy variable x(t). For the reconstructed time series, the L max takes the value of ∼4, while for the surrogate data L max takes the value of ∼15. On the other hand, for the V 1 -component and the surrogate data no difference was found, as the value of L max is ∼1 for both of them. Finally, the significance of the statistical tests of the L max for the noisy, its reconstructed, and V 1 -component time series are shown in Fig. 5d. We observe that the discriminating statistics approaches the value of ∼50 sigmas for the reconstructed time series, obtaining larger values than those of the noisy variable x(t), which is fluctuated in the region of 20 sigmas. These results permit us to reject the null hypothesis with a confidence larger than 95% for both cases, while for the V 1 -component it is not possible, as the quantity S takes values smaller than 2 sigmas.
Studying the previous results, related to the estimation of L max and the test for the surrogate data, it turns out that for a chaotic system perturbed with coloured noise the reconstructed time series preserves its main characteristics, while the V 1 -component exhibits the characteristics of the external coloured noise perturbation. Similar results were found from the estimation of geometrical characteristics of the Lorenz system (Pavlos et al., 2003;Athanasiu and Pavlos, 2001). At this point we must emphasize the strong similarity between the results obtained from the L max estimation for the energetic ions' time series and the Lorenz noisy variable x(t), as well as for the reconstructed V 2−10 and the V 1 -component, as shown in Figs. 3 and 5. These results are in agreement with similar results related to the AE index (Athanasiu and Pavlos, 2001) and support the theoretical concept that different data reveal two different physical processes: the one process which shows low dimensional and chaotic behavior corresponding to the internal magnetospheric dynamics and the other which looks like a linear stochastic process corresponds to an external coloured noise perturbation of magnetosperic dynamics.

Mutual information
Mutual information for the energetic ions time series and the surrogate data has been estimated by implementing the algorithm of Eq. (7) in Sect. 2.2. Figure 6a shows the mutual information estimated for the energetic ions' time series and the surrogate data as a function of the lag time τ . The mutual information for the energetic ions time series is slightly larger than that for the surrogate data. This difference is significant enough for the discrimination between the energetic ions time series and the surrogate data (Fig. 6b). In particular, for small τ , e.g. the first 15 lags, with the exception of τ = 1, the values of the statistical test were significant in the range of 2-4.5 sigmas. For larger lags the significance remains smaller than two sigmas. The above results permit us to reject the null hypothesis and support the nonlinearity of the energetic ions time series. Figure 6c shows the mutual information for the reconstructed time series that has smaller values than those of the energetic ions' time series for the same lag time. However, the difference with the surrogate data is large enough, as shown in Fig. 6d, where the significance of statistics approaches the value of ∼40 sigmas for 12 lags (∼2.5 min), permitting the rejection of the null hypothesis at a level of confidence larger than 95%. Figure 6e is similar to Fig. 6a but corresponds to the V 1 -component and the surrogate data. The values of mutual information for the V 1 -component decay slower than those of the energetic ions and the reconstructed time series. This result is expected if we take into account the high linear autocorrelations of the V 1 -component. The significance of the statistical test for the mutual information of the V 1 -component remains lower than two sigmas (see Fig. 6f), except in a narrow region of 2-7 lags, where it approaches the value of ∼4 sigmas. This result is in contrast to several results coming from the other tests of chaotic analysis for the V 1 -component, suggesting that a part of the V 1component is partially explained by nonlinear terms, something which is not excluded by the method of its construction.  Eventually, the study of mutual information exhibits weak nonlinearity in the original time series of energetic ions, possibly because of the existence of a strong linear component (see Figs. 6b, f). On the other hand, the reconstructed time series after the removal of the first component, reveals nonlinear behaviour, as shown in Fig. 6d. The above results suggest that the SVD method can be a very useful tool in the filtering and detection of nonlinearity in observable data.

Parametric local linear prediction
In the following we present the results of the hypothesis test using the prediction and modeling methods described in Sect. 2.3. In particular, we use the local linear prediction (LLP) model, as well as global linear and nonlinear polynomial fitting. In the first case we estimate the local linear maps using ordinary least-square fitting (OLS). Figure 7a presents the results of LLP with the surrogate data. The correlation coefficient (CC) for T = 1 − 50 step ahead predictions obtained iteratively for the energetic particles' time series and the 40 surrogate data are presented. In order to reduce the noise of the estimations we used a low embedding dimension m = 6, the best reconstruction delay time τ = 30, and k = 2(m + 1) = 16 neighbours (Farmer and Sidorowich, 1987). The training set consisted of the first 18 000 values and the test set consisted of 1800 values, which is 10% of the training set, of the time series. The significance of the above statistics is shown in Fig. 7b. The small values of the significance (0.0 − 1.4) show that there is no significant difference between the original time series and the surrogate data. This result is in agreement with the observation of Theiler and Prichard, that the prediction error generally has a large variance and provides a poor discriminating statistic (Theiler and Prichard, 1997).  It is known for a chaotic time series that the predictability depends on the width of the training set. The prediction error is decreased when the width of the training set is increased (Farmer and Sidorowich, 1987;Smith, 1992). If the energetic ions' time series includes chaotic determinism, then we expect an improvement of the predictability with the increase of state vectors in the training set. Also, a similar behaviour is expected for the surrogate data because of their reconstruction (they have the same autocorrelation function and amplitude distribution with the original time series). However, this improvement must be smaller than that of the original time series, because of the stochastic character of the surrogate data by construction (in this case the phases have been destroyed). In the following paragraph we are going to check this crucial hypothesis on the energetic ions' time series and the surrogate data. The correlation coefficient between real and predicted values has been used as a prediction estimator.
In Fig. 8a the correlation coefficient between the real and predicted values is shown as a function of prediction time when the training set has the first 3000, 9000, 18 000 data points of the time series. In every case the width of the test set is chosen to be equal to 10% of the corresponding training set. In addition, the parameters of the prediction, embedding dimension, delay time, and number of neighbour state vectors are equal to those of Fig. 7a. As we can see, the values of the correlation coefficient are increased when the number of data in the training set is increased. Figure 8b is the same as Fig. 8a but corresponds to the mean value of correlation coefficient for the surrogate data. Comparing the two figures, we observe that the mean predictability of the surrogate data is better than the original when the training region is the first 3000 points. Surprisingly, in the case of the original time series the difference between the correlation coefficients corresponding to 3000 and 18 000 data points of the training set is approximately twice that of the surrogate data. This result by itself constitutes evidence of determinism and chaoticity in the energetic ions' time series. In order to es-  timate the significance of the alteration as the width of the training set is increased, we calculated the absolute difference between the correlation coefficients corresponding to a training set of 3000 and 18 000 data points for the energetic ions' time series and the surrogate data (see Fig. 8c). The quantity S shown in Fig. 8d tends to be significant in the region of 2-5 prediction steps as it obtains values larger than two sigmas. Figure 8e is similar to Fig. 8c but with a training set of the first 3000 and 9000 data points. The corresponding significance of the statistics is shown in Fig. 8f. In this  case the small values of the significance indicate that there is no significant difference between the original time series and the surrogate data. The above results concerning the predictability in relation to the width of the training set support further the existence of nonlinearity and determinism in the magnetosperic dynamics.

Global linear and nonlinear polynomial fitting
The modeling with global polynomials gives good discrimination between the original and surrogate data, especially when we define the discrimination statistic to be the change in the modeling error as we go from linear to nonlinear polynomial terms. Here, we use the NRMSE to quantify the modeling error as presented in Sect. 2.3.2. In Fig. 9a the NRMSE for energetic ions and 40 surrogate data is shown as a function of the polynomial terms of the Volterra-Wiener series, using m = 6 and q = 2 and 5 steps ahead mappings for all the polynomials of the Volterra-Wiener series (see Eq. 8). The first 7 polynomial terms are linear (the first term is the constant a 0 ) and the rest are nonlinear interactions of the 30 delays. The existence of noticeable nonlinearity for the energetic ions' time series is proved by the abrupt reduction of the energetic ions' error, as we pass from the linear to the nonlinear fitting. This significant result is clearly depicted in Fig. 9b, that shows the significance of the reduction error statistic corresponding to the difference of the polynomial fitting of the first seven linear terms from the polynomial fitting containing all the possible linear-nonlinear terms. For the difference corresponding to linear terms, the significance of the statistic reveals values lower than 2 sigmas, while for nonlinear terms the significance varies at higher values in the range of 4.0 − 5.0 sigmas. Figure 9c shows the discriminating statistic of the reduction error for the energetic ions (tall line) and the surrogate data (short line) corresponding to the difference of the polynomial fitting of the first seven linear terms from the polynomial fitting containing all the linear-nonlinear terms. This figure shows that the difference is clearly significant permitting the rejection of the null hy-  pothesis with a confidence larger than 95%.
In Fig. 10a, the NRMSE for the reconstructed time series of energetic ions and 40 surrogate data is shown as a function of the polynomial terms of the Volterra-Wiener series, using parameters m = 6, q = 2, τ = 7, and T = 1 step ahead. As in the case of the energetic ions' time series, an abrupt decay is observed again for the NRMSE as we pass from linear (the first seven polynomial terms) to the nonlinear fitting. Figure 10b is similar to Fig. 9b but corresponds to a reconstructed time series. The significance of statistics remains lower than two sigmas for linear models, while it is increased with respect to the number of the nonlinear terms included in the model approaching the value of 5.25 sigmas. Moreover, in Fig. 10a, a reduction of NRMSE is observed for the reconstructed time series as nonlinear terms are added in Volterra-Wiener series. This result means that for the reconstructed time series, nonlinearity does not only appear in the first nonlinear terms as in the case of the original time series (see Fig. 9a). On the other hand, the NRMSE of the surrogate data remains unchanged when the first nonlinear polynomial  term is included in the model. This result is quantified in Fig. 10c, which shows the significance of the reduction error statistics for the difference between the polynomial fitting of the first 8 terms from the polynomial fitting containing all the possible linear-nonlinear terms of polynomial fitting. For the difference corresponding to linear terms, the significance of the statistics reveals values with almost zero sigmas, while for nonlinear terms, the significance is increased, approaching a saturation value of ∼10 sigmas. The discriminating statistic of the reduction error of the polynomial fitting with the first 8 terms from the polynomial fitting with all 28 terms is shown in Fig. 10d, using tall lines for the reconstructed time series and short lines for the surrogate data.
Finally, Fig. 11 presents the NRMSE for the V 1component and the 40 surrogate data as a function of the polynomial terms of the Volterra-Wiener series, using m = 6, τ = 70, q = 2 and T = 10 steps ahead. Obviously, the addition of nonlinear terms does not cause a reduction of NRMSE, as it does in the original and reconstructed time series. This result is depicted in Fig. 11b, which shows the significance of the reduction error statistic corresponding to the difference between the polynomial fitting of the first seven linear terms from the polynomial fitting containing all the possible linear-nonlinear terms. As is expected from the previous tests, the significance obtains values lower than 2 sigmas and thus, the null hypothesis cannot be rejected.
The above results obtained by using a polynomial fitting, especially for the reconstructed time series reject the null hypothesis and strongly support the nonlinearity of internal magnetospher a dynamics with a confidence larger than 95%.

Summary and discussion
The dynamical characteristics of the magnetospheric ions' time series that were studied in this work reveal strong nonlinear and chaotic characteristics of the underlying magnetospheric dynamics which produces the energetic ions' signal. The null hypothesis was tested for the original time series, as well as for the first SVD component and the reconstructed signal by adding the other SVD components. This test was applied for the largest Lyapunov exponent, the mutual information, the local linear prediction and the Voltera-Wiener modeling. For the application of the test we have used surrogate data constructed according to Schreiber's algorithm, in order to mimic the amplitude distribution and the power spectrum of the magnetosperic signal and its SVD components.
In Table 1 we summarize the significance of the estimated statistics for the case of the largest Lyapunov exponent, the mutual information, the local linear prediction and the Voltera-Wiener modeling, for the ion time series, the first SVD component and the reconstructed signal. The significance of the statistics for the original signal obtains values slightly higher than two sigmas, except in the case of the local linear prediction, while the significance of the first SVD component reveal values lower than two sigmas, except in   The discriminating statistic of the reduction error for the energetic ions (tall line) and the surrogate data (short line), corresponding to the difference between the polynomial fitting including the first nonlinear term (8th term) from the polynomial fitting containing all the linear-nonlinear terms.  the case of mutual information. On the other hand, the significance of the statistics for the reconstructed signal by using the higher SVD components obtains values in the region of 10-40 sigmas, except in the case of local linear prediction and the Voltera-Wiener linear modeling. Moreover, the largest Lyapunov exponent was found to be positive, with values ∼0.07 s −1 for the original and the reconstructed time series, while for the first SVD component the value is significantly lower ∼0.01 s −1 . These results concerning the dynamical characteristics of the magnetospheric ions' time series are in good agreement with the results obtained in Pavlos et al. (2003), corresponding to the geometrical characteristics (correlation dimension, false nearest neighbors and singular values). Also, the results of chaotic analysis applied for the magnetospheric ions' time series are quite similar to the results of the chaotic analysis of the magnetospheric AE index (Pavlos et al., 1999a, b, c). As we know the magnetospheric energetic ions are observed during periods of strong substorms of the Earth's magnetosphere (Pavlos et al., 1985(Pavlos et al., , 1989. Therefore, it is acceptable from the physical point of view to assume that for both time series the AE index and the energetic ions correspond to the same dynamical system. This hypothesis is strongly supported by the new results of chaotic analysis obtained for the AE index and the energetic ions signal, and for both time series the chaotic analysis reveals a low dimensional process. As we have shown in Pavlos et al. (2003) the correlation dimension for the energetic ions was found to be ∼3.5 and the independent dynamical degrees of freedom estimated by false nearest neighbours and singular spectrum analysis was found to be ∼7. Finally, the low dimensional chaotic dynamics of the magnetospheric system is now better proved to be two different magnetospheric signals, one measured at the Earth (AE index) and the other at the distant magnetosphere, exhibiting common geometrical and dynamical characteristics in accordance with the low dimensional chaos hypothesis. This concept of low dimensional chaotic dynamics of the magnetospheric system is also supported strongly by the similar behaviour of the Lorenz dynamic system perturbed by external noise.
In the previous analysis the low dimensional chaotic behaviour of energetic ions' time series was mainly supported by the SVD reconstructed signal after rejection of the first SVD component. The comparison of the results of the magnetospheric signals with the results produced by the Lorenz system after its perturbation with a strong colored noise, indicates that the internal low dimensional chaotic dynamic of the magnetospheric system is perturbed by an external coloured noise signal that could be related to the stochastic dynamics of the solar wind system. From this point of view, in order to test the above hypothesis concerning the coloured noise, the extended chaotic analysis used in this paper must be applied simultaneously to the solar wind and magnetospheric system.