Integrability and chaos: determinism without predictability

Idea. Complete integrability is \(n\) functions, pairwise commuting, independent. The counting is the easy half; the involution is where systems fail, and the standard failure is using several components of an angular momentum at once.

1. The anisotropic oscillator. \(n = 2\). Take \(F_1 = \tfrac12(p_1^2+\omega_1^2q_1^2)\), \(F_2 = \tfrac12(p_2^2+\omega_2^2q_2^2)\), so \(H = F_1 + F_2\).

Brackets: \(F_1\) depends only on \((q_1,p_1)\) and \(F_2\) only on \((q_2,p_2)\), so in \(\{F_1,F_2\} = \sum_i(\partial_{q_i}F_1\partial_{p_i}F_2 - \partial_{p_i}F_1\partial_{q_i}F_2)\) every term has one factor from a variable the other function does not contain. Term by term: \(\partial_{q_1}F_1\cdot\partial_{p_1}F_2 = \omega_1^2q_1\cdot0 = 0\), \(\partial_{p_1}F_1\cdot\partial_{q_1}F_2 = p_1\cdot0 = 0\), and the \(i=2\) pair vanishes the same way. So \(\{F_1,F_2\} = 0\) ✓, and \(\{F_i,H\} = \{F_i,F_1+F_2\} = 0\) ✓.

Independence: \(\mathrm dF_1 = (\omega_1^2q_1, 0, p_1, 0)\) and \(\mathrm dF_2 = (0,\omega_2^2q_2,0,p_2)\) are linearly independent unless one of them vanishes, which happens only when \((q_i,p_i) = (0,0)\) — the two planes of pure motion in one coordinate. On the open set where both \(F_i\gt0\), completely integrable.

2. The planar central force. \(n = 2\). \(\theta\) is cyclic, so \(F_1 = H\) and \(F_2 = p_\theta\) are both conserved. Bracket: \(H\) does not contain \(\theta\), so \(\partial_\theta H = 0\), and \(p_\theta\) does not contain \(p_r\) or \(r\), so \(\{H,p_\theta\} = \partial_rH\cdot0 - \partial_{p_r}H\cdot0 + \partial_\theta H\cdot1 - \partial_{p_\theta}H\cdot0 = 0\) ✓. Independence: in the order \((r,\theta,p_r,p_\theta)\), \(\mathrm dH = (\partial_rH,\,0,\,p_r/m,\,p_\theta/mr^2)\) and \(\mathrm dp_\theta = (0,0,0,1)\). These are dependent exactly when the first three entries of \(\mathrm dH\) vanish: \(p_r = 0\) and \(\partial_rH = V'(r) - p_\theta^2/mr^3 = 0\) — the circular orbits, where the level set of \((H,p_\theta)\) degenerates to a single circle. That set is closed, so its complement in \(\{r\gt0\}\) is open, and there the system is completely integrable. On all of \(\{r\gt0\}\) it is not, by the letter of the definition — a first taste of where the definition bites.

3. Three dimensions. \(n = 3\), so three integrals are needed — but the proposed triple fails, because unit 14 computed \(\{L_1,L_2\} = L_3\), which is not identically zero. The repair: use \(H\), \(L_3\) and \(\lVert\vec L\rVert^2\). Then \(\{L_3,\lVert\vec L\rVert^2\} = \sum_i 2L_i\{L_3,L_i\} = 2\sum_i L_i\epsilon_{3ik}L_k = 0\), because \(\epsilon_{3ik}\) is antisymmetric under \(i\leftrightarrow k\) while \(L_iL_k\) is symmetric, so the sum cancels in pairs ✓. Both commute with \(H\) as conserved quantities. Their differentials are independent away from a closed set — the circular orbits again, and the states with \(\vec L\) along the third axis, where \(\mathrm d\lVert\vec L\rVert^2 = 2L_3\,\mathrm dL_3\) — so on its open complement the system is completely integrable.

Which of the seven Kepler integrals may be used together: \(H\), one component of \(\vec L\), and \(\lVert\vec L\rVert^2\). You may not add a second component of \(\vec L\), nor any component of the Runge–Lenz vector \(\vec A\), since unit 14 found \(\{L_i,A_j\}\ne0\) in general.

The last part. No. Kepler has seven conserved quantities, five of them independent — at most \(2n-1 = 5\) can be (P3), and five is more than \(n = 3\) — and it is integrable; yet the number alone proves nothing, because what is required is \(n\) that commute. The counterexample is exactly part 3: \(H, L_1, L_2\) are three independent integrals of a three-degree-of-freedom system and do not establish integrability.

Check. Dimensions: each \(F_i\) in part 1 is an energy ✓; \(p_\theta\) in part 2 is an angular momentum and \(H\) an energy, and involution is a statement about a bracket, which has dimensions energy \(\times\) angular momentum / action \(=\) energy — nonzero dimensions do not stop it vanishing ✓. Special case in part 1: \(\omega_1 = \omega_2\) gives the isotropic oscillator, which has the extra integrals of the hidden \(\mathrm{SU}(2)\) symmetry; our two still work, so "integrable" is not disturbed by having more ✓.

Idea. Each degree of freedom is an independent one-dimensional oscillator, so each action is the area of an ellipse over \(2\pi\). Everything else follows by differentiation and one piece of arithmetic about integers.

1. The actions. Fix \(F_i\). The level curve \(\tfrac12p_i^2 + \tfrac12\omega_i^2q_i^2 = F_i\) is the ellipse \(\frac{q_i^2}{(\sqrt{2F_i}/\omega_i)^2} + \frac{p_i^2}{(\sqrt{2F_i})^2} = 1\), with semi-axes \(a = \sqrt{2F_i}/\omega_i\) and \(b = \sqrt{2F_i}\). Its area is \(\pi ab = \pi\cdot\frac{\sqrt{2F_i}}{\omega_i}\cdot\sqrt{2F_i} = \frac{2\pi F_i}{\omega_i}\). Hence \[I_i \;=\; \frac1{2\pi}\oint p_i\,\mathrm dq_i \;=\; \frac{F_i}{\omega_i}.\]

2. The Hamiltonian and the frequencies. \(F_i = \omega_iI_i\), so \[H = F_1 + F_2 = \omega_1I_1 + \omega_2I_2,\qquad \frac{\partial H}{\partial I_i} = \omega_i.\] The frequencies come out equal to the \(\omega_i\) already in the Hamiltonian, which is the check: the oscillator's angular frequency is \(\omega\), known since unit 00.

3. \(\omega = (1,\sqrt2)\). Suppose \(k_1 + k_2\sqrt2 = 0\) with \((k_1,k_2)\in\mathbb Z^2\setminus\{0\}\). If \(k_2 = 0\) then \(k_1 = 0\), excluded. So \(k_2\ne0\) and \(\sqrt2 = -k_1/k_2\) is rational, which it is not. Hence \(\vec\omega\) is non-resonant, and by the Liouville–Arnold discussion every orbit on such a torus is dense in it: a single trajectory eventually passes arbitrarily near every state of that energy pair. Its Lissajous figure never closes.

4. \(\omega = (2,3)\). Take \(\vec k = (3,-2)\): \(3\cdot2 - 2\cdot3 = 0\) ✓. The torus is resonant, and every orbit on it is periodic. The two separate periods are \(2\pi/\omega_1 = \pi\) and \(2\pi/\omega_2 = 2\pi/3\); the orbit closes at their least common multiple, \(T = 2\pi\), for which \(\omega_1T = 4\pi\) and \(\omega_2T = 6\pi\) are both multiples of \(2\pi\) ✓ — two oscillations of \(q_1\) against three of \(q_2\). The orbit is a closed Lissajous curve, a one-dimensional subset of the two-dimensional torus.

5. Non-degeneracy. \(H = \omega_1I_1 + \omega_2I_2\) is linear in the actions, so \(\partial^2H/\partial I_i\partial I_j = 0\) identically and the determinant is zero. The condition fails. What follows is a warning: the frequencies do not vary from torus to torus, so the resonant and non-resonant tori are not interleaved — either all tori are resonant or none are, depending on a single ratio. Standard KAM does not apply, and a perturbation of a degenerate system can behave quite differently from the generic picture. This is exactly the situation of the Hénon–Heiles system's integrable limit, where \(\omega_1 = \omega_2 = 1\) for every torus.

Check. Dimensions: \(I = F/\omega\) is energy over frequency, an action ✓, and \(\partial H/\partial I\) is energy over action, a frequency ✓. Limit: set \(\omega_2\to\omega_1\) and the system becomes isotropic; then \(H = \omega_1(I_1+I_2)\) depends on the sum only, every orbit is closed, and the extra degeneracy shows up as the vanishing of \(\omega_1-\omega_2\) ✓. A number: with \(\omega = (1,\sqrt2)\) and both actions equal to \(1\), \(H = 1+\sqrt2 = 2.414\).

Idea. Involution is a statement that a set of gradients spans a subspace on which an antisymmetric nondegenerate form vanishes identically. Such subspaces cannot be large, and the bound is exactly half the dimension — which is why the definition of integrability asks for \(n\) integrals and not more.

1. The matrix form. With \(z = (q_1,\dots,q_n,p_1,\dots,p_n)\), \(\nabla F = (\partial_{q}F, \partial_pF)\) and \[(\nabla F)^{\!\top}J\nabla G = (\partial_qF,\partial_pF)\begin{pmatrix}0&I\\-I&0\end{pmatrix} \begin{pmatrix}\partial_qG\\ \partial_pG\end{pmatrix} = (\partial_qF,\partial_pF)\begin{pmatrix}\partial_pG\\ -\partial_qG\end{pmatrix} = \partial_qF\cdot\partial_pG - \partial_pF\cdot\partial_qG,\] which is the coordinate formula ✓. \(J^{\!\top} = -J\) by inspection, and \(J^2 = -I_{2n}\) so \(J^{-1} = -J\): invertible ✓.

2. Isotropic subspaces have \(\dim\le n\). Let \(W\) have basis \(w_1,\dots,w_k\) and define \(\Phi:\mathbb R^{2n}\to\mathbb R^k\) by \(\Phi(u)_i = w_i^{\!\top}Ju\). By construction \(\ker\Phi = W^{\perp J}\).

\(\Phi\) is surjective. The \(i\)-th component of \(\Phi\) is the linear functional \(u\mapsto(J^{\!\top}w_i)^{\!\top}u\), so the rows of the matrix of \(\Phi\) are the vectors \(J^{\!\top}w_i = -Jw_i\). Since \(J\) is invertible and the \(w_i\) are independent, the \(Jw_i\) are independent, so the matrix has rank \(k\) and \(\Phi\) is onto.

Rank–nullity gives \(\dim W^{\perp J} = 2n - k\). Isotropy says precisely \(W\subseteq W^{\perp J}\), so \(k\le 2n-k\), that is \[\boxed{\ \dim W\le n.\ }\qquad\square\]

3. The deduction. Let \(v_i = \nabla F_i(z_0)\). Involution says \(v_i^{\!\top}Jv_j = \{F_i,F_j\}(z_0) = 0\) for all \(i,j\), so by bilinearity every pair of vectors in \(W = \operatorname{span}\{v_1,\dots,v_k\}\) pairs to zero: \(W\) is isotropic. Independence of the differentials says \(\dim W = k\). By part 2, \(k\le n\). \(\square\)

4. Without involution. The gradients then span a subspace on which \(J\) need not vanish, and no bound below \(2n\) applies — except that the \(F_i\) are constant along the flow, so their common level set contains the orbit and has dimension at least \(1\), giving \(k\le2n-1\). That bound is attained: the two-dimensional isotropic oscillator with \(\omega_1 = \omega_2\) has \(2n-1 = 3\) independent integrals, namely \(F_1\), \(F_2\) and the angular momentum \(L = q_1p_2 - q_2p_1\), since every orbit is a closed ellipse and three conditions pin it down. They are not all in involution: \(\{F_1,F_2\} = 0\), but \(\{F_1, L\} = q_1\cdot(-q_2) - p_1\cdot p_2 = -(q_1q_2 + p_1p_2) \ne 0\).

Check. Dimensions: part 2 is dimensionless linear algebra ✓. Special case \(n = 1\): the bound says at most one integral can be in involution with itself and others, and indeed \(H\) alone is all a one-degree-of-freedom system has — a second independent commuting integral would force \(\dim W = 2 \gt 1\) ✓. A number: for \(n = 3\), the bound is \(3\), which is why part 3 of P1 needed exactly three and not four.

Another way

Avoid \(W^{\perp J}\) entirely: suppose \(k\ge n+1\) and consider the \(2n\times k\) matrix \(A\) whose columns are the \(v_i\). Isotropy says \(A^{\!\top}JA = 0\), so the column space of \(JA\) annihilates the column space of \(A\). But \(\operatorname{rank}(JA) = \operatorname{rank}A = k\) because \(J\) is invertible, and two subspaces of \(\mathbb R^{2n}\) of dimension \(k\) each, annihilating one another under the standard pairing, force \(2k\le2n\). Same conclusion, fewer words, but the first route explains what \(W^{\perp J}\) is, which recurs whenever "Lagrangian submanifold" is said.

Idea. Infinitely many disjoint copies of a set of positive area cannot fit into a region of finite area. That is the whole proof, and it uses no dynamics beyond area preservation — which is why it applies to every Hamiltonian system without exception.

1. Some point returns. Write \(a = \operatorname{area}(U)\gt0\), which is positive because \(U\) is open and nonempty. Since \(P\) is an area-preserving bijection of \(D\), so is \(P^{-1}\), and therefore \(\operatorname{area}(P^{-i}(U)) = a\) for every \(i\ge0\).

Suppose, for contradiction, that the sets \(U, P^{-1}(U), P^{-2}(U),\dots\) are pairwise disjoint. Then for every \(N\), \[A \;=\; \operatorname{area}(D)\;\ge\;\operatorname{area}\Big(\bigcup_{i=0}^{N}P^{-i}(U)\Big) \;=\; (N+1)\,a,\] so \(A\ge(N+1)a\) for all \(N\) — impossible, since \(a\gt0\) and \(A\) is finite. Take \(N\gt A/a\) to see it explicitly.

So there are \(0\le i\lt j\) with \(P^{-i}(U)\cap P^{-j}(U)\ne\emptyset\). Apply the bijection \(P^{i}\) to both sides: \(U\cap P^{-(j-i)}(U)\ne\emptyset\). Hence there is \(u\in U\) with \(P^{m}(u)\in U\) for \(m = j - i\ge1\). \(\square\)

2. Infinitely many such \(m\). Fix any \(N\ge1\). The map \(P^N\) is again an area-preserving bijection of \(D\), so part 1 applied to \(P^N\) gives an \(r\ge1\) with \(P^{Nr}(U)\cap U\ne\emptyset\): a return to \(U\) at time \(m = Nr\ge N\). So the set of return times contains an element \(\ge N\) for every \(N\) — it is unbounded, hence infinite. \(\square\)

3. On a section. The section of a bounded two-degree-of-freedom Hamiltonian system has finite area and its return map preserves area, so: every open patch of the section, however small, contains a point whose orbit comes back to that patch, infinitely often. In particular chaotic motion cannot simply drain away from a region of the section: whatever patch it passes through, orbits through that patch come back. (That the same orbit returns is the almost-every-point version of part 4, and needs measure theory.)

4. What it does not say. It does not say that every point returns — only that some point of every open set does. (The stronger "almost every point returns" is true and needs measure theory.) And it puts no bound at all on the return time. A crude estimate: an orbit that visits the section uniformly enters a target of area \(a\) about once every \(A/a\) crossings, so for a disc of radius \(\varepsilon\), \[\langle m\rangle\;\sim\;\frac{A}{\pi\varepsilon^2}.\] With \(A = 1\) and \(\varepsilon = 10^{-3}\), that is \(1/(\pi\times10^{-6}) = 3.2\times10^{5}\) crossings. Halving \(\varepsilon\) quadruples it; asking for \(\varepsilon = 10^{-9}\) — the accuracy of a decent measurement — gives \(3.2\times10^{17}\) crossings, which at a few time units each is longer than any computation. Recurrence is a theorem about eternity, not a prediction.

Check. Dimensions: \(A/(\pi\varepsilon^2)\) is area over area, a pure count ✓. Limit \(\varepsilon\to\sqrt{A/\pi}\), a target as big as the whole section: the estimate gives \(\langle m\rangle\to1\), return on the first step ✓. Special case: if \(P\) is the identity, every point returns at \(m = 1\) and the theorem is trivially true — the proof never claimed the return was interesting ✓.

Idea. Area preservation is one determinant; the type of a fixed point of an area-preserving map is one trace. Both are mechanical once written down, and the second says everything about the local picture.

1. The standard map. \(DP = \begin{pmatrix}\partial\theta'/\partial\theta&\partial\theta'/\partial p\\ \partial p'/\partial\theta&\partial p'/\partial p\end{pmatrix} = \begin{pmatrix}1+K\cos\theta&1\\ K\cos\theta&1\end{pmatrix}\), so \(\det = (1+K\cos\theta)\cdot1 - 1\cdot K\cos\theta = \mathbf{1}\) for every \(K\) and \(\theta\). Yes, it can be such a return map — and it is: the kicked rotor's.

2. The Hénon map. \(DP = \begin{pmatrix}-2ax&1\\ b&0\end{pmatrix}\), \(\det = 0\cdot(-2ax) - 1\cdot b = -b\). Area-preserving exactly when \(|b| = 1\); orientation-preserving as well only for \(b = -1\). For the usual \(b = 0.3\) it contracts area by a factor \(0.3\) per step, so no: it is a dissipative map, and its strange attractor is a phenomenon this unit's systems cannot have.

3. The damped kicked rotor. \(\partial p'/\partial p = 1-\gamma\), \(\partial p'/\partial\theta = K\cos\theta\), \(\theta' = \theta + p'\) so \(\partial\theta'/\partial\theta = 1 + K\cos\theta\), \(\partial\theta'/\partial p = 1-\gamma\). Then \(\det = (1+K\cos\theta)(1-\gamma) - (1-\gamma)K\cos\theta = \mathbf{1-\gamma}\lt1\). No: area shrinks, and attractors become possible.

Fixed points of the standard map. \(p' = p\) needs \(K\sin\theta = 0\), so \(\theta\in\{0,\pi\}\); then \(\theta' = \theta + p'\equiv\theta\ (\mathrm{mod}\,2\pi)\) needs \(p = 0\). The fixed points are \((0,0)\) and \((\pi,0)\).

• At \((0,0)\): \(\cos\theta = 1\), so \(DP = \begin{pmatrix}1+K&1\\ K&1\end{pmatrix}\), \(\tau = 2+K\). For \(K\gt0\), \(\tau\gt2\): hyperbolic, always. Eigenvalues \(\mu_\pm = \big(2+K\pm\sqrt{(2+K)^2-4}\big)/2\).
• At \((\pi,0)\): \(\cos\theta = -1\), so \(DP = \begin{pmatrix}1-K&1\\ -K&1\end{pmatrix}\), \(\tau = 2-K\). Elliptic when \(|2-K|\lt2\), that is \(\mathbf{0\lt K\lt4}\); hyperbolic for \(K\gt4\).

At \(K = 4\) the trace reaches \(-2\), both eigenvalues equal \(-1\), and the elliptic point loses stability: a period-doubling, after which a stable orbit of period two appears nearby. This is the standard map's first bifurcation, and it happens well above \(K_c = 0.9716\), so the last spanning torus is long gone by then.

Numbers. \(K = 0.5\): at \((\pi,0)\), \(\tau = 1.5\), \(\mu_\pm = (1.5\pm\mathrm i\sqrt{4-2.25})/2 = 0.75\pm0.661\mathrm i\), modulus \(\sqrt{0.5625+0.4371} = 1.000\) ✓ — on the unit circle, elliptic. \(K = 5\): at \((\pi,0)\), \(\tau = -3\), \(\mu_\pm = (-3\pm\sqrt5)/2 = -0.382,\ -2.618\), real reciprocals (\(0.382\times2.618 = 1.000\) ✓) — hyperbolic, with a reflection.

Check. Dimensions: all three maps are dimensionless, and a determinant of a Jacobian between like variables is dimensionless ✓. Limit \(K\to0\) in the standard map: \(DP\) becomes \(\begin{pmatrix}1&1\\0&1\end{pmatrix}\), a shear, \(\tau = 2\) — parabolic, the boundary case, and correctly so, since at \(K = 0\) the whole line \(p = 0\) is fixed ✓. Limit \(\gamma\to0\) in part 3 recovers part 1 and \(\det\to1\) ✓.

Idea. At large \(K\) each step multiplies the \(\theta\)-component of a tangent vector by \(K\cos\theta\), and \(\theta\) wanders uniformly. So the exponent is the average of the logarithm of that factor, and the whole problem reduces to one classical integral.

1. The variational relations. Differentiate the map. From \(p_{n+1} = p_n + K\sin\theta_n\), \(\delta p_{n+1} = \delta p_n + K\cos\theta_n\,\delta\theta_n = \delta p_n + c_n\delta\theta_n\) ✓. From \(\theta_{n+1} = \theta_n + p_{n+1}\), \(\delta\theta_{n+1} = \delta\theta_n + \delta p_{n+1} = \delta\theta_n + \delta p_n + c_n\delta\theta_n = (1+c_n)\delta\theta_n + \delta p_n\) ✓.

2. The large-\(K\) reduction. If \(|c_n|\gg1\) and \(|\delta p_n|\lesssim|c_n\delta\theta_n|\), then \(\delta\theta_{n+1} = (1+c_n)\delta\theta_n + \delta p_n \approx c_n\delta\theta_n\). The assumption is self-consistent: with that approximation \(\delta p_{n+1} \approx c_n\delta\theta_n \approx \delta\theta_{n+1}\), so \(|\delta p_{n+1}|\approx|\delta\theta_{n+1}|\), which is smaller than \(|c_{n+1}\delta\theta_{n+1}|\) whenever \(|c_{n+1}|\gg1\) ✓. Hence \(\lVert v_N\rVert \approx \lVert v_0\rVert\prod_{n=0}^{N-1}|c_n|\) and \[\lambda_1 = \lim_N\frac1N\sum_{n=0}^{N-1}\ln|c_n| = \big\langle\ln|K\cos\theta|\big\rangle,\] the last step using that, in the chaotic regime at large \(K\), \(\theta_n\) is uniformly distributed on the circle.

3. The integral. Put \(I = \int_0^{\pi/2}\ln\cos x\,\mathrm dx\). Substituting \(x\mapsto\pi/2-x\) gives \(I = \int_0^{\pi/2}\ln\sin x\,\mathrm dx\). Adding, \[2I = \int_0^{\pi/2}\ln(\sin x\cos x)\,\mathrm dx = \int_0^{\pi/2}\ln\big(\tfrac12\sin2x\big)\,\mathrm dx = \int_0^{\pi/2}\ln\sin2x\,\mathrm dx - \frac\pi2\ln2 .\] In the remaining integral put \(u = 2x\): \(\int_0^{\pi/2}\ln\sin2x\,\mathrm dx = \tfrac12\int_0^{\pi}\ln\sin u\,\mathrm du = \tfrac12\cdot2\int_0^{\pi/2}\ln\sin u\,\mathrm du = I\), using the symmetry of \(\ln\sin u\) about \(u = \pi/2\). Hence \(2I = I - \frac\pi2\ln2\), so \(I = -\frac\pi2\ln2\) and \[\frac1{2\pi}\int_0^{2\pi}\ln|\cos\theta|\,\mathrm d\theta = \frac{1}{2\pi}\cdot4I = \frac{2}{\pi}I = -\ln2,\] where the factor uses that \(|\cos\theta|\) repeats four times over \([0,2\pi)\) on intervals of length \(\pi/2\). Therefore \[\lambda_1 \approx \ln K + \langle\ln|\cos\theta|\rangle = \ln K - \ln2 = \boxed{\ \ln\frac K2\ }\] (The integral converges despite the logarithmic singularities at \(\theta = \pi/2, 3\pi/2\): near a zero of \(\cos\), \(\int_0\ln s\,\mathrm ds\) is finite — unit 00's integrable singularity.)

4. Against the measurements.

• \(K = 6\): \(\ln3 = 1.0986\) against \(1.136\). Low by \(3.3\%\).
• \(K = 10\): \(\ln5 = 1.6094\) against \(1.620\). Low by \(0.7\%\).
• \(K = 2\): \(\ln1 = 0\) against \(\approx0.44\). Wrong by everything.

The error shrinks as \(K\) grows, as an asymptotic formula should. At \(K = 2\) the derivation's premise — \(|K\cos\theta|\gg1\) for most \(\theta\) — is simply false: \(|2\cos\theta|\lt1\) on a third of the circle, and there the neglected \(\delta p_n\) dominates. Worse, at \(K = 2\) a large elliptic island survives around \((\pi,0)\) (P5: elliptic for \(K\lt4\)), so \(\theta_n\) is not uniformly distributed and the average is over the wrong measure. An asymptotic formula quoted outside its range is not a check on anything.

Check. Dimensions: \(\lambda_1\) for a map is per iteration, dimensionless, and so is \(\ln(K/2)\) ✓. Limit: \(K\to2^+\) gives \(\lambda_1\to0\), which would mean the map is regular at \(K = 2\) — it is not, and that is exactly the failure just diagnosed; the formula's zero is at the wrong place ✓ as a warning. A further number: \(K = 100\) predicts \(\ln50 = 3.912\), and the true value, averaged over three starting points, is within \(0.2\%\) of it.

Idea. The apparent paradox — determinism plus continuous dependence, yet unpredictability — dissolves the moment the two statements are written with their quantifiers in order. Everything else here is arithmetic in support of that.

1. The bound. Both solutions satisfy the integral equation of unit 01, \(y(t) = y(0) + \int_0^tf(y(s))\,\mathrm ds\). Subtracting and setting \(u(t) = \lVert y(t)-\tilde y(t)\rVert\), \[u(t) \;\le\; u(0) + \int_0^t\big\lVert f(y(s)) - f(\tilde y(s))\big\rVert\,\mathrm ds \;\le\; u(0) + L\int_0^tu(s)\,\mathrm ds,\] using the triangle inequality for integrals and then the Lipschitz condition. Gronwall's inequality in the form if \(u(t)\le a + L\int_0^tu\) with \(u\ge0\) continuous and \(a,L\ge0\), then \(u(t)\le a\,e^{Lt}\) — proved in unit 01 — gives \(u(t)\le\lVert\delta_0\rVert e^{Lt}\). \(\square\)

2. Equality. Take \(m = 1\) and \(f(y) = Ly\). This is Lipschitz with constant exactly \(L\), since \(|f(y)-f(\tilde y)| = L|y-\tilde y|\). The solutions are \(y = y_0e^{Lt}\), so \(u(t) = |\delta_0|e^{Lt}\) — equality for every \(t\ge0\). Hence no smaller constant works for all Lipschitz \(f\): the bound is sharp as stated.

3. The factor of 53. \(L = \sup\lVert DX\rVert_2\) is a pointwise supremum over the whole accessible region: it bounds the instantaneous stretching rate at the worst point, and Gronwall assumes that worst point is met at every instant. \(\lambda_1\) is a time average: from the variational equation, \(\frac{\mathrm d}{\mathrm dt}\ln\lVert\delta\rVert = \frac{\delta^{\!\top}DX\,\delta}{\lVert\delta\rVert^2} \le\lVert DX\rVert_2\le L\), so \[\lambda_1 = \lim_{T\to\infty}\frac1T\int_0^T\frac{\mathrm d}{\mathrm dt}\ln\lVert\delta\rVert\,\mathrm dt \;\le\; L,\] a mean of quantities each at most \(L\). The measured value cannot exceed \(L\), and in a bounded system it is far below it, because the orbit spends most of its time where the local stretching is small or negative — it must, since the motion is confined and area is preserved. The factor \(53\) is the price of a bound that has to hold for every orbit and every instant.

4. The two propositions.

Continuous dependence (unit 01). \[\forall\varepsilon\gt0\ \ \forall T\gt0\ \ \exists\delta\gt0:\quad \lVert\delta_0\rVert\lt\delta \;\Longrightarrow\; \sup_{0\le t\le T}\lVert y(t)-\tilde y(t)\rVert\lt\varepsilon.\]

Positive exponent. \[\exists\lambda_1\gt0:\quad \lVert\delta(t)\rVert \approx \lVert\delta_0\rVert e^{\lambda_1t} \ \text{ for typical }\delta_0,\ t\ \text{large}.\]

These are consistent, and the proof of consistency is the construction: given \(\varepsilon\) and \(T\), take \(\delta = \varepsilon e^{-LT}\), which is legitimate because \(\delta\) is chosen after \(T\) is given. The quantifier that does the work is the \(\forall T\) standing before the \(\exists\delta\). A contradiction would require \(\exists\delta\ \forall T\), that is, continuity uniform in \(t\) — and unit 01 never asserts that. Indeed part 2's example shows it is false even for the linear equation, which nobody calls chaotic.

5. The numbers. Require \(\varepsilon_0e^{\lambda_1T} = 10^{-2}\) at \(T = 1000\): \[\varepsilon_0 = 10^{-2}e^{-0.044\times1000} = 10^{-2}\times e^{-44} = 10^{-2}\times7.78\times10^{-20} = \mathbf{7.8\times10^{-22}}.\] Double precision carries about \(2.2\times10^{-16}\), nearly \(3\times10^{5}\) times larger: the computation is impossible in double precision, and would be impossible with exact arithmetic too, since no measurement delivers \(8\times10^{-22}\).

The largest \(t\) for which double-precision data suffice for one part in \(100\): \[t_{\max} = \frac1{\lambda_1}\ln\frac{10^{-2}}{2.2\times10^{-16}} = \frac{\ln(4.55\times10^{13})}{0.044} = \frac{31.45}{0.044} = \mathbf{715}\] time units — about \(114\) orbital periods. Each additional decimal digit of initial accuracy buys \(\ln10/\lambda_1 = 52\) further time units, or about eight orbits.

Check. Dimensions: \(L\) and \(\lambda_1\) are both inverse times ✓; \(\varepsilon_0\) and \(\varepsilon\) are the same kind of quantity, so their ratio is dimensionless inside the logarithm ✓. Limit \(\lambda_1\to0\): \(t_{\max}\to\infty\), an integrable system is predictable for ever, as it should be ✓. Special case \(T\to0\) in part 4: \(\delta\to\varepsilon\), the statement becomes trivial, which is right — nothing has had time to go wrong ✓.

Idea. \(\varphi\) is a root of \(x^2-x-1\), so the quantity \(p^2-pq-q^2\) — the norm of \(p - q\varphi\) in \(\mathbb Z[\varphi]\) — is an integer. A nonzero integer is at least \(1\), and that single fact bounds how well \(\varphi\) can be approximated.

1. The norm form. \[(q\varphi-p)(q\bar\varphi-p) = q^2\varphi\bar\varphi - pq(\varphi+\bar\varphi) + p^2 .\] From \(\varphi,\bar\varphi\) being the roots of \(x^2-x-1\): \(\varphi+\bar\varphi = 1\) and \(\varphi\bar\varphi = -1\). Hence the product equals \(p^2 - pq - q^2\), an integer ✓. It vanishes only if \(\varphi = p/q\) or \(\bar\varphi = p/q\) — both impossible for \((p,q)\ne(0,0)\), since \(\varphi\) and \(\bar\varphi\) are irrational (if \(\varphi\) were rational so would \(\sqrt5\) be). For \(q = 0\) it is \(p^2\ne0\). So \(|p^2-pq-q^2|\ge1\).

2. The lower bound. Write \(x = q\varphi - p\). Then \(q\bar\varphi - p = x - q(\varphi-\bar\varphi) = x - q\sqrt5\), so part 1 gives \[1 \;\le\; |x|\cdot|x - q\sqrt5| \;\le\; |x|\big(|x| + q\sqrt5\big).\] If \(|x|\le1\) then \(|x|+q\sqrt5 \le 1 + q\sqrt5 \le q(1+\sqrt5)\) for \(q\ge1\), and therefore \[|q\varphi-p| = |x| \;\ge\; \frac1{q(1+\sqrt5)} . \qquad\square\] (If \(|x|\gt1\) the claimed bound is immediate, since the right side is at most \(1/(1+\sqrt5)\lt1\).)

3. The Diophantine constant. For \(\vec k = (k_1,k_2)\ne0\), \(\vec k\cdot\vec\omega = k_1\varphi + k_2\). Put \(q = |k_1|\), \(p = -k_2\,\mathrm{sgn}(k_1)\) when \(k_1\ne0\); then \(|\vec k\cdot\vec\omega| = |q\varphi-p|\ge\frac1{q(1+\sqrt5)}\) and \(q\le\lVert\vec k\rVert_\infty\), so \[|\vec k\cdot\vec\omega| \;\ge\; \frac{\gamma}{\lVert\vec k\rVert_\infty},\qquad \gamma = \frac1{1+\sqrt5} = \mathbf{0.309017}.\] If \(k_1 = 0\) then \(|\vec k\cdot\vec\omega| = |k_2|\ge1\ge\gamma/\lVert\vec k\rVert_\infty\) ✓. So \(\tau = 1\), the smallest exponent for which any \(\vec\omega\) in two dimensions can work, with a constant of order a third.

4. A Liouville number. Let \(p_N/q_N\) be the truncation of \(\alpha = \sum_{j\ge1}10^{-j!}\) at \(j = N\), so \(q_N = 10^{N!}\) and \(p_N = q_N\sum_{j\le N}10^{-j!}\), an integer. Then \[|q_N\alpha - p_N| = q_N\sum_{j\gt N}10^{-j!} \le q_N\cdot2\cdot10^{-(N+1)!} = 2\cdot10^{N! - (N+1)!} = 2\,q_N^{-N},\] using \((N+1)! - N! = N\cdot N!\). So \(|q_N\alpha - p_N|\le2q_N^{-N}\), which for any fixed \(\tau\) and \(\gamma\gt0\) eventually beats \(\gamma q_N^{-\tau}\), because \(N\to\infty\). Hence no Diophantine condition holds, for any exponent. Prediction: the torus with frequency ratio \(\alpha\) is excluded by KAM for every \(\varepsilon\gt0\), however small — it is as fragile as a torus can be, and should break first.

5. Why the golden circle is last. KAM guarantees survival for tori whose frequencies satisfy the Diophantine condition with \(\gamma\) above a threshold that grows like \(\sqrt\varepsilon\). As the perturbation grows, the threshold rises, and each torus is abandoned by the theorem when its own \(\gamma\) falls below it: the tori are destroyed in increasing order of how badly their frequency ratio is approximated by rationals. Part 3 shows the golden ratio has the largest possible \(\gamma\) at \(\tau = 1\) — Hurwitz's theorem says no irrational does better — and part 4 shows what the other extreme looks like. So the golden circle is the one whose guarantee survives longest, and the numerical fact that the last spanning invariant circle of the standard map has golden rotation number, breaking at \(K_c = 0.9716\), is the arithmetic ordering made visible. The logic is one-way: KAM gives sufficiency, not necessity, so this explains the observation rather than proving it.

Check. Dimensions: everything here is dimensionless, as a frequency ratio must be ✓. Numerical check of part 2 at \(q = 8\): the best rational is \(13/8 = 1.625\) against \(\varphi = 1.618034\), so \(|8\varphi - 13| = 0.05573\); the bound gives \(1/(8\times3.23607) = 0.03863\) ✓, correctly below. At \(q = 13\): \(|13\varphi - 21| = 0.03444\) against the bound \(1/(13\times3.23607) = 0.02377\) ✓. The ratio of truth to bound tends to \((1+\sqrt5)/\sqrt5 = 1.447\), so the constant \(0.309\) is within \(45\%\) of the best possible \(1/\sqrt5 = 0.447\) ✓.

Idea. Two estimates of the same threshold, from opposite directions: a hand criterion that is systematically high, and an experiment that is systematically high for a different reason but whose bias is quantified by a scaling law and can therefore be removed.

Part A. Over one step with \(p\) and \(\Delta\theta\) small, the map is the Euler step of \(\dot\theta = p\), \(\dot p = K\sin\theta\), which is Hamiltonian with \(H_{\text{eff}} = \tfrac12p^2 + K\cos\theta\), since \(-\partial_\theta H_{\text{eff}} = K\sin\theta\). The unstable fixed point of that pendulum is at \(\theta = 0\), \(p = 0\), where \(H_{\text{eff}} = +K\) — the hyperbolic point P5 finds. Its separatrix is the level set \(H_{\text{eff}} = K\); at the stable point \(\theta = \pi\), \[\tfrac12p^2 - K = K \;\Longrightarrow\; p = \pm2\sqrt K,\qquad \Delta p = 2\sqrt K .\] The primary resonances sit at \(p = 2\pi m\), so the spacing is \(\delta p = 2\pi\). Overlap when the two half-widths together span the gap, \(2\Delta p\ge\delta p\): \[4\sqrt K\ge2\pi \;\Longrightarrow\; K \ge \frac{\pi^2}{4} = \mathbf{2.4674}.\] Against \(K_c = 0.9716\) this is high by a factor \(2.54\). It errs high because it counts only the primary resonances: between them lie higher-order ones, each with its own width and its own chaotic layer, and those help bridge the gap at a smaller \(K\). The criterion therefore gives an upper bound in practice, not a two-sided estimate.

Part B1. Below \(K_c\) an invariant circle spans the cylinder and the flux is exactly zero, so no experiment of any length reports transport: the recorded threshold can never be below \(K_c\). Above \(K_c\) the circle becomes a cantorus with tiny gaps, and the flux is positive but vanishes as \(K\downarrow K_c\); a run of length \(N\) sees transport only where the escape time is under \(N\). Lengthening the run reaches smaller fluxes and therefore smaller \(K\), so \(K_N\) decreases towards \(K_c\) — from above, always.

Part B2 — the extrapolation. Model the first-escape time as \(N = C(K_N-K_c)^{-3.01}\). With \(N_1 = 10^{5}\), \(K_1 = 1.020\) and \(N_3 = 10^{7}\), \(K_3 = 0.983\), dividing kills \(C\): \[\frac{K_1-K_c}{K_3-K_c} = \left(\frac{N_3}{N_1}\right)^{1/3.01} = 100^{1/3.01} = 10^{0.66445} = 4.6180 .\] Solve: \(K_1 - K_c = 4.6180(K_3-K_c)\), so \(1.020 - K_c = 4.5395 - 4.6180K_c\), giving \(3.6180\,K_c = 3.5195\) and \[K_c = \mathbf{0.9728}.\] Greene's value is \(0.971635\). The extrapolation is high by \(0.0011\), a relative error of \(1.2\times10^{-3}\) — three digits recovered from an experiment whose raw numbers were wrong in the second.

Part B3 — the spread. Same computation on the other pairs. \(N = 10^{5},10^{6}\): ratio \(10^{1/3.01} = 2.1489\), giving \(K_c = (1.020 - 2.1489\times0.994)/(1-2.1489) = 0.9714\). \(N = 10^{6},10^{7}\): ratio \(2.1489\) again, giving \(K_c = (0.994 - 2.1489\times0.983)/(1-2.1489) = 0.9734\). So the three estimates are \(0.971\), \(0.973\), \(0.973\). Honest quotation: \[K_c = 0.972 \pm 0.001,\] which contains Greene's \(0.9716\). The widest-separated pair is the most reliable, as usual for an extrapolation, and the spread is dominated by the fact that the raw thresholds are known only to three decimals, and that "first escape among \(400\) orbits" is itself a noisy statistic.

Check. Dimensions: \(K\), \(p\) and \(N\) are all dimensionless for a map ✓, and the exponent \(3.01\) is a pure number. Limit: as \(N\to\infty\) the model gives \(K_N\to K_c\) ✓. A consistency check on the exponent: the raw data give \((K_1-K_c)/(K_3-K_c) = (1.020-0.9716)/(0.983-0.9716) = 4.26\) using Greene's value, against the predicted \(4.62\) — agreement to \(8\%\), so the scaling law is doing real work and is not being fitted to itself ✓.

Idea. A resonance is a statement about periods; Kepler's third law turns it into a statement about distances, which is what a telescope measures. No chaos enters the arithmetic — chaos enters only in explaining why those particular distances should be empty.

1. The formula. A \(p\!:\!q\) mean-motion resonance means the asteroid completes \(p\) orbits in the time Jupiter completes \(q\): \(pT_a = qT_{\mathrm J}\), so \(T_a/T_{\mathrm J} = q/p\). Kepler's third law (unit 08), \(T^2\propto a^3\) for bodies orbiting the same central mass, gives \((a/a_{\mathrm J})^{3} = (T_a/T_{\mathrm J})^{2}\), hence \[a = a_{\mathrm J}\left(\frac qp\right)^{2/3}.\] The law holds to the extent that the asteroid's and Jupiter's masses are negligible against the Sun's, which is the approximation being made.

2. The \(3\!:\!1\) and \(5\!:\!2\).

• \(3\!:\!1\): \((1/3)^{2/3} = e^{-\frac23\ln3} = e^{-0.73241} = 0.480750\), so \(a = 5.2044\times0.480750 = \mathbf{2.5020}\) AU, and \(T_a = 11.862/3 = \mathbf{3.954}\) yr.
• \(5\!:\!2\): \((2/5)^{2/3} = e^{-\frac23\ln2.5} = e^{-0.61086} = 0.542884\), so \(a = 5.2044\times0.542884 = \mathbf{2.8254}\) AU, and \(T_a = 11.862\times2/5 = \mathbf{4.745}\) yr.

Against the observed gap centres \(2.50\) and \(2.82\) AU, the agreement is to three significant figures. What this does and does not show: it shows the gaps sit at the resonances, which was Kirkwood's 1866 observation and is arithmetic. It does not show why a resonance is empty; that needs the chaotic growth of eccentricity at the resonance, which is Wisdom's computation and no part of this calculation.

3. Two more.

• \(2\!:\!1\): \((1/2)^{2/3} = e^{-0.46210} = 0.629961\), \(a = \mathbf{3.2786}\) AU, against the observed \(3.27\) ✓, \(T_a = 5.931\) yr.
• \(3\!:\!2\): \((2/3)^{2/3} = e^{-0.27031} = 0.763143\), \(a = \mathbf{3.9717}\) AU, \(T_a = 7.908\) yr.

4. The Hildas. Being at a resonance means the resonant angle is nearly stationary; whether that is dangerous depends on where in the resonance the body sits. A Hilda sits deep inside the resonance's elliptic island, where the angle librates about a value that puts the asteroid far from Jupiter at every conjunction — protected, not endangered. The bodies that are cleared are those near the resonance's hyperbolic separatrix, where the chaotic layer lives and the eccentricity wanders; in the language of this unit, an island centre and a separatrix layer are opposite fates of one resonance.

5. The self-check. In units of AU and years with the Sun's mass, Kepler's third law is \(T = a^{3/2}\). Then \(a_{\mathrm J}^{3/2} = 5.2044^{1.5} = 11.8729\) yr against the observed \(T_{\mathrm J} = 11.862\) yr: a residual of \(0.011\) yr, or \(0.09\%\). That residual measures exactly what the law neglects — chiefly Jupiter's own mass, which enters as \(T^2 = a^3/(1+m_{\mathrm J}/M_\odot)\); with \(m_{\mathrm J}/M_\odot = 9.55\times10^{-4}\) that predicts \(T\) low by a factor \((1+9.55\times10^{-4})^{-1/2}\), i.e. by \(0.048\%\), which accounts for about half of it, the rest being the difference between mean and osculating elements.

Check. Dimensions: \((q/p)^{2/3}\) is a pure number and \(a_{\mathrm J}\) a length ✓; \(T_a = (q/p)T_{\mathrm J}\) is a time ✓. Limit \(p\!:\!q\to1\!:\!1\) gives \(a\to a_{\mathrm J}\) — the co-orbital case, and there sit the Trojans, which are a crowd and not a gap ✓, consistent with part 4. Ordering check: the five resonances \(4\!:\!1, 3\!:\!1, 5\!:\!2, 7\!:\!3, 2\!:\!1\) must give increasing \(a\) as \(q/p\) increases, and \(0.25, 0.333, 0.4, 0.4286, 0.5\) do ✓.

Idea. Every one of these is about the direction of an inference or the budget of a measurement. A section refutes an integral; it never establishes one. An exponent is an asymptotic limit of a linearised quantity; anything measured over a short time, or after saturation, is a different number.

A — curves everywhere. May be concluded: nothing about existence. The implication runs one way only: a second integral forces every orbit onto a curve, so seeing an orbit cover an area refutes it; seeing curves is consistent with an integral and equally consistent with chaotic layers narrower than the plotting resolution, which is exactly the situation at \(E = 1/12\) in this unit's running example, where \(0.1\%\) of the area is chaotic and invisible. May not: "twenty orbits" is a sample, and the surviving tori have empty interior as a set, so no finite sample is representative. Settle it with: a Lyapunov exponent for each of the twenty, run long enough that the regular baseline \(\ln T/T\) is well below the threshold — and, better, a scan on a grid, which is what turns the statement into the percentage.

B — a scattered orbit with \(3\%\) energy drift. May be concluded: nothing whatever. A \(3\%\) drift means the integration is not solving the intended system; the "orbit" is wandering across energy surfaces, and the scatter may be entirely an artefact of that wandering, since the section's accessible region itself depends on \(E\). Settle it with: re-integrate with a smaller step or a symplectic method until the relative drift is \(10^{-7}\) or better, and check that the picture is unchanged; only then is the scatter evidence. This unit's runs stayed below \(10^{-7}\) relative, even over \(T = 10^6\).

C — \(40\%\) above \(0.1\) at \(T = 50\). May be concluded: that \(40\%\) of the section has a finite-time exponent above \(0.1\) at \(T = 50\), which is a real and useful quantity but is not \(\lambda_1\). May not: that \(40\%\) is chaotic. A regular orbit's measured exponent decays like \(\ln T/T\), which at \(T = 50\) is \(\ln50/50 = 0.078\) — within a factor \(1.3\) of the threshold, so regular orbits are certainly being miscounted; and finite-time exponents fluctuate strongly along a single chaotic orbit, so the classification is unstable. Settle it with: the same scan at \(T = 2000\), where the baseline is \(0.0038\), and a check that the answer is unchanged when the threshold is moved by a factor of two.

D — two nearby orbits to separation \(0.4\). May be concluded: that the separation grew by a factor \(4\times10^{7}\) in \(430\) time units, which is genuine evidence of exponential growth. May not: that \(\lambda_1 = 0.0407\) to three digits. Two objections. The final separation \(0.4\) is a large fraction of the accessible region, so the last part of the growth is saturated and nonlinear, which biases the estimate down; and a single stretch of \(430\) units is one sample of a fluctuating quantity — the running example's own \(\lambda_1(T)\) wanders between \(0.037\) and \(0.047\) over \(T\in[4000,20000]\), and settles at \(0.044\) only by \(T = 10^6\). Settle it with: Benettin's algorithm on the variational equation, which never leaves the linear regime, run to \(T\sim10^{6}\), and quote the plateau rather than the final value. As it happens \(0.0407\) is \(7\%\) below this unit's \(0.044\) — the direction saturation predicts — so the method is not useless; it is merely unquotable.

The protocol. (1) Integrate with a step small enough that the relative energy drift over the whole run is below \(10^{-7}\), and report that number. (2) Produce the section from a single orbit per region and report how many crossings it contains, so that "fills an area" is a claim with a count attached. (3) Compute \(\lambda_1\) by Benettin's algorithm on the variational equation, report the plateau value with the horizon \(T\) it was measured over, and check that a regular orbit in the same system returns \(\lambda_1(T)\approx\ln T/T\). Three steps, three numbers: a drift, a count, and an exponent with a budget.

Check. Dimensions: exponents are inverse times throughout; the drift and the fractions are dimensionless ✓. Limit: apply the protocol to an integrable system — step 2 gives curves, step 3 gives \(\lambda_1(T)\) falling like \(\ln T/T\), and the verdict is correctly "no evidence of chaos", which is a different statement from "proved regular" ✓. A number: in this unit's running example the three outputs at \(E = 1/8\) are a drift of \(9\times10^{-8}\) over \(T = 10^6\), \(850\) crossings in the figure, and \(\lambda_1 = 0.044\) over \(T = 10^6\).