Fluctuations of finite-time stability exponents in the standard map and the detection of small islands

Fluctuations of finite-time stability exponents in the standard map and the detection of small islands

Steven Tomsovic111Martin Gutzwiller Fellow, 2006/2007; permanent address: Department of Physics and Astronomy, Washington State University, Pullman, WA 99164-2814    Arul Lakshminarayan222Permanent address: Department of Physics, Indian Institute of Technology Madras, Chennai, 600036, India. Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Strae 38, D-01187 Dresden, Germany
July 20, 2019

Some statistical properties of finite-time stability exponents in the standard map can be estimated analytically. The mean exponent averaged over the entire phase space behaves quite differently from all the other cumulants. Whereas the mean carries information about the strength of the interaction, and only indirect information about dynamical correlations, the higher cumulants carry information about dynamical correlations and essentially no information about the interaction strength. In particular, the variance and higher cumulants of the exponent are very sensitive to dynamical correlations and easily detect the presence of very small islands of regular motion via their anomalous time-scalings. The average of the stability matrix’ inverse trace is even more sensitive to the presence of small islands and has a seemingly fractal behavior in the standard map parameter. The usual accelerator modes and the small islands created through double saddle node bifurcations, which come halfway between the positions in interaction strength of the usual accelerator modes, are clearly visible in the variance, whose time scaling is capable of detecting the presence of islands as small as of the phase space. We study these quantities with a local approximation to the trace of the stability matrix which significantly simplifies the numerical calculations as well as allows for generalization of these methods to higher dimensions. We also discuss the nature of this local approximation in some detail.

05.45.-a, 05.45.Ac, 45.10.-b,45.05.+x

I Introduction

A central and widely used measure of instability in dynamical systems is the Lyapunov exponent Lichtenberg and Lieberman (1992); Ott (1997). If a positive Lyapunov exponent is found for a continuous family of trajectories, then it is taken as a measure of their dynamical chaos. It is natural that the Lyapunov exponent, which is defined as a time-approaching-infinity limiting value, when calculated numerically, would be approximated by considering finite time segments of the infinite history of a trajectory. This has led to the study of the so-called finite time or local Lyapunov exponents (FTLE) Ott (1997) in a variety of contexts. For example, FTLE in one-dimensional maps Fujisaka (1983); Theiler and Smith (1995); Prasad and Ramaswamy (1999); Anteneodo (2004), Hamiltonian systems Sepúlveda et al. (1989); Amitrano and Berry (1992, 1993), advection and turbulent flows Beigie et al. (1993); Lapeyre (2002); Arratia and Gollub (2005), form just a part of a rather large literature. Additionally in most experimental situations one encounters finite time quantities. Such FTLE have been the subject of numerous studies and have been experimentally measured in many contexts such as for instance in the study of passively advected particles Arratia and Gollub (2005). There is also a very intimate relationship between the theory, methods, and measures found in the subject of Anderson localizationCrisanti et al. (1993) and FTLE. Although, typically dynamical correlations are assumed vanishing or if correlations are accounted forRechester and White (1980), then they are not the kind which is found in Hamiltonian flows where all the complexity of boundaries between chaotic and regular motion may exist. In other words these correlations are then calculated assuming that a typical trajectory covers phase space uniformly, whereas the presence of islands, however small, ensures that this is not true.

A very closely related, but less directly studied, quantity is the stability exponent which is defined via the eigenvalues of the stability matrix. A detailed study of this measure vis-a-vis the Lyapunov exponent was carried out in Goldhirsch et al. (1987) where it was conjectured that the Lyapunov and stability exponents were equal under a broad range of dynamical assumptions. It turns out that some of the studies that purported to be about the FTLE have actually been studies of the stability exponents instead. In semiclassical mechanics, such as in the Gutzwiller trace formula Gutzwiller (1990), one of the natural classical quantities that appears is, in fact, the stability exponent and not the Lyapunov exponent. In higher dimensional systems (larger than two-degrees-of-freedom), it is the sum of the positive stability exponents that appears and in the infinite-time limit this tends to the Kolmogorov-Sinai (K-S) entropy. Since it has long been appreciated that the small wavelength limit and the large time limit in quantum mechanics are not interchangeable, we may expect that for fixed wavelength, it is the finite-time stability exponents, to be defined below more precisely, which may be relevant; see for example Silvestrov et al. (2003). In fact, many of the interesting ways in which the fluctuations in the stability exponents might potentially show up in quantum mechanics and semiclassical theories remains to be investigated. This forms our primary motivation for studying the stability exponents, their distribution and their scaling properties. Nevertheless, in this paper we restrict our attention to purely classical and low-dimensional systems (two-dimensional area preserving maps) to establish basic properties and approximations, as well as to demonstrate their usefulness in detecting small stable dynamical structures. Naturally there are many parallels in the large literature on FTLE both of a quantitative and qualitative character given that for a broad class of dynamical systems the infinite-time stability exponents are identical to the Lyapunov exponent. Our results presented here generally for stability exponents are thus of relevance to the literature on the FTLE as well.

In the following section we first introduce the model in use throughout (the standard map) and an approximation for the stability exponent which provides a practical, calculational simplification that can be generalized to higher dimensions. Its ergodic value has corrections to the well-known formula of Chirikov Chirikov (1979). In section III, the local properties of the stability exponents are studied. We image the finite-time stability exponents in phase space and show their close relationship to the stable and unstable manifolds of the map, and show their connection to locating small non-hyperbolic regions. In section IV, the fluctuations of the stability exponents are investigated, and a formula for the ergodic average of the variance is derived. The scaling of the variance is shown to be a sensitive sensor for small islands. It is also shown that quantities that incorporate higher cumulants and that are of natural semiclassical relevance are even better sensors of such phase-space structures. In the appendices, formulas are collected for an expansion of the exact trace of the stability matrix as well as results for the time scaling of the stability exponent higher cumulants.

Ii Background, an approximation, and ergodic average

In this section we introduce the kicked rotor, provide definitions of the quantities of interest, explore an approximation of great practical utility, and use it to generate an improved analytic estimate of the Lyapunov exponent for the standard map.

ii.1 The Hamiltonian

Consider the 2-D standard map Lichtenberg and Lieberman (1992); Ott (1997) as a well-suited vehicle for our investigations. The methods and conclusions we reach will have wider applicability as this map represents the generic, local dynamics of any integrable system with shear in action-angle coordinates and shows how that structure gets altered by perturbation. A general kicked rotor is a mechanical-type particle constrained to move on a ring that is kicked instantaneously every multiple of a unit time, . Supposing the radius of the ring to be and , the Hamiltonian takes the form


where is a function periodic on the interval . From , mapping equations can be given:


where we have made the choice to apply the potential kick before the free motion. The notation indicates the derivative of with respect to . We have also restricted the phase space to the unit torus () for convenience; this is of no consequence for the results presented in this paper. The simplest periodic function on a ring is just the lowest harmonic and leads to the standard map.


Many results have long been known for the standard mapChirikov (1979). At the map is integrable and is essentially a stroboscopic map of a freely rotating particle. There are both rational and irrational tori, depending on whether the frequency of rotation is commensurate or not with the frequency of the strobe. As the kicking strength is increased from zero with exactly the same frequency of stroboscopic observation, the incommensurate tori (irrational) survive small perturbations in accordance with the KAM theorem, while the commensurate (rational) ones break up into a pair of stable and unstable orbits in accordance with the Poincare -Birkhoff theorem. The phase space becomes mixed with stable and chaotic orbits for increasing . At around the last rotational irrational KAM tori breaks and this leads to global diffusion. Up to a stable fixed point persists. Beyond the standard map is considered to be largely chaotic, although it is also not proven to be completely chaotic for any value of as far as we know. At an infinity of values of stable fixed points are known to appear in the line for the map on the torus that are accelerator modes for the map on the cylinder Lichtenberg and Lieberman (1992). These typically occupy regions in the phase space whose areas scale as Chirikov (1979). Otherwise, the Lyapunov exponent of the map, to a very good approximation (Chirikov’s result), increases with as .

ii.2 The finite time stability exponents

For the map of Eq. (II.1), the stability matrix at time is given by


where the are local one-step stability matrices. The finite time Lyapunov exponents are defined via the singular values of , which are the eigenvalues of and turn out to be non-negative real numbers. The eigenvalues of stability matrices on the other hand can in general be complex. To be technically correct we refer to the positive finite-time “stability exponents” (FTSE), , where the time is the integer number of applications of the mapping equations and is the initial condition. Consider only initial conditions for which the orbits are unstable, in which case for area preserving maps such as we are presently studying the eigenvalues have to be real. The can be expressed in terms of the trace of the stability matrix as follows


where denotes the trace operation. The approximation applies for large enough. In that case, the approximation error is exponentially small in the time except for stable trajectories. The generalization to higher dimensions is through the exponential growth rate of the magnitudes of the eigenvalues of the stability matrix (the Jacobian) at time .

The natural quantity, which occurs in semiclassical expressions, is whose large time behavior is determined by the sum of the positive stability exponents. For 2D area-preserving maps one has , and the maximum stability exponent coincides with the K-S entropy in the limit. Thus for the kicked rotor, studying , , or the K-S entropy is equivalent and our interest is in these quantities’ (or their logarithms) finite-time fluctuations over phase space. Therefore our study is naturally close to the well-explored field of FTLE. There are strong reasons to believe that in the infinite time limit the definitions of FTLE and FTSE are equivalent in the chaotic cases that we are concerned with, although non-generic examples are known for which this is not true Goldhirsch et al. (1987). Since we are dealing with finite time measures we wish to be technically correct and refer to FTSE, although these typically deviate from the FTLE by terms of the order of . Recall that the FTLE can also be defined via the way initial error vectors grow in length:


This definition of FTLE has the added complication that they will depend on the initial error vector as well. However, in the infinite time limit this corresponds to the exponent coming from the singular values of the limit of . Thus, is the Lyapunov exponent and it turns out not to have initial condition dependence or fluctuations. There are reliable methods of evaluating each of these quantities Lichtenberg and Lieberman (1992); Greene and Kim (1987) (whether FTLE or FTSE) and they give essentially the same information, especially regarding the fluctuations of the exponents that we are primarily concerned with here. Higher dimensional generalizations, which are currently under study Tomsovic and Lakshminarayan (2007), involve the fluctuations of and therefore we prefer to concentrate on FTSE in this paper on 2D maps.

ii.3 Approximating the trace

Before turning to a more detailed analysis of , its average, variance and higher order moments and cumulants, we investigate approximations to the trace. These approximations make the trace calculations both efficient and simple, and extend to higher-dimensional systems, where they provide significant advantages. For any stability matrix of the form of Eq. (4), it turns out that can be expressed in an exact, closed form of products of the individual traces ; the full expressions to all orders are found in Appendix A. We display here the zeroth order approximation, indicated by , and its first order correction, indicated by :


These expressions represent well the large- behavior of the exact trace , even though, a priori, one might suspect that estimating this trace by Eq. (7) is rather inadequate. After all, consider the fact that the leading correction, shown in Eq. (8), contains terms of inverse products of pseudo-random factors which have the possibility of nearly vanishing. Crudely speaking, as multiple trajectories are considered or as increases for an individual trajectory, it becomes more and more likely that a nearly vanishing factor will be encountered thus creating a correction term that dominates the leading term. Worse, all the successive corrections involve combinations of the same products of inverse factors (Appendix A). It is possible that no particular correction term dominates either.

Fortunately and oddly enough, Eq. (7) (the term) is an excellent approximation, and in fact, is more than adequate for analytically estimating the leading asymptotic behaviors of the finite time stability exponents. It is even good enough to capture the power-laws for both local averages and variances of FTSE, which we investigate ahead.

Figure 1: Comparison of the two approximations, Eqs. (7,8), for the trace of with the exact value. In panel (a), the three curves are for K=8 with the initial condition , , the exact trace is the solid line, the approximation is the dotted line (lowest), and the approximation is the dashed line. In panel (b), . For the solid line, is evaluated to the approximation, and for the dashed line, is evaluated to the approximation. Similarly to panel (b), in panel (d), the ratio of the two approximations is given. Panel (c) is similar to panel (b), except that with the initial condition (which is a stable trajectory). The solid line is for the approximation whereas the lower curve is for the approximation.

Figure (1) compares for a single trajectory the natural logarithm of the exact trace, the product form estimation of Eq. (7), and the estimation of Eq. (8) where the mean exponential increase has been subtracted from all three in the first panel. In these panels several general features are readily seen. As a function of time, the leading approximation follows the exact trace time step by time step except for sudden vertical displacements at time steps where the factor nearly vanishes. In this example, the corrected approximation follows the exact trace longer in time before departing; although, this does not have to be the case. A second point is that the exponential growth dwarfs the errors of the approximation. Even though the approximation of Eq. (7) is off by a factor of roughly before , this is miniscule compared with the average exponential growth of the trace (in this example it is and so it alters the stability exponent very little. From the form of the leading correction, a Levy-flight nature might be expected for the corrections. This can be seen in the natural logarithm of the ratio of Eqs. (7,8) respectively where only a couple of jumps dominate the differences between the zeroth and first approximation. Although not shown, as K increases, a corresponding curve seems to be rather “quiet” punctuated by relatively fewer jumps.

These approximations naturally are worst for segments of trajectories that are stable or nearly stable. For them is repeatedly close to zero. The large value of thus does not order corrections. An example is given of a nearly stable trajectory in a region near a double saddle point bifurcation. The value of is larger than used previously to illustrate the behavior of a highly unstable orbit. Nevertheless, the approximation is far worse. On the other hand however, it gives a stability (or Lyapunov) exponent of about 0.067, a small number compared with . As long as the percentage of stable and nearly-stable orbits is very small, the errors in using Eq. (7) for all the orbits regardless of their nature are small after averaging over initial conditions. For the usual Lyapunov exponent there exists a “local” expression which is in some sense the equivalent of the expression for the trace we have just examined Fujisaka (1983).

ii.4 The average finite-time stability exponents

Much of what is known applies readily if a product form can be found for . This is precisely what the kicked rotor has for sufficiently strong potentials and times not too long. There the term would be sufficient, and from Eqs. (5,7), the are approximately given by


The superscript indicates that it is the zeroth order approximation. Assuming that the system is fully chaotic, the expectation value is given by integrating (averaging) over all initial conditions and


For the standard map, the result is ; the integral is evaluated exactly. Not surprisingly, this coincides exactly with the Lyapunov exponent approximation made by Chirikov Chirikov (1979) even though this result is for the stability exponent. The phase-space mean of the finite-time stability exponent is independent of integration time and time is dropped from its notation. This indicates that the mean is to within approximations made so far equal to . In contrast, the entire phase space average of the FTLE is not equal to its infinite time ergodic average, but has a transient that decays as a power law in time and in the absence of significant islands of regular motion goes as .

We plot the mean stability exponent in Fig. (2).

Figure 2: The exact mean stability exponent calculated as a function of kicking strength with the zeroth order approximation, , subtracted. Other than small oscillations, the difference is seen to be almost exactly equal to the first correction, , given in Eq. (16). This is plotted as the dashed curve .

It is interesting that there is a simple systematic correction to the well-known formula of Chirikov above that is of the order of , and that other than the value of (excluding sample size noise), there is essentially no indication of any other kind of dynamical information.

The origin of the corrections lie in the terms with in the approximations to the exact trace, however its rigorous derivation seems fraught with difficulties. A somewhat heuristic derivation of the -correction to Chirikov’s expression can be derived from Eq. (8). The product term gives an independent contribution and is exactly as already noted. Subtracting that term, one is left with the ergodic average of at least the term to determine the correction;


We are not aware of a way to evaluate this integral without making a few simplifications. Consider first that the correction is essentially being expanded about or in other words the expansion is about . In this limit at some fixed value , the probability is vanishing that more than one factor of for a given trajectory is close enough to zero that it can meaningfully contribute. Let denote the value for the minimum of the set of . Although that involves two terms in the series, only one will dominate. However, since there are -terms in the cyclic sum, there are ways this can happen; together this is consistent with what happens in the subject of extreme statistics. The basic expression to be evaluated is


Due to the momentum not entering the expression, can be integrated as two independent variables each uniformly distributed on the interval . The next step follows from the discontinuity in the expression for the integral


It is straightforward to verify that to the extent that the upper form of the integral applies, the expression for vanishes. Therefore it turns out that only values of for which


contribute a correction to the average stability exponent and the limits of integration on are restricted to that domain (actually there are two such domains). Switching variables to , one has (accounting for both domains)


If the leading square root factor is expanded in a series of powers of , the integrals can be analytically evaluated. This gives,


As seen in Fig. (2), the first term accounts completely for the overall correction to the Chirikov result to the accuracy level of the numerical calculations. The trace correction appears to have given properly the stability (Lyapunov) exponent to . This is interesting because it seems that it is not necessary to consider the and higher corrections to the trace whereas one might have naively expected the trace correction to give an contribution. Chirikov’s accurate estimate for the Lyapunov exponent also has the same corrections derived above. This may be expected as at infinite times the Lyapunov and stability exponents indeed converge for the standard map. However, numerical calculations of the FTLE, even averaged over the entire phase space, mask this with transients that are of the order of .

Iii Local properties of finite-time stability exponents

The FTSE is also a local stability exponent as emphasized by its dependence on initial conditions. The mean over a local “area” that is much smaller than the entire phase-space () includes transients that are power-laws in time. Its interesting variations over phase space can be displayed by assuming the following particular form:


The area is ideally a very small local region of the available phase space centered at and the final two terms are responsible for removing , thus leaving just the fluctuating part. Here (in practice ) and the time is larger than the logarithmically short mixing time scale, . This time scale arises as the “log-time” during which the initial region of volume through propagation essentially elongates to the extent that it enters each distinct local region of area once on average. The averaging in the LHS of the above expression is over the uniform Lebesgue measure, which is the invariant measure of the map. Note that when the entire phase-space, i.e. the unit square for the kicked rotor, is used for averaging the FTSE, the transients in time seem to be exponentially small and essentially .

Expressing the finite time stability exponents in this way allows us also to draw attention to anomalous features that may be present for various values of , especially those with small, but significant, stable islands in phase space that survive for large , such as the accelerator modes. The quantity fluctuates both in time and phase space point around which the averaging is done, it can also have a strong dependency on . The index serves to remove the principal non-fluctuating time scaling part of the FTSE and determines the rate at which it approaches the ergodic average. For example, in the case when there are no apparent small islands (see text ahead), while in the presence of stable regions . The area serves as a region of coarse-graining and it must be emphasized that unless we choose to be the entire phase space (), the regions will exclude regular islands and trajectories. We present results for an exponent similar to when studying the variance below. Previous studies of the FTLE have also presented such forms of convergence to the infinite time Lyapunov exponents Grassberger et al. (1988); Abarbanel et al. (1991); Goldhirsch et al. (1987), but typically in dissipative systems.

To get an idea of the complexity of the phase space distribution of the FTSE as we move around, always averaging over a small area of phase space centered at this point, we display density plots in Fig. (3).

Figure 3: A density plot of the average (top, bottom respectively) as a function of . High to low density is indicated as a color transition from purple to red. The area over which the averages are calculated are squares of area . The case (center column) has no apparent stable islands, whereas at (left column) there are two small islands each of period four where the density is enhanced and purple (see the next figure for a magnification of the dynamical structure underlying the purple zone). Note the prominent linear structures for a large value of =40 in the right column.

The figures show forward and backward in time evolution for three cases of , namely , , and respectively. The intricate structures that arise highlight areas of reduced hyperbolicity and hence smaller FTSE. For the case of forward evolution in time, trajectories which are trapped in more stable areas or propagate into those areas have lowest stability and appear purple. Those trajectories propagating into more stable zones must do so in a direction aligned along the greatest exponential compression. This direction is determined by the stable manifolds of the short periodic orbits. Thus, one sees a structure that mirrors the foliation of phase space by these stable manifolds in the top row of the figure. Similarly, for reverse evolution in time, one sees a structure that mirrors the corresponding unstable manifolds in the lower row. In the first column, the phase space has small islands, while the center column, , does not have such prominent islands. The absence of a detectable island of regular motion results in a less contrasted variation of the FTSE. At (right column), which is very strongly chaotic, there are linear regions of low FTSE that are clearly visible. However, any other regions have more or less merged into the background. These linear regions coincide with the independent borders of the non-hyperbolic regions for a one-step forward (the vertical lines at and ) and one-step backward iteration (the lines and ). The fact that these lines are not dense in the phase space and survive for long times at large underlies the coexistence of strongly chaotic dynamics with the fact that the map may not possess a single value for which it becomes mathematically fully hyperbolic. Similar plots for the FTLE have been plotted earlier in the contexts of area-preserving maps for smaller values of Eckhardt and Yao (1993) and passive scalar advection Lapeyre (2002).

It is possible to use this plot to help locate small islands. Investigating the purple zones seen in the left column of Fig. (3) more closely reveals the interesting structure found in Fig. (4).

Figure 4: Expanded plot of two islands of period-four regular motion and the stable and unstable manifolds associated with two unstable orbits born at the same double saddle node bifurcation (K=9.20667769750) as the islands themselves. Only one of this structure’s four images is plotted. Transport in and out of the inner region is severely limited by the small turnstile fluxes visible.

It arises for values greater than approximately at which point there is a double saddle node bifurcation. As increases from the bifurcation, the structure grows to a maximum size near . The full structure includes not only the small stable islands, but also the short unstable periodic orbits along with their homoclinic tangles, i.e. crossings of stable and unstable manifolds. The turnstiles formed by these crossings limit the rate trajectories can move in and out of the region surrounding the islands of stable motion MacKay et al. (1987). Therefore in this structure, it is not just the immediate interface between the chaotic and regular motion which is generating trapping, although that, of course, is happening additionally just beyond the outermost tori pictured. In fact, chaotic trajectories may remain trapped inside the tangle for hundreds or more time steps in the roughly 2/3 portion outside the islands themselves. For this reason, the entire structure was visible in Fig. (3), not just the islands, and multiple trapping time scales would be expected. One of the more unusual features is a highly asymmetric instability (factor of 1000 or so) in two of the stable and unstable manifolds pictured. The interior middle line is actually a stable and unstable manifold whose turnstile is so tiny as to appear as a single line dividing the interior, thus making transport between the upper and lower region nearly non-existent.

Iv Fluctuations of the finite-time stability exponent

The statistical properties of can be studied via the construction of a probability density , i.e. the probability of finding in the interval centered at . For a large class of chaotic systems, a great deal is known about general properties of , although strictly speaking just for the FTLE. The approach to a unique Lyapunov exponent is connected to approaching a -function limit as . Through the study of the density’s cumulants for fully hyperbolic systems, one sees that there is a slow approach to an ever narrower, roughly Gaussian density as Ott (1997). A detailed study with analytic results has also been performed for “white-noise” random media Schomerus and Titov (2002). Interestingly enough, it is also known that this limiting process is too slow to understand the behavior of the generating function for the cumulants


which is also a quantity of direct interest in semiclassical studies or the theory of Anderson localization. This relation should be thought of as an expansion in the real variable .

If one were to calculate the full probability density for the kicked rotor, one would see that it qualitatively takes on one of two forms. If the kicking strength is set to a value at which there is no discernible evidence for any stable orbits (no detectable non-exponentially decaying dynamical correlations), the density has a well defined peak near the Lyapunov exponent, and is not too far off a Gaussian form. If the kicking strength is set to a value at which both stable and strongly chaotic trajectories co-exist in the phase space, then it takes on a bi-modal form with a peak near the mean Lyapunov or stability exponent and another representing algebraic behavior near zero Sepúlveda et al. (1989). The density is continuous between the two peaks. Nevertheless, as , a unique Lyapunov exponent emerges albeit much more slowly. This has a profound effect on the time scaling properties of the density or more specifically its cumulants (moments) with the exception of its mean. The algebraic peak and the ‘neck’ connecting the two peaks disappears in the limit very differently than the rate at which the width of the main peak narrows. Time scaling therefore gives a great deal of information about non-hyperbolicities.

iv.1 Time scaling of the variance

We proceed next to a more detailed analysis of the variance (second cumulant) by first deriving the contributions that neglect dynamical correlations, which leads to an explicit formula for the standard map; the first cumulant (the mean) is already treated above. It is shown that while the general time scaling of the variance is , this is significantly modified by the presence of even very small stable islands, thus providing a practical way of detecting a ‘barely’ mixed phase space. We anticipate that this would be especially useful in its higher-dimensional generalizations Tomsovic and Lakshminarayan (2007). A first study appears to support our expectations Beims et al. (2007).

The variance of the FTSE over can be expressed simply using the approximation of Eq. (7) for the exact trace; the exact one, , differs from its approximation above, , by very small amounts which can be ignored. One finds:


The variance may be be separated into a term that is explicitly independent of dynamical correlations and a series of terms that contain the contributions of dynamical correlations, which respectively are:


The contribution of the first term above is


Here we have introduced as a term that does not carry correlation information.

Depending on the system, strong behavior connected to dynamical correlations but not reflected in the mean, is possible in the variance, which the summation above makes clear. Note that the dynamical correlation function is bounded for a broad class of potentials. Furthermore, as long as the correlations decay faster than , they cannot alter the basic time-dependence of the variance. If they decay more slowly, they could well become the leading contribution in time. As seen in Fig. (5) the variance is strongly dependent on the value of the kicking strength and shows marked peaks at various values.

Figure 5: The order calculation of the variance as a function of kicking strength. It is multiplied by and averaged over . It illustrates that the variance is very sensitive to dynamical correlations. The dashed line at is the value of the integral from Eq. (21).

These include conventional accelerator modes as well as other stable islands, such as those illustrated in Fig.(4) that appear near .

iv.2 The variance as a sensor for small elliptic islands

As just noted, the variance scales with time differently from the law above as a result of correlations decaying as slow or more slowly than . For the kicked rotor, consider ensembles of trajectories that originate from small regions of the phase space rather than the whole space; the reduced averaging region is motivated by the interest in a process that more practically generalizes to systems with greater numbers of degrees of freedom. The average behavior may obscure the scaling by the presence of the infinite-time asymptotic value, whereas the fluctuations die out asymptotically and are more stable for evaluating power laws. Although, note that “wiggles” in the FTLE were attributed to stability islands in an earlier work on kicked tops Constantoudis and Theodorakopoulos (1997). The variance can be written


where . After accounting properly for the dependence, one sees that fluctuates as a function of time similarly to a diffusive variable. Generally speaking, the greater the time range over which is calculated and the greater the number of trajectories run, the greater precision with which is determined. In this regard, the simplicity of the approximation is quite helpful because it makes it possible to increase both in a practical sense. One suspects that the smaller and less influential the island of regular motion, the closer approaches unity and thus -accuracy is directly linked to how small an island one can detect in this way. Earlier works such as Grassberger and Kantz (1985); Horita et al. (1990) have studied scaling laws that arise in the fluctuations of the FTLE and have already emphasized the role of sticky regions and islands. Ahead, we show islands of measure roughly of that were found with the FTSE time scalings relatively easily.

In Fig. (6) we show the scaling of the FTSE variance for a few values of the kicking parameter compared with .

Figure 6: The scaling of the variance for an ensemble of initial conditions uniformly distributed in the square .

The regular island structures of Fig.(4) at have a very significant effect on the time scaling. One also sees the appearance of more than one time scale, i.e. slope, depending on the time regime. The turnstiles previously pictured would at first lead to only small numbers of trajectories entering into the interior region of the homoclinic tangle, and one would expect to begin close to unity. Later, a certain proportion would get trapped inside for up to hundreds or thousands of time steps. During this period, the greatest proportion would resemble nearly stable orbits and would deviate most from unity. Further in time, only the trapping right at the regular island-chaos interface would contribute to deviations from unity and the slope increases somewhat again. It is known that different trapping mechanisms can co-exist in phase space and lead to a “multifractal” process Zaslavsky (2002). This particular example is instructive because the turnstile time scale is fairly well separated from the others. The smaller the stable islands, the more the power laws recover the uncorrelated form. At there are only small detectable deviations from unity, and in both cases it is possible to locate very small islands; see Fig. (7).

Figure 7: Small regular islands found at and their surrounding homoclinic tangles. Not only are the islands much smaller than for , but the turnstile fluxes are larger relative to the measure of the enclosed area, making the flow much less restricted in the neighborhoods of the islands. The measure of the regular region for is roughly of counting all the unique iterations of the structure.

At we could not locate islands in this way and the exponent is as close to unity as the precision of the calculation.

It is worthwhile remarking that while accelerator modes give rise to anomalous diffusion (in momentum) in the standard map on the cylinder Ishizaki et al. (1991), the different scaling regimes discussed above hold whether the stable regions are accelerator modes or not. However, the suppression of the rate of the variance’s approach to zero and the average Lyapunov or stability exponent to the ergodic average in the presence of small islands has similar origins. Namely, they are due to trajectories spending intermittently long intervals of time near sticky islands or regions of almost marginal stability and nonhyperbolicity.

We point out that although power-laws such as that given in Eq. (22) are susceptible to variations in the initial position , to changes of the trajectory sampling, and to the time window in which they are under study, the deviations from unity are robust structures and can be correlated to phase space structures.

Figure 8: The exponent in the decay of the variance of the FTSE as a function of the kicking strength . The ensemble is the same as that used in Fig. (6). A region around is shown at a resolution of in .

In Fig. (8) we show the variation of as is tuned through a range which has small stable regions. The “normal ” exponent is realized to a very good approximation for most values of indicating the general absence of significantly large, low instability traps. The prominent deviation around is due to the structure of Fig. (4) born out of a double saddle node bifurcation around and it persists to near . It grows to its largest size around where the variance decays the most slowly, . In terms of the phase space area, the islands occupy approximately of at and at (for which ). Notice robust, though small, deviations from unity at several other values of . For example, even for an island as small as that found at (, there is a range of values where is significantly and consistently different from unity; the inset shows a finer scan of this region. The fine oscillatory structure that is apparent here could be reflecting the fractal behavior of such exponents themselves which in turn are due to bifurcations in the map that produce or destroy very tiny islands almost continuously. The Newhouse phenomenon wherein homoclinic tangencies produce an infinite number of sinks is replaced in conservative scenarios by a large number of elliptic islands Duarte (1994). Knowing that there is possibly an island somewhere in phase space does not locate it, and a blind search could well be futile. Some direction is important. The islands are found more easily by searching around regions of the lowest stability exponent, regions that are highlighted in figures such as Fig. (3), and correspond to zones of low hyperbolicity. Note one last observation, for the three islands constructed, , the degree of ’s deviation from unity and the measure of the islands are in corresponding order, smallest to largest.

iv.3 The higher order cumulants

Similar more complicated dynamical correlation expressions enter for all the higher cumulants. They naturally weight more the extremes of the probability density, and thus, one might expect them to be even more sensitive to stable islands. To be useful, they also need a standard time scaling for their correlation independent components as found for the variance. In fact they do have a simple time dependence and it involves generalizations of . By defining:


the correlation-independent component of the higher order cumulants are expressible as polynomials in these quantities with a single overall power-law decay in time and no other parametric dependencies. Details are presented in Appendix B, but note here an interesting consequence of these expressions that imply a special feature found in any of the 2D kicked rotors. Substitution gives


which shows that to the extent that dynamical correlations can be ignored, every cumulant is independent of the overall kicking strength. The probability density thus has a single functional form for all kicking strengths. Therefore, although the Lyapunov exponent grows logarithmically with , the probability density of FTSE retains the same shape and width (it just translates with as increases). Assuming the kicking strength is large enough that the system is strongly chaotic, any dependence of the probability density, other than where it is centered, is due to variation of dynamical correlations and oscillatory.

Even stronger dynamical correlations and more sharply defined fluctuations show up in the higher cumulants. In Fig. (9), the fourth and second cumulants are compared.

Figure 9: The fourth and second cumulants as a function of using the approximation. The solid line is the fourth cumulant and the dashed line is the second.

The fourth cumulant displays much sharper features than the second. On the other hand, the features are present in the same regions of , so the connection to stable islands is the same as for the variance. We did not carry out the corresponding time scaling calculations as we have no reason to doubt that they reflect the stable structures for the same reasons as for the variance.

iv.4 The fluctuations in

Since higher cumulants carry information and have increased sensitivity to stable structures, a quantity that contains information from all the cumulants has the potential to be roughly optimal for its sensitivity to stable structures. Through relations such as Eq. (18), one sees that fits into this category. Furthermore, it is a quantity of considerable interest in semiclassical theory, where the square root occurs as the amplitude along classical paths. We therefore define


In a uniformly hyperbolic system, is the Lyapunov exponent. However, if the hyperbolicity is not uniform, as is generically the case, then is less than the corresponding stability exponent. In the presence of stable motion, it must tend toward vanishing. The averaging above can be done over any ensemble of phase space trajectories over time ; we use a uniform density over . Shown in Fig. (10) is the variation of this quantity

Figure 10: The quantity as a function of at . The averaging involved initial conditions spread uniformly across phase space. The solid line is the exact trace. The relatively smoothly behaving dashed line is what the result would be for a normal variable, i.e. . The dotted line derives from the approximation. In principle, the approximation’s errors are much worse for a quantity such as relative to one of the low order cumulants, .

across a range of values. Its extreme sensitivity to small islands and sticky regions of phase space is quite apparent. This may be expected as the denominator can come arbitrarily close to zero in the neighborhood of points undergoing bifurcations. We show two calculations in the figure, one of which is an “exact” evaluation based on finding the trace of at large times, using a determinant method used initially for periodic orbits Bountis and Helleman (1981); Greene (1979), and the other is by using the product () approximation in Eq.(7). The latter procedure is a stringent test for this approximation and it is interesting that it picks out almost all of the “hotspots”. It is however much faster and easier to implement than the exact method. The smooth line shown is what one may expect from a normally distributed variable.

On comparison with and with the variance itself, we find that the ranges of where deviates from that expected of a normal random variable, is also reflected in these quantities. Note the consistency in all these measures with respect to structures due to stable, if small, islands in the phase space. For example, the structures around reflect those already illustrated in Fig. (8), while those around correspond to similar islands. We again expect that higher dimensional generalizations of the approximation studied here will provide a method of ascertaining higher dimensional, hard to detect, ‘weakly’ mixed phase spaces.

V Discussion

The consideration of finite-time stability exponents has some advantages vis-a-vis finite-time Lyapunov exponents. For one, they are more directly related to quantities that appear naturally in quantum mechanics. For another, their generalization to many-degree-of-freedom systems tend to the Kolmogorov-Sinai entropy in the limit. There is a local approximation just as there is for FTLE, but the time-dependent transient may disappear exponentially fast in many cases for FTSE and not FTLE, depending on averaging techniques employed. Corrections to the local approximation are also relatively simple and the first one is sufficient to calculate improved estimates of Lyapunov exponents. We showed that this gives a leading correction to of exactly for the standard map. In fact, although not explicitly mentioned, we explored a wide variety of techniques all of which failed to give this correction properly (most vanished or diverged). The corrections themselves are interesting as well. One would guess that they should lead to a Levy-process and numerically this is borne out. At least for low-order cumulants in the large kicking strength regime, the local approximation along with the first correction gives the full picture, and for many quantities the local approximation is more than sufficient.

For the standard map, the mean (first cumulant) FTSE converges rapidly to the Lyapunov exponent. This is understood to be a consequence of tending exponentially quickly to zero with time in Eq. (17). The local approximation gives a mean FTSE that is independent of direct dynamical correlation contributions. Indirectly, in the mostly chaotic regime, the exclusion of chaotic trajectories from small regular regions (islands) would slightly alter the ergodic averaging, but never more than the relative measure of those regions. Consistent with this picture, other than small deviations, one sees very little effect of small islands on the mean FTSE.

The fluctuations behave quite the opposite. They can be investigated through the creation of a probability density for the FTSE. In the absence of small regular regions, the probability density takes on a single form, independent of system parameters (such as the kicking strength), which slowly shifts towards a Gaussian form as of width . Therefore, all deviations from this form must be due to dynamical correlations, and in particular, correlations introduced by the existence of small stable islands of motion. Since very strong dependence is seen, the higher cumulants are very sensitive sensors of the presence of such non-hyperbolicity in the system. As such, they can be used as ‘detectors’. An important quantity which shows up in semiclassical theories, , in principle, contains information carried in all of the cumulants. In fact, it may be the most sensitive detector of dynamical correlations induced by tiny regions of non-hyperbolicity in the dynamics. All of the techniques and main results discussed have their generalizations to many degree-of-freedom systems. Those studies are currently underway Tomsovic and Lakshminarayan (2007).

Appendix A

The relation is


where is the integer part of and the odd arrangements of the product of factors is specified just ahead. Also note that if is even, the final term (the one for which ) is .

There are two key properties required for the above identity. First, there must exist at least one column or row in the set of both of whose elements are the same constants for all members of the set. Second, the determinants of the set of must all be equal to the same value (if this value is not equal to unity, then the terms in the equation above must be multiplied by its power) .

It turns out there is a fairly simple rule for which products can contribute to a particular “-term”. Imagine the integers uniformly spaced on a ring. Those are the “locations” that each factor in a product can occupy, so-to-speak, determined by its own subscript. For a given time with factors in the product, each arrangement that has exclusively odd-unit-spacings between occupied locations contributes once to that -term. All cyclically connected arrangements are automatically included in this construction, and thus the cyclic invariance of the trace is easy to see. The number of such odd spacing arrangements turns out to be


which is related to the summation and recursion relation


The combinatorials are related to an expansion of the functional form of eigenvalues raised to various powers of a two-by-two matrix. The term is simply the product of all factors


One could also label this term indicating that there are unit spacings. Similarly, the term is always the cyclic sum of which there are terms. Thus,


where the last term has indices . Pulling the product of all factors out front emphasizes that if the are generally much larger than , then for times not too long the terms can be thought of as corrections, each successive one roughly proportional to weaker that the previous term. In some sense, the first non-trivial -term comes for and . One expects double products of factors. The allowed odd-spacing arrangements are and . There are terms related to the cyclic permutations of the product () and only terms for the product () since it repeats itself after half a cycle, and thus , the correct count of terms, is consistent with Eq. (27). The above expression is the expansion of the determinant given by Bountis and Helleman Bountis and Helleman (1981) nearly 30 years ago. Greene developed its small- expansion in his study of the behavior of KAM (Kolmogorov-Arnol’d-Moser) surfaces Greene (1979).

Appendix B

It is more convenient to focus on the cumulants than the central moments as it simplifies the time dependence. There is also the advantage of interpretation of the probability density. For a Gaussian density, only the first two cumulants are non-zero. For other densities, the higher cumulants describe shape deviations. There is a third advantage, which comes into play when considering the relationship with quantum mechanics, but we do not elaborate here. If the dynamical correlations vanish or are small enough to be ignored, then the expression for the cumulants can be simplified. With these quantities, the cumulants have the same polynomial form, , as with the central moments and a very simple dependence on .

Figure 11: The first four cumulants as a function of time. The cumulants are multiplied by the time factors needed to make them independent of time according to the expressions leaving out dynamical correlations. The kicking strength is 8.9633781407, which gives the prediction that the Lyapunov exponent (first cumulant) is equal to .

The result is


For example,


An amusing, known consequence of the time dependence in Eq. (31) is that there is a sense in which the distribution collapses to a -function density as as expected (every trajectory has the same Lyapunov exponent), and the time-dependence of the reduced cumulants goes as


which for as . In other words, the density slowly approaches a Gaussian form. Nevertheless, the -function limit is approached slowly enough that the width of the probability density for grows without bound as increases, and it is not well approximated by a lognormal density. An example of the time dependence of the first four cumulants for a kicked rotor can be seen in Fig. (11). It illustrates that the time dependence is to within fluctuations as predicted for values at which dynamical correlations are not observed.


The authors gratefully acknowledge important discussions with Denis Ullmo, Alfredo Ozorio de Almeida, and Shmuel Fishman. They also gratefully acknowledge Roland Ketzmerick for finding kicking strength of the double saddle-node bifurcation, and S. T. gratefully acknowledges support from the US National Science Foundation grant PHY-0555301.


  • Lichtenberg and Lieberman (1992) A. J. Lichtenberg and M. A. Lieberman, Regular and Chaotic Dynamics (Springer, New York, 1992).
  • Ott (1997) E. Ott, Chaos in Dynamical Systems (Cambridge University Press, New York, 1997).
  • Fujisaka (1983) H. Fujisaka, Prog. Theor. Phys. 70, 1264 (1983).
  • Theiler and Smith (1995) J. Theiler and L. A. Smith, Phys. Rev. E 51, 3738 (1995).
  • Prasad and Ramaswamy (1999) A. Prasad and R. Ramaswamy, Phys. Rev. E 60, 2761 (1999).
  • Anteneodo (2004) C. Anteneodo, Phys. Rev. E 69, 016207 (2004).
  • Sepúlveda et al. (1989) M. A. Sepúlveda, R. Badii, and E. Pollak, Phys. Rev. Lett. 63, 1226 (1989).
  • Amitrano and Berry (1992) C. Amitrano and R. S. Berry, Phys. Rev. Lett. 68, 729 (1992).
  • Amitrano and Berry (1993) C. Amitrano and R. S. Berry, Phys. Rev. E 47, 3158 (1993).
  • Beigie et al. (1993) D. Beigie, A. Leonard, and S. Wiggins, Phys. Rev. Lett. 70, 275 (1993).
  • Lapeyre (2002) G. Lapeyre, Chaos 12, 688 (2002).
  • Arratia and Gollub (2005) P. E. Arratia and J. P. Gollub, J. Stat. Phys. 121, 805 (2005).
  • Crisanti et al. (1993) A. Crisanti, G. Paladin, and A. Vulpiani, Products of Random Matrices in Statistical Physics (Springer-Verlag, Berlin, 1993).
  • Rechester and White (1980) A. B. Rechester and R. B. White, Phys. Rev. Lett. 44, 1586 (1980).
  • Goldhirsch et al. (1987) I. Goldhirsch, P.-L. Sulem, and S. A. Orszag, Physica 27D, 311 (1987).
  • Gutzwiller (1990) M. C. Gutzwiller, Chaos in Classical and Quantum Mechanics (Springer-Verlag, New York, 1990).
  • Silvestrov et al. (2003) P. G. Silvestrov, J. Tworzydlo, and C. W. J. Beenakker, Phys. Rev. E 67, 025204 (R) (2003).
  • Chirikov (1979) B. V. Chirikov, Phys. Rep. 52, 263 (1979).
  • Greene and Kim (1987) J. M. Greene and J.-S. Kim, Physica 24D, 213 (1987).
  • Tomsovic and Lakshminarayan (2007) S. Tomsovic and A. Lakshminarayan, in preparation (2007).
  • Grassberger et al. (1988) P. Grassberger, R. Badii, and A. Politi, J. Stat. Phys. 51, 135 (1988).
  • Abarbanel et al. (1991) H. D. I. Abarbanel, R. Brown, and M. B. Kennel, J. Nonlinear. Sci. 1, 175 (1991).
  • Eckhardt and Yao (1993) B. Eckhardt and D. Yao, Physica 65D, 100 (1993).
  • MacKay et al. (1987) R. S. MacKay, J. D. Meiss, and I. C. Percival, Physica 27D, 1 (1987).
  • Schomerus and Titov (2002) H. Schomerus and M. Titov, Phys. Rev. E 66, 066207 (2002).
  • Beims et al. (2007) M. W. Beims, C. Manchein, and J.-M. Rost, Phys. Rev. E (2007), submitted.
  • Constantoudis and Theodorakopoulos (1997) V. Constantoudis and N. Theodorakopoulos, Phys. Rev. E 56, 5189 (1997).
  • Grassberger and Kantz (1985) P. Grassberger and H. Kantz, Phys. Lett. 113A, 167 (1985).
  • Horita et al. (1990) T. Horita, H. Hata, R. Ishizaki, and H. Mori, Prog. Theor. Phys. 83, 1065 (1990).
  • Zaslavsky (2002) G. M. Zaslavsky, Phys. Rep. 371, 461 (2002).
  • Ishizaki et al. (1991) R. Ishizaki, T. Horita, T. Kobayashi, and H. Mori, Prog. Theor. Phys. 85, 1013 (1991).
  • Duarte (1994) P. Duarte, Ann. Inst. Henri Poincare 11, 359 (1994).
  • Bountis and Helleman (1981) T. Bountis and R. H. G. Helleman, J. Math. Phys. 22, 1867 (1981).
  • Greene (1979) J. M. Greene, J. Math. Phys. 20, 1183 (1979).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description