Publishers imprint. Printed in Malaysia. A Preliminary Nonlinear Analysis of the

The Chandler wobble (CW) is a resonant response of the Earth rotational pole wandering around its figure axis whose excitation mechanism is still uncertain. It appears as polar motion oscillations with an average period of about 433 days and a slowly varying amplitude in the range (0–300) milliarcsec (mas). We here perform a nonlinear analysis of the CW via a time-delay coordinate embedding of its measured X and Y components and show that the CW can be interpreted as a low dimensional unstable deterministic process.


INTRODUCTION
The Earth rotation rate and spin axis position relative to the Earth's figure axis, experience fluctuations at all time scales.Variations in the length of day (LOD) and polar motion (PM) are observed in response to applied external (e.g.tidal) and internal (e.g.topographic coupling between Corresponding author.E-mail: fred@hpvlbi.obspm.fr.core and mantle) torques, to angular momentum exchanges between the solid and fluid parts of the Earth and to global mass redistributions primarily associated to geophysical fluid motions which per- turb the Earth inertia tensor (Lambeck, 1980).Now- adays, the high accuracy geodetic measurements of these fluctuations provide useful constraints on the time evolution of the Earth-Moon system and on the modelling of the earth's core and mantle as well as hydrosphere and atmosphere dynamics (Hide and Dickey, 1991;Dehant et al., 1997).
The Chandler wobble (CW) is a resonant response of the PM.The corresponding free oscillation was already predicted by Euler in 1765 with a period of about 306 solar days for an axi- symmetric flattened rigid Earth model.But it was only observed by Chandler in 1891 from data of the variation of latitudes, with a period of about 14 months ( 436 days).This increase of the period is well accounted for by the elastic yielding of the Earth.Several processes in the Earth's global dynamics being dissipative (e.g. the wind and oceanic current drag on the crust), the CW must be sustained by an excitation source which has resisted any reliable identification until today.
Several candidates have been proposed in the past but none of them appear to be appropriate.The energy input by earthquakes is probably between two or five orders of magnitude too small to excite the wobble (Souriau, 1986; Chao and  Gross, 1987).Other possible excitation mechanisms were searched for in pre-or post-seismic deforma- tions (Sabadini et al., 1984) and in convective processes of the core external layers (Le Mouel  et al., 1985) without reaching conclusive results.
Early investigations about the possible excitation by atmospheric or oceanic variabilities (Ooe, 1978;  Wahr, 1983; Chao, 1984) are somewhat supported by recent advances in global ocean-atmosphere modelling (Vondrfik, 1990; Chao, 1993) but with rather mitigated results.Significant departures between the predicted and observed wobble are still unexplained so that other mechanisms might at least partly sustain the Chandler oscillation (Eubanks, 1993).
Other studies have been dedicated to the char- acterization of the CW from the data series.Variations in amplitude and departures from a purely periodic motion are reported since a long time (see e.g.Jeffreys, 1940).Frequency modula- tions were also sought with diverging conclusions that should result from differences in the way data were analyzed and from data contaminations by a time-dependent noise.Indeed phase varia- tions of the CW could be an artifact of observa- tional errors or result from random unknown excitations (Okubo, 1982).Two or more free frequencies could also coexist near the Chandler frequency supposed constant and modulate the signal (Colombo and Shapiro, 1968).A single Chandler frequency might be amplitude-dependent due to nonlinear response to the oceans and change with time (Carter, 1981).
Major information on the CW energy sources and sinks can be gained from the estimation of its quality factor Qcw.This factor, which is inversely proportional to the ratio of the energy loss -AEcw in one wobble cycle to the energy Ecw of an unsustained wobble, can be retrieved from the width of the CW spectral line (Stacey, 1992).Recent estimates of Qcw range from 50 to 180 (Furuya  and Chao, 1996) that correspond to decays of an unmaintained CW with typical time scales of 20- 70 years respectively.
Until now almost all analyses aiming at the identification of excitation sources or at extracting CW parameters from the data were conceptually confined to the "stochastic process paradigm".
However the resonant response of a damped oscillator is stronger when the forcing has a rich spectral content about its free period.Following this idea, quasi-periodic angular momentum fluc- tuations in the atmosphere near the 14 month period were recently retrieved by Furuya et al. (1996).Unfortunately their study relies on linear analyses (Fourier spectra, coherence functions) though it is known that the dynamics of the various geophysical fluids (core, mantle, ocean, ground- water, atmosphere) are driven by nonlinear pro- cesses at sub-daily to climate time scales.
Results presented here suggest that the CW is sustained by a low dimensional nonlinear determin- istic process.The nonlinear data time series analysis of the rapid fluctuations of the LOD and PM in the frequency band (2-100) days already supported a similar conclusion (Frede and Mazzega, 1999a).In particular nonlinear signatures of El Niffo Southern Oscillation events (Harrison and Larkin,  1998) were recovered in the time-dependent local Lyapunov exponents, indicating that at these periods the unstable modes of the Earth's rotation are prompted by the ocean-atmosphere coupled dynamics (Frede and Mazzega, 1999b).
In the next section we present the data time series and explain how the CW signal is routinely separated from the annual wobble, trend and residues.In Section 3, the embedding De and local DL dimensions are determined from global and local neighbourhood analyses of the embedded data vectors.The Lyapunov exponents and dimensions are then computed from the data time series as explained Section 4 and the results discussed in Section 5. However this analysis should be con- sidered as preliminary.Indeed we do not tried here to test the robustness of these results.Similar computations should also be undertaken with PM data time series published by other Earth's rota- tion center.The algorithms and sensitivity tests (to noise contaminations, limited data series, etc.) are detailed in the context of Earth's rotation studies in the twin papers by Frede and Mazzega (1999a,b).
From our experience, we think that the results presented below are robust and ascertain the non- linear determinist origin of the CW excitation mechanism.
TIME SERIES OF THE CHANDLER  WOBBLE (1846-1997)   The CW is extracted from data series of the raw PM components X (along the Greenwich meridian) and Y (along the 90 E meridian) as edited by the International Earth Rotation Service (IERS, Observatory of Paris).The longest available series (the C01 series, see IERS, 1997) extends from 1846 up to 1997.During this one and half century of survey the techniques of PM observation have considerably evolved and the resulting C01 series is a concatenation or, when possible, an optimal combination of partial data series (Vicente and  Wilson, 1997).For example in the early period (1846-1890) the PM is deduced from astrometric optical measurements of absolute declinations with a frequency band covering almost only the Chandler and annual wobbles and an accuracy probably not better than a few tens of milliarcsec (mas).The latest two decades of data combine geodetic observations from Very Long Base Inter- ferometry (VLBI), Global Positioning System (GPS) and Satellite and Lunar Laser Ranging (SLR and LLR respectively, IERS, 1997) with an accuracy of 0.3 mas.The sampling rate is 10 data/ year from 1846 to 1889 and 20 data/year from 1890 up to now.
The PM is dominated by two main components: the annual wobble which is presumably excited by the ocean-atmosphere system and the CW.They can be quite well separated from each other with a 7 year-long record of PM as during such a time interval 6 full cycles of the CW are completed.The splitting of a time series into trend, seasonal and irregular components is a practical issue in financial time series analysis.The Census X-11 filter (Shiskin  et al., 1967) achieves this kind of decomposition on the basis of local trend estimators for non-seasonal time series using finite moving averages (Gray and  Thomson, 1997).It is routinely implemented in an iterative data analysis scheme at IERS in order to separate the annual and Chandler wobbles (see IERS, 1991, for details).
In a first step, the X and Y time series of the PM have been re-interpolated with a homogeneous sampling rate of 100 data/year all over the (1846-1997) interval with a cubic spline.Then we applied the iterative decomposition with the finite moving averages being estimated over moving 7-year sub- intervals in order to separate the annual and Chandler wobble without smothing too much the time variability of these signals.We retain a central period of 433 solar days for the CW in good agreement with the recent determinations of this parameter (see e.g.Vicente and Wilson, 1997).It should be noted that the frequency resolution of the moving averages over 7-year time span is quite low so that weak frequency modulations of the CW about its central frequency are implicitly permitted in the decomposition process.The annual wobble is about 60 mas RMS, the CW about 100 mas RMS and the irregular part about 40 mas RMS on both components of the PM (means and RMS are summarized in Table I).The statistics of the wobbles and irregular parts of the X and Y series are similar.The main difference between the two raw series is a consequence of the strong trend in the Y component of the PM (4mas/year since the beginning of this century).We will not be interested in this secular drift that might be caused by the postglacial rebound or changes in polar ice covering (McCarthy and Luzum, 1996).
The raw series, trend, wobbles and irregular part are drawn in Fig. for the Y component (similar curves are obtained with the X component of PM).
The quite high noise level in the years (1846-1890) is responsible for 40 mas RMS of the irregular (or residual) part of PM.It would rather be of the order of 10 mas RMS when considering the data obtained after 1890.The noise is also spoiling the annual and CW recovery before 1890 as can be seen on Fig. 1.Then the most striking phenomenon is the low amplitude of the CW in the years (1925)(1926)(1927)(1928)(1929)(1930)(1931)(1932)(1933)(1934)(1935)(1936)(1937)(1938)(1939)(1940).This strong modulation has been interpreted as a phase reversal of the wobble (see e.g.Vondrfik, 1990)  whose origin is still unknown (this modulation occurs on the X component as well).In fact the time variations of the excitation mechanism should be searched in these slow modulations of the resonant response of the Earth's rotation (and may be on related frequency modulations that are more difficult to track).
A detailed examination of the annual and Chandler wobbles also reveals that these two components might not be perfectly separated by the Census X-11 algorithm.This is not so surprising as we operate the decomposition over the smallest possible moving window of 7 years.In the com- putations described below, the first 10 years of data are skipped in order to limit the effect of the data noise on the results.Moreover a weak smoothing via local parabolae fitting to data included in a (moving) window of 150 days, is applied to the Chandler series before entering the nonlinear analysis.The instantaneous wobble is taken as the value of the fitted local parabola middle point.This filtering only reduces the X and Y CW series by 2mas RMS though removing most of the residual noise contaminations in the data series.

DIMENSIONAL ANALYSIS OF THE CHANDLER WOBBLE
Dimensional and dynamical invariants of a com- plex system can be extracted from the nonlinear analysis of an univariate data time series issued from that system.Data vectors are constructed in a DE,-dimensional phase space E* in a frame of time delayed coordinates.In E* the time series of data vectors draws an orbit whose statistical and geometrical properties, thanks to the embedding theorem (Takens, 1981), we know to be directly inherited from the unknown source dynamics (Eckmann and Ruelle, 1985).The deep theoretical background of this approach will not be addressed at all in this study.The interested reader will refer to the papers cited above or to the recent overview given by Abarbanel (1996).
Here the complex system is the Earth with its various interacting envelopes (solid parts, geophysi- cal fluids, etc.) that we study from the point of view

FIGURE
The raw time series of the Y-PM (upper panel) is splitted into a secular trend (scaled by a factor 2), the Chandler (period 433 days) and annual (period 365.25 days) wobbles and an irregular or residual part (lower panel).Units are in mas.Very similar plots are obtained with the X component of PM.
of its rotational dynamics.In the computations below, we consider separately the Y and Y com- ponents of the CW.Being extracted from the same dynamics they should lead us to the same estimates of the system invariants.The series of the CW being limited in length and contaminated by noise, the coherence of the results extracted from both series is a crude way to give some confidence in the characterizations so obtained.
The dimensional analysis starts with the determi- nation of the smallest dimension of the space E* required to unfold the reconstructed orbits.Note that this dimension Du, might depend on the data series extracted from the system and so is not an invariant of the dynamics.From the data time series at hand x(t) we construct/)trial-dimensional vectors y(t) using a time delay r by: Ix(t), + + In Erial, the time series of the vectors y(t) draws some parametric curve with a lot of neighbour points (or even intersections).These neighbour- hoods, which can be easily counted, result either from the system passing again similar states in phase space during its evolution or from undesir- able projections of the reconstructed orbit in a too low dimensional space.These "false" neighbours are progressively eliminated when increasing the trial embedding dimension/)trial Only true neigh- e* bours remain as soon as the smallest dimension required to completely unfold the orbit, say De,, is reached.Stochastic processes being infinite dimen- sional, their data contamination tends to give an overestimated embedding dimension, and even sometimes compromises its determination (Frede and Mazzega, 1999a).
In Eq. ( 1), the definition of a system of statistically independent coordinates relies on the choice of the time-delay r.In the context of nonlinear dynamics, it has been shown that the best way to estimate a reliable delay is to calculate the Average Mutual Information (AMI) function (Fraser and Swinney, 1986) which, shortly speak- ing, is based on the probability distribution of events in the time series.The joint probability distribution between the data series x(n) and its delayed version x(n + T), say P[x(n), x(n + T)] is obtained from the joint histogram in the [x(n), x(n + t)] variables (for simplicity we here assume that the data are given at discrete times n 1, N).
The AMI for any delay T is defined by: AMI(T) The first minimum of the AMI function is the sought delay.We obtain time delays of 105 and 115 days respectily for the X and Y components of the CW.It can be seen from the AMI function for I/ plotted in Fig. 2, that the minimum is surrounded by weak gradients.So the estimation of r is not very acute and could easily vary within + 10 days or so.
Anyway, the various orbits that would be recon- structed on the basis of slightly different time delays are diffeomorphic and equally lead to the sought invariant estimates.
With the optimal delay for a given series, we build data vectors as in Eq. (1) for a trial dimension/)trial For each vector, we count the number of neighbour vectors and compare it to the number of neighbours optimal delay: 115 days 8 o ;' ',;,' ',;' '2-5' ';  in dimension /)trial we, / 1.Those neighbourhood of dynamical origin are not affected while those resulting from a projection are gradually eliminated by the unfolding of the embedding space.In Fig. 3 we plot the global percentage of false neighbours as a function of the trial embedding dimension/)trial E, for the Y series with a delay r 115 days (once again similar plots are systematically obtained for the X series).Almost all the neighbours in dimen- sion are false neighbours (globally 100% of false neighbours) as the data are artificially pro- jected on a line.This percentage decreases to zero for a trial dimension of 4. At this stage it is difficult to know if the very low percentage in dimension 3 is not a remanent of data contamination by a weak noise.However, exactly the same behaviour being observed for the X component, we retain an embedding dimension De, 4 in the following.It can be shown that these results are not sensitive to the value of the radius used to define any neighbourhood as long as it is chosen as a small fraction of the data series RMS.
From the knowledge of the time-delay r and embedding dimension De,, the orbit is recon- structed in pseudo-phase space E*.A 3D projection of the Y-Chandler orbit is drawn in Fig. 4. It looks like interweaving "circles" more or less confined in a bounded disk-like sub-space.Much will be learnt about the dynamics of the CW from the geometry of such a reconstructed orbit.
The next dimension to be determined is the local dimension.Whereas the embedding dimension De, is a global dimension, the dynamical dimension DL is related to the exact number of degrees of freedom activated during the system evolution.It is an invariant of the system dynamics that can be also extracted from a time series of observations.The broad lines of the algorithm are summarized as follows (see Abarbanel and Kennel, 1993 for details): the orbit is reconstructed in a working space with high enough dimension Dw > De, to ascertain that true neighbours only are kept.We then consider a given point of the orbit, say y(n), select its Nb nearest neighbours [yi nghb (f/); i--1, Nb] and compute the corresponding covariance matrix.The eigenvectors of this matrix give the principal local directions along which are distributed the selected data vectors.Keeping only the D ial local directions, we build a maping which carries the nghb nghb (/7 / m) after m vectors v (n) to their images .J/i time steps of the ambient flow.This mapping is finally applied to the vector y(n) and the result Ypred(n / m) is compared to its true image y(n + m).
The prediction is considered successful if the euclidian distance Ily(n+m)-yprd(n +rn)l is below a predefined fraction of the data series RMS (typically 10%).These operations are repeated for all the points y(n) of the orbit and the percentage of bad predictions is estimated.All the predictions are successful only when the local mappings are built in the right dynamical dimension DL.In this case the performance of the prediction does not depend on the number of neighbours Nb nor on the success criterion.
The choice of Nb for a limited time series in the presence of noise contamination is discussed in Frede and Mazzega (1999b).Indeed we found this algorithm to be quite sensitive to noise because it basically works with neighbourhoods so that the local evolution of the flow is partially hidden by noisy features of the reconstructed orbit.The FIGURE 4 Y component of the CW reconstructed as an orbit in 4D pseudo-phase space E* (one of its 3D projection is drawn here).The delay for embedding is 115 days.The radius of the torus-like structure is about twice the data series RMS ( 100 mas).algorithm applied to the complete Chandler series spanning the period (1847-1997) gives no results: the percentage of bad predictions remains about 25-40% whatever the trial local dimension is.It converges to DI 3 when skipping the first 10 years of observations of the Y component (Fig. 5).Thirty years of data (1847-1876) have to be skipped for observing a convergence to the 0% level of bad prediction for the X component of the CW.Once again the corresponding local dimension is DL 3.
A few tens of years of very noisy observations (see the residues in Fig. 1) can spoil the whole determination of DL,,because such a data sub-series is whirling several times around the attractor draw in Fig. 4 and contaminates the decomposition of all the local data covariance matrices.Indeed in phase space E* the set of the nearest neighbours to some recent accurate data vector often includes old and  The 0% prediction error level reached in dimension DL 3 does not depend on the number of nearest neighbours fixed to build the local covariance matrices (see text) and indicated in the window.noisy measurements.Anyway, the removal of the too noisy early data ensures the convergence of the local dimension determination to the value DL 3 for both series.

STABILITY ANALYSIS OF THE CHANDLER WOBBLE
Information related to the stability of the observed CW can also be gained from the nonlinear data analysis.The first point is to understand what can be the fate of a small perturbation of the wobble orbit, say 6y: will it fade away, remain unchanged or grow with time?These three possibilities, encapsulated in the following formula: @(t) By(0)exp[At], (3) are clearly determined by the sign of the Lyapunov exponent A (here @(0) is the perturbation of the reconstructed orbit at an arbitrary time considered as "initial" and @(t) the perturbation at time t).Any infinitesimal perturbation taken along the direction locally tangent to the orbit, cannot grow nor fade with time as the perturbated state y / @ belongs to the same orbit and so actually corresponds to a near past or near future state of the system.As a consequence the Lyapunov exponent associated to this neutral direction is zero, a value that should be retrieved from the data (up to noise contamina- tion).The CW being modellized by a three degrees of freedom system (DL 3), two more directions must be considered that can be associated to stable (negative exponent) or unstable (positive expo- nent) manifolds.The wobble resulting from the interaction of dissipative processes, the sum of the exponents must he negative, a condition that ex- presses the overall contraction of volumes in phase space and the confinement of orbits in a bounded domain.
Without entering technical details, we briefly report how the Lyapunov spectrum can be ex- tracted from a data series.Consider the CW as a dynamical system observed at discrete times indexed by integers i: yi+ F(yi). (4) After K time steps, the perturbation 5; applied to y; is transformed by the composed action of the local jacobian matrices DF of the dynamics along the orbit, say: 6i+K DF(yi+K-1) DF(yi+K-2) DF(yi) (i DFK(yi) (i. It can be shown that the sought Lyapunov exponents are related to the eigenvalues of the Osedelec matrix S(yi, K) (see e.g.Eckmann and Ruelle, 1985) defined by: S(yi, K) [DFK(yi) T DFK(yi)] /(2K), (6) where the superscript T is for transpose.Besides those computational difficulties involved in the decompostion of S because of its intrinsic bad conditioning (see Abarbanel, 1996 for algorithmic details), the important point to be noted here is that the DF's can be estimated from the data time series, without prior knowledge of the dynamics in closed formulae.They are obtained from the partial derivatives of local maps Y Yv;i+l (7) developed on basis functions so to fit the time evolution of the embedded data vectors (yv;;+l is the vth neighbour of y; at time i+ 1).The global exponents are obtained from the Oseledec matrix which spans all the data series (starting from yl and K being the length of the vector series).Theoreti- cally the Lyapunov spectrum is invariant with regard to changes in the initial conditions, as long as the data series is infinite (K-+ oc).The possible effects of series finiteness and noise contaminations on the estimation and reliability of Lyapunov exponents are empirically studied and discussed in a similar context in Frede and Mazzega (1999b).We here just mention that our series seem to be long enough and not too much spoiled by noises so that significative estimates of the three exponents can be extracted.The coverage of the phase space sub- volumes visited by the CW system is dense enough and allows the proper computations to be robust.This statement will be further documented below when considering time-reversed data series.However, the genuine "global" Lyapunov exponents of the CW are definitely beyond the observational grasp.The exponents we actually compute are average values over the measurement time span.
The Lyapunov exponents for the X and Y components of the CW have been computed.For both series the principal exponent '1 is positive (-(X) 0.046 and ,-(Y) 0.086 respectively, see Table II).The exponent associated with the neutral manifold is smaller than 0.01 in absolute value for X and Y.These departures from the exact zero value result from the series finiteness, noise contaminations and numerical approximations.principal Lyapunov exponent Al +/--; sum of the Lyapunov exponents A; + Lyapunov dimension DLyap and horizon of the prediction Hp.The '+' and '-' refer to the forward and backward time series respectively.
However it is well below those estimates of the A1 exponents which are, from a computational point of view, better approximated.As expected, both third exponents are associated with stable mani- folds (-0.115 for X and -0.179 for Y), the wobble process being dissipative and the sum of the exponents being negative (Table II).
Apart from the exact value of ,1, the key information is about its sign.The presence of a positive exponent in the CW series would imply that nonlinear instabilities can grow on a Du= 1D manifold so putting an intrinsic limitation to the possibility to predict its time evolution over a finite time interval (the so-called horizon of prediction), whatever might be the nature and skill of the algorithm designed for this task.Our confidence in the positiveness of the ,)k exponents can be reinforced by considering now the CW series with reversed time (starting from 1997 and ending on 1857).Any direction expanding (contracting) for- ward in time will contract (expand) backward so that the corresponding exponent will change sign (note that sign corrections are applied in the discussion below and Table II for comparisons).
With time reversed, we find again the principal exponents to be positive but smaller than before ()q-(X) 0.032 and A-((Y) 0.018, the minus superscript indicating time reversal).The exponents associated with the neutral manifolds are ,(X)=-0.004and A-(Y)=-0.001and the sums of the exponents always negative as it should be.Moreover these characterizations of the Xand Y components of the CW processed independently from each other are coherent.
All these results strongly suggest that the CW is an unstable (or equivalently "chaotic") low dimen- sional deterministic process.From our past experi- ence, we know that the nonlinear analysis of a finite time series dominated by a stochastic component sometimes converges to a finite embedding De, and (much less convincingly) local DL dimensions.
But until now we found not a single example of such a series providing forward and backward negative exponent sums aud sign-invariant princi- pal exponent.
On the basis of the average principal exponents (,1-[,i --Xi-]/2), we estimate the horizon of prediction Hp for both series by: Hp C. ln(50)/A. ( This somewhat arbitrary definition corresponds to the time necessary for a 1% RMS error on the system state to grow up to 50% RMS (for the CW a mas error will grow to about 50 mas) along the unstable manifold.The constant C in ( 8) is required to scale the horizon in days.We obtain 368 days and 276 days for the X and Y wobble components respectively.If these values are to be confirmed, this would mean that reliable predictions of the CW cannot be given beyond about one year.The observed departure between the horizons for X and Y probably reflects uncertainties in the esti- mates of the principal Lyapunov exponents.It might also result from some noise or physical process being predominantly "polarized" on the Y component of the wobble.Such a possibility is not unrealistic as for example we noticed before that the trend component of PM is preferentially polarized on that component (see Section 2).
Another fundamental invariant of the CW dynamics, the Lyapunov dimension DLyap is extracted from the data series.It provides a characterization about the way the orbit (as drawn in Fig. 4) fulfills the phase space.It is conjectured that DLyap equals the attractor dimension D A.
However we found its estimation to be much more robust than the one of DA based on the correlation integral (Frede and Mazzega, 1999b).The Lyapu- nov dimension is define from the exponent spec- trum by (Kaplan and Yorke, 1979): DLyap K + -JK'-I Aj with }-_ ,j > 0 and jK__, ,j < 0. A striking result is that the estimates of DLyap is very robust when considering the forward and backward time series though the Lyapunov spectra exhibit serious (but not sign affecting) departures as reported above.
Indeed we find DLyap+(el() 2.34 to be compared to Dyap(X 2.52 and, for the Y series, DLyap+ (Y)   2.43 to be compared to D[yap(Y --2.41.The average dimensions of the two series are even in better agreement with DLyap(X)=2.43 and DLyap(Y)=2.42 (these are the values quoted in Table II).It should be noted that these fractal values do not characterize the orbits themselves (which are of course differentiable) but the way the bounded domain of phase space is fulfilled by the system state whirling on the attractor.Such non- integer dimension is often associated with chaotic time series (strange attractor).

DISCUSSION
At this point of our study what have we learned?Probably the most promising break-through is about the usefulness of considering the CW as a deterministically driven process of low dimension.
The assumption of a stochastic origin has always been at the fundament of the various analyses of wobble observations.Such a partiality is well explained from an historical perspective: first, it is only recently that a set of nonlinear concepts and analysis tools are available for an alternative approach; second, the complex behaviours of the forcing geophysical fluids were themselves tradi- tionally recasted in the paradigm of stochastic processes.
The nonlinear data analysis by phase space reconstruction also provides us with qualitative and quantitative results.As illustrations of the first one, we have shown the CW to evolve in time with a 1D unstable manifold embedded in its dynamics and that any model of the wobble given in form of a system of coupled ordinary differential equations requires three degrees of freedom (DL=3).A quantitative bound of about year is found beyond which any attempt to reliably predict the CW is doomed to failure.
Have we identified the forcing mechanism of the CW?Not yet but this could be hardly done from the analysis of PM observations only.However the various invariants we have extracted from the wobble data series are related to those of the sought forcing.Indeed the dynamics of a linear dissipative pendulum (modellized by Liouville equations, see Lambeck, 1980) does not exhibit any chaotic time evolution.Starting from a non-equilibrium state, the orbit is spiraling during the transcient response, down to the fixed point equilibrium.The corre- sponding local and Lyapunov dimensions are D ans-2 and Fltrans_ 0 (no unstable manifold) Lyap fP for the transient and D 0 and DLyap 0 for the fixed point.All departures from these values observed in the CW are inherited from the time evolution of the forcing(s).
A similar comment might apply as well to the quality factor Q. Is the Q factor estimated from Earth Orientation Parameters relevant to the way the solid Earth relaxes external energy inputs or is it a measure of the dissipative processes at work in the forcing dynamics?We have no direct answer to this question.However the rate of volume shrinking in phase space to a lower dimensional structure (the attractor set) as estimated from the sum of the Lyapunov exponents should give some clue.In particular the relashionship between the working definition of the quality factor Q and the exponent sum should be investigated.
A further step can be done with the nonlinear data analysis.Usually the identification of a forcing field of Earth rotation is done by comparing the corresponding time series (projected in terms of angular momentum or inertia tensor changes) with the measured variation of the relevant rotation parameter.A complementary approach would consist in comparing the variation in stability of both series.In the previous section we derived the global Lyapunov exponents from the Osedelec matrix S(y, K) built with all the jacobians of the time series (see Eq. ( 6)).In the same way, local exponents can be computed from local Osedelec matrices starting at sample y; with k jacobians (k integer).The exponents so obtained give a measure of the average stability (or, more precisely, unstabil- ity if we consider a positive exponent) over the k following time steps in the data series.Time varia- tions of the stability with state evolution is a well Year FIGURE 6 Time series of the local principal Lyapunov exponent for the X (upper panel) and Y (lower panel) components of the CW.The plotted exponents are averaged over running windows of 6 and 24 years.High picks in the series correspond to relative losses of stability and reduced horizons of prediction (note that the global exponents are 0.046 and 0.086 for the X and Y series respectively).
known feature of dynamical systems (Abarbanel  et al., 1991).In Fig. 6 are plotted the time series of the local principal exponents for the X and Y components of the CW.In each case, we computed the average exponents over 6 year and 24 year intervals.The choice for a 6 year window is related to the choice of the same moving window used to separate the annual and CW (Section 2).Siginficant time variations in stability can probably not be retrieved with a higher resolution.The 24 year window is within the accepted range of dissipation rate (deduced from the quality factor; see Section 1) at the CW period.When increasing the window length, the variations in the local exponent are progressively smoothed out.We notice first that the principal Lyapunov exponents are always positive and exhibit large variations, up to a multiplicative factor 3 when compared to the global exponents of Table II.Then the corresponding horizons of prediction are divided by the same factor (see Eq. (8)).For example, in the decades 1860-1870 and 1940-1950, predictions of the future state of the CW attempted from past observations were intrin- sically limited to 4 or 6 months rather than year as stated in the previous Section.Stability variations are even more sensible in the Y component of the wobble.The 1920-1930 decade of weak wobble (see Fig. for Y) is also quite stable (small exponents) though still chaotic (positive exponents).
Once again we stress on the fact that the observed stability variations are inherited from the unknown forcing of the CW.The time series of a candidate forcing should not only match as much as possible the observed wobble but it should also present the same time evolution of the local Lyapunov exponent.Such a comparison will bring a strong constraint for the identification of the Chandler forcing(s).Such analyses revealed quite success- ful in recognizing the ocean-atmosphere system coupled in strong El Nio events as a major source of instability in the rapid fluctuations of the Earth rotation (Frede and Mazzega, 1999b).
The various results presented in this study are derived from the data series of the International Earth Rotation Service (IERS, Paris).Similar computations should be performed with PM series estimated by other institutes and the results com- pared.These "challenging" series cannot be statis- tically independent because only a few historical records going back to the middle of the previous century exist.However the data selection, pre- processing and combination can be performed in different ways and lead to more or less diverging results.The global characteristics (time delay, embedding and local dimensions, positivity of the principal Lyapunov exponent) extracted from the IERS data will probably not be affected by such a change in the reference series.But the local characteristics, in particular the local exponents, are likely to change significantly.This kind of result robustness should be checked before proposing detailed interpretations of the time variations of the CW stability.

CONCLUSION
We give strong evidences that the CW is driven by a deterministic low dimensional process.The X and Y components of the wobble extracted from PM data series, are embedded in multi-dimensional pseudo-phase spaces with a system of delayed coordinates.The optimal time delay as deduced from the average mutual information function is 105 and 115 days respectively.The embedding dimension is De, 4. Both series can be locally approximated by a three degree of freedom system (local dimension DI 3).Then the CW is shown to be a chaotic process with a Du= 1D unstable manifold.The Lyapunov dimension is probably in the range 2.40-2.50.The nonlinear analysis of both components of the wobble converges to the same kind of conclusions, even when considering the data series with reversed time.The average horizon of prediction is about year.Nevertheless strong fluctuations in the wobble stability can be seen from the time series of the local Lyapunov exponents.These stability crises are inherited from the unknown wobble excitation source and might contribute to its identification.
FIGURE 2 AMI as a function of the time-delay T. The delay r for data embedding corresponds to the first minimum of the AMI function.

FIGURE 3
FIGURE 3 Global percentage of false neighbours as a function of the trial embedding dimension /)trial for the Y-Chandler data time series.The 0% in dimension De, 4 corresponds to the complete unfolding of the reconstructed orbit in pseudo-phase space E*.

FIGURE 5
FIGURE 5 Percentage of bad predictions as a function of the trial local dimension D ial for the Y-Chandler data time series.

TABLE Mean and
RMS of the raw series of the X-Y components of the PM and of their respective trend, CW, annual wobble and irregular part (unit: mas)