Nonlinear time-series analysis of Hyperion's lightcurves

Hyperion is a satellite of Saturn that was predicted to remain in a chaotic rotational state. This was confirmed to some extent by Voyager 2 and Cassini series of images and some ground-based photometric observations. The aim of this aticle is to explore conditions for potential observations to meet in order to estimate a maximal Lyapunov Exponent (mLE), which being positive is an indicator of chaos and allows to characterise it quantitatively. Lightcurves existing in literature as well as numerical simulations are examined using standard tools of theory of chaos. It is found that existing datasets are too short and undersampled to detect a positive mLE, although its presence is not rejected. Analysis of simulated lightcurves leads to an assertion that observations from one site should be performed over a year-long period to detect a positive mLE, if present, in a reliable way. Another approach would be to use 2---3 telescopes spread over the world to have observations distributed more uniformly. This may be achieved without disrupting other observational projects being conducted. The necessity of time-series to be stationary is highly stressed.


Introduction
Saturn's seventh moon, Hyperion, was discovered in the XIX century by Bond (1848) and Lassel (1848), but it took more than a century to obtain its images due to Voyager 2 (Smith et al. 1982) and Cassini (Thomas 2010) missions. Its shape is highly elongated (360×266×205 km), making it the biggest known highly aspherical celestial body in the Solar System. Wisdom et al. (1984) predicted Hyperion to remain in a chaotic rotational state due to its high oblateness and relatively high eccentricity, e = 0.1. In dynamical system theory, a chaotic behaviour is recognised through a positive maximal Lyapunov Exponent (mLE), which describes the rate of divergence (or convergence in the negative case) of initially nearby phase-space trajectories. The Lyapunov spectrum is relatively easy to calculate in the case when the differential equations are known (Benettin et al. 1985;Wolf et al. 1985;Sandri 1996;Baker & Gollub 1998;Ott 2002). On the other hand, there exist algorithms allowing to obtain an mLE from an experimental or observational time-series (Wolf et al. 1985;Kantz 1994), although they are to be used with carefullness, using at least a few hundred data points (Rosenstein et al. 1993;Katsev & L'Hereux 2003) for nonlinear analysis. It is a hard task in astronomy to obtain long-term, well-sampled lightcurves. Although, despite this difficulty, it has been efficiently shown that pulsar spin-down rates exhibit chaotic dynamics (Seymour & Lorimer 2013): by re-sampling the original measurements, artificial time-series were produced, equivalent to the original ones, containing such a number of data points that the calculation of the correlation dimension and the mLE of the attractor, reconstructed via Takens time delay embedding method, was possible.
Hyperion's long-term observations were carried out twice in the post Voyager 2 era. In 1987, Klavetter (1989 (hereinafter, K89) performed photometric R band observations over a timespan of more than 50 days, resulting in 38 high-quality data points. In 1999 and is excluded from this work. The brightness was reported to be constant over a time period of 6 hours at 0.01 mag level. Each night resulted in multiple, independent observations, such that some uncertainties were smaller than 0.01 mag. The data was corrected to mean opposition distances and to zero solar phase angle (for further details, see K89). For the purpose of this work it is useful to resample the unequally spaced data to form an artificial, uniformly spaced lightcurve. In order to do that a cubic spline was formed, which was next sampled with a step equal to the mean of the original dataset to retrieve a set consisting of the same number of points as the original. The statistical properties are gathered in Table 1 and the result of this procedure is displayed in Fig. 1a. This was to verify whether the cubic spline introduces or cancels some structures in the original lightcurve. The inspection of the Lomb-Scargle periodogram (Scargle 1982) in Fig. 1b shows that both time-series posses roughly the same modes. Next, the cubic spline was sampled with a step chosen so to form a 5000-points dataset. Fig. 1b presents also the periodogram for this final time-series, that will be used in the estimation of the mLE. Devyatkin et al. (2002) performed observations broken into two parts by a long interval: from September 1999 to March 2000 and from September to October 2000. Each of the CBV R bands in the respective period will be therefore referred to as C1, C2, and so on.
The intrinsic accuracy was obtained by means of the standard deviation of the average brightness relative to each of the comparison stars in the frame. These values have an average of 0.12 mag in the C band, 0.06 mag for the B band and 0.10 mag for both V and R bands (for further details, see D02). The B1 and V 1 datasets contain only 5 points, and therefore are excluded from the analysis in this work; the R1 and R2 ones contain 11 points each, though the first observation in R1 occured more than a month before the next one, therefore is excluded from the calculations. The second period is better sampled, with mean spacings from 3.21 to 4.39 days, and standard deviations about 80% of the mean. The C1 sample from the first period is spread over a longer interval, but with a mean spacing equal -6to 7.86 days (which is roughly 25% of the minimal Lyapunov time (Shevchenko 2002)) with a standard deviation of 9.50 days. Each sample was used to form a cubic spline to be spaced with a step equal to the mean of the original dataset. The power spectrum, constructed as in Fig. 1b, shows this procedure not to be as valid as in the K89 case, although with an undersampled time-series it is difficult to formulate unambiguous conclusions. On the other hand, it is visible in Fig. 1a that the cubic spline and sampling reproduce the original observations well and do not introduce any auxiliary peaks that would not follow the trend of the inspected time-series, what is also the case in D02 data, therefore we proceed with the analysis. Next, each cubic spline was sampled with a step chosen so to obtain artificial datasets consisting of 5000 points. These, as well as the one obtained from K89, are the subject of the current work.

Takens reconstruction
Before evaluating the mLE, it is insightful to reconstruct the phase-space trajectories via Takens time delay embedding method (Takens 1981). Having a series of scalar measurements x(t) uniformly distributed at times t one can form an m-dimensional location vector using only the values of x(t) at different times given by the classical formula where m is the embedding dimension, and τ is the time delay chosen so that the components of S(t) are independent or uncorrelated. These parameters are estimated via the False Nearest Neighbor (FNN) algorithm (i.e., embedding dimension m) and as the first minimum of the Mutual Information (MI) method (i.e., time delay τ ), using the programs described in (Kodba et al. 2005;Perc 2005aPerc ,b, 2006. These .exe programs (fnn, mutual and others) are available from the website 1 . The TISEAN software package (Hegger et al. 1999) is also widely used throughout this article 2 . The MI is an alternative for the commonly used autocorrelation function, where the time delay τ used to be estimated as the delay at which autocorrelation drops to 1/e (to ensure e-folding, or, rarely, to 0, in order to form uncorrelated location vectors). This approach will be also explored in the subsequent sections. 3 The MI gives the amount of information one can obtain about x t+τ given x t . In short, the absolute difference between x max and x min is binned into j bins, j being large enough, and for each bin the MI, denoted by I(τ ), is constructed from the probabilities that the variable lies in the h-th and k-th bins (P h and P k respectively and h, k = 1, . . . , j) and the joint probability P h,k that x t and x t+τ are in bins h-th and k-th, respectively: The FNN fraction is calculated under the assumption that the trajectory folds and unforlds smoothly. Roughly speaking, when two initially nearby points diverge under forward iteration (typically not longer than for a time τ ), they do so not faster than R tr , where is the initial separation (not greater than the standard deviation of the data) and a value significantly close to zero (in practical implementations to the threshold R tr ) for the embedding dimension considered to be a correct estimate of the proper one. This is achieved by calculating the FNN fraction for nearby points S(i), S(t) such that ||S(i) − S(t)|| < .
For a detailed description of these algorithms, see (Kodba et al. 2005;Perc 2005aPerc ,b, 2006 and (Hegger et al. 1999) for another common implementation. The parameters m and τ are both required for the mLE calculation. The embeddings using time delays from the autocorrelation function were also performed and are displayed in the Appendix.

Correlation dimension
The correlation dimension d C is a measure of how much space is covered by a set. For usual 1D, 2D or 3D cases the d C is equal to 1, 2 and 3, respectively. However, there are sets called fractals, which posses a fractional correlation dimension. For dissipative dynamical systems, e.g. the Lorenz system (Lorenz 1963), chaotic motion manifests itself through a limitting trajectory called the strange attractor, due to its fractal properties. Herein we examine the fractal dimension of the reconstructed phase-space trajectories in order to further constrain the proper embedding dimension. For Hamiltonian systems, i.e. volume preserving, this is not an indicator of chaoticity (Greiner 2010).
The correlation dimension is defined as with the estimate for the correlation function C(r) as where the Heaviside step function H adds to C(r) only points x i in a distance smaller than r from x j and vice versa. The limit in Eq. (4) is attained by using multiple values of r and fitting a straight line to the linear part of the obtained dependencies. The correlation dimension is estimated as the slope of the linear regression. The calculations were performed with a MATLAB code with the Theiler window equal to zero (Seymour, priv. comm.) and the error calculations are described in (Seymour & Lorimer 2013). The results in a graphical form are shown in Fig. 3; the numerical values are gathered in Table 2.
Following the reasoning of Seymour & Lorimer (2013) one could infer, because the fractal dimensions are not much greater than unity, that there would be a total of two dynamical variables governing Hyperion's rotation at the time of observations. However, it is well known that chaos can not occur in a two dimensional continuous dynamical system (see Poincaré-Bendixson theorem, e.g. (Alligood et al. 2000)); therefore, there must be at least three variables to consider the rotation being chaotic. This requirement is met for example in the simplified model (with 1.5 dof) in which the axis of rotation is fixed and perpendicular to the orbit plane (Wisdom et al. 1984;Celletii & Chierchia 2000). This leads to the suspicion that the underlying dynamics are not governed by the chaotic zone, i.e. Hyperion remained in a regular (quasi-periodic) state, possibly influenced by noise.
On the other hand, the datasets with m = 2 are undersampled and the chaotic behaviour may not be visible. Datasets with m = 3 consist of more observations, therefore the FNN algorithm may have caught the occurence of nonlinear phenomena. There might be also a possibility that Hyperion remained in the state of temporary low-dimensional chaos. For possible rotational-lightcurve models see (Hicks et al. 2008).

Maximal Lyapunov Exponent
The mLE, denoted λ in general, was calculated using two distinct algorithms: the Wolf et al. (1985) and Kantz (1994) methods, incorporated in the programs lyapmax and lyapmaxk (Kodba et al. 2005;Perc 2005aPerc ,b, 2006, and lyap k from the TISEAN package.
Herein the results of those investigations are presented.
The time delay τ was calculated in Section 3.1. The embedding dimension m was not set according to the FNN results, but was varied from m = 2 to m = 10 and first the mLE was computed using the approach of Wolf et al. The algorithm finds a nearest neighbor to an initial point and evolves them both until the separation becomes too big; next, the distance is being rescaled in order to stay in the small-scale structure, and this repeats to the end of the time-series. Then, the average of the logarithms of the displacement ratios is the estimate of the mLE. All of the numerical results are gathered in Table 3 The exclusion leads to lessening the standard deviation of theλ max corresponding to C1 data and changing the sign of K89 and R1 datasets' mLE to negative with significant diminishing their standard deviation. Because, as mentioned, chaos cannot be present in a two-dimensional continuous phase-space, the results from Fig. 4b seem to be more realistic.
The FNN convergence to m = 2, as shown in Table 3 require the data to fulfill the stationarity assumption, i.e. that the statistical properties of a time-series (e.g., mean and standard deviation) are constant in time. We therefore perform a test in a following way (Perc 2006).

Stationarity test
Let us consider a point p(t) as an event and find all similar events, i.e. those points p(i) that lie not further than from p(t). We average all values of x i and call this a prediction of a future observation based on the value of x t . The key now is to use cross-prediction, i.e.
to partition the whole dataset into non-overlapping segments and use the j-th segment to make a prediction of a k-th segment. We quantify the correctness by calculating the error δ jk via square root of mean square deviations from the mean in segment k and repeat this for all j and k: wherex is the prediction and x is the true value in a k-th segment. If δ jk is significantly larger than the average, then either the dynamics are not conserved from one segment to another, or the data is undersampled. Both cases yield a conclusion that the data is non-stationary.
The program stationarity was applied to sampled lightcurves under consideration and the results are gathered in Fig. 7.
The immediate denouement is that the lightcurves from K89 and D02 are too short, undersampled, or both. Therefore, it is justified to ask a question: how long and how dense should photometric observations be in order to reveal a positive mLE in a lightcurve?

Numerical experiments
To answer the last question, we examine simulated lightcurves of Hyperion for chaotic and regular solutions of the Euler system of equations. These lightcurves, as well as the LEs spectra, were obtained from (Mel'nikov, priv. comm.) and were computed using a procedure described in D02. The algorithm gives Hyperion's stellar magnitude m in time (JD), corrected to zero solar phase angle and mean opposition magnitude, as well as the time evolution of the Euler angles (θ, ϕ, ψ) and the corresponding angular velocities (ω 1 , ω 2 , ω 3 ).
Parameters and initial conditions are presented in the Appendix. Full spectra of LEs were calculated using the HQRB method (von Bremen et al. 1997) realized as a software complex in (Shevchenko & Kouprianov 2002;Kouprianov & Shevchenko 2003). The described data are shown graphically in Fig. 8 together with the output of the stationarity test. The system is Hamiltonian, therefore the six LEs are paired so that λ i + λ j = 0 and the plots show only three positive LEs. The Lyapunov time T L for the chaotic solution is equal to 44 d.
Since the lightcurves have a constant time step ∆t = 0.1 d, in order to produce time-series more astronomically realistic only three first points during each day were left and averaged (see Fig. 8). Then a cubic spline and sampling were performed to produce datasets consisting of 5000 points. From these sampled lightcurves, intervals of lengths: 2 months, 6 months and 1 year were chosen randomly; each had ten realisations both for the chaotic and regular solution. The whole 3 year lightcurves were taken as single realisations.
Next, the routines false nearest and autocor from the TISEAN package were applied for obtaining the time delay τ and lyap k for embedding dimension from 2-10 to extract the mLE. In the same way time dependencies of dynamical variables (θ, ϕ, ψ, ω 1 , ω 2 , ω 3 ) were investigated. It may be surprising that there are clearly linear parts in the stretching factor S(n) plots in the regular case, however, a closer look at the stationarity tests show that in any of the S(n) plots, as expected, but due to non-stationarity of the sampled data, in some of the stretching factors a linear rise is more evident than for the simulated chaotic case. This proves that stationarity is a necessary condition for a dependable detection of chaotic phenomena.

Discussion and conlusions
The aim of this paper was to verify whether it is possible to infer the value of the mLE from photometric observations of Hyperion. Firstly, existing datasets (K89 and D02) were investigated using Takens phase-space reconstruction and its correlation dimension, stationarity tests and finally the mLE was estimated using two algorithms: Wolf et al. (1985) and using stretching factors (Kantz 1994). Wolf et al. method yielded a positive detection for K89 observations, nevertheless the S(n) plots showed no linear region, implying lack of chaoticity. As elaborated, the Wolf et al. method is likely to yield spurious detections, especially for short datasets. We therefore conclude that existing datasets are too short and undersampled to detect chaotic rotation using the mLE.
In order to list conditions allowing to obtain the mLE from potential ground-based observations, simulated lightcurves spanning 3 years (with a constant time step of 0.1 d) were examined. As suggested from previous considerations, two-month long subsamples appeared to be too short to yield a sign of chaos. Half-year data retained a clearly visible linear rise in some of the S(n) plots. On the other hand, to make these samples more atronomically realistic, only magnitudes spanning ≈ 7.5 h each day (more precisely, night observations) were left and averaged. This led to a conclusion that only one year subsets are long enough to reveal the presence of chaos in the stretching factor.
Additionally, a false detection of chaos was observed in the case of lightcurves based on regular solutions of the Euler equations. To explain this extraordinary behaviour a careful inspection of the stationarity test outputs was conducted. It was found that the time-series underlying the simulated lightcurve are non-stationary (as well as the lightcurve itself) which violates the assumptions underlying the mLE calculation algorithm. Therefore, obtaining a positive mLE is by itself not sufficient to claim detection. A necessary stationarity condition must be also fulfilled.
Based on computations described herein we assert than to reliably estimate presence of chaos in ground-based photometric observations of Hyperion via mLE, these observations should be performed over a time period of approximately one year. A way to shorten this period is to obtain well-sampled photometry, e.g. by observing with 2-3 telescopes spread over the world. As was noted, in case of data points distributed uniformly with a time step equal to 0.1 d, timespan may be shortened to half year. We remind that the resulting time-series should be stationary.
The author is grateful to Andrew Seymour for discussions and sharing the MATLAB code, and especially to Alexander Mel'nikov for helpful discussions and providing useful information as well as the simulated lighcurves.

A. Appendix
Embedding delays τ for phase-space reconstruction via Takens method were calculated as the delay at which the autocorrelation function drops to 1/e, leading to values being an order of magnitude greater than those from the MI algorithm. The reconstructions are displayed in Fig. 11, while correlation dimensions of these embeddings are uninsightful and therefore not presented herein. Note that all trajectories are characterized by correlation dimension not much greater than unity; furthermore, the 3D embeddings in both cases show no intersections, which is a premise that m = 3 is sufficient. Figs. 2 and 11 reveal intersections in 2D plots, what is naturally a projection effect.
Average values of mLEs computed with the Wolf et al. method are presented in Fig. 12 and in many instances give results contrary to those from Fig. 4 obtained using the MI for estimating time delays. Moreover, the mLE convergence plots frequently shows behaviour oscillating around zero, preventing to see a tendency for a certain sign. As argued in the main text, even when the Wolf et al. method shows convergence to a positive mLE, this is an ambiguous detection, which may be spurious due to an assumption of exponential divergence of initially nearby trajectories, not necessarily to be met in actual time-series. Table 4 gathers parameters and initial conditions necessary to run simulations of lightcurves as described in D02 and obtained from (Mel'nikov, priv. comm.).