Semi-analytic computations of the speed of Arnold diffusion along single resonances in a priori stable Hamiltonian systems

TL;DR

This paper develops a semi-analytic model to predict the extremely slow evolution of action variables along resonances in a priori stable Hamiltonian systems, aligning well with numerical simulations and advancing understanding of Arnold diffusion.

Contribution

It introduces a novel semi-analytic approach combining phase analysis and Melnikov integrals to accurately estimate Arnold diffusion speeds in multi-degree-of-freedom Hamiltonian systems.

Findings

Model accurately predicts Arnold diffusion speed over several orders of magnitude.

Only a small fraction of terms in the Nekhoroshev normal form significantly influence action evolution.

Good agreement between semi-analytic predictions and numerical simulations.

Abstract

Cornerstone models of Physics, from the semi-classical mechanics in atomic and molecular physics to planetary systems, are represented by quasi-integrable Hamiltonian systems. Since Arnold's example, the long-term diffusion in Hamiltonian systems with more than two degrees of freedom has been represented as a slow diffusion within the `Arnold web', an intricate web formed by chaotic trajectories. With modern computers it became possible to perform numerical integrations which reveal this phenomenon for moderately small perturbations. Here we provide a semi-analytic model which predicts the extremely slow-time evolution of the action variables along the resonances of multiplicity one. We base our model on two concepts: (i) By considering a (quasi-)stationary phase approach to the analysis of the Nekhoroshev normal form, we demonstrate that only a small fraction of the terms of the associated optimal remainder provide meaningful contributions to the evolution of the action variables. (ii) We provide rigorous analytical approximations to the Melnikov integrals of terms with stationary or quasi-stationary phase. Applying our model to an example of three degrees of freedom steep Hamiltonian provides the speed of Arnold diffusion, as well as a precise representation of the evolution of the action variables, in very good agreement (over several orders of magnitude) with the numerically computed one.

Tables

Table 1. Table 1: Summary of the numerical experiments on Hamiltonian ( 66 ). The upper table reports the parameters of the Hamiltonian normalizing algorithm and some of the informations that we can extract from its output: Δ ​ S Δ 𝑆 \Delta S denotes the amplitude of the domain in the resonant action S 𝑆 S , 𝒥 𝒥 {\cal J} the truncation order, 𝒥 N subscript 𝒥 𝑁 {\cal J}_{N} the optimal normalization order, ‖ r N ‖ norm superscript 𝑟 𝑁 \left\|r^{N}\right\| the norm of the remainder expansion ( 6 ) close to I ∗ = ( 0.664887 , 0.955495 , 1 ) subscript 𝐼 0.664887 0.955495 1 I_{*}=(0.664887,0.955495,1) , T α subscript 𝑇 𝛼 T_{\alpha} the period of the resonant variables computed using ( 62 ); | Δ ​ F 1 | N ​ e ​ k ​ h subscript Δ subscript 𝐹 1 𝑁 𝑒 𝑘 ℎ |\Delta F_{1}|_{Nekh} represents the a priori upper bound of the maximum variation of F 1 subscript 𝐹 1 F_{1} over a period T α subscript 𝑇 𝛼 T_{\alpha} forced by the remainder r N superscript 𝑟 𝑁 r^{N} . The lower table concerns the numerical computation and the analytic estimates about the variations of the normalized adiabatic action F 1 subscript 𝐹 1 F_{1} during a resonant period: | Δ ​ F 1 | M ​ a ​ x subscript Δ subscript 𝐹 1 𝑀 𝑎 𝑥 |\Delta F_{1}|_{Max} denotes the maximum variation of F 1 subscript 𝐹 1 F_{1} after a full resonant period for a swarm of 100 orbits with initial actions close to I ∗ subscript 𝐼 I_{*} obtained from numerical integrations of the Hamilton equations; | Δ ​ F 1 | N ​ P subscript Δ subscript 𝐹 1 𝑁 𝑃 |\Delta F_{1}|_{NP} denotes the semi-analytic estimate of the maximum variation obtained by computing numerically the Melnikov integral whose phase is stationary or quasi-stationary (since the numerical computation of the Melnikov integrals is more precise than the linear approximation, we include for safety a larger number of terms in the category (III), by checking directly the value δ c subscript 𝛿 𝑐 \delta_{c} for which the integrals are negligible with respect to δ = 0 𝛿 0 \delta=0 ; the number of terms, reported in the column (II)+(III) 1 , is still in a ratio of 1 ∼ 1000 similar-to 1 1000 1\sim 1000 of the total number); | Δ ​ F 1 | P subscript Δ subscript 𝐹 1 𝑃 |\Delta F_{1}|_{P} is analogous to | Δ ​ F 1 | N subscript Δ subscript 𝐹 1 𝑁 |\Delta F_{1}|_{N} , but obtained with the ’patched’ formula ( 68 ); | Δ ​ F 1 | s ​ a subscript Δ subscript 𝐹 1 𝑠 𝑎 |\Delta F_{1}|_{sa} is the value obtained using the asymptotic expansions of Lemmas 1 or 2 ((II)+(III) 2 represents the number of terms included in this computation).

— $\epsilon$ —	$\Delta S$	$𝒥$	${𝒥}_N$	$\Vert r^N\Vert$	$T_{\alpha }$	${\left|\Delta F_1\right|}_{Nekh}$
$\mathrm{\hspace{0.33em}0.08}$	$\mathrm{\hspace{0.33em}0.114}$	$\mathrm{\hspace{0.33em}9}$	$\mathrm{\hspace{0.33em}6}$	$\mathrm{\hspace{0.33em}1.179}\times {10}^{[}-4]$	$\mathrm{\hspace{0.33em}165.0}$	$\mathrm{\hspace{0.33em}1.95}\times {10}^{[}-2]$
$0.05$	$0.090$	$9$	$6$	$3.01\times {10}^{[}-5]$	$265.9$	$7.99\times {10}^{[}-3]$
$0.02$	$0.057$	$10$	$7$	$2.13\times {10}^{[}-6]$	$519.6$	$1.11\times {10}^{[}-3]$
$0.01$	$0.040$	$12$	$9$	$2.07\times {10}^{[}-7]$	$864.2$	$1.78\times {10}^{[}-4]$
$0.008$	$0.036$	$13$	$10$	$8.43\times {10}^{[}-8]$	$1028.4$	$8.67\times {10}^{[}-5]$
$0.005$	$0.029$	$13$	$10$	$1.24\times {10}^{[}-8]$	$1468.5$	$1.82\times {10}^{[}-5]$
$0.002$	$0.018$	$13$	$10$	$2.85\times {10}^{[}-10]$	$2705.0$	$7.70\times {10}^{[}-7]$
$0.001$	$0.013$	$13$	$10$	$3.05\times {10}^{[}-11]$	$4216.4$	$1.29\times {10}^{[}-7]$
$0.0005$	$0.009$	$13$	$10$	$4.01\times {10}^{[}-12]$	$6442.7$	$2.58\times {10}^{[}-8]$
— $\epsilon$ —	${\left|\Delta F_1\right|}_{Max}$	(II)+(III)1	${\left|\Delta F_1\right|}_P$	${\left|\Delta F_1\right|}_{NP}$	(II)+(III)2	${\left|\Delta F_1\right|}_{sa}$
$0.08$	$\mathrm{\hspace{0.33em}5.22}\times {10}^{[}-4]$	$\mathrm{\hspace{0.33em}1334}$	$\mathrm{\hspace{0.33em}3.45}\times {10}^{[}-4]$	$\mathrm{\hspace{0.33em}3.35}\times {10}^{[}-4]$	$\mathrm{\hspace{0.33em}648}\left(83\right)$	$\mathrm{\hspace{0.33em}5.42}\times {10}^{[}-4]$
$0.05$	$1.75\times {10}^{[}-4]$	$1124$	$1.56\times {10}^{[}-4]$	$1.32\times {10}^{[}-4]$	$513\left(48\right)$	$2.28\times {10}^{[}-4]$
$0.02$	$1.0\times {10}^{[}-5]$	$1696$	$1.2\times {10}^{[}-5]$	$5.2\times {10}^{[}-6]$	$601\left(50\right)$	$8.34\times {10}^{[}-5]$
$0.01$	$2.08\times {10}^{[}-7]$	$3838$	$3.6\times {10}^{[}-7]$	$2.1\times {10}^{[}-7]$	$1300\left(202\right)$	$2.91\times {10}^{[}-6]$
$0.008$	$7.0\times {10}^{[}-8]$	$5470$	$4.8\times {10}^{[}-8]$	$2.5\times {10}^{[}-8]$	$1703\left(180\right)$	$4.36\times {10}^{[}-7]$
$0.005$	$1.0\times {10}^{[}-8]$	$6476$	$1.2\times {10}^{[}-8]$	$6.7\times {10}^{[}-9]$	$634\left(201\right)$	$6.79\times {10}^{[}-9]$
$0.002$	$2.36\times {10}^{[}-9]$	$7346$	$4.38\times {10}^{[}-9]$	$2.79\times {10}^{[}-9]$	$916\left(7\right)$	$6.66\times {10}^{[}-9]$
$0.001$	$1.08\times {10}^{[}-9]$	$8350$	$1.32\times {10}^{[}-9]$	$3.87\times {10}^{[}-10]$	$885\left(12\right)$	$7.68\times {10}^{[}-10]$
$0.0005$	$2.12\times {10}^{[}-10]$	$9364$	$1.87\times {10}^{[}-10]$	$1.86\times {10}^{[}-10]$	$836\left(26\right)$	$2.04\times {10}^{[}-10]$

Equations

$$H_{\epsilon }(I,\phi )=H_0(I)+\epsilon f(I,\phi )$$

$$|I(t)-I(0)|\leq\varepsilon^{b}\ \ {\rm for}\ \ |t|\ \leq\ T_{N}:=\frac{1}{\varepsilon}\ \exp\Big{(}{\frac{1}{\varepsilon^{a}}}\Big{)}.$$

$$K_{\epsilon }(I,\phi )=H_{\epsilon }(I,\phi )+{\tau }^{\nu }W(I,\phi ;\epsilon )\text{leavevmode},$$

$${\mathcal{R}}_{\Lambda }=\{I\in A:\ell \cdot \nabla H_0(I)=0,\forall \ell \in \Lambda \},$$

$$(S,F,\sigma ,\varphi )=\mathcal{C}({\Gamma }^{-T}I,\Gamma \phi )$$

$$(\tilde{S},\tilde{F},\tilde{\sigma },\tilde{\varphi })=({\Gamma }^{-T}I,\Gamma \phi )$$

$$\tilde{\mathcal{R}}=\{(\tilde{S},\tilde{F})\in {\Gamma }^{-T}A:\frac{\partial }{\partial {\tilde{S}}_j}h(\tilde{S},\tilde{F})=0,\forall j=1,\dots ,d\}\text{leavevmode}.$$

$$H_{\epsilon ,\Lambda }=h(S,F)+\epsilon f_{\Lambda }(S,F,\sigma )+r_{\Lambda }(S,F,\sigma ,\varphi )\text{leavevmode},$$

$$H^N=h(F,S)+\epsilon f^N(F,S,\sigma )+r^N(F,S,\sigma ,\varphi )\text{leavevmode},$$

$$r^N=\sum _{m\in {\mathbb{N}}^d}\sum _{\nu \in {\mathbb{Z}}^d}\sum _{k\in {\mathbb{Z}}^{n-d}}{r}_{\nu ,k}^m(F-F_{\ast })(S_1-S_{\ast ,1}{)}^{m_1}\cdots (S_d-S_{\ast ,d}{)}^{m_d}{e}^{i\nu \cdot \sigma +ik\cdot \varphi }$$

$$\left|F_j(t)-F_j(0)\right|\le \left|t\right|\Vert \frac{\partial r_{\Lambda }}{\partial {\phi }_j}\Vert \text{leavevmode}.$$

$${\overline{H}}^N=h(F,S)+\epsilon f^N(F,S,\sigma ).$$

$$\frac{\partial h}{\partial S}(S_{\ast },F_{\ast })=0.$$

$${\overline{H}}^N={\overline{H}}_0+...,{\overline{H}}_0={\omega }_{\ast }\cdot \hat{F}+\frac{A}{2}{\hat{S}}^2+\hat{S}B\cdot \hat{F}+\frac{1}{2}C\hat{F}\cdot \hat{F}+\epsilon v(\sigma )\text{leavevmode},$$

$$\hat{S}=\sqrt{\epsilon }{s}_{\alpha }(\sigma ):=\pm \sqrt{\epsilon }\sqrt{\frac{2}{\left|A\right|}(M(1+\alpha )-v(\sigma ))}\text{leavevmode},$$

$$M=\underset{\sigma \in [0,2\pi ]}{\max }v(\sigma ),$$

$$\bar{M}=\underset{\sigma \in [0,2\pi ]}{\min }v(\sigma )\text{leavevmode},$$

$$\Delta F_j(T)=-\sum _{m,\nu ,k}{\int }_0^{T}i{k}_j{r}_{\nu ,k}^m(F(t))\hat{S}(t{)}^m{e}^{i\nu \sigma (t)+ik\cdot \varphi (t)}dt:=\sum _{m,\nu ,k}\Delta F_{j,T}^{m,\nu ,k}\text{leavevmode}\text{leavevmode}.$$

$$-\sum _{m,\nu ,k}{\int }_0^{T}i{k}_j{r}_{\nu ,k}^m(F_{\ast }){{\hat{S}}^0(t)}^m{e}^{i\nu {\sigma }^0(t)+ik\cdot {\varphi }^0(t)}dt\text{leavevmode}\text{leavevmode},$$

$$(F_{\ast },S^0(t),{\sigma }^0(t),{\varphi }^0(t))=(F_{\ast },S_{\ast },0,0)+(0,{\hat{S}}^0(t),{\sigma }^0(t),{\varphi }^0(t))$$

$$\Delta F_{j,T}^{m,\nu ,k}\simeq {\Delta }^0{F}_{j,T}^{m,\nu ,k}:=-i{k}_j\frac{r_{\nu ,k}^m(F_{\ast }){\epsilon }^{\frac{m-1}{2}}}{A}{e}^{ik\cdot \varphi (0)}{\int }_0^{{\sigma }^0(T)}[s_{\alpha }(\sigma ){]}^{m-1}{e}^{i\theta (\sigma )}d\sigma$$

$$\theta (\sigma )=\mathcal{N}\sigma +\frac{\Omega }{A\sqrt{\epsilon }}{\int }_0^{\sigma }\frac{dx}{s_{\alpha }(x)}$$

$$\mathcal{N}=\nu +k\cdot B/A,\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\Omega =k\cdot {\omega }_{\ast }\text{leavevmode}\text{leavevmode}\text{leavevmode}.$$

$$I_{\lambda }={\int }_a^b\eta (\sigma )e^{i\lambda \Phi (\sigma )}d\sigma \text{leavevmode},$$

$$I_{\lambda}\sim{\eta(\sigma)\over\lambda\Phi^{\prime}(\sigma)}e^{i(\lambda\Phi(\sigma))}{\Bigg{|}}_{a}^{b},$$

$$I_{\lambda }\sim \eta ({\sigma }_c)e^{i\left(\lambda \Phi ({\sigma }_c)\pm \frac{\pi }{4}\right)}\sqrt{\frac{2\pi }{\lambda \left|{\Phi }^{\prime \prime }({\sigma }_c)\right|}}\text{leavevmode},$$

$$I_{\lambda }\sim \eta ({\sigma }_c)e^{i\lambda \Phi ({\sigma }_c)}\sqrt{3}\Gamma (4/3){\left(\frac{6}{\lambda \left|{\Phi }^{\prime \prime \prime }({\sigma }_c)\right|}\right)}^{\frac{1}{3}}.$$

$$\mathcal{I}={\int }_0^{{\sigma }^0(T)}[s_{\alpha }(\sigma ){]}^{m-1}{e}^{i\theta (\sigma )}d\sigma$$

$$\frac{d\theta }{d\sigma }=\mathcal{N}+\frac{\Omega }{A\sqrt{\epsilon }}\frac{1}{s_{\alpha }(\sigma )}$$

$$\mathcal{I}({\sigma }_c)=[s_{\alpha }({\sigma }_c){]}^{m-1}{e}^{i\theta ({\sigma }_c)\pm i\frac{\pi }{4}}\frac{\sqrt{2\pi }}{\sqrt{\frac{|\Omega |}{A^2\sqrt{\epsilon }}\frac{|v^{\prime }({\sigma }_c)|}{{|s_{\alpha }({\sigma }_c)|}^3}}}\text{leavevmode},$$

Full text

1 Introduction

Fundamental problems of Physics are often modeled with small Hamiltonian perturbations of integrable systems. For example, the problem of the stability and long–term evolution of the Solar System can be modeled in terms of perturbations to Kepler’s motion of each planet under the gravity of the Sun. Similar perturbative approaches are employed in some of the most classical problems of Mechanics appearing from the microscopic scale (e.g. the semi-classical treatment of atomic and molecular dynamics) up to the astronomical one (e.g. solar systems and galaxies). The above and other important applications, as e.g. in plasma and accelerator physics, or statistical mechanics, have rendered Hamiltonian near-integrable systems a fundamental topic in physics (see [54] for a collection of basic papers and reviews in this field).

One of the most interesting questions in near-integrable systems is the long–term fate of trajectories which belong to the so-called ‘Arnold web’. Following the pioneering work of V.I. Arnold [1], the Arnold web is understood as an intricate in shape and connected set in phase–space which contains chaotic trajectories. The Arnold web is tightly related to the existence, in phase–space, of a corresponding ‘web of resonances’, i.e., domains where the trajectories undergo near-oscillatory motions with a commensurable set of frequencies. ‘Arnold diffusion’ is a theoretically predicted phenomenon, according to which a trajectory with initial conditions within the distorted separatrices of the resonances undergoes slow chaotic diffusion. When the number $n$ of the degrees of freedom is equal to 3 or larger, such diffusion renders possible, in principle, to connect every part of the Arnold web within sufficiently long times. Let us therefore consider a $n$-degree of freedom Hamiltonian of the form:

[TABLE]

where $(I,\varphi)\in A\times{\mathbb{T}}^{n}$ are action-angle variables, $A\subseteq{\mathbb{R}}^{n}$ is open bounded, the integrable approximation $H_{0}$ and the perturbation $f$ are real analytic, $\varepsilon$ is a small parameter. The problem we address in this paper is the following:

PROBLEM 1: For given $H_{0},f$, small $\varepsilon>0$, and $I_{*}\in A$ such that $\ell\cdot\nabla H_{0}(I_{*})=0$ for a unique $\ell\in{\mathbb{Z}}^{n}\backslash 0$ (with its multiples), provide a formula which gives the maximum speed of the drift along the resonance $\ell\cdot\nabla H_{0}(I)=0$ (averaged on time intervals $T$ longer than $1/\varepsilon$) among all the solutions of Hamilton’s equations with initial conditions $I(0),\varphi(0)$ with $\varphi(0)\in{\mathbb{T}}^{n}$ and $I(0)$ in the ball $B(I_{*},C\sqrt{\varepsilon})\subseteq A$ of center $I_{*}$ and radius $C\sqrt{\varepsilon}$, with some $C>0$.

Remarks:

• (i)

For special choices of $H_{0}$, $f$, and of the resonant vector $\ell\in{\mathbb{Z}}^{n}\backslash 0$, the previous problem has a simple solution. In fact, Nekhoroshev provided a class of quasi-integrable Hamiltonian systems with variations of the actions of order 1 already on times of order $1/\varepsilon$ which can be explicitly computed with a simple quadrature.

• (ii)

For $n=2$ and $H_{0}$ iso–energetically non–degenerate, the KAM theorem provides a topological obstruction to the drift along the resonances of the system, so the previous problem is not interesting.

• (iii)

For $H_{0}$ satisfying a transversality condition, called by Nekhoroshev “steepness”, the stability time of the action variables improves dramatically to an exponential order in $1/\varepsilon$ [55, 56]: precisely, there exist positive constants $a$, $b$ and $\varepsilon_{0}$ such that for any $0\leq\varepsilon<\varepsilon_{0}$ the solutions $(I(t),\varphi(t))$ of the Hamilton equations of $H_{\varepsilon}(I,\varphi)$ satisfy

[TABLE]

According to Nekhoroshev’s theorem, any large drift of the action variables needs time intervals longer than the exponentially long–time $T_{N}$. For systems with $n\geq 3$ satisfying the hypotheses of Nekhoroshev’s theorem, proving the existence of orbits with variations of the actions of order 1 in some suitable long time for any small value of $\left|\varepsilon\right|$, is highly non trivial, and these are the conditions under which Problem 1 is interesting and, up to now, unsolved.

• (iv)

We do not provide here a rigorous solution to Problem 1, but we provide formulas which match the very slow drifts observed in numerical experiments. These formulas are obtained by combining: – a semi-analytic argument including the computer assisted computation of normal forms, whose coefficients are provided in floating point arithmetics, – a rigorous approximation of the Melnikov integrals using methods of asymptotic analysis based on the so called stationary-phase approximation, – a random–phase assumption used in the Melnikov approximation. We call our approach ’semi-analytical’, since it combines rigorous results (in the stationary-phase method, see Section 3) with ones based on the numerical (computer-assisted) computations of the Nekhoroshev normal forms. Whether these ideas can be transformed into a fully rigorous argument is a question beyond the purpose of this paper.

• (v)

The estimated speed of drift along a resonance expected from the solution of Problem 1 should depend on $\varepsilon$, $I_{*}$ and $\ell$. Of course, close to $I_{*}$ the resonance $\ell\cdot\nabla H_{0}(I)=0$ may intersect an infinite number of other resonances $\tilde{\ell}\cdot\nabla H_{0}(I)=0$ with $\tilde{\ell}\in{\mathbb{Z}}^{n}\backslash 0$ independent on $\ell$, and solutions with initial conditions close to $I_{*}$ may leave the resonance $\ell\cdot\nabla H_{0}(I)=0$ and drift along different resonances, possibly of different multiplicities. Problem 1 concerns only the orbits which drift along the fixed resonance $\ell\cdot\nabla H_{0}(I)=0$.

• (vi)

For systems of $n=3$ degrees of freedom a solution of Problem 1 gives the opportunity to compare the time needed to diffuse along the resonances of multiplicity 1 with the stability time $T_{N}$ of the Nekhoroshev theorem. In fact for $n=3$ distant points of the action–space on the same energy level are connected through paths of the Arnold web which are mostly contained in resonances of multiplicity one, where Problem 1 is applicable. The resonances of multiplicity two are just at the points of intersection of the resonances of multiplicity one. The transit of the orbits through resonances of multiplicity two, the so–called ’large gap problem’, is one of the hardest theoretical difficulties in rigorously proving the existence of Arnold diffusion. Numerical studies, instead, provide overwhelming evidence for the existence of such transits, see [43, 44], while the key question regarding the quantification of Arnold diffusion is Problem 1.

• (vii)

Throughout this paper, and mostly in Section 5, we compare our semi-analytic solution of Problem 1 with numerical experiments. The long-term behaviour of Hamiltonian systems, including Arnold diffusion, can be numerically investigated with symplectic integrators (see [5, 57, 32]). In fact, depending on the order of the integration scheme, every step $\phi^{\tau}$ of the integrator is exponentially close, with exponential factor $-1/\tau$, to the exact Hamiltonian flow of a modified Hamiltonian

[TABLE]

where the integer $\nu$ and $W$ both depend on the integration scheme. Therefore, for suitably small $\tau$ the spurious term $\tau^{\nu}W(I,\varphi;\varepsilon)$ is just a perturbation of the original Hamiltonian, and the exponential factor becomes negligible with respect to any observed diffusion.

• (viii)

While the KAM and Nekhoroshev theorems, as well as the examples of Arnold diffusion (and also Problem 1), are usually formulated for quasi-integrable hamiltonians (1), many quasi-integrable systems of interest for Physics and Celestial Mechanics are characterized by degeneracies and singularities of the action variables which introduce additional complications. Many researches recovered the proofs of the KAM and Nekhoroshev theorems, as well as of the existence of hyperbolic tori, also for the cases of interest for Celestial Mechanics, see for example, [15, 3, 35, 37, 16, 58, 23].

Problem 1 is an applicative spin-off of the problem of Arnold diffusion, which started with the fundamental paper published by Arnold in 1964 [1], first providing a quasi–integrable Hamiltonian system with non trivial long–term instability. Since Arnold’s pioneering paper, a rich literature has appeared on attempts to prove of existence of Arnold diffusion for more general quasi–integrable Hamiltonian systems, called, in the context of Arnold diffusion, a priori stable systems. A simpler, albeit still highly non trivial, case is the one of a priori unstable systems. In the latter case, the existence of diffusing motions has been proved using different models and techniques, including Mather’s variational methods, geometrical methods and the so-called separatrix and scattering maps (among the rich literature see [15, 7, 20, 62, 18, 6, 33, 47, 46, 21, 11] and references therein).

Due to the long timescales involved, also numerical or experimental observations of Arnold diffusion are hard to achieve. Already few years after the first numerical detection of chaotic motions [45], the long-term instability in Hamiltonian systems was discussed from both an analytical and numerical point of view in [17]. However, only modern computers rendered possible to simulate the phenomenon in simple physical models. In the last decades, diffusion through the resonances has been clearly detected [49, 24, 48, 34, 50, 43, 30, 40, 34, 26, 27, 61, 41, 42]. Then, in the series of papers [50, 43, 30, 40, 42], diffusion of orbits has been also detected for values of the perturbation parameters so small that the set of resonant motions has the structure of the Arnold web embedded in a large volume of invariant tori (the distributions of resonances and tori being computed numerically with chaos indicators [29, 50, 41]). In these experiments, the instability was characterized by diffusion coefficients decreasing faster than power laws in $\varepsilon$, compatibly with the exponential stability result of Nekhoroshev’s theorem. This was confirmed by a direct comparison of the numerical diffusion coefficient with the size of the optimal remainder of the Nekhoroshev normal form in [26].

In this paper we propose a semi-analytic solution to Problem 1 which is obtained through the following steps:

(a) Given $\varepsilon$ and $I_{*}$ we construct a computer assisted normal form adapted to the local resonance properties at $I_{*}$ up to an optimal normalization order, by following the construction of normal forms which appears within the proof of the Nekhoroshev theorem. The computer assisted construction of normal forms is mandatory, since our purpose is to compare the predicted values of the drifts with the numerically observed ones. As it is well known (see [12, 13] for the KAM theorem, and [14, 37] for the Nekhoroshev theorem) purely analytic estimates which do not use computer assisted methods are highly unrealistic. (b) We represent the variation of the actions along the resonance with Melnikov integrals defined from the normal form constructed as indicated in (a). Since the remainders of these normal forms are represented as expansions of millions of very small terms, the variation of the actions is represented as a sum of millions of Melnikov integrals, which have to be computed in order to solve Problem 1. The problem becomes prohibitive if the goal is to maximize the result with respect to some variables in order to compute the orbits with largest instability. To overcome this difficulty, we require an analytic method that allows to descriminate between the terms of the remainder associated with large contribution to the variations of the actions and those of negligible contribution, therefore reducing the total amount of terms to consider.

(c) We represent the Melnikov integrals with a method from asymptotic analysis, the so-called method of the stationary-phase (see [10]). In fact, for quasi-integrable systems, the Melnikov integrals can be reformulated as integrals with a rapidly oscillating phase, and the computation of the critical points of this phase provides an estimate of the integral. We find that only the Melnikov integrals whose phase either 1) has critical points, or 2) the derivative of the phase with respect to the slow angle variable of the resonance is suitably small, provide major contributions to the Arnold diffusion. We call the corresponding terms in the remainder stationary or quasi-stationary, respectively. The Melnikov integrals whose phase is neither stationary nor quasi-stationary represent the large majority of terms, and their cumulative contribution to the Arnold diffusion is negligible with respect to the cumulative contribution of the stationary or quasi-stationary terms. Therefore, we provide a rigorous criterion to select, from the millions of harmonics of the remainder of the Nekhoroshev normal form, a few thousand ones. All the relevant integrals of Melnikov theory can be explicitly represented with an asymptotic formula or directly computed by quadratures. The asymptotic formula, providing the variation of the actions during a resonant libration, depends on the initial phases $\varphi(0)$. For all possible values of these phases the formula represents closely the spread of the actions which is observed with numerical integrations. (d) Finally, by maximizing with respect to $\varphi(0)$ the variations of the actions obtained from the Melnikov integrals we obtain the orbits with largest variation of the actions at each homoclinic loop, as well as the rare initial conditions whose orbit, in a sequence of homoclinic loops, have a systematic variation of the action variables. Thus we predict which orbits undergo the ’fastest’ Arnold diffusion which we can observe.

From (a), (b), (c), (d) we have a qualitative and quantitative description of the drift along the resonances of multiplicity one. The qualitative picture of the diffusion is in agreement with the idea having its roots in Chirikov’s fundamental paper [17] and recently recovered e.g. in [19], namely that the diffusion along a resonance is not uniform in time, but it is produced by impulsive ‘kicks’ or ‘jumps’ at every homoclinic loop, see [2, 22, 60]. The new quantitative analysis allows us to determine the frequency of occurrence and amplitude of these jumps as the resonant angle becomes critical for some Melnikov integrals; also, we are able to select the initial conditions whose orbits have the fastest Arnold diffusion. Therefore, for given values of $\varepsilon$, we are able to predict the minimum timescales needed to observe long–term diffusion along any single resonance.

The paper is organized as follows. In Section 2 we define the Melnikov integrals from the normal forms of Nekhoroshev theory. In Section 3 we provide rigorous asymptotic representations of the Melnikov integrals using the method of stationary phase approximation. In Section 4 we present a semi-analytic solution to Problem 1. Section 5 is devoted to a numerical demonstration of the theory presented in Sections 2, 3, 4.

2 Nekhoroshev normal forms and Melnikov integrals

The long–term dynamics of the quasi–integrable Hamiltonian (1) is traditionally studied using the averaging method. In the refined version of the method defined within the proof of Nekhoroshev’s theorem, for a $d$–dimensional lattice $\Lambda\subseteq{\mathbb{Z}}^{n}$ defining the resonance

[TABLE]

one constructs a canonical transformation

[TABLE]

defined in a suitable resonant domain $D_{\Lambda}\times{\mathbb{T}}^{n}$, where

$\Gamma$ is a matrix with $\Gamma_{ij}\in{\mathbb{Z}}$ and $\det\Gamma=1$, that defines a linear canonical map (see [4])

[TABLE]

and conjugates $H_{0}(I)$ to $h(\tilde{S},\tilde{F})$ such that the resonance ${\cal R}_{\Lambda}$ is transformed into

[TABLE]

Equivalently, in the new variables the resonant lattice $\Lambda$ is transformed into the lattice generated by $e_{1},\ldots,e_{d}$, ($e_{1},\ldots,e_{n}$ denotes the canonical basis of ${\mathbb{R}}^{n}$).

• -

$\sigma\in{\mathbb{T}}^{d},\phi\in{\mathbb{T}}^{n-d}$ are angles conjugate to the actions $S\in{\mathbb{R}}^{d},F\in{\mathbb{R}}^{n-d}$.

• -

$D_{\Lambda}$ is a domain whose definition depends both on the resonant lattice $\Lambda$ and on $\varepsilon$.

• -

${\cal C}$ is a near to the identity transformation and, when composed with (3), conjugates $H_{\varepsilon}$ to the Nekhoroshev normal form Hamiltonian

[TABLE]

where the remainder $r_{\Lambda}$ has norm bounded by a factor exponentially small with respect to $-1/\varepsilon^{a}$ (see [55, 59, 52, 53, 39] for precise definitions and statements).

The integer $d\in{1,\ldots,n-1}$ is called the multiplicity of the resonance.

Although the proof of Nekhoroshev’s theorem grants the existence of normal forms (4) for suitably small values $0\leq\varepsilon\leq\varepsilon_{0}$ (see (2)), the precise value of the threshold parameter $\varepsilon_{0}$ is believed to be largely underestimated by the general proofs, while the Fourier coefficients of $f_{\Lambda},r_{\Lambda}$ are estimated in norm, but not explicitly provided. Both problems can be overcomed by constructing the normal forms (4) with computer assisted methods [35, 14, 36, 26]. We call Hamiltonian normalizing algorithm (HNA) a computer-algebraic implementation which provides the coefficients of the Nekhoroshev normal form (4). We use the HNA introduced in [26], which normalizes quasi-integrable Hamiltonians $H_{\varepsilon}$. The HNA is constructed by composing $N$ elementary transformations; the input of the HNA is the Hamilton function, a resonance lattice $\Lambda$, a domain $D\times{\mathbb{T}}^{n}$ where the transformation if defined; the output of the HNA is a canonical transformation $(F,S,\sigma,\phi)={\cal C}^{N}(I,\varphi)$ and a normal form Hamiltonian

[TABLE]

conjugate to $H_{\varepsilon}$ by ${\cal C}^{N}$. The remainder $r^{N}$ is provided as a Taylor-Fourier series

[TABLE]

expanded at a suitable $(F_{*},S_{*})$, with computer-evaluated truncations involving a large number (typically $\sim 10^{7},10^{8}$) of terms. $N$ is chosen so that $\left\|r^{1}\right\|>\ldots>\left\|r^{N}\right\|$ and $\left\|r^{N+1}\right\|>\left\|r^{N}\right\|$, thus the normal form is called optimal, $r^{N}$ the optimal remainder and $N$ the optimal normalization order (the norm definitions in the selected domain are as in [26]).

If we artificially suppress the remainder $r_{\Lambda}$ in Eq. (4), or $r^{N}$ in (5), we obtain an exponentially small perturbation of the original Hamiltonian which possibly exhibits chaotic motions due to homoclinic and heteroclinic phenomena (for $d>1$), but in which the actions $F_{j}$, which we call ’adiabatic’, remain constant in time. Therefore, in the flow of the complete Hamiltonian, any long–term evolution of the adiabatic actions is due to the accumulation of the effects of the very small remainder on very large times. In particular, the adiabatic actions $F_{j}$ have a long-term variation, representing the drift along the resonance, bounded for an exponentially long time by

[TABLE]

The a priori estimate (7) obtained from the Nekhoroshev normal form (4) provides an upper bound to the average variation of the adiabatic actions; establishing lower bounds to $\left|F_{j}(t)-F_{j}(0)\right|$ is the fundamental brick in the theory of Arnold diffusion.

From now on we focus our discussion on resonances of multiplicity $d=1$. Denoting by $(S,\sigma)$ the resonant action-angle pair, we first consider the dynamics of the approximated normal form which is obtained from (5) just by dropping the remainder $r^{N}(F,S,\sigma,\phi)$:

[TABLE]

Since the corresponding Hamiltonian $\overline{H}^{N}$ depends only on one angle, it is integrable, and we represent its motions as follows. Following [4], we first expand $\overline{H}^{N}$ at $(F_{*},S_{*})$ identifying the center of the resonance, precisely such that

[TABLE]

We obtain

[TABLE]

where $\hat{F}=F-F_{*}$, $\hat{S}=S-S_{*}$, and the Hamiltonian $\overline{H}_{0}$ is represented using a number $A\in{\mathbb{R}}$, two vectors $\omega_{*},B\in{\mathbb{R}}^{n-1}$, a square matrix $C$ and a function $v(\sigma)$, all these quantities depending parametrically on $S_{*},F_{*}$.

The actions $\hat{F}$ are constants of motion for the Hamiltonian flow of $\overline{H}^{N}$ as well as of $\overline{H}_{0}$. We parameterize by

[TABLE]

the level curves of $\overline{H}_{0}$ for $\hat{F}=0$, where

[TABLE]

and $\alpha$ is a convenient label for the energy levels of $\overline{H}_{0}$. For any $\alpha\neq 0$, we denote by $T_{\alpha}$ the period of the corresponding solutions of Hamilton’s equations under $\overline{H}_{0}$ (for $\hat{F}=0$). We also define

[TABLE]

and, without loss of generality, we assume $\bar{M}\leq 0$, $M>0$.

When considering the solutions $(F(t),S(t),\sigma(t),\phi(t))$ of Hamilton’s equations under the complete Hamiltonian (5), the adiabatic actions $F_{j}$ can have a slow evolution forced by the remainder $r^{N}$, whose variation $\Delta F_{j}(T)=F_{j}(t)-F_{j}(0)$ in a time interval $[0,T]$ is given by

[TABLE]

According to the well known Melnikov approach (see [17] for a review) we approximate $\Delta F_{j}(T)$ with

[TABLE]

obtained by replacing the solution $(F(t),S(t),\sigma(t),\phi(t))$ in the integrals with

[TABLE]

where $(0,\hat{S}^{0}(t),\sigma^{0}(t),\phi^{0}(t))$ is a fixed solution of Hamilton’s equations under $\overline{H}_{0}$ (see remark (ix) below). Finally, by changing the integration variable from $t$ to $\sigma=\sigma^{0}(t)$, we have

[TABLE]

where the phase $\theta(\sigma)$ is defined by:

[TABLE]

with

[TABLE]

Remark (ix). Usually, Melnikov approximations are introduced to compute the splittings of stable-unstable manifolds, so that integrals like (12) are approximated by choosing $(\hat{S}^{0}(t),\sigma^{0}(t),\phi^{0}(t))$ to be the solution of the approximate normal form $\overline{H}_{0}$ corresponding to a separatrix homoclinic loop ($\alpha=0$ in our notation). Our method exposed below differs from the usual Melnikov approach since, in order to find the orbits which diffuse in shorter time along the resonance, we evaluate the integrals for a solution $(\hat{S}^{0}(t),\sigma^{0}(t),\phi^{0}(t))$ which is suitably close to, but not exactly on the separatrix, precisely $\alpha\sim\left\|r^{N}\right\|$, with finite period $T_{\alpha}$.

3 Stationary phase approximation of Melnikov integrals

3.1 The principle of stationary phase

While formulas (12) or (13) allow, in principle, to compute the time evolution of the adiabatic actions during consecutive homoclinic loops along the resonance, the evaluation of the sums over millions of remainder terms $r^{m}_{\nu,k}$ is hardly tractable in practice. We find that most of these terms (including some of the largest in norm) contribute very little to the sum (13). This fact can be explained by invoking methods of asymptotic analysis inspired by the so–called principle of the stationary phase (PSP hereafter). In its classical formulation (e.g., see [10]) the principle concerns the asymptotic behaviour of the parametric integrals

[TABLE]

when the parameter $\lambda$ is large. With mild conditions on the amplitude function $\eta(\sigma)$, we have the following cases:

• (A)

The phase $\Phi(\sigma)$ has no stationary points, i.e. $\Phi(\sigma)\neq 0$ for all $\sigma\in[a,b]$, for large $\lambda$ we have

[TABLE]

where the neglected contributions are of order smaller than $1/\lambda$.

• (B)

The phase $\Phi(\sigma)$ has a non–degenerate stationary point $\sigma_{c}\in(a,b)$, i.e. $\Phi^{\prime}(\sigma_{c})=0$ and $\Phi^{\prime\prime}(\sigma_{c})\neq 0$. Then, if $\eta(\sigma_{c})\neq 0$, for large $\lambda$ we have

[TABLE]

where the $\pm$ is chosen according to the sign of $\Phi^{\prime\prime}(\sigma_{c})$. If there are more stationary points in $(a,b)$, we must sum all the corresponding terms.

• (C)

The phase $\Phi(\sigma)$ has a degenerate stationary point $\sigma_{c}\in(a,b)$, i.e. $\Phi^{\prime}(\sigma_{c})=0,\Phi^{\prime\prime}(\sigma_{c})=0$ and we assume $\Phi^{\prime\prime\prime}(\sigma_{c})\neq 0$. Then, if $\eta(\sigma_{c})\neq 0$, for large $\lambda$ we have

[TABLE]

3.2 Heuristic discussion of the PSP for Melnikov integrals

Let us consider the Melnikov integrals (14)

[TABLE]

and identify the oscillating phase with $\theta(\sigma):=\lambda\Phi(\sigma)$. Since the derivative of $\theta(\sigma)$

[TABLE]

is, in principle, divergent for $\varepsilon$ going to zero, we will justify in the next sub-section the use of the method of the stationary phase by showing how to identify the large parameter $\lambda$ in terms of specific parameters entering into the calculus of (20). According to this idea, the Melnikov integrals whose oscillating phase $\theta(\sigma)$ has critical points are expected to be dominant over those whose oscillating phase has no critical points. However, we find that the asymptotic behaviour of the Melnikov integrals is more complicated than the behaviour of the integrals (16) thus rendering necessary to use a refinement of the stationary phase method: Specifically, we need to consider also a case which is intermediate between (A) and (B), called hereafter the quasi-stationary case, produced by the disappearance of couples of non-degenerate critical points $\sigma^{1}_{c},\sigma^{2}_{c}$ after they merge into a degenerate critical point. The quasi-stationary case represents a transition between the stationary and the non-stationary case, which is not considered in the usual formulations of the PSP method. To be more precise, depending on the values of $m,\nu,k$, we will consider three cases:

• (I)

the phase $\theta(\sigma)$ is ’fast’ for all $\sigma\in[0,2\pi]$, i.e. we have $\left|\theta^{\prime}(\sigma)\right|>\gamma$, with $\gamma$ a large parameter to be defined later.

In this case the integral in (14) is estimated smaller than order $\lambda^{-1}$, see Eq. (17), and we will assume that the contribution of $\Delta F_{j,T}^{m,\nu,k}(T)$ to the series expansion in Eq. (12) can be neglected.

• (II)

the phase $\theta(\sigma)$ is ’slow’ only close to non-degenerate critical points $\sigma_{c}$, provided they are distant enough with respect to $1/\lambda^{1\over 3}$.

In this case, by invoking the PSP (see (18)), the integrals $\int_{0}^{\sigma}[s_{\alpha}(\sigma)]^{m-1}e^{i\theta(\sigma)}d\sigma$ are approximated by the sum on all the stationary points $\sigma_{c}\in[0,\sigma]$ of contributions ${\cal I}(\sigma_{c})$ defined by

[TABLE]

where the $\pm$ depends on the sign of $v^{\prime}(\sigma_{c})$. The above formula is valid if $\theta^{\prime\prime}(\sigma_{c})\neq 0$. Consequently, we obtain (see Lemma 1):

[TABLE]

• (III)

the phase $\theta(\sigma)$ is quasi-stationary, i.e. $\left|\theta^{\prime}(\sigma)\right|\leq\gamma$ in an interval of size of order $1/\lambda^{1\over 3}$ (notice that this condition can occur in absence of critical points, or in presence of two very close non-degenerate critical points or of one degenerate critical point). While in this case we cannot directly apply (A), (B), (C), we will obtain an asymptotic formula for the integrals stemming from formula (C) (see Lemma 2).

Remarks:

• (x)

from estimates (A) and (B), any individual integral estimated using (18) is of order $1/\sqrt{\lambda}$, larger with respect to the integrals estimated using (17) which are of order $1/\lambda$. Moreover, in the non–stationary case (I) the integrals in (17) are estimated only by the difference of a function computed at the border values $a,b$. Since Arnold diffusion is produced by the variations of $F_{j}$ through a sequence of circulations or librations, the border values of the sequence may cancel (as a matter of fact they may cancel only partially, since from a libration/circulation to the next one there can be small variations of $\alpha$ and a change of the fast phases $\phi(0)$ in the factor multiplying the integral in (14)).

• (xi)

The practical classification of the integrals in one of the categories (I), (II), or (III) will be done by a fast algorithmic criterion (see below), based only on each term’s integer labels $m,\nu,k$. Since for the large majority of $\nu,k$ the phase satisfies (I) (see Table 1 of the numerical examples), we have a criterion to select the few harmonics ($\sim$few in 1000) belonging to (II) and (III), hence, producing the dominant terms in the time evolution of $F_{j}$.

• (xii)

Despite being more complicated than (II), the inclusion in the computation of the quasi–stationary terms (III) is essential, since these individual contributions can be as large as those of (II) and quite often we find algebraic near-cancellations between terms of the groups (II) and (III) leaving residuals of order only few percent of the absolute values of the corresponding Melnikov integrals.

3.3 Rigorous discussion of the PSP for Melnikov integrals

In this subsection we state the results which allow us to assign all the terms in Eq. (14), labeled by the integers $m,\nu\in{\mathbb{Z}}$ and $k\in{\mathbb{Z}}^{n-1}$, into the categories (I), (II) or (III). We then provide rigorous asymptotic representations for the integrals of the terms in (II) and (III). To simplify the discussion, we assume $\alpha>0$; the modifications needed to represent the case $\alpha<0$ are straightforward.

We first analyze the conditions on $m,\nu,k$ which imply that the corresponding term has a phase with stationary points. To simplify the analysis we assume mild conditions on the potential $v(\sigma)$ in the normal form Hamiltonian (10). The first hypothesis is that $v(\sigma)$ has only one point of maximum. Up to a translation of the angle $\sigma$, we assume it at $\sigma=0$, so that with the notations introduced in Section 2, we have $v(0)=M$, implying

[TABLE]

with the sign $\pm$ chosen according to the signs of $s_{\alpha}(\sigma)$ and $A$. We now have the following: Lemma 1. Consider a phase $\theta(\sigma)$ defined by the labels $m,\nu,k$. Then:

If ${\cal N}\cdot{\cal W}>0$ the phase $\theta(\sigma)$ has no stationary points;

• -

If ${\cal N}\cdot{\cal W}<0$, the phase $\theta(\sigma)$ has stationary points if and only if

[TABLE]

Furthermore, suppose that $v(\sigma)$ has only one non-degenerate local maximum at $\sigma=0$, one non-degenerate local minimum at $\sigma=\bar{\sigma}$, and $v^{\prime}(\sigma)\neq 0$ elsewhere. For any given $\Delta_{Max}>\sqrt{1-{\bar{M}\over M}}$, consider $\varepsilon$ suitably small and

[TABLE]

Then

• (i)

If the inequality (24) is strictly satisfied, the phase $\theta(\sigma)$ has two non–degenerate critical points $\sigma^{1}_{c},\sigma^{2}_{c}$ and for all ${\cal N}=-{\cal W}\Delta$ with $\Delta\in\left({\sqrt{1-{\bar{M}\over M}}\over\sqrt{1+\alpha-{\bar{M}\over M}}},\Delta_{Max}\right)$, we have

[TABLE]

where ${\cal I}(\sigma)$ has been defined in (21), and the functions $a_{1}$, $a_{2}$ satisfy

[TABLE]

• (ii)

If ${|{\cal N}|\over{|{\cal W}|}}$ is equal to the lowermost bound of Eq.(24), the two non-degenerate critical points merge into one degenerate critical point $\sigma_{c}=\bar{\sigma}$.

Proof of Lemma 1. Since we have

[TABLE]

if ${\cal N}\cdot{\cal W}>0$ there are no stationary points. If ${\cal N}\cdot{\cal W}<0$ the stationary points are the solutions of the equation

[TABLE]

which exist if and only if ${|{\cal N}|/{|{\cal W}|}}$ satisfies (24). By assumption, $v(\sigma)$ has only one local maximum $\sigma=0$ and one local minimum $\bar{\sigma}$. If (24) is strictly satisfied, the function

[TABLE]

has a strict maximum at $\sigma=\bar{\sigma}$ with ${\tilde{\theta}}(\bar{\sigma})>0$ and converges to $-1\over\sqrt{\alpha}$ for $\sigma$ tending to [math] or $2\pi$. Therefore there are two values $\sigma^{1}_{c},\sigma^{2}_{c}\in(0,2\pi)\backslash\{\bar{\sigma}\}$ such that ${\tilde{\theta}}^{\prime}(\sigma_{c}^{i})=0$. Then, from

[TABLE]

we have $\theta^{\prime\prime}(\sigma^{i}_{c})\neq 0$, and consequently the two critical points are non-degenerate. Instead, when the inequality (24) is satisfied at its lower extremum, we have only one critical point $\sigma_{c}=\bar{\sigma}$, which is degenerate.

It remains to prove Eq. (26). With no loss of generality, we consider the case ${\cal N}>0,{\cal W}<0$. Setting ${\cal N}$ within the phase $\theta(\sigma)$ with ${\cal N}=-{\cal W}\Delta=\Delta\left|\cal W\right|$, we obtain $\theta(\sigma)=\left|\cal W\right|\Phi(\sigma)$ with

[TABLE]

Since, depending on the values of $m$, the integral

[TABLE]

may not be smooth at $\sigma=0,2\pi$, we use the technique called neutralization of the extremals. Precisely, choose a small $\mu>0$, depending possibly on the given $\Delta_{Max}$, but independent of $\left|\cal W\right|$ and $\Delta$. In the following we denote by $k_{1},k_{2},\ldots$ suitable constants which do not depend on $m,{\cal W},\Delta,\varepsilon$, while they may depend on $\mu,\Delta_{Max}$.

We first prove that in the hypothesis of the Lemma there exists a small $\mu$ such that both critical points $\sigma^{j}_{c}$ are in $(\mu,2\pi-\mu)$, and for any $\sigma\in[0,\mu]$, we have

[TABLE]

In fact, for the given $\Delta,{\cal W}$, the critical points $\sigma_{c}^{i}$ satisfy

[TABLE]

or equivalently

[TABLE]

Choosing $\mu$ to satisfy

[TABLE]

and using (25), we obtain

[TABLE]

and therefore we have $\sigma^{i}_{c}\in(\mu,2\pi-\mu)$. Then, for all $\sigma\in[0,\mu]$, we have

[TABLE]

as soon as

[TABLE]

which is satisfied if

[TABLE]

From (25) and (29) we obtain

[TABLE]

Therefore, for $\mu$ satisfying (29), for all $\sigma\in[0,\mu]$ we have

[TABLE]

Analogously, we have

[TABLE]

[TABLE]

Let us now consider an infinitely differentiable function $\rho(x,\mu)$ such that $\rho(x,\mu)=0$ for $x\leq\mu/2$ and $x\geq 2\pi-\mu/2$, $\rho(x,\mu)=1$ for $x\in[\mu,2\pi-\mu]$, ${\partial^{j}\over\partial x^{j}}\rho(\mu/2,\mu)={\partial^{j}\over\partial x^{j}}\rho(2\pi-\mu/2,\mu)$ for all $j\geq 0$. Then, we define

[TABLE]

and

[TABLE]

The integrand of $\hat{I}(\left|\cal W\right|,\Delta)$ is smooth and bounded for all $m\geq 0$, and vanishes at the extrema together with all its derivatives. Therefore it has the form suitable for the application of the rigorous version of PSP (see, for example, [10]). As a consequence, taking into account that the phase $\Phi(\sigma)$ has two non-degenerate critical points $\sigma_{c}^{1},\sigma_{c}^{2}$, $\hat{I}(\left|\cal W\right|,\Delta)$ is represented by (26) with $a_{1},a_{2}$ satisfying the limits (27). It remains to estimate the integral

[TABLE]

[TABLE]

We prove that there exists a constant $\kappa$ independent on $\left|\cal W\right|$, $\Delta$ and $m$, such that

[TABLE]

and therefore also the integral $I(\left|W\right|,\Delta)$ has the representation (26) with $a_{1},a_{2}$ satisfying the limits (27).

Since the phase $\Phi(\sigma)$ has no stationary points in $[0,\mu]$, and if $m\geq 1$, integrating by parts we obtain

[TABLE]

We consider the following cases:

If $m\geq 3$ or $m=1$, using (28) we obtain

[TABLE]

Then we estimate

[TABLE]

[TABLE]

Using again (28) we obtain

[TABLE]

Since $1/\Phi^{\prime}(\sigma)$ is strictly monotone in $[0,\mu]$, its derivative has the same sign in $[0,\mu]$, and therefore we have

[TABLE]

[TABLE]

By collecting all these estimates, we obtain

[TABLE]

with $k_{5}>0$ independent on $\left|\cal W\right|$, $\Delta$ and $m$.

• -

If $m=2$, we estimate the border contribution in (32), and first and the third integral in (33) as in the case $m\geq 3$. It remains to estimate the second integral, whose denominator is only apparently divergent, since because of (28) we have $\left|\Phi^{\prime}(\sigma)s_{\alpha}(\sigma)\right|\geq k_{1}$. Therefore we have

[TABLE]

• -

If $m=0$, for any arbitrary small $\xi\in(0,\mu)$, we have

[TABLE]

[TABLE]

Using (28) we obtain, uniformly on $\xi$,

[TABLE]

Since $1/(\Phi^{\prime}(\sigma)s_{\alpha}(\sigma))$, is strictly monotone in $[0,\mu]$, by proceeding as in the cases $m\geq 1$, we obtain

[TABLE]

uniformly in $\xi$, so that

[TABLE]

By repeating the argument to estimate the integral on $[2\pi-\mu,2\pi]$ we obtain (31). $\square$

Lemma 1 provides an explicit criterion allowing to classify Melnikov integrals as belonging to the categories (I) or (II) depending on the values of of ${\cal N},{\cal W}$. Precisely, if ${\cal N}\cdot{\cal W}>0$ the phase is considered in the category (I), and the contribution of the corresponding Melnikov integral to the Arnold diffusion will be considered negligible. Instead, if ${\cal N}\cdot{\cal W}<0$, and $|{\cal N}|/|{\cal W}|$ satisfies strictly (24), the phase has two non-degenerate critical points. Then we distinguish two subcases:

• • $|{\cal N}|/|{\cal W}|$ is not too close to its lower extremum, according to to a criterion specified by Lemma 2 below. Then, the phase is considered in the category (II) and the contribution of the corresponding Melnikov integral to the Arnold diffusion can be estimated analytically (26).
• • ${\cal N}\cdot{\cal W}<0$ and ${|{\cal N}|\over|{\cal W}|}$ suitably close to its lower extremum. We find that such a term, while formally ’stationary’, contributes to the Melnikov integral similarly as ’quasi-stationary’ terms satisfying

[TABLE]

In fact, a careful investigation of the transition of ${|{\cal N}|\over|{\cal W}|}$ from values higher than $Q_{\alpha}$ to smaller ones, reveals that the transition corresponds to a degenerate critical point for the phase. The corresponding integral blows to values of order $1/\left|\cal W\right|^{1\over 3}$, at the transition, and depending linearly on the distance $\delta=\frac{{\cal N}}{{\cal W}}-Q_{\alpha}$, for small $|\delta|$. Thus, these intermediate cases will be considered in the quasi-stationary category (III).

The values of the Melnikov integrals for terms in category (III) are estimated according to the following: Lemma 2. Let the potential $v(\sigma)$ have only one local non–degenerate maximum at $\sigma=0$ and one local non–degenerate minimum at $\sigma=\bar{\sigma}$. Let us consider $\varepsilon$ suitably small and $\alpha\in(0,\left\|r^{N}\right\|)$. For any phase $\theta(\sigma)$ defined by the labels $m,\nu,k$ such that ${\cal N}>0$, ${\cal W}<0$ and

[TABLE]

with some $\delta>0$, and defining

[TABLE]

we have

[TABLE]

where

[TABLE]

with constants $\kappa$ and $c_{j}$ independent on $\left|W\right|$ and $\delta$, and $\theta_{*}=\theta(\bar{\sigma})_{|\delta=0}$. In particular, we have

[TABLE]

and, for all $j\geq 2$,

[TABLE]

*with $\alpha_{j}\in\{0,1,\sqrt{3},2\}$ depending on $j$. * Proof of Lemma 2. Let us choose $\mu>0$ small enough, but independent on $\left|\cal W\right|$, $\delta$ and $\alpha$; let us define an infinitely differentiable function $\rho(x,\mu)$ such that $\rho(x,\mu)=0$ for $x\leq 0$ and $x\geq 2\pi$; $\rho(x,\mu)=1$ for $x\in[\mu,2\pi-\mu]$; ${\partial^{j}\over\partial x^{j}}\rho(0,\mu)={\partial^{j}\over\partial x^{j}}\rho(2\pi,\mu)$ for all $j\geq 0$. Then, we define

[TABLE]

and

[TABLE]

Let us preliminarly write the phase $\theta(\sigma)$ and its derivative by replacing ${\cal N}$ using (35):

[TABLE]

as well as the expansions at $\sigma=\bar{\sigma}$

[TABLE]

In the following we denote by $k_{1},k_{2},\ldots$ suitable constants which do not depend on $m,{\cal W},\delta,\varepsilon,\alpha$, while they may depend on $\mu$. Since (for small $\mu$) $\bar{\sigma}\in(\mu,2\pi-\mu)$, for all $\sigma\in[0,\mu]$ we have (see (46))

[TABLE]

as well as

[TABLE]

To estimate

[TABLE]

[TABLE]

we first notice that the phase $\theta(\sigma)$ has no stationary points in $[0,\mu]$ and therefore integrating by parts we obtain

[TABLE]

We consider the following cases:

If $m\geq 3$ or $m=1$, using (48) we obtain

[TABLE]

Then we can estimate

[TABLE]

[TABLE]

Using (48) we obtain

[TABLE]

Since $1/\theta^{\prime}(\sigma)$ is strictly monotone in $[0,\mu]$, its derivative has the same sign in $[0,\mu]$, and therefore we have

[TABLE]

[TABLE]

By collecting all these estimates, we obtain

[TABLE]

• -

If $m=2$, we estimate the border contribution in (50), and first and the third integral in (51) as in the case $m\geq 3$. It remains to estimate the second integral, whose denominator is only apparently divergent, since because of (49) we have $\left|\theta^{\prime}(\sigma)s_{\alpha}(\sigma)\right|\geq\left|\cal W\right|k_{1}$. Therefore we have

[TABLE]

• -

If $m=0$, for any arbitrary small $\xi\in(0,\mu)$, we have

[TABLE]

[TABLE]

Using (49) we obtain uniformly on $\xi$

[TABLE]

Since $1/(\Phi^{\prime}(\sigma)s_{\alpha}(\sigma))$, is strictly monotone in $[0,\mu]$, by proceeding as in the cases $m\geq 1$, we obtain

[TABLE]

uniformly in $\xi$, so that

[TABLE]

By repeating the argument to estimate the integral on $[2\pi-\mu,2\pi]$ we obtain that there exists a constant $\kappa$ independent on $\left|\cal W\right|$ and $\delta$ such that $\left|b(\left|\cal W\right|,\delta)\right|\leq\kappa$.

Then, since the integral $\hat{I}(\left|\cal W\right|,\delta)$ is smooth with respect to $\delta$ at $\delta=0$, we have

[TABLE]

where

[TABLE]

Therefore, by using the degenerate version of the principle of stationary-phase (see [28]), by identifying in $\left|\cal W\right|$ the large parameter, and by considering that $\eta(\sigma)(\sigma-\bar{\sigma})^{j}$ vanishes with all its derivatives at $\sigma=0,2\pi$, we obtain (38), (40), (43), (44). $\square$

Lemma 2 allows us to study the transition in the representations of the Melnikov integrals from the regime of stationary phases to the regime of non stationary phases. Hence:

In the case $\delta\geq 0$, a non-stationary phase is considered quasi-stationary, and the corresponding Melnikov integral is approximated by

[TABLE]

if

[TABLE]

otherwise the phase is considered non-stationary and the corresponding Melnikov integral is neglected.

• -

In the case $\delta<0$, since $I(\left|\cal W\right|,\delta)$ is smooth in $\delta$, at $\delta=0$, we compare two estimates: one coming from Lemma 1 (stationary phase approximation) and another coming from the extension of the linear law (52) (quasi–stationary phase approximation) to small negative $\delta$. We find that, for negative $\delta$ suitably close to [math] the quasi-stationary phase approximation provides a better estimate with respect to the stationary phase approximation. To determine a threshold to decide which one to use, we compared the numerical computation of the integrals with the estimates provided by both Lemmas. Let us, for example, consider $v(\sigma)=\cos\sigma$; for $\alpha=0$, we have

[TABLE]

and the integrals

[TABLE]

Figure 1 shows the values of the integrals $\Delta{\cal I}$, computed numerically, for several values of ${\cal W}<0$, and fixed $\delta$ (left panel), or for fixed ${\cal W}$ and several values of $\delta$ (right panel). The left panel shows that the value of the integrals as computed numerically by solving Eq. (55) (blue dots) is well approximated by the corresponding asymptotic law $1/\left|\cal W\right|^{1\over 3}$ for the values of $|{\cal W}|$ considered (blue line). In the right-panel we compare the numerical computations of (55) (blue dots), with the corresponding estimate provided by the stationary-phase approximation (green curve) and with the linear law (52) (red line), for the sample value $\left|W\right|=15$ (very similar pictures are obtained for different values). We see that the stationary phase estimates reproduce well the values of the integrals for $\delta\leq-\delta_{c}/2$. For $\delta\in[-\delta_{c}/2,0]$ we have a divergence of the stationary phase approximation formula, indicating that the approximation is no more valid since we are entering the regime of quasi-stationary phase. In fact, we observe that the linear law (52) represents much better the value of the integral for both positive and negative $\delta$ in the interval $-\delta_{c}/2<\delta<\delta_{c}$, and therefore, we use Eq. (52) with $c_{0}$, $c_{1}$ given by Lemma 2 also to estimate those integrals. By using the formula down to $\delta=\delta_{c}$ we introduce some errors, which could be reduced by considering the non linear corrections (see formula (56) below). On the other hand, $\delta_{c}$ as computed from (52) (represented by the point at which the linear law crosses the x–axis) is clearly underestimated, since the non-linear contributions determine that $\Delta{\cal I}$ has a tail extending only asymptotically to zero (see remark (xiii)). Remarks:

• (xiii)

The non-linear terms of the expansion (52) provide corrections to the critical value $\delta_{c}$ necessary to discriminate between quasi-stationary or non stationary phase. Since $c_{2}=0$ and $c_{3}<0$, a more careful analysis of the non-linear terms provides

[TABLE]

so that, for smaller values of $\left|\cal W\right|$ the quadratic and the cubic terms in $\delta$ can produce a small variation of $\delta_{c}$.

• (xiv)

The above analysis applies to the terms with suitably large values of

[TABLE]

where $\left|k\cdot\omega_{*}\right|$ are non-resonant divisors. In Nekhoroshev theorem, also in the more generic steep case, if $\left|k\right|\leq N$ these divisors are estimated by (see for example, the discussion in the introduction of [39])

[TABLE]

so that, correspondingly, $\left|\cal W\right|$ is large. For $\left|k\right|\geq N$, we have no theoretical lower bounds on $\left|W\right|$. However, by analyticity, the corresponding Fourier harmonics in the original perturbation are exponentially small in $\left|k\right|$, and hence also exponentially small in $N$. Thus, after the normalization procedure we do not expect to find such harmonics in the dominant terms of $r^{N}$.

• (xv)

Lemma 1 and 2 are derived by considering the upper branch of the separatrix solution $\theta(\sigma)$ and $\alpha>0$. Equivalent results are found for the lower branch, and for $\alpha<0$, after some obvious modifications in the formulas.

4 A semi-analytic solution to Problem 1

Following the analytical results of section 3, the semi-analytic solution to Problem 1 that we provide in this paper goes through the following steps (we assume here that, for the integer vector $\ell$, the function $v(\sigma)$ satisfies the hypotheses of Lemmas 1 and 2, otherwise obvious modifications apply).

Semi-analytic representation of $\Delta F_{j}$ during a resonant libration. On the basis of Lemma 1 and Lemma 2 we first define the algorithm which approximates the variation of the adiabatic actions $\Delta F_{j}$ in the time interval $[0,T_{\alpha}]$ ($t=0$ is chosen so that $\sigma(0)=0$ and $T_{\alpha}$ is the circulation period of the dynamics of the resonant normal forms, see Section 2), with the function (see (12))

[TABLE]

[TABLE]

where the coefficients $f_{j,m,\nu,k}$ are provided as floating point numbers obtained by replacing the integral in

[TABLE]

according to the following fast algorithm:

For any $m,\nu,k$ such that $r^{m}_{\nu,k}(F_{*})\neq 0$ compute $\Omega$, ${\cal N}$ and ${\cal W}$ (see (15), (54)).

• -

If ${\cal N}\cdot{\cal W}>0$, or if $|\Omega|>1$, set $f_{j,m,\nu,k}=0$. If ${\cal N}\cdot{\cal W}<0$ and if $|\Omega|\leq 1$, check if condition (24) and $\delta<-\delta_{c}/2$ are satisfied. In such a case, compute $f_{j,m,\nu,k}$ by replacing the integral with its asymptotic expression as indicated in Section 3, Lemma 1.

• -

If $-\delta_{c}/2\leq\delta<\delta_{c}$ compute $f_{j,m,\nu,k}$ by replacing the integral with its asymptotic expression as indicated in Eq. (52), with the coefficients $c_{0}$ and $c_{1}$ given in Section 3, Lemma 2. Otherwise, if $\delta\geq\delta_{c}$ set $f_{j,m,\nu,k}=0$.

Randomization of the phases, a refinement of equation (58). The above representation is obtained by first approximating the integrals in (12) with Melnikov integrals, and then by computing the Melnikov integrals using Lemmas 1 and 2. Then, we describe the long–term diffusion of the actions caused by a sequence of resonant circulations by applying iteratively formula (58) and by updating the values of the phases $\phi(0)$ at the beginning of each circulation, assuming a random variation. In fact, during any resonant circulation, the angles $\phi(t)$ deviate from the approximation considered in the Melnikov integrals, i.e. $k\cdot\phi(t)\sim k\cdot\phi(0)+k\cdot\Omega t$. Since the dynamics is chaotic and the phases $k\cdot\phi(t)$ are fast, we expect a random deviation from the linear approximation $k\cdot\Omega t$ to appear during a circulation period. The small errors introduced by the randomization of the phases during a circulation period is reduced if we split the period $[0,T_{\alpha}]$ into two time intervals:

[TABLE]

where $\sigma(\bar{T}_{\alpha})=\bar{\sigma}$, and we compute two semi-analytic formulas:

[TABLE]

representing the change of the adiabatic actions in the first part of the resonant circulation ($\sigma(t)\in[0,\bar{\sigma}]$) and in the second part ($\sigma(t)\in[\bar{\sigma},2\pi]$) respectively. The value of the phases $\phi$ are then updated also when $\sigma=\bar{\sigma}$.

The reason for this improvement is the symmetry of the distribution of the critical points with respect to the minimum $\bar{\sigma}$, so that the change of the actions is really split into two well differenciated parts, the first one taking place before $\bar{\sigma}$ and the second one taking place after $\bar{\sigma}$.

Largest $\Delta F_{j}$ during a resonant libration. From the numerical values of the coefficients $f_{j,m,\nu,k},f^{i}_{j,m,\nu,k}$, one can estimate the maximum variation $\Delta F_{j}$ during a resonant libration by computing the maximum of the function

[TABLE]

with respect to all the possible initial values of the phases $\phi(0)$. We remark that to obtain results which compare to the numerical experiments we are not allowed to replace the maximum of the series (60) with the value of the majorant series $\sum_{m,\nu,k}\left|f_{j,m,\nu,k}\right|$. In fact, there are examples (see Section 5) where the majorant series is one order of magnitude larger than (60).

Orbits with the fastest long-term instability, ballistic diffusion. The long-term instability of an orbit may arise from a sequence of circulations/librations of $S,\sigma$, which produce very small jumps of $F_{j}$ and $\alpha$, while the phases $\phi$ are treated as random variables. Since $\Delta F_{1},\Delta\alpha$ are very small at each step, their variations along several circulations/librations are mainly determined by the values of the phases $\phi$ at the beginning of each circulation/libration. Random variation of the phases yields a random walk of $F_{1}$ and, by selecting an initial condition such that the values of the phases at each circulation/libration produce the maximum $\Delta F_{1}$, we obtain a monotonic ballistic motion along the resonance. The conditions to observe these ballistic motions from swarms of ${\cal K}$ diffusive orbits are determined as follows: by assuming a randomization of the phases occurring at each resonant libration (random phase approximation), half of the orbits will have $\Delta F_{1}\geq 0$ and the other half $\Delta F_{1}<0$, within the range determined by $\left|\Delta F_{1}(T_{\alpha})\right|$ computed as indicated in the previous step. Therefore, we observe orbits with $\Delta F_{1}$ of the same sign for a number of ${\cal M}$ randomizations as soon as ${\cal K}/2^{\cal M}\geq 1$. Correspondingly, given ${\cal K}$, we observe orbits with $F_{1}$ which increases (or decreases) almost monotonically in time for a time interval $(\log{\cal K}/\log 2)T_{\alpha}$. By denoting with $10^{-p}$ the precision of the numerical integration, this time interval is bounded by $(p/\log 2)T_{\alpha}$.

The speed of the ballistic diffusion in the sequence of resonant librations is represented by

[TABLE]

where we estimate the variation $\left|\Delta F^{(i)}_{j}\right|$ occurring at the $i$–th step by the maximum value computed using the semi-analytic theory previously indicated. The sum of the libration periods $T_{\alpha^{(i)}}$ is instead estimated by assuming an average period needed by any libration

[TABLE]

with $\alpha=\left\|r^{N}\right\|$. In the case $v(\sigma)=\varepsilon M\cos\sigma$, we find

[TABLE]

obtained when the energy of the libration differs from the separatrix value by the norm of the remainder. Therefore, we obtain the formula for the average speed of the ballistic diffusion as

[TABLE]

Numerical computation of $\Delta F_{j}(T)$ during a resonant libration. The analytic formulas provided above allow us to compute the maximum speed of diffusion along the resonance. If we are also interested in following, for any given value of the phase $\phi(0)$ at the beginning of the resonant libration, the individual variation $\Delta F_{j}(T)$ for all $T\in[0,T_{\alpha}]$, it is possible to compute numerically the function

[TABLE]

where the coefficients

[TABLE]

are obtained, only for the terms in the category (II) or (III), by evaluating numerically the integrals in (65). Since the terms in the category (II) or (III) are just  1/1000 of terms of the remainder, the numerical computation of all these integrals is well within the possibility of modern computers.

Remark: (xvi) Formula (63) has been obtained from the analysis of the optimal normal form constructed as indicated in the Nekhoroshev theorem. Therefore it represents an improvement of the a priori estimate obtained from the same normal form, for the diffusion along the resonances of multiplicity 1. In the examples of Section 5 the improvement is of some orders of magnitude.

5 Numerical demonstrations on a three degrees of freedom steep

Hamiltonian model

We illustrate our theory for the 3-degrees of freedom Hamiltonian (introduced following [40, 61])

[TABLE]

satisfying the hypotheses of the Nekhoroshev theorem ($H_{0}$ is steep and the perturbation is analytic), for the resonance defined by the integer vector $\ell=(1,1,0)$, and by $I_{*}=(21\pi/100,3\pi/10,1)$. The upper value of $\varepsilon$. Using the method of the Fast Lyapunov Indicator (see [29, 41], FLI hereafter) we preliminary checked that for the largest value of $\varepsilon$ that we considered in our experiments, i.e. $\varepsilon=0.08$, the resonance ${\cal R}_{\ell}$ close to $I_{*}$ is embedded in a domain dominated by regular motions, with the other resonances forming a web, a circumstance ensuring that the diffusion occurs mainly along the resonance ${\cal R}_{\ell}$. Evidently, the condition persists for smaller values of $\varepsilon$.

Computation of the normal form using a HNA. For a sample of values of $\varepsilon$ we computed the normal form of Hamiltonian (66) by implementing the HNA described in [26]. Following the notations of Lemma 2, for the selected resonance determined by $\ell=(1,1,0)$, we preliminary define the canonical transformation

[TABLE]

[TABLE]

with

[TABLE]

and then, implementing the HNA, we obtain a canonical transformation

[TABLE]

conjugating the Hamiltonian to the normal form (5)

[TABLE]

with optimal normalization order $N$ depending on the specific value of $\varepsilon$. For all the details about the HNA we refer the reader to [26]. Nevertheless we provide below some details about the output of the algorithm for the case treated in this paper.

Truncation order, optimal normalization order, optimal reminder. Since the HNA is implemented on a computer algebra system, any function $Z(S,F,\sigma,\phi)$ is stored in the memory of the computer as a Taylor–Fourier expansion defined by its series of terms111Recall that for the specific Hamiltonian (66) the action $F_{2}=I_{3}$ appears only as an isolated linear term.

[TABLE]

truncated to some suitably large truncation order. To define the truncation order, as well as other orders within the algorithm, the series is modified by multiplying each term by

[TABLE]

where $\xi$ is a formal parameter (which at the end of the computation will be set equal to 1), and $\mu$ is defined so that the perturbation is analytic in the complex domain $\{\varphi:\ \left|\Im\varphi_{j}\right|\leq\mu\}$. Then we represent the modified series ${\cal Z}(S,F,\sigma,\phi,\xi)$ obtained in this way as a Taylor expansion with respect to the parameter $\xi$

[TABLE]

truncated at some suitable order ${\cal J}$. The truncation order of $Z$ is decided as the truncation order of the Taylor expansion of ${\cal Z}$ with respect to the parameter $\xi$.

The expansions in the formal parameter $\xi$ are used also to define the optimal normalization order. In fact, if we consider all the intermediate Hamiltonians which are constructed within the algorithm

[TABLE]

any remainder $r^{i}$ has a truncated Taylor expansion in the formal parameter

[TABLE]

starting from a minimum order ${\cal J}_{i}$ such that ${\cal J}_{i+1}={\cal J}_{i}+1$. The optimal value of $N$ is then chosen so that

[TABLE]

A necessary condition for the correct execution of the algorithm is that the truncation order ${\cal J}$ is larger than ${\cal J}_{N}$. Therefore, the practical limitation for its implementation is due to the limited memory of the computer to store all the series expansions required by the HNA to work within the truncation order. Since the optimal normalization order increases as $\varepsilon$ decreases, for any given computer memory we have a lower bound on the value of $\varepsilon$ such that we are able to construct the normal form Hamiltonian. For the practical pourpose of this work, we considered a lower bound of $\varepsilon=0.0005$.

• -

Domains of the normal forms. In order to solve Problem 1 we need to provide an estimate of the norms of the normal form remainders (computed as the series of the absolute values of the Taylor-Fourier coefficients) in a domain of the actions $S,F$ which is bounded, in principle, by order $\sqrt{\varepsilon}$. The numerical bounds of the action variables are chosen, for each value of $\varepsilon$, according to the amplitude of the separatrices of the resonant motions.

• -

Estimates on the canonical transformation. The canonical transformation $\tilde{\cal C}$ (see (67)) is near to the identity, and in particular the difference $\left|F_{j}-\tilde{F}_{j}\right|$ can be uniformly bounded by $\varepsilon^{b}$ (with some $b>0$ defined as in (2)) which is a quantity much larger than the norm of the optimal remainder $r^{N}$. As a consequence, even if we suppress from the normal form the remainder $r^{N}$, so that the normalized actions $F_{j}$ are constants of motion, the non-normalized actions have a variation of order $\varepsilon^{b}$ which cannot be ascribed to the Arnold diffusion which is instead produced by a variation of the normalized actions $\tilde{F}_{j}$ (the so–called ’deformation’ in Nekhoroshev theory).

In Table 1 we summarize the values of the orders of truncation and of optimal normalization, as well as the norm of the optimal remainders, for a sample of values of $\varepsilon$ from $\varepsilon=0.08$ down to $\varepsilon=0.0005$. The computations were performed with double floating point precision for the largest values of $\varepsilon$, and with quadruple floating point precision for the smaller ones. The CPU time required by the execution of the HNA on a modern fast multi-processor workstation ranges from few minutes for $\varepsilon=0.08$ to some hours for $\varepsilon=0.0005$. We notice that the norm of the optimal remainder spans 9 orders of magnitude in this range of variation of $\varepsilon$.

To provide an idea of the efficiency of the normalizing transformations, in Fig. 2 we compare the time evolution of $F_{1}$ with the time evolution of $\tilde{F}_{1}$ for a swarm of solutions with initial conditions in a small neighborhood of the separatrix of the resonant normal form, for $\varepsilon=0.01$. The solutions $(\tilde{S}(t),\tilde{F}(t),\tilde{\sigma}(t),\tilde{\phi}(t))$ have been obtained from a numerical integration of Hamilton’s equations of the original Hamiltonian (66); the evolution of the adiabatic action $F_{1}(t)$ has been obtained by transforming the numerical solution with the canonical transformation $\tilde{\cal C}$: $(S(t),F(t),\sigma(t),\phi(t))=\tilde{\cal C}(\tilde{S}(t),\tilde{F}(t),\tilde{\sigma}(t),\tilde{\phi}(t))$. In the left panel we see that the variation of $\tilde{F}_{1}$ produces a swarm of points rapidly oscillating in a band of width $6\times 10^{-3}$, which is due to the terms of order $\varepsilon^{b}$ which bound $\left|F_{1}-\tilde{F}_{1}\right|$. A totally different picture appears in the right panel, where the variation of the normalized action $F_{1}$ is represented: in this case the slow time evolution is well defined, characterized by jumps of order $10^{-7}$ (typical values of long–term diffusion of the action variables for this value of $\varepsilon$), that are detectable on such time intervals only thanks to the implementation of the normalizing transformation.

Before applying the theory developed in Sections 2, 3, 4, we provide as in the proof of Nekhoroshev theorem an upper bound to the variation of the action variables by computing the right-hand side of inequality (7). The upper bound computed for a period (62) is reported in the column $\left|\Delta F_{1}\right|_{Nekh}$ and is larger up to two order of magnitudes with the numerically computed variations $\left|\Delta F_{1}\right|_{Max}$. We therefore proceed by estimating these variations with the Melnikov integrals.

Estimate of $\Delta F_{1}$ during a resonant libration. Let us analyze more in detail the variation of the adiabatic action $F_{1}$. In Fig. 3, as before, we represent the time evolution of $F_{1}(t)$ during a circulation of the variables $\sigma,S$ obtained from a numerical integration of Hamilton’s equations of (66) for a swarm of 100 orbits, for two sample values of $\varepsilon$. The spread of $F_{1}(t)$ after the circulation is due to the different values of $\phi(0)$. We are now able to predict the time evolution of all these orbits by using the semi-analytic theory developed in Section 3.

Since, due to the discrimination between phases, the number of Melnikov integrals to take into account is now small, we have the opportunity to compute these integrals also numerically for all the intermediate times $t\in[0,T_{\alpha}]$. For these computations, we can safely extend the value of $\delta_{c}$ computed from the linear approximation as soon as the phases with $\delta>0$ provide non negligible contributions, and still have a small number of terms (see Table 1, column (II)+(III)1).

The red curves of Fig. 3 represent the orbits yielding the maximum negative jump obtained for the numerical integration of the Hamilton equations, while the black and blue curves represent the Melnikov approximations (without and with the patched formula (68), respectively). One sees that for both values of $\varepsilon$ all the curves are sticked up to a time corresponding approximately to half a period of a complete homoclinic loop. In the middle of the homoclinic loop, we distinguish two cases. In the first case the jump is due mostly to remainder terms which become locally stationary at angles $\sigma_{c}$ sufficiently far from $\bar{\sigma}=\pi$, while the slope $d\theta/d\sigma$ is substantially larger than unity at $\sigma=\pi$. In such cases, the jumps are localized around the two stationary values symmetric with respect to the middle of the loop, while the associated remainder terms yield a rapid oscillatory evolution of the actions $F_{1}$ in between the two jumps. Since the motion is in reality chaotic, the orbits during the rapid oscillations undergo also a randomization of the phases, implying that the predictions obtained by computing (14) may introduce an error. This can be remedied using both representations (59): precisely, the blue curve represents the ‘patched’ evolution given by:

[TABLE]

On the other hand, in cases where important quasi-stationary terms enter into play, the phase $\theta(\sigma)$ remains at small values over a large interval around $\sigma=\pi$. Then no rapid oscillations of the fast variables are observed, and the variations become predictable along the whole homoclinic loop using the original estimate (57). In fact, these are cases where the method illustrates its full power, as it is able to capture large cancellations taking place between stationary terms (II), which, however, exhibit near-stationarity in the whole interval between the two (symmetric with respect to $\pi$) critical values $\sigma_{c}$, and true quasi-stationary terms (III). An example is provided in Fig 4: the terms of groups (II) and (III) independently produce jumps of order $10^{-5}$, which nearly cancel, leaving a residual of order $10^{-7}$ which fits exactly the numerical evolution of the action $F_{1}$. Since no rapid oscillations are observed in the middle of the homoclinic loop, the non patched estimate is more precise than the patched estimate, as also shown in the right panel of Fig. 3.

In Table 1, for several different values of $\varepsilon$, we report the values $|\Delta F_{1}|_{Max}$ representing the maximum variation of $F_{1}$ in a full resonant libration, obtained from the numerical integration of the Hamilton equations for a swarm of 100 orbits with initial actions close to $I_{*}$; $|\Delta F_{1}|_{NP}$ denotes the semi-analytic estimate of the maximum variation obtained by computing numerically the Melnikov integral whose phase is stationary or quasi-stationary; $|\Delta F_{1}|_{P}$ is analogous to $|\Delta F_{1}|_{N}$, but obtained with the patched formula (68); $|\Delta F_{1}|_{sa}$ is the value obtained using the asymptotic expansions of Lemmas 1 or 2. By comparing $\left|\Delta F_{1}\right|_{Max}$ with $|\Delta F_{1}|_{N},|\Delta F_{1}|_{NP}$ we have a good agreement between the numerical integrations and the predictive model for all the values of $\varepsilon$ (we notice that for a given $\varepsilon$, only one of the two values $|\Delta F_{1}|_{N},|\Delta F_{1}|_{NP}$ is applicable), to within a factor 2 in variations over 6 orders of magnitude as $\varepsilon$ varies between $0.0005$ and $0.08$). The values $|\Delta F_{1}|_{sa}$ are expected to be slightly less precise than $|\Delta F_{1}|_{N},|\Delta F_{1}|_{NP}$, since they rely on the linear law (52) for the quasi-stationary cases, and do not take into account the patched formula (68); we expect that the errors can be more important for larger values of $\varepsilon$. Here we have an agreement within a factor 3 as $\varepsilon$ varies between 0.0005 and 0.08, except in the interval $0.008\leq\varepsilon\leq 0.02$, where the cancellations (as in Fig. 4) become important.

Diffusion and ballistic orbits. As discussed in Section 4, the long-term instability of an orbit may arise from a sequence of circulations/librations of $S,\sigma$, which produce very small jumps of $F_{1}$ and $\alpha$, while the phases $\phi$ are treated as random variables. The random variation of the phases determines the random walk along the resonance in jumps of maximum amplitude estimated according to the theory of Section 3; for special initial conditions the sequence of jumps has the same sign, so that we have the orbits which move along the resonance with the largest speed (ballistic orbits). An illustration of this phenomenon is represented in Fig. 5 where we represent a ballistic orbit through a sequence of 14 circulations, which is the limit of the quadruple precision. The speed of the ballistic orbits numerically measured is in agreement with formula (63) (see Table 1). Note also the overall random walk nature of the jumps $\Delta F_{1}$ for most other orbits nearby to the ballistic one. Since estimates on $\Delta F_{j}$ can be regarded as providing the one–step size in the random walk, they are crucial in modelling the diffusion process for a large measure of trajectories over times of practical interest in the applications.

Acknowledgments. This research has been supported by: ERC project 677793 Stable and Chaotic Motions in the Planetary Problem and the Research Committee of the Academy of Athens Grant 200/854, during non-overlapping periods (R.I.P.); project CPDA149421/14 of the University of Padova “New Asymptotic Aspects of Hamiltonian Perturbation Theory” (M.G.). C. Efthymiopoulos acknowledges the hospitality of the University of Padova, under the program Visiting Scientist.