Lyapunov Exponents

The previous article ended with a number it never properly defined. Sensitive dependence says neighbouring trajectories diverge like $e^{\lambda t}$ — but which neighbours, in which direction, measured how, and why should the wildly fluctuating local stretching of a chaotic orbit average out to any well-defined rate at all? This article builds the answer: the theory of Lyapunov exponents, the spectrum of numbers that has become the working currency of chaos — the quantity researchers actually compute, report, and argue about when they claim a system is chaotic. The right mental picture is an infinitesimal ball of initial conditions riding the flow. The dynamics shears it into an ellipsoid, and each principal axis grows or shrinks at its own asymptotic exponential rate: one rate per phase-space dimension, ordered $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_n$ . That full spectrum, not just its famous leading member, is a portrait of the dynamics. Its sign pattern classifies the attractor at a glance — $(-,-,-)$ a fixed point, $(0,-,-)$ a limit cycle, $(+,0,-)$ chaos. Its sum equals the average phase-volume divergence, tying it to the dissipation bookkeeping of the last three articles. In Hamiltonian systems it is exactly symmetric, $\pm$ pairs, symplectic geometry made numerical. And weighted partial sums of it estimate the fractal dimension of the attractor the dynamics builds — the Kaplan–Yorke bridge to the next article. Two theorems and one algorithm make all this respectable: Oseledets' multiplicative ergodic theorem, which guarantees the limits exist; the symplectic pairing rule, which disciplines the Hamiltonian case; and the Benettin grow-log-renormalize recipe, which turns the definitions into measurements. By the end, "this system is chaotic" becomes a statement with error bars.

Key takeaways

The Lyapunov spectrum $\lambda_1 \geq \cdots \geq \lambda_n$ is the set of asymptotic exponential growth rates of the principal axes of an infinitesimal phase-space ball carried by the flow; $\lambda_1 > 0$ on a bounded set defines chaos.
Perturbations obey the variational equation $\dot{\delta\mathbf{x}} = D\mathbf{f}(\mathbf{x}(t))\,\delta\mathbf{x}$ ; the propagator $M(t)$ accumulates the stretching, and Oseledets' theorem guarantees the limiting rates exist and are trajectory-generic.
Every bounded non-equilibrium orbit of an autonomous flow has one exponent exactly zero — the direction along the flow.
Sum rule: $\sum_i \lambda_i = \langle\nabla\cdot\mathbf{f}\rangle$ — negative for dissipative flows (attractors), exactly zero for Hamiltonian ones, where the spectrum pairs as $\lambda_i = -\lambda_{n+1-i}$ .
The spectrum's sign signature classifies attractors: $(-,-,-)$ point, $(0,-,-)$ cycle, $(0,0,-)$ torus, $(+,0,-)$ strange attractor.
Measurement: Benettin's grow-log-renormalize algorithm (with Gram–Schmidt for the full spectrum) in simulations; delay-embedding methods (Wolf, Rosenstein) for experimental time series.
Key idea: an infinitesimal ball of initial conditions is deformed by the flow into an ellipsoid whose principal axes grow like $e^{\lambda_i t}$ ; the ordered rates $\lambda_1 \geq \cdots \geq \lambda_n$ — existence guaranteed by Oseledets' theorem — classify the dynamics, sum to the mean volume divergence, and certify chaos when $\lambda_1 > 0$ on a bounded invariant set.

The core insight

An infinitesimal ball of initial conditions is deformed by the flow into an ellipsoid whose principal axes grow like $e^{\lambda_i t}$ ; the ordered rates $\lambda_1 \geq \cdots \geq \lambda_n$ classify the dynamics, sum to the mean phase-volume divergence, and certify chaos when $\lambda_1 > 0$ on a bounded invariant set.

Where this lives today

Lyapunov analysis is now instrumentation across an absurd range of scales. Solar-system dynamicists compute exponent maps of asteroid belt regions to chart chaotic transport routes (the Kirkwood gaps are Lyapunov-positive zones); accelerator physicists tune storage-ring lattices by tracking tangent dynamics so that stored beams sit in Lyapunov-quiet regions of phase space; and climate-adjacent fluid dynamics has upgraded from single exponents to covariant Lyapunov vectors — the intrinsic directions attached to each exponent — to dissect which physical structures carry the instability in turbulence. The same mathematics surfaced, half unexpectedly, in machine learning: a recurrent neural network processing a sequence is an iterated map, its exploding/vanishing gradient problem is literally a positive/negative Lyapunov exponent of the training dynamics, and architectures like LSTMs are engineering solutions for parking a system near $\lambda \approx 0$ — the edge of chaos as a design specification. Within the chapter, this article is the metrology station. The sensitive-dependence article supplied the concept; here it became a spectrum with theorems (Oseledets, pairing) and an algorithm (Benettin) attached; and the Kaplan–Yorke formula has already leaked the next article's subject — the fractal geometry of the sets that chaotic dissipative dynamics builds. The handoff is direct: the spectrum says the Lorenz attractor has dimension $2.06$ ; the strange attractors article shows what a $2.06$ -dimensional object is.

The physical picture

Drop a tiny ball of uncertainty — all the states compatible with one imperfect measurement — into phase space and let the dynamics carry it. Locally the flow is linear (that is what the variational equation says), and a linear map takes balls to ellipsoids. Watch the ellipsoid over a long stretch of orbit: each principal axis settles into its own exponential habit, growing or shrinking like $e^{\lambda_i t}$ . The rates, ordered, are the Lyapunov spectrum. The longest axis is the star — its rate $\lambda_1$ is what a generic error grows at, since any realistic uncertainty has some component along it — but the supporting cast matters just as much: the zero exponent along the flow direction, the negative exponents that fold uncertainty onto the attractor, and the sum of all of them, which must reproduce the volume budget the flow's divergence dictates.

Figure

A ball of uncertainty becomes an ellipsoid — one exponent per principal axis

The spectrum is a fingerprint. A stable fixed point contracts every axis: all exponents negative. A limit cycle carries one marker of time-translation freedom — perturbations along the orbit neither grow nor die — so its signature is $(0, -, \ldots)$ . A quasiperiodic torus has two neutral directions, $(0, 0, -, \ldots)$ . And chaos wears the signature this chapter has been building toward: $(+, 0, -)$ in three dimensions — stretch, slide, squeeze — with the sum negative for a dissipative system (a strange attractor of zero volume) or exactly zero for a Hamiltonian one (a chaotic sea, exponents in $\pm$ pairs). Read the signs and you have read the phase portrait. The honest subtlety is that no single glance at the Jacobian reveals any of this. Along a chaotic orbit the local stretching fluctuates violently — the trajectory passes through regions that compress the very direction that elsewhere expands — so the exponent must be defined as an infinite-time average, and it is a genuine theorem (Oseledets, 1968) that the average exists and comes out the same for almost every starting point in an ergodic component. That theorem is what licenses the entire enterprise: it turns "chaotic-looking" into a measurable number attached to the system, not to the whim of one trajectory. Measured across a family of systems, the exponent becomes a narrative device. Below is $\lambda$ for the logistic map as a function of $r$ — the entire story of the last five articles compressed into one computed curve. Through the periodic cascade the exponent is negative, plunging toward $-\infty$ at each superstable orbit and touching zero at every doubling threshold (where the multiplier sits at $-1$ and perturbations neither grow nor decay). At the Feigenbaum accumulation point it arrives at zero with nowhere left to dip — and crosses. Beyond, the positive plateau is chaos, gashed by the negative dips of the periodic windows, and at $r = 4$ the exponent reaches its ceiling $\ln 2$ : one binary digit of the initial condition consumed per iteration.

Figure

The Lyapunov exponent of the logistic map versus r, computed — negative through the cascade, zero at each doubling, positive past the accumulation point, dipping in every window

And it is measurable in practice, by an algorithm of almost embarrassing simplicity. Track one tangent vector alongside the trajectory; it aligns itself with the most-expanding direction (everything does); every so often, record how much it grew and rescale it to unit length before the computer overflows. The average of the logged growth rates converges to $\lambda_1$ . Carry a full orthonormal frame instead, re-orthogonalizing as you go, and the successive frame vectors — each forced to live orthogonal to the faster ones — report the entire spectrum. That is Benettin's method, and it is how every Lyapunov number quoted in this chapter was actually obtained.

Figure

The Benettin recipe — grow, log, renormalize, repeat: how Lyapunov spectra are actually measured

The formal core

Everything begins where the previous article ended: linearize about a moving trajectory. Perturbations obey the variational equation

\dot{\delta\mathbf{x}}(t) = D\mathbf{f}\big(\mathbf{x}(t)\big)\,\delta\mathbf{x}(t) \qquad (1)

where: $D\mathbf{f}$ — the Jacobian along the reference orbit, a time-dependent matrix even for autonomous systems. Its solutions are organized by the propagator (fundamental matrix)

\delta\mathbf{x}(t) = M(t)\,\delta\mathbf{x}(0), \qquad \dot M = D\mathbf{f}(\mathbf{x}(t))\,M, \quad M(0) = \mathbb{I} \qquad (2)

where: $M(t)$ — the time-ordered accumulation of all the local stretchings and rotations the orbit has visited; because the Jacobians at different times do not commute, $M$ is not the exponential of anything simple, and its long-time structure is a nontrivial question. The exponent of a chosen initial direction is

\lambda(\delta\mathbf{x}_0) = \lim_{t\to\infty}\frac1t \ln\frac{\lVert M(t)\,\delta\mathbf{x}_0\rVert}{\lVert\delta\mathbf{x}_0\rVert} \qquad (3)

where: the limit is independent of the norm chosen — a factor bounded above and below dies under $\tfrac1t\ln$ . What makes (3) more than wishful notation is Oseledets' multiplicative ergodic theorem: for an ergodic invariant measure, the limit

\Lambda = \lim_{t\to\infty}\big(M(t)^{\mathsf T} M(t)\big)^{1/2t} \qquad (4)

where: $M^\top M$ — the symmetric matrix whose eigenvalues are the squared singular values of $M$ , i.e. the squared stretch factors of the ellipsoid axes — exists, and its eigenvalue logarithms are the Lyapunov exponents $\lambda_1 \geq \cdots \geq \lambda_n$ , the same for almost every initial condition. A generic direction picks up the top exponent, $\lambda(\delta\mathbf{x}_0) = \lambda_1$ almost surely: components along slower directions are outrun by factors $e^{(\lambda_1 - \lambda_i)t}$ , so only measure-zero conspiracies see anything but the maximum. Two exact structural rules discipline the spectrum. First, the sum rule: by Jacobi's formula, $\det M$ — the volume stretch factor of the ellipsoid — evolves as $\tfrac{\mathrm d}{\mathrm dt}\ln\det M = \operatorname{tr} D\mathbf{f} = \nabla\cdot\mathbf{f}$ , so

\sum_{i=1}^{n}\lambda_i = \big\langle \nabla\cdot\mathbf{f} \big\rangle \qquad (5)

where: $\langle\cdot\rangle$ — the time average along the orbit. This is the bridge to the volume bookkeeping of the last three articles: for the driven pendulum the right side was the constant $-2\beta$ ; for any Hamiltonian flow it vanishes identically. Second, for Hamiltonian systems the symplectic structure sharpens "sums to zero" into a term-by-term symmetry (Derivation C):

\lambda_{n+1-i} = -\lambda_i \qquad \text{(symplectic flows)} \qquad (6)

where: every stretching direction has a mirror-image contracting partner at exactly the opposite rate — the double pendulum's $\{+\lambda, 0, 0, -\lambda\}$ was this rule in action. Two derived quantities convert the spectrum into forecasts and geometry. The Lyapunov time and prediction horizon (imported from the sensitive-dependence article):

t_\lambda = \frac{1}{\lambda_1}, \qquad t_{\text{pred}} \sim \frac{1}{\lambda_1}\ln\frac{\Delta}{\delta_0} \qquad (7)

where: $t_\lambda$ — the e-folding time of errors, the natural clock of unpredictability. And the Kaplan–Yorke dimension:

D_{KY} = k + \frac{1}{|\lambda_{k+1}|}\sum_{i=1}^{k}\lambda_i, \qquad k = \max\Big\{j : \textstyle\sum_{i=1}^{j}\lambda_i \geq 0\Big\} \qquad (8)

where: $k$ — the number of directions whose combined stretching is still non-negative; the fractional part interpolates: the attractor is "as many dimensions as the dynamics can keep from contracting." For the Lorenz attractor, spectrum $(0.906, 0, -14.57)$ : $k = 2$ and $D_{KY} = 2 + 0.906/14.57 = 2.062$ — a fractal barely thicker than a surface, and the subject of the next article.

Which tool answers which question

you want	use	notes
does an error grow, and how fast?	$\lambda_1$ via (3)	generic directions inherit the top rate
the full portrait	spectrum signature	$(-,-,-)$ point · $(0,-,-)$ cycle · $(0,0,-)$ torus · $(+,0,-)$ chaos
consistency check on numerics	sum rule (5)	must match $\langle\nabla\cdot\mathbf{f}\rangle$ exactly
Hamiltonian bookkeeping	pairing (6)	spectrum symmetric; no attractors possible
attractor dimension	Kaplan–Yorke (8)	dynamics → geometry, no box counting
exponents from data	delay embedding + Rosenstein/Wolf	only $\lambda_1$ (and 0) are robust from experiments

Derivation AThe exponent of a one-dimensional map — the ergodic average of the log-slope

We derive the workhorse formula of computational chaos: for a map $x_{n+1} = g(x_n)$ , the Lyapunov exponent is the orbit average of $\ln\lvert g'\rvert$ . The strategy is to chain the one-step linearization into a product and let the logarithm convert products into averages — three lines of algebra hiding one deep idea, the ergodic exchange of time averages for space averages. The step to watch is the final one, where a statement about one orbit becomes a statement about an invariant density: that exchange is what makes the exponent a property of the system.

Step 1. Linearize a single step: a small error $\delta_n$ at $x_n$ becomes $\delta_{n+1} = g'(x_n)\,\delta_n$ — the map's local slope is the one-dimensional Jacobian, stretching where $\lvert g'\rvert > 1$ and compressing where $\lvert g'\rvert < 1$ .

Step 2. Chain $N$ steps: $\delta_N = \delta_0\prod_{n=0}^{N-1} g'(x_n)$ — the 1D avatar of the propagator $M(t)$ , with the time-ordering trivial because scalars commute.

Step 3. Define the rate and take logarithms:

\lambda = \lim_{N\to\infty}\frac1N \ln\Big|\frac{\delta_N}{\delta_0}\Big| = \lim_{N\to\infty}\frac1N \sum_{n=0}^{N-1}\ln\big|g'(x_n)\big|

Why the logarithm is essential: the growth factors multiply, and multiplicative fluctuations (slopes ranging from 0 to 4 for the logistic map) average meaningfully only in the log — the exponent is a geometric mean rate, not an arithmetic one.

Step 4. Exchange time for space. For an ergodic invariant density $\rho(x)$ , the time average equals the space average:

\lambda = \int \rho(x)\,\ln\big|g'(x)\big|\,\mathrm{d}x

Why this matters: the right side makes no reference to any particular orbit — the exponent is revealed as a property of the map plus its invariant measure, computable analytically whenever $\rho$ is known.

Step 5. Cash the formula on the benchmark: for $g(x) = 4x(1-x)$ , the invariant density is $\rho(x) = 1/[\pi\sqrt{x(1-x)}]$ (from the tent-map conjugacy of the previous article), and the integral evaluates exactly to $\lambda = \ln 2 = 0.693$ per iteration. The numerically computed curve of Figure 2 lands on this value at $r = 4$ to three decimals — and the same formula, evaluated on a stable 2-cycle, returns $\lambda = \tfrac12\ln\lvert(f^2)'\rvert < 0$ , negative through the periodic cascade and zero exactly at each doubling threshold, reproducing the sawtooth of the figure. Verified: $\ln 2$ integral per Ott §2.3 and Strogatz §10.5; the Figure 2 computation gives $\lambda(4.0) = 0.6927$ vs $\ln 2 = 0.6931$ ; 2-cycle value at $r = 3.2$ : $\tfrac12\ln(0.16) = -0.92$ per iteration, matching direct simulation.

Go deeper · Derivation B · 1 minThe sum rule — spectrum meets phase-volume budget

Derivation BThe sum rule — spectrum meets phase-volume budget

The goal is equation (5): the Lyapunov exponents sum to the time-averaged divergence of the flow. This one identity welds the spectrum to everything the course knows about volumes — Liouville's theorem, dissipative contraction, the existence or impossibility of attractors. The strategy is to track the volume of the evolving ellipsoid two ways: as the product of its axis stretches (giving $\sum\lambda_i$ ) and via Jacobi's formula for $\det M$ (giving the divergence). Watch the logarithm once again do the work of converting a product of fluctuating factors into a clean time average.

Step 1. The ellipsoid's volume factor after time $t$ is $\det M(t)$ — the product of the singular values $\mu_i(t)$ , i.e. of the axis stretch factors: $\det M = \prod_i \mu_i$ .

Step 2. Differentiate the determinant along the flow with Jacobi's formula:

\frac{\mathrm{d}}{\mathrm{d}t}\ln\det M = \operatorname{tr}\big(M^{-1}\dot M\big) = \operatorname{tr}\big(M^{-1} D\mathbf{f}\, M\big) = \operatorname{tr} D\mathbf{f} = \nabla\cdot\mathbf{f}

Why the trace collapses: the trace is invariant under conjugation, so the propagator sandwich drops out — the volume growth rate at each instant is just the local divergence, regardless of how the ellipsoid is currently oriented.

Step 3. Integrate and average: $\ln\det M(t) = \int_0^t \nabla\cdot\mathbf{f}\,\mathrm{d}t'$ , so

\lim_{t\to\infty}\frac1t\ln\det M(t) = \big\langle\nabla\cdot\mathbf{f}\big\rangle

Step 4. Compute the same limit from the axes: $\tfrac1t\ln\det M = \tfrac1t\sum_i\ln\mu_i(t) \to \sum_i \lambda_i$ by the definition of the exponents as axis growth rates. Equate:

\sum_{i=1}^n \lambda_i = \big\langle\nabla\cdot\mathbf{f}\big\rangle

Step 5. Harvest the consequences. Dissipative flow ( $\langle\nabla\cdot\mathbf{f}\rangle < 0$ ): the sum is negative, so even with $\lambda_1 > 0$ the attractor has zero volume — strangeness is forced, and the contraction directions must overpower the stretching. Hamiltonian flow: sum exactly zero (Liouville), so a positive exponent requires an equal-magnitude negative partner — no net contraction, no attractor, chaotic seas instead. And numerically, the rule is a free integrity test: for Lorenz, the computed spectrum $(0.906, 0, -14.57)$ must sum to the constant divergence $-(\sigma + 1 + \beta) = -13.67$ — it does, to the third digit (Worked Example 3). Verified: Jacobi's formula and the transport theorem per Ott §4.2; Lorenz check $0.906 + 0 - 14.57 = -13.66$ vs $-41/3 = -13.667$ ; driven-pendulum case reproduces the exact $-2\beta$ of the previous article.

Go deeper · Derivation C · 2 minSymplectic pairing — why Hamiltonian spectra are mirror-symmetric

Derivation CSymplectic pairing — why Hamiltonian spectra are mirror-symmetric

We prove the pairing rule (6): for Hamiltonian flows the exponents come in exact $\pm$ pairs. This is stronger than the zero-sum rule and is pure symplectic geometry — the same structure that powered the canonical-transformation articles now disciplining chaos. The strategy: show the tangent propagator is a symplectic matrix at all times, then convert the symplectic condition into a statement about singular values. The step to watch is the similarity trick in Step 3 — one line of matrix algebra that forces the entire mirror symmetry.

Step 1. Write Hamilton's equations as a flow: $\dot{\mathbf z} = J\,\nabla H(\mathbf z)$ with $\mathbf z = (\mathbf q, \mathbf p)$ and the symplectic unit $J = \begin{pmatrix} 0 & \mathbb{I} \\ -\mathbb{I} & 0 \end{pmatrix}$ , $J^{\mathsf T} = -J = J^{-1}$ . The variational equation is then $\dot M = J\,(\nabla^2 H)\,M$ , with $\nabla^2 H$ the (symmetric) Hessian along the orbit.

Step 2. Show the propagator preserves the symplectic form:

\frac{\mathrm{d}}{\mathrm{d}t}\big(M^{\mathsf T} J M\big) = M^{\mathsf T}(\nabla^2 H)^{\mathsf T} J^{\mathsf T} J M + M^{\mathsf T} J J (\nabla^2 H) M = M^{\mathsf T}\big(\nabla^2 H\big)\big(J^{\mathsf T}J + JJ\big)...

carried through carefully: $\dot M^{\mathsf T} J M = M^{\mathsf T}(\nabla^2 H) J^{\mathsf T} J M = M^{\mathsf T}(\nabla^2 H) M$ and $M^{\mathsf T} J \dot M = M^{\mathsf T} J J (\nabla^2 H) M = -M^{\mathsf T}(\nabla^2 H) M$ ; the two terms cancel. Hence $M^{\mathsf T} J M = J$ for all $t$ : the tangent flow is symplectic. Why this is the linearized Liouville-and-more: preservation of $J$ preserves all the Poincaré invariants of the W4 articles, volume among them — but the matrix identity carries finer information than the determinant alone.

Step 3. Convert to a spectral statement. From $M^{\mathsf T} J M = J$ : $M^{-1} = J^{-1} M^{\mathsf T} J$ — the inverse propagator is similar to the transpose. Therefore $M^{-1}$ and $M^{\mathsf T}$ (hence $M$ ) share eigenvalue spectra, and more usefully, $M^{\mathsf T} M$ is similar to $(M M^{\mathsf T})^{-1}$ : if $\mu^2$ is an eigenvalue of $M^{\mathsf T}M$ , so is $\mu^{-2}$ . Why singular values and not eigenvalues matter: the exponents are growth rates of lengths — singular value data — and the similarity chain above forces the singular values into reciprocal pairs $\mu, 1/\mu$ up to bounded factors, which is what survives the $\tfrac1t\ln$ limit.

Step 4. Take the limit: reciprocal singular-value pairs give exponent pairs $\lambda, -\lambda$ . Ordered, the spectrum satisfies $\lambda_{n+1-i} = -\lambda_i$ , and the sum vanishes term by term — Liouville recovered as the shadow of a sharper symmetry.

Step 5. Add the two structural zeros for bounded motion: the flow direction always contributes $\lambda = 0$ (perturbations along the orbit merely reparametrize time), and its symplectic partner — the direction across energy surfaces — contributes the paired zero. For the double pendulum: $\{+\lambda, 0, 0, -\lambda\}$ , exactly as measured in the previous article. Verified: symplecticity of the tangent map and reciprocal pairing per Arnold §44 and Goldstein §9.4; the flow-direction zero follows since $\delta\mathbf x = \mathbf f(\mathbf x(t))$ solves the variational equation identically (differentiate $\dot{\mathbf x} = \mathbf f$ ).

Limiting cases

Case	Spectrum signature	Meaning
stable fixed point	$(-, -, \ldots, -)$	all directions contract; no zero (nothing flows)
stable limit cycle	$(0, -, \ldots, -)$	one neutral direction along the orbit
quasiperiodic 2-torus	$(0, 0, -, \ldots)$	two angle directions, both neutral
strange attractor (dissipative)	$(+, 0, -, \ldots)$ , $\sum < 0$	stretch + slide + overpowering squeeze; needs $n \geq 3$ for flows
Hamiltonian chaos	$\lambda_i = -\lambda_{n+1-i}$ , $\sum = 0$	mirror pairs; chaotic sea, no attractor
conservative regular motion (KAM torus)	all zeros ( $n = 2$ DOF: $0,0,0,0$ )	polynomial separation only
1D map	single $\lambda = \langle\ln\lvert g'\rvert\rangle$	per-iteration rate; $\ln 2$ at full logistic chaos

Common trapA positive local stretching rate is not a positive Lyapunov exponent — and naive twin-trajectory measurements produce neither

Two failure modes account for most wrong Lyapunov claims. First, localism: a positive Jacobian eigenvalue at a point (or even along a stretch of orbit) can be cancelled by contraction elsewhere — the exponent is the infinite-time average, and only Oseledets-style limits count. The logistic map at $r = 3.2$ has local slopes up to $3.2$ (log-positive) yet $\lambda = -0.92$ : stable. Second, saturation: integrating two nearby full trajectories and fitting $\ln$ (separation) works only until the separation approaches the attractor size, after which the "exponent" decays toward zero and contaminates everything; and using a finite separation from the start measures a mixture of exponents, not $\lambda_1$ . The professional cures are exactly the algorithmic points of this article: evolve tangent vectors (infinitesimal by construction) with the variational equation, renormalize periodically (Benettin), and for the full spectrum re-orthogonalize a frame — otherwise every vector collapses onto the leading direction and the subleading exponents are garbage. Always report the sum-rule check: a spectrum that fails $\sum\lambda_i = \langle\nabla\cdot\mathbf f\rangle$ is not a measurement, it is a bug.

Worked examples

1 · Reading a measured spectrum (confidence builder)

A 3D dissipative flow yields the spectrum $(0.5, 0, -2.0)\ \mathrm{s^{-1}}$ . Classify the dynamics, compute the Kaplan–Yorke dimension, the Lyapunov time, and the horizon for initial error $10^{-6}$ with tolerance $1$ (attractor units).

Step 1. Signature $(+, 0, -)$ on a bounded set: chaotic. Sum $= -1.5 < 0$ : dissipative — a strange attractor. The zero confirms a genuine flow (not a map) away from equilibrium.

Step 2. Kaplan–Yorke (8): partial sums $0.5 \geq 0$ , $0.5 + 0 = 0.5 \geq 0$ , $0.5 + 0 - 2.0 = -1.5 < 0$ , so $k = 2$ and $D_{KY} = 2 + 0.5/2.0 = \boxed{2.25}$ .

Step 3. Lyapunov time $t_\lambda = 1/0.5 = \boxed{2\ \mathrm{s}}$ ; horizon $t_{\text{pred}} = t_\lambda \ln(10^6) = 2 \times 13.8 = \boxed{27.6\ \mathrm{s}}$ — about 14 e-foldings, the standard budget.

2 · The exponent as cascade bookkeeper (exam level)

For the logistic map: (i) verify $\lambda$ at the superstable period-2 parameter $r = 1+\sqrt5$ is $-\infty$ in principle and explain why; (ii) compute $\lambda$ on the 2-cycle at $r = 3.2$ ; (iii) state $\lambda$ at $r_1 = 3$ , at $r_\infty$ , and at $r = 4$ , and connect each to Figure 2.

Step 1. At $r = 1+\sqrt5 = 3.236$ , the 2-cycle contains the map's maximum $x = \tfrac12$ where $g'(\tfrac12) = 0$ : the orbit average of $\ln\lvert g'\rvert$ contains $\ln 0 = -\infty$ . Perturbations are annihilated (not just shrunk) once per period — superstability, the deep negative spikes of Figure 2.

Step 2. At $r = 3.2$ the 2-cycle multiplier (previous article, Derivation A) is $(f^2)' = 4 + 2r - r^2 = 0.16$ , so per iteration $\lambda = \tfrac12\ln\lvert 0.16\rvert = \boxed{-0.92}$ — strongly stable, matching the deep valley between $r_1$ and the spike.

Step 3. At $r_1 = 3$ : multiplier $-1$ , so $\lambda = \ln\lvert{-1}\rvert = \boxed{0}$ — the exponent touches zero at every doubling threshold, the sawtooth teeth of Figure 2. At $r_\infty = 3.5699$ : $\lambda = 0$ arriving from below with no room left — the marginal Feigenbaum attractor. At $r = 4$ : $\lambda = \ln 2 = \boxed{0.693}$ , the ceiling — one bit of initial-condition information consumed per iteration ( $h_{KS} = \ln 2$ by Pesin). The whole cascade is thus visible as one curve: negative sawtooth, zero at the accumulation, positive plateau with window dips.

3 · The Lorenz spectrum audited (applied)

The Lorenz system ( $\sigma = 10$ , $\rho = 28$ , $\beta = 8/3$ ) has measured spectrum $(0.906, 0, -14.57)$ . Audit it against the sum rule, and give $D_{KY}$ .

Step 1. The divergence is constant: $\nabla\cdot\mathbf f = -\sigma - 1 - \beta = -10 - 1 - 2.667 = -13.667$ — no averaging needed, the cleanest possible test case.

Step 2. Spectrum sum: $0.906 + 0 - 14.57 = -13.664$ . Agreement to $0.02\%$ : the measurement passes its integrity audit. (Any proposed spectrum failing this test — e.g. one omitting the zero exponent — is instantly refuted.)

Step 3. Dimension: $k = 2$ , $D_{KY} = 2 + 0.906/14.57 = \boxed{2.062}$ — the Lorenz attractor is two-dimensional plus six percent, in excellent agreement with direct fractal-dimension measurements. The next article takes exactly this handoff: what that fractional $0.062$ looks like.

Origins · 1 minAleksandr Lyapunov and how this was discovered

Aleksandr Lyapunov

1857–1918 · St. Petersburg · Kharkov

Lyapunov built the theory of stability of motion two generations before anyone could compute with it. His 1892 doctoral dissertation, The General Problem of the Stability of Motion, written in Kharkov under the long shadow of his teacher Chebyshev, introduced both pillars of the modern subject: the "first method," which extracts stability from the growth exponents of linearized solutions — the direct ancestor of the numbers bearing his name — and the "second method" of Lyapunov functions, which certifies stability with no linearization at all and now runs through control theory, from power grids to rocket guidance. The motivating problem was celestial and concrete: whether rotating self-gravitating fluid masses (model planets) hold their shape, a question on which he duelled politely with Poincaré for decades, deciding delicate stability cases the Frenchman had left open. His exponents waited three-quarters of a century for their apotheosis: Oseledets' 1968 theorem secured their existence for generic trajectories, and the chaos revolution of the 1970s made them the field's standard measure. Lyapunov did not live to see any of it — he shot himself in Odessa in 1918, three days after his wife's death from tuberculosis, ending a life of severe and solitary rigor whose definitions turned out to be exactly what a science he never imagined would need.

History

1892

— Lyapunov's dissertation founds stability theory: characteristic exponents of linearized motion and the direct method of auxiliary functions.

1968

— Oseledets proves the multiplicative ergodic theorem: the exponents exist as almost-sure limits — the field's existence license.

1976

— Pesin ties the positive exponents to entropy, equating information production with stretching for smooth chaotic systems.

1980

— Benettin, Galgani, Giorgilli and Strelcyn publish the grow-log-renormalize algorithm with Gram–Schmidt, making full spectra computable.

1985

— Wolf, Swift, Swinney and Vastano extract exponents from experimental time series via delay embedding — chaos becomes measurable in the lab.

Going further

References. Strogatz, Nonlinear Dynamics and Chaos, 2nd ed., §9.3, §10.5, §11 · Ott, Chaos in Dynamical Systems, 2nd ed., Ch. 4 · Goldstein, Poole & Safko, Classical Mechanics, 3rd ed., §11.6–11.8 · Landau & Lifshitz, Mechanics, Vol. 1 (3rd ed.), discussion of phase-space flow · Benettin, Galgani, Giorgilli & Strelcyn, Meccanica 15, 9 (1980) · Wolf, Swift, Swinney & Vastano, Physica D 16, 285 (1985)

Provenance. Variational equation, fundamental matrix and the Oseledec definition follow Strogatz §9.3 and Ott Ch. 4; verify the tangent-equation derivation by direct Taylor expansion. The sum rule via Jacobi's formula and Liouville's theorem is standard (Ott §4.2). Symplectic pairing checked against Arnold §44 and Goldstein §9.4. Numbers: logistic-map exponent $\ln2\approx0.693$ at $r=4$ and Feigenbaum point $r_\infty\approx3.5699$ (Strogatz §10.5, §10.6); Lorenz spectrum $(0.906,0,-14.57)$ and $D_{KY}\approx2.06$ (Strogatz §9.3, Wolf et al. 1985) — recompute $D_{KY}=2+0.906/14.57$ . Solar-system Lyapunov time $\sim5$ Myr from Laskar's integrations. All exponents ordered $\lambda_1\ge\dots\ge\lambda_n$ .

Builds on: Linear Stability Analysis and the Jacobian · Phase Portraits and the Phase Plane · Classification of Fixed Points: Nodes, Saddles, Spirals, Centers · Conservative Systems and Hamiltonian Phase Flow Leads to: Period-Doubling Cascade and the Route to Chaos · Bifurcations and the Loss of Stability · Limit Cycles and Self-Sustained Oscillations · Hopf Bifurcation and Birth of Limit Cycles

Why do people integrate the variational equation instead of just tracking two nearby trajectories?

Because the definition demands the limit of infinitesimal separation, and finite twins violate it in two ways. First, saturation: a real pair separates until it spans the attractor, after which the log-separation curve flattens and the fitted "exponent" sinks — you measure attractor diameter, not dynamics. Second, contamination: a finite separation vector has weight on all the Oseledets directions and samples nonlinear terms; its growth mixes exponents unpredictably. The variational equation solves both at a stroke: the tangent vector is infinitesimal by construction (it lives in the tangent space, evolved by the exact linearization along the true orbit), so it never saturates, and periodic renormalization handles the numerical overflow that would otherwise occur after $\sim 700/\lambda$ digits of growth. The twin-trajectory method survives only as a data-analysis tool (Rosenstein/Kantz) where no equations are available — and there it is used carefully, fitting only the pre-saturation window at small scales.

Can the full spectrum be extracted from an experimental time series?

The top exponent and the zero, yes, with care; the negative tail, essentially no. From a single measured signal one reconstructs the attractor by delay embedding — Takens' theorem guarantees the reconstruction is diffeomorphic to the true attractor for generic measurements — and then tracks how small neighbourhoods separate (Rosenstein/Kantz for $\lambda_1$ , checking stability against embedding dimension and delay). The zero exponent shows up as a reliable second plateau. But negative exponents describe how off-attractor perturbations decay, and the data contains no off-attractor points: everything measured has already collapsed onto the attractor, so the contracting directions are sampled only through their fractal residue, and estimates of $\lambda_3$ and beyond are noise-dominated. Consequently experimental Kaplan–Yorke dimensions are usually computed the other way — dimension measured geometrically (next article), exponents inferred — or the negative exponents are taken from a fitted model. Simulation and experiment thus divide the spectrum between them: models give you all of it, data gives you the top.

Are Lyapunov exponents coordinate-independent? I could measure separation in angles or in energies and get different numbers.

The exponents are invariant under any smooth change of coordinates with bounded derivatives — precisely because a diffeomorphism multiplies separations by factors bounded above and below, which contribute $\tfrac1t\ln(\text{bounded})\to 0$ . Angles, Cartesian positions, canonical pairs: same spectrum. This is what licenses delay-embedding measurements (the embedding is a coordinate change) and what makes the exponent a property of the dynamics rather than the chart. The fine print: the transformation must be smooth and non-degenerate on the invariant set (singular changes can distort exponents), the finite-time exponents and the covariant directions are NOT invariant (only the asymptotic rates are), and for flows the exponents carry units of inverse time, so rescaling time rescales the spectrum — per-second and per-orbit conventions differ by the mean period, a factor worth pinning down before comparing literature values.

What do the fluctuations of finite-time exponents mean — is the infinite-time limit hiding physics?

Yes, and modern usage increasingly mines it. The finite-time (local) exponent $\lambda_1(T)$ — the growth rate over a window $T$ — fluctuates around the asymptotic value with a distribution whose width shrinks like $1/\sqrt T$ but whose shape encodes real structure: intermittency shows up as heavy tails (episodes of anomalous stability or instability), sticky KAM-island boundaries in Hamiltonian systems produce long low- $\lambda$ episodes, and in forecasting, today's local exponent is tomorrow's error bar — weather predictability genuinely varies from day to day, and ensemble spread tracks the local, not global, stretching. Large-deviation theory formalizes this: the distribution of finite-time exponents obeys a rate function whose Legendre transform generates the generalized dimensions and entropies of multifractal analysis. So the single number $\lambda_1$ is the headline; the fluctuation spectrum around it is the full story, and for mixed phase spaces (double pendulum at moderate energy) a map of local exponents is the honest portrait.

How exactly does Gram–Schmidt pull out the subleading exponents rather than k copies of the top one?

By forcing each successive vector to live in the orthogonal complement of the faster ones, where the leading growth has been projected away. Evolve an orthonormal frame $\{u_1, \ldots, u_k\}$ with the variational equation: every vector, left alone, would rotate toward the top Oseledets direction and grow at $\lambda_1$ . At each renormalization, Gram–Schmidt keeps $u_1$ 's direction (it legitimately tracks the fastest direction), then strips from $u_2$ its component along $u_1$ before renormalizing — so the logged growth factor of $u_2$ measures expansion in the 2D subspace modulo the 1D one, which grows like $e^{(\lambda_1+\lambda_2)t}$ divided by $e^{\lambda_1 t}$ : rate $\lambda_2$ . Inductively, the $j$ -th vector's logged factors average to $\lambda_j$ because the $j$ -volume of the frame grows at $\sum_{i\le j}\lambda_i$ and Gram–Schmidt factorizes that volume growth into per-vector pieces (a running QR decomposition of $M(t)$ , in linear-algebra language). This is also why the sums $\sum_{i \le j}\lambda_i$ converge fastest and are the natural objects — the Kaplan–Yorke formula is built from exactly those sums.

Is a positive Lyapunov exponent by itself proof of chaos?

It is the strongest single indicator, but the honest checklist has three lines. (1) $\lambda_1 > 0$ — exponential separation. (2) Boundedness/recurrence — the motion stays in a compact invariant set, so the stretching is recycled rather than a one-shot escape (a linear saddle has positive local exponents and no chaos; escape to infinity can also fake positivity over finite windows). (3) Aperiodicity of the sampled orbit — ruling out the measure-zero pathology of sitting exactly on an unstable periodic orbit, whose exponent is positive but whose motion is periodic. In practice a fourth line matters most: robustness — the estimate must be stable under longer integration, different initial conditions in the same basin, refined tolerances, and it must pass the sum-rule audit. For experimental data, add stability under embedding parameters and surrogate-data tests against colored noise. When all lines check, "chaotic" is as certified as empirical dynamics gets; any single line alone has famous counterexamples.

Lyapunov Exponents

Where this lives today

The physical picture

The formal core

Which tool answers which question

Limiting cases

Worked examples

1 · Reading a measured spectrum (confidence builder)

2 · The exponent as cascade bookkeeper (exam level)

3 · The Lorenz spectrum audited (applied)

Going further

Before this

Where this leads