On the chaotic diffusion in multidimensional Hamiltonian systems

We present numerical evidence that diffusion in the herein studied multidimensional near-integrable Hamiltonian systems departs from a normal process, at least for realistic timescales. Therefore, the derivation of a diffusion coefficient from a linear fit on the variance evolution of the unperturbed integrals fails. We review some topics on diffusion in the Arnold Hamiltonian and yield numerical and theoretical arguments to show that in the examples we considered, a standard coefficient would not provide a good estimation of the speed of diffusion. However, numerical experiments concerning diffusion would provide reliable information about the stability of the motion within chaotic regions of the phase space. In this direction, we present an extension of previous results concerning the dynamical structure of the Laplace resonance in Gliese-876 planetary system considering variations of the orbital parameters accordingly to the error introduced by the radial velocity determination. We found that a slight variation of the eccentricity of planet c would destabilize the inner region of the resonance that, though chaotic, shows stable when adopting the best fit values for the parameters.


Introduction
An issue of relevance in several fields of astronomy, such as planetary and galactic dynamics among many others, is whether chaos could lead to large variations of the orbital parameters. Indeed, the Solar System as well as almost all the discovered exoplanet systems are chaotic as revealed in Laskar (1990) and more recently in Gayon et al. (2008), Rivera et al. (2010), Deck et al. (2012), Martí et al. (2013), Barnes et al. (2015), among others. Also in most models of elliptical and disk galaxies, the amount of chaos turns out to be significant as already shown many years ago in for instance Merritt and Friedman (1996), Merritt and Valluri (1996), Papaphilippou and Laskar (1998), Contopoulos and Grosbol (1986) and Contopoulos (2009). Though chaos is often associated with large instabilities, it should be noticed that this is not necessarily the case, as the so-called stable chaos observed in Milani and Nobili (1992) illustrates.
The drift of the unperturbed integrals or orbital elements due to a chaotic dynamics is known as chaotic diffusion, and in general, this instability is a global property of nearintegrable Hamiltonian systems with more than two degrees of freedom. The latter problem, particularly for large instabilities, is still open as it was pointed out in Lochak (1999) almost two decades ago. However, progress in this direction can be found for instance in Delshams and Shaefer (2017) and references therein. Thus, we briefly summarize below what is actually well known about these intriguing phenomena.
A first point to be stressed is that local exponential divergence of nearby orbits (i.e., a positive Lyapunov characteristic number) does not necessarily imply chaotic diffusion. (We refer again the example of the observed stable chaos.) In general, "fast diffusion" occurs when an overlap of several resonances takes place, as shown in Chirikov (1979) for lowdimensional systems like the standard map for a large value of the perturbation parameter [this issue was previously discussed in Izrailev and Chirikov (1968)] or in Wisdom (1980) for the planar-circular restricted three-body problem in asteroidal dynamics. The overlap criterion requires, in general, that the perturbation exceeds some critical value; probably, the first report of this instability condition is given in Chirikov (1966). So much so that usually a system is said to be in Chirikov's regime when most of the invariant tori are destroyed by the overlap of resonances, thus giving rise to large chaotic domains. In such a case, diffusion is assumed to be fast. Conversely, a system is said to be in Nekhoroshev's regime when chaos (or any instability) is confined to the thin layers surrounding resonances (Nekhoroshev 1977). KAM theory is then required, i.e., the size of the perturbation should be small enough, and from Nekhoroshev theorem the timescale for any instability is exponentially large. Let us notice that though both KAM theory and Nekhoroshev estimates are rigorous, they only yield upper bounds for stability conditions and for the speed of the rather slow diffusion along the narrow chaotic layers, the so-called Arnold diffusion [see for instance Froeschlé et al. (2005) or (2006) for a detailed discussion]. On the other hand, Chirikov's approach (Chirikov 1979) (Ch79 hereafter) though heuristic provides constructive ways to compute local and global diffusion coefficients in both scenarios, characterized either by fast or slow diffusion like in the standard map for large values of the parameter or along the very thin chaotic layer of a single resonance, respectively. Therefore, strictly speaking, it seems not completely correct to refer to Chirikov's regime as a strongly chaotic one.
Studies on diffusion in multidimensional near-integrable Hamiltonian systems were conducted by many authors, most of them concerning the slow diffusion regime in the limit of weak chaos where KAM and Nekhoroshev theories apply; for example, we refer to Lega et al. (2003Lega et al. ( , 2008, Froeschlé et al. (2005Froeschlé et al. ( , 2006, Guzzo et al. (2005), Efthymiopoulos and Harsoula (2013) and Cincotta et al. (2014). There, the effort was devoted to show both local and global diffusion for very large timescales along the narrow chaotic layers of simple 3 degree of freedom Hamiltonian systems or 4D-maps and to establish the relation between the obtained diffusion coefficient with the perturbation parameter of the model. Only a few works deal with diffusion in the strong chaotic regime for multidimensional realistic systems. In Tsiganis et al. (2005) and Tsiganis (2008), the general diffusion theory is revisited, with applications to the chaotic dynamics of asteroids and Jupiter's Trojans, while in Cachucho et al. (2010) Chirikov's theory is used to investigate diffusion in the asteroidal three-body resonance (5, − 2, − 2). More recently in Batygin et al. (2015) and Martí et al. (2016), the dynamics and diffusion in the Laplace resonance in the GJ-876 planetary system is discussed. Mathematical studies about the Arnold mechanism and diffusion, by means of geometrical methods in different scenarios of the three-body problem, are presented in Delshams et al. (2016), Capiński et al. (2017) and Féjoz et al. (2016), while by recourse of variational methods are addressed in for instance Bessi (1997), Bernard et al. (2016) and references therein.
Let us recall that the derivation of a diffusion coefficient in near-integrable Hamiltonian systems (or in systems with divided phase space) is still an open problem. In fact, in all the above-mentioned works the Brownian motion or normal diffusion approximation is assumed and, in spite of our systematic and profuse search, we failed to find any other approach to estimate a measure of diffusion in the literature. Thus, the aim of the present work is to review the past and present efforts in this direction and study the limits of application of normal diffusion approximation in a simple model.
In the present work, we revisit and discuss the nature of diffusion in a well-known 2 1 2 degrees of freedom Hamiltonian, the Arnold Hamiltonian introduced in Arnold (1964), which is the paradigmatic model for the so-called Arnold diffusion (actually the Arnold mechanism that leads to an instability) but for perturbation parameters sufficiently large so that no analytical estimates apply. We present evidence that for the adopted values of the parameters and considered timescales, the diffusion is far from normal and thus the derivation of a diffusion coefficient from the variance evolution of the actions does not seem to be well sustained in case of large instabilities. Anyway, we show that diffusion experiments would serve to discriminate stable and unstable regions within a chaotic domain. Further, an extension of the results presented in Martí et al. (2016) is addressed, varying the main parameters of the system within the error tolerance.

Diffusion in the Arnold Model
The classical Arnold Hamiltonian (Arnold 1964) [already discussed in Ch79, Giorgilli (1990), Simó (2001), Lega et al. (2008), Efthymiopoulos (2012) among others] has the form 1 B(θ 2 , t) = sin θ 2 + cos t, with I 1 , I 2 ∈ R, θ 1 , θ 2 , t ∈ S 1 , where μ should be exponentially small with respect to ε, so that εμ ε 1. (The implications of the latter assumption will be discussed later on.) For ε = 0 the system allows two global integrals of motion, namely I 1 and I 2 , which determine the invariant tori supporting the quasiperiodic motion with frequencies ω 1 = I 1 , ω 2 = I 2 . Therefore, we have a very simple dynamical system consisting of two uncoupled free rotators, so that, θ 1 (t) = I 1 t + θ 0 1 , θ 2 (t) = I 2 t + θ 0 2 . For the case in which ε = 0, μ = 0 we still have two integrals, and the unperturbed Hamiltonian could be written as Notice that H ε 1 is the pendulum model for the resonance ω 1 = 0; H ε 1 ≡ h 1 = −2ε corresponds to the exact resonance or stable equilibrium point at (I 1 , θ 1 ) = (0, π) while h 1 = 0 to the separatrix, and thus, (I 1 , θ 1 ) = (0, 0) is the unstable point or whiskered torus. 2 The frequencies are now where ω p (h 1 , ε), the pendulum frequency is given by are the small oscillation and the half-rotation frequencies, respectively, being K (κ) the complete elliptical integral of the first kind. For rotations, the second line in (4) provides the half-rotation frequency in order to keep similar values of the frequency at both sides of the separatrix.
The half-width of the resonance ω 1 = 0 in action space is ( I 1 ) r = 2 √ ε, so within this resonance the variation of I 1 is bounded by | I 1 | ≤ 2 √ ε while I 2 remains constant.
However, due to the dependence of V ε on θ 2 , also the integral I 2 changes, and therefore, motion along the stochastic layer could proceed. Due to the chaotic character of the dynamics inside this narrow layer, the variation of I 2 should also be chaotic, giving rise to diffusion in I 2 . Thus, as I 2 might change unboundedly, a large instability could exist. This is the way in which Arnold diffusion is discussed in Ch79. In this particular model, since the perturbation V ε vanishes at (I 1 , θ 1 ) = (0, 0), for μ ε it is possible to build up a transition chain (Arnold 1964;Giorgilli 1990) such that if ω 2 is irrational, then all tori defined by I 1 = 0, I 2 = ω 2 are transition tori, so, when t → ∞, |I 2 (t) − I 2 (0)| = O(1), independently of ε and μ. Therefore, a "large variation" of I 2 could take place.
In (6) ω 1 = 0 is just one of the six first-order resonances involved. Indeed, it is easy to show that the primary resonances at order ε and εμ are Regarding the resonances in action (or energy) space, we should consider (7) with ω 1 = ω p (h 1 , ε) for h 1 < 0 and ω 1 = 2ω p (h 1 , ε) for h 1 > 0.
The resonances intersect at seven fixed different points in frequency space, namely (ω 1 , ω 2 ) = (0, 0), (0, ±1), (±1, ±1). 3 Hence, as pointed out in Ch79, the diffusion would spread over all this resonance set. Notice, however, that for εμ ε 1 the diffusion rate should be negligible along all resonances except for ω 1 = 0, since the latter has the main strength, its amplitude being ε, while all the remaining resonances have amplitudes εμ ε. Indeed, it is known that the theoretical diffusion rate depends exponentially on −1/ √ V mn , where V mn stands for the amplitude of the above-considered resonances [see for instance Ch79 and Cincotta (2002)].
Considering the fully perturbed motion, besides the ones given in (7), the full set of resonances (in action-energy space) is given by a linear combination of the form Figure 1(left) shows the resonance web in action space for ε = 0.25, for |m 1 | + |m 2 | + |m 3 | ≤ 6 and setting θ 1 = π so that I 2 1 = 2h 1 + 4ε. We clearly observe vertical resonances corresponding to m 2 = 0, horizontal ones to m 1 = 0 and an infinite but countable set of curves for m 1 , m 2 = 0 that accumulate toward the separatrix at I 1 = 2 √ ε.
For the sake of illustration, in Fig. 1 (right) we present the result of a numerical experiment for ε = 0.25 and μ = 0.025. This figure displays the actual resonances while plotting just the  (8) with ε = 0.25, taking θ 1 = π so I 2 1 = 2h 1 + 4 and |m 1 | + |m 2 | + |m 3 | ≤ 6 (left). Actual resonances in the Arnold model according to a MEGNO mapping for I 1 > 0, I 2 > 0, θ 1 = π, θ 2 = t = 0, ε = 0.25, μ = 0.025. Region in black corresponds to chaotic domains, while those in white correspond to periodic or quasiperiodic motion (right). For the given value of ε, the separatrix lies at I 1 = 1 MEGNO (Cincotta and Simó 2000;Cincotta et al. 2003;Cincotta and Giordano 2016) values for 10 6 initial conditions along a given region in the (I 1 , I 2 ) space, fixing θ 1 = π, θ 2 = t = 0 and after a total motion time t f = 10 4 . The actual resonance web is well revealed by the MEGNO contour plot; in particular, take notice of the complex chaotic structure around the separatrix of the resonance ω 1 = ω p = 0. The strongest resonances, besides the main one, can be distinguished in the figure with their corresponding chaotic layers surrounding the separatrices: ω 1 = 1, ω 2 = ω 1 and ω 2 = 0. Clearly the largest width of the chaotic layers, w s , corresponds to the resonance of amplitude ε (ω 1 = 0). A similar plot for ε = μ is given in Simó (2001), while the one for small values of the parameters and θ 1 = θ 2 = t = 0, obtained by recourse to the FLI [introduced first in Froeschlé and Lega (2000)], is presented in Lega et al. (2008).

Numerical experiments
Let us numerically investigate the diffusion in this model. To this end, we fix the value of ε and take two different values of μ. We will be concerned with diffusion on the chaotic layer of the main resonance, ω 1 = 0, and for different values of ω 2 = I 2 . According to Chirikov's theoretical results Ch79, when ε 1 and μ ε, the local diffusion rate for |ω 2 | < 1 should be larger than for |ω 2 | > 1. The estimates of the width of the chaotic layer of the main resonance and the diffusion coefficient (assuming a normal local character, see below) along the latter are given by (see Ch79 eqs. (7.17)-(7.20) and the concomitant discussion) where T a is the mean period of motion within the chaotic layer of the resonance ω 1 = 0: that in terms of the perturbation parameters and ω 2 becomes The above derivation of the diffusion coefficient rests on the assumption that for |ω 2 | < 1, the strongest perturbation to the motion close to the separatrix of the guiding resonance ω 1 = 0 corresponds to the closest resonance to the guiding one in frequency space, the layer resonances ω 2 = ±ω 1 (corresponding to the harmonics θ 1 ± θ 2 ) that give rise to the chaotic layer around the separatrix of the guiding resonance, and thus the diffusion along this chaotic layer is due to the smaller perturbing terms (the more distant resonances to the guiding one in frequency space), the driving resonances ω 1 = ±1. For |ω 2 | > 1, the layer resonances become ω 1 = ±1, while the driving resonances are ω 2 = ±ω 1 .
Around |ω 2 | ≈ 1, all these arguments and estimates no longer applies, as well as when ε 1. Indeed, only when ε 1, the Melnikov integral (see Ch79) could be approximated by an exponential of argument 1/ √ ε and thus, the relative amplitude In fact, Chirikov's estimates are given for diffusion along the thin chaotic layer of the guiding resonance, i.e., Arnold diffusion.
On adopting not too small values of the parameters, ε = 0.25 and μ = 0.025, 0.01 and |ω 2 | ≤ 2, the above estimates for the diffusion coefficient are not applicable and numerical experiments are definitively required. For the given values of the parameters, the center of the chaotic layer (the unperturbed separatrix) of the guiding resonance lies at I 1 = 1 and typical timescales for the motion inside the chaotic layer are T a 20. We focus on the diffusion within a range of |ω 2 | in which the guiding and secondary resonances exhibit overlapping and crossings as Fig. 1 (right) shows. We will only discuss the results corresponding to μ = 0.025 since those for μ = 0.010 are rather similar. From that figure, we observe that for |ω 2 | > 2 resonances interaction decreases significantly with |ω 2 | and thus any instability is expected to be almost negligible over any physical timescale. In fact, for |ω 2 | 1, the relative amplitude v is small and thus Chirikov's estimates become a plausible assumption.
For the numerical experiments, we selected action values on the chaotic layer, I 1 = 1, I 2 = ω 2 , and adopted θ 1 = π, θ 2 = t = 0 as initial values for the angle variables (the same initial values of the angles as those used to compute the MEGNO in Fig. 1). In rigor, we considered a set of six ensembles of 100 initial values of the actions around (1, ω 2 ) of size ∼ 10 −7 . The center of each ensemble lies at ω 2 = − 1.69, − 1.13, − 0.173, 0.053, 0.55, 1.55, respectively. We integrated the equations of motion with a time-step of size 0.01, for a total motion time 4 × 10 6 ≈ 2 × 10 5 T a , which is large enough in order to look for any significant diffusion at moderate times. Indeed, we are not focusing on instabilities over timescales that exceed any possible physical time as we shall see in the next section. The results are reported in Fig. 2. There, the wandering of the actions for the initial ensemble (depicted in white) are pursued and superimposed on the MEGNO contour plot (for |I 1 | ≤ 1.5, |I 2 | ≤ 2), each color Six initial ensembles (indicated as white circles) are followed onto the MEGNO contour plot for μ = 0.025, the concomitant trajectories that intersect the section |θ 1 (t) − π | + |θ 2 (t)| < 0.01 are depicted with different colors identifying a different set of 100 trajectories. In order to visualize the diffusion in the 2D plane, (I 1 , I 2 ), we considered the section defined by |θ 1 (t) − π| + |θ 2 (t)| < 0.01. From this set of experiments, we clearly observe that despite the initial location of the ensemble, the diffusion spreads almost uniformly over the range |I 2 | ≤ 2 and in all cases the same secondary resonances are visited after a motion time t ≈ 2 × 10 5 T a . Then, for these values of the adopted parameters, the instability is large with an expected nearly uniform mean diffusion rate of ∼ 10 −6 in the domain |I 2 | 2.
In order to envision the diffusion in the 3D space (I 1 , θ 1 , I 2 ), we regarded the section |θ 2 | < 10 −5 . The plot in Fig. 3 respects the same palette for individualizing the first four ensembles, where the normal form or pendulum model for the guiding resonance shows up. All the ensembles appear almost completely overlapped for |I 2 | 2, only for I 2 ≈ 2 the ensemble that starts with the larger value of ω 2 being distinguishable.

Numerical derivation of a diffusion coefficient
In dynamical systems, the diffusion coefficient is introduced through a generalization of the statistical properties of Brownian motion (i.e., Varvoglis 2005). In this direction, Chandrasekhar review (1943) is one of the most outstanding references in the field. In such systems, where the motion is ergodic, the diffusion coefficient D is just the constant rate at which the variance of any coordinate q evolves with time. This is usually called normal diffusion in the sense that being D independent of the position, direction and time. In other words, after certain initial transient T m , the system reaches the steady state and turns out completely mixed, diffusion becoming homogeneous and isotropic. Under this assumption, time and space averages are equivalent and thus the variance σ q could be defined in two ways. In the first case, a single orbit is followed along a large time span T and the average is performed over time intervals t T . For the second average, an ensemble of several trajectories in a small neighborhood around the same initial orbit is considered and the average is done over the whole ensemble at different times (for instance, t k = kδt with δt t). The space average naturally turns out to be much less sensitive to the selected orbit than the time average and of course, in any system exhibiting a divided phase space, both averages in general do not coincide [see for instance Meiss (1992) or Tsiganis (2008)].
Therefore, following the above-mentioned arguments, the numerical derivation of the diffusion coefficient should rest on the computation of the variance of a given ensemble of some appropriate variables. Just for the sake of simplicity, we discuss three different ways of this calculation for the results presented in Fig. 2.
Let N p be the number of initial orbits in a small neighborhood of (I 1 (0), I 2 (0)) and define, in case of a single resonance, (I r , I f ) as the corresponding resonant and fast actions and their conjugate angles (θ r , θ f ) (in the present example I r = I 1 , I f = I 2 ) 4 and take t j = t 0 + jδt, j ∈ Z, δt being for instance the integration time step (or a given multiple of 4 For a multiple resonance or an overlapped domain of resonances where there is not any specific direction it); thus, the ensemble variance is defined as In general σ 2 e (t) could be rather noisy and for small perturbations, its time evolution may hide any slow secular growth. However, it could be significantly improved if instead of using of the original variables (I r , I f ), new "optimal" actions are considered by a normal form construction to eliminate the deformation effect due to oscillations, as for instance it is discussed in Giorgilli (1990), Efthymiopoulos (2012) and Cincotta et al. (2014). In the present example, this procedure is not necessary since we have already shown that the Arnold model is in fact a normal form for the resonance ω 1 = 0, as Fig. 3 shows.
Another way to eliminate oscillations is to compute the variance over a double section [see for example Lega et al. (2003)], like that used to obtain Fig. 2. Formally, let S be defined as the section |θ r | + |θ f | < η 1, and take a fixed time interval t δt. Let t l = l t and consider motion times t l−1 < t ≤ t l , l ∈ Z. Assume that the N p initial trajectories have N l 1 intersections with S at different times τ l ∈ (t l−1 , t l ]; thus, the variance over the section S is defined as A variant is the cumulative variance over the section S. Let now t 0 < t ≤ t l , and assume that the N p initial trajectories have M l 1 intersections with S at different times τ l ∈ (t 0 , t l ], then It is simple to show that if N l ≈ N 0 1, nearly constant, then and for any power law σ 2 In Fig. 4 we present as an example, the evolution of these three variances for the ensemble given at the top right of Fig. 2. Since σ 2 s looks nearly linear with time and N l ≈ N 0 , we observe that 2σ 2 c behaves very similar to σ 2 s but in a smoother way. The ensemble variance σ 2 e also evolves, besides oscillations, as σ 2 s . Since in all cases we observe the very same behavior of the three variances, in what follows we shall present the results just for σ 2 e . Let us mention that when this linear trend is observed, the diffusion is called normal, while if σ e ∝ t b , for b < 1 we speak about subdiffusion and when b > 1 the process is called superdiffusion as Cordeiro (2006) and Cordeiro and Mendes de Souza (2005) show for the dynamics in the asteroid belt or Venegeroles (2008) for the standard map.
Taking into account the nearly linear character of the variances for the ensemble located at ω 2 = 1.55 and considering that for the six ensembles presented in Fig. 2  In the numerical investigations about diffusion we have reviewed, the estimation of a constant coefficient is performed by a linear fit on the evolution of any of these variances (e.g., Lega et al. 2008;Froeschlé et al. 2005;Cincotta et al. 2014). For instance, in Cincotta et al. (2014) the ensemble variance is computed after a normal form construction in order to compare the numerical estimates of the diffusion coefficient with the theoretical ones provided by Ch79 and also with Nekhoroshev estimates. In this work, the local diffusion coefficient is derived by a linear fit on the evolution of σ 2 e for different values of the perturbation parameter in a fixed ensemble of orbits on the very same location in the chaotic layer of a single resonance. On the other hand, in Cachucho et al. (2010) Chirikov's formalism is applied to an asteroidal three-body resonance, where the diffusion coefficient is computed locally for a grid of initial conditions along the chaotic layer of the guiding resonance by means of a time average on a single orbit. In Guzzo et al. (2005), Lega et al. (2008) and Guzzo et al. (2011) the local diffusion and global diffusion in the limit of weak chaos are studied, where the diffusion coefficient is estimated by a linear fit on σ 2 s , since the normal diffusion approximation seems to work in such a case. Figure 5 shows the evolution of σ 2 e for the six ensembles adopting the same color palette revealing very different evolution, in all cases σ 2 e /t depends on both, the initial position and time. Any attempt to derive a mean diffusion coefficient through a linear fit on σ 2 e would provide quite different values depending on the considered ensemble. It could be possible that for extremely large motion times they reveal a similar behavior; however, we are interested in the character of the diffusion on moderate or physical times (see next section).

the diffusion
All the numerical experiments considering different values of the parameters exhibit the same behavior. Therefore, we observe that when a large instability sets up, it could not be characterized by a standard diffusion coefficient obtained from the variance evolution of the actions. Indeed, the diffusion coefficient is a local property, as for instance Chirikov's estimates (10) show for the Arnold model. For a generic multidimensional Hamiltonian, Ch79 derives an expression for the diffusion coefficient which is clearly local and it is only valid in a small neighborhood of the initial condition on the guiding resonance [see also Cincotta et al. (2014)]. Thus, if the diffusion is fast and it spreads far away from the starting point P after a certain time span (as shown in Fig. 2), then any value of the diffusion coefficient at P would be determined by the dynamical structure of regions of the phase space very distant from P. Even though the layer of the guiding resonance within the domain |I 1 | 2 is highly chaotic and unstable, its motion is far from ergodic. Stickiness among other effects seriously prevents free diffusion, 5 and thus, it would not be right to measure the diffusion at P following the evolution of the ensemble over long time intervals. Let us then consider the evolution of the variance in much shorter times, for instance t 50T a 1000, which are presumably large enough for the ensemble to move all over across the chaotic layer but sufficiently short in order to locally characterize the diffusion at ω 2 . Let us define the statistical moments of the actions and energies as where h 1 and h 2 are the energies of each degree of freedom of the Hamiltonian H ε 0 (I 1 , I 2 , θ 1 ) given in (3); and we compute where I 2 j is a slight variation of σ 2 e and h j is the mean deviation of the energies of the unperturbed motion, which is quadratic in the actions. Then, we follow the evolution of both magnitudes in (16). Figure 6 shows I 2 j (t) for an ensemble of 2000 orbits located at ω 2 = 1.55. After an initial transient, both variances increase, I 2 1 exhibiting fast oscillations and reaching an asymptotic value three orders of magnitude larger than I 2 2 . About t ≈ 500 the ensemble fills completely the width of the chaotic layer and thus I 2 1 cannot increase any longer. Approximately at the same time, the oscillations in I 2 2 decay and this quantity starts to increase almost linearly with a slope D ≈ 4 × 10 −7 . Though such a value matches Chirikov's estimates (10), this does not mean that they are applicable in our experiments, as it has already been discussed. The evolution of h j is given in Fig. 7. Since we are computing energies at ω 2 1, h 1 and h 2 are of the same order of magnitude. For t 500, the perturbation μV ε to H ε 0 allows both degrees of freedom to vary (and exchange) its energy, but when the motion reaches the borders of the chaotic layer, h 1 is bounded by εw s and thus any increase of the total energy of the system could only be due to changes in h 2 or I 2 . The fact that h j is quadratic in I j explains the similarities of Figs. 6 and 7.
For ensembles at different locations, the results are in a broad sense similar: for small values of ω 2 , motion times of the order of 100T a are needed in order to observe a nearly linear increase of h 2 or I 2 2 . However, in some experiments, for instance for ω 2 = 0.55, we did not find any similar evolution as those presented in Figs. 6 and 7. This could be explained by the fact that at ω 2 ≈ 0.5 we observe the crossing of the guiding resonance with the coupling resonance ω 2 = ω 1 (see Fig. 2) and thus I 1 , I 2 are not the resonant and fast actions, respectively, the width of the chaotic layer around this value of ω 2 is not well defined, and thus, the arguments on the expected behavior given in the above paragraph do not apply in the case of a multiple resonance.
From the above discussion, it seems that the normal character of diffusion could be assumed locally, for a relatively short timescale, the latter being determined by the time required for the ensemble to reach the steady state in its motion across the layer. But for larger timescales, nothing could be said about the character of the diffusion. In this direction, we expect that in a generic near-integrable Hamiltonian system diffusion would be time dependent and inhomogeneous. Thus, a diffusion coefficient derived from the time evolution of any action variance over a large time span might not be a meaningful measure for large instabilities. In our experiments, a realistic measure of the diffusion over an extended region of phase space could be provided, for instance, by a mean or macroscopical coefficient defined as the size or diameter of the region swept by diffusion over the considered time span. However, in Giordano and Cincotta (2017) an alternative measure of diffusion is provided by means of an entropy-like estimator.
In conclusion, though it seems that a diffusion coefficient derived from the variance evolution of the actions could not be well defined in case of a macroscopic instability, we suggest that numerical experiments like those presented in Fig. 2 are actually very useful to investigate stability/instability domains within chaotic regions of phase space over physical timescales. Moreover, the evolution of the variance would provide information about the character of the diffusion as we shall show in the next section.

The Gliese-876 planetary system
In this section, we analyze how diffusion works in the case of an observed planetary system. Possibly one of the most interesting cases is Gliese-876, consisting of four known planets orbiting a low-mass star (M = 0.32M Marcy et al. 1998). The innermost planet, known as GJ-876 d, is very small, extremely close to the central star, and dynamically detached from the rest of the system. Consequently, we will restrict our dynamical analysis to the three other bodies (i.e., GJ-876 b, GJ-876 c, GJ-876 e-hereafter referred to as planets 1, 2 and 3). Two independent studies of radial velocity data, Rivera et al. (2010), Baluev (2011, have indicated that these planets lie in the vicinity of a Laplace-type multiplanet resonance such that the mean motion ratios of the consecutive bodies satisfy the relation n 1 /n 2 n 2 /n 3 2/1. Though the inherent uncertainties of the detection method imply that the planetary masses and orbital elements are not determined precisely, it is generally accepted that the dynamics of the system is dominated by the Laplace resonance, whose critical angle librates with a moderate amplitude (e.g., Baluev 2011;Nelson et al. 2015). Some resonant angles associated with the individual 2/1 two-planet commensurabilities exhibit librational behavior, indicating that the system is under the effects of several distinct resonant perturbation terms.
In Batygin et al. (2015), the authors showed that for all the best fits the dynamics is strongly chaotic with Lyapunov times of only a few decades although the orbits seem stable for timescales of the order of the system's lifespan. They also estimated diffusion timescales using a simple 2D analytical model assuming normal diffusion, finding an analytical estimation in excellent agreement with numerical simulations. Under the same assumptions, a crude estimation of a diffusion coefficient in the context of a non-adiabatic regime is given, yielding a slightly poor estimation compared to the coefficient derived via numerical simulations. Although the authors claim that this discrepancy is due to the coarseness of the model, a different interpretation is given in Martí et al. (2016) (from now on MCB16). In fact, therein a more detailed analysis of the resonant domain associated with the Laplace resonance of the GJ-876 system is provided. Namely, from dynamical maps drawn in the (a 3 , e 3 ) plane, we found the existence of two distinct regions, both within the librational domain: an outer one, characterized by a strong instability similar to that pointed out in Batygin et al. (2015), and also a chaotic but more stable inner resonance domain. The outer region seems to be associated with the overlap of different two-body mean motion resonant terms, while the inner one looks much more stable. Let us recall that even though MEGNO calculations indicated a chaotic dynamics within the whole resonance domain, none of our N-body simulations led to disruptions of the orbits within the integration timescales.
Also in MCB16 diffusion studies were performed for several initial conditions selected onto the (a 3 , e 3 ) plane. Nearly normal diffusion was detected only outside the libration domain, in highly chaotic, unstable regions of the phase space. Meanwhile, the many other considered initial conditions, particularly those in the inner resonant domain, led to diffusion rates (i.e., variances evolution) much smaller than those predicted by analytical models and then inconsistent with Batygin's conjectures (Batygin et al. 2015).
The studies presented in MCB16 considered variations in the values of the outer (less massive) planet's semimajor axis and eccentricity, while the remaining parameters of the system (masses included) were taken equal to the best fit values given in Rivera et al. (2010). Herein instead we will address the issue of the robustness of the diffusion rates estimates regarding the uncertainty in the remaining parameters of the planetary system. This is an important issue to be addressed when the times involved in the simulation are relatively short in comparison with most studies on diffusion. In other words, we inquire how representative of the global dynamics of Gliese-876 are the results delivered in MCB16.

Numerical experiments
To address this issue, we performed a new series of numerical experiments regarding the evolution of ensembles of initial conditions for different values of the parameters of the system.
In the first place, we considered a set of 8 ensembles, each consisting of 256 initial conditions taken within a range of 10 −3 in e 3 and 2 × 10 −4 AU in a 3 , which from now on will be referred to as our "nominal set." Such initial conditions were integrated for a total time span of 2 × 10 5 years, adopting for m i , a j , e j , i = 1, . . . , 3, j = 1, 2 the values of the best fit given in Table 1, which displays also the errors in their concomitant determination. 6 We performed three different sets of simulations for the 8 ensembles, varying in turn the parameters e 1 , e 2 and m 3 (which from now on will be termed e 1 , e 2 and m 3 sets), taking the values: e 1 = 0.25498, m 3 = 13.75M Earth ≈ 0.0433M Jup and e 2 = 0.0311, respectively, which are close to the extreme values within the error of the corresponding best fit. (We refer again to Table 1.) During the evolution of each initial condition in the 8 ensembles, we kept track of each time the system was close enough to its initial (or reference) plane, in order to account for the diffusion in the parameters a 3 and e 3 . As in MCB16, we considered that the system crossed the reference plane (section S) whenever the following three conditions were satisfied during the integration: ang , e and a predefined values, i denotes the longitude of the perihelium of the planets and 0 1 = 0 2 = 0, 0 3 = 180 • . In particular, for our experiments we adopted ang = 6 • , a = 0.005AU and e = 0.005. Let us recall that the values of e 0 i and a 0 i for i = 1, 2 correspond to the initial values of the adopted eccentricities and semimajor axes of the two innermost planets of the system. Since we are dealing with a six degrees of freedom Hamiltonian, the above-defined conditions confine the phase space motion to the nearly twodimensional section S: (a 3 , e 3 ). As a consequence, the number of intersections of a given These values were extracted from the best orbital coplanar fit from Rivera et al. (2010), with their respective uncertainties trajectory with S strongly depends on the stability of the motion. In case of stable regular motion, the section S is just a slice of the six-dimensional torus and thus many crossings would occur. Instead, for a unstable chaotic orbit, since the tori structure does not exist, only a few intersections with S are expected. A qualitative comparison between the diffusion of each ensemble of the different sets of simulations with the nominal ensemble of MCB16 reveals that the only significant differences are observed in the set of experiments corresponding to the modified value of e 2 , i.e., the e 2 set. The results concerning any of the other sets are analogous to those displayed for the nominal set in MCB16, in particular they confirm the existence of dynamically detached regions in the (a 3 , e 3 ) plane.

Ensembles for the e 2 set
It is worth noting that, as Table 1 shows, e 2 is the parameter with the greatest observational indetermination. The results of the integration of the ensembles in the e 2 set are shown in Fig. 8. Let us state clear that the contour plot serving as background corresponds actually to a e 3 dynamical map performed with the nominal set of orbital parameters. The general overview of the evolution of the 8 ensembles discloses a significant decay in the relative number of crossings to the number of initial conditions. As already mentioned, though this is not a rigorous dynamical indicator, it should be expected that less unstable trajectories would have much more intersections with the reference plane.
The differences between the set of ensembles corresponding to a modified value of e 2 and any of the others are more than evident. In the first place, there is a significant drop in the number of crossings with the reference plane for this set of simulations, which barely amounts to one order of magnitude lower than that obtained for any of the other three sets. The main differences are noticed in panels corresponding to ensembles 5, 7 and 9, which starting inside the inner domain eventually reach the outer region.
Though none of the initial conditions of these highly chaotic ensembles becomes unstable at least up to 2 × 10 5 years, since the orbits can reach the outer resonant domain from within the inner region, such domains are no longer dynamically detached and these trajectories are most likely to lead to disruption for a sufficiently large time span.
In order to gain some insight on this issue, we performed the numerical integration of ensembles 7 and 9 from the set corresponding to the modified e 2 value, but now taking 900 initial conditions for each ensemble so that the number of crossings with the reference plane is increased.
The results concerning ensemble 9 are presented in Fig. 9, where we have zoomed out the map range in order to clearly distinguish every crossing with the reference plane. Therein a few intersections surrounding the initial condition and inside the stable region are observed.  Moreover a non-negligible number of crossings were obtained in the outer region where the system is supposed to become unstable. It is remarkable that all these crossings are widely spread along the semimajor axis of the third planet, though a much more restricted diffusion is observed in e 3 . Similar results were obtained for ensemble 7, for which diffusion in the a 3 parameter is also much more effective than the one in e 3 .
Finally in Fig. 10 we display the values of a 3 and e 3 corresponding to the crossing points, in terms of the time at which each crossing occurs, both, for the nominal ensemble 9 (in blue), and for the ensemble of 900 initial conditions corresponding to the e 2 set (in red). Despite the far lower amount of crossings for the latter ensemble, we can see that diffusion is much more noticeable in the a 3 parameter. For e 3 instead, diffusion seems to be as limited as for the nominal ensemble.

Diffusion rates
In order to characterize the diffusion in the parameters a 3 and e 3 , we followed their variance evolution for each particular ensemble. (This procedure is the same as that carried out in MCB16.) As already discussed herein, any estimation of a diffusion coefficient by just adjusting a simple Brownian model does not work in general and, further, it fails inside the resonant region, as already shown in MCB16. Nonetheless, it could well serve to compare the different situations described in the previous subsection. Since the resonant domain is dominated by a chaotic dynamics due to the overlap of several resonances (namely twobody and three-body mean motion resonances), the ensemble variances could hide any slow secular change of the orbital parameters, in particular when diffusion is confined to a small region of the resonance. Thus, we decided to use the cumulative variance over the section to characterize the diffusion evinced in Fig. 8 (see, however, the end of this subsection).
Therefore, in each panel of Fig. 11 we display the cumulative variances over the section S attained for a 3 corresponding to the 8 ensembles for each different set: the nominal, the e 1 , e 2 and m 3 sets. Notice that in every plot we have also depicted the curve corresponding to normal diffusion for the sake of comparison and further analysis. Indeed, curves with slopes close to that corresponding to the normal case would imply that the orbits are dominated by an almost fully chaotic, diffusive dynamics, while variances evolving with slopes smaller than the normal one are representative of a chaotic but much more stable dynamics. Aside from a considerable lack of crossing points, which makes the sample statistically poor, a substantial different variances' evolution can be appreciated in the case of the e 2 set from those corresponding to the remaining sets. Indeed, almost every ensemble in the e 2 set evolves with a slope much closer to the normal one than those ensembles in any of the other scenarios (except maybe for some particular ensemble in the other sets, such as ensemble 4). For nearly all the ensembles in the nominal, the e 1 and m 3 sets, the curves describing the variances are significantly different from those corresponding to free diffusion. Figure 12 exhibits similar results as those presented in Fig. 11, but for e 3 . There we also observe, as in the case of a 3 , that the variances' evolution for the e 2 set is much closer to the normal behavior than those corresponding to ensembles in any of the other sets.
The results allow the conclusion that the rate in which the variances evolve is, in general, much smaller than what would be expected for a Brownian motion. Indeed, this is the case not only for the nominal set, which has been discussed extensively in MCB16, but also for the e 1 and m 3 sets, which should present qualitatively the same diffusion rates. Meanwhile, even though the e 2 set exhibits variances' evolution sensibly different from those corresponding to the other sets of ensembles, their diffusion rates depart appreciably from those associated with the normal case.
Aside from ensemble 4, let us notice that the variances' evolution for every ensemble in a particular set is quite similar. Such a behavior clearly shows that the associated diffusion rates are, not only lower than the expected one for normal diffusion, but they also appear to be independent on the initial conditions. As long as the initial location of the ensembles is taken inside the resonant region, the expected diffusion rates of the different ensembles become In each of these plots, the curve corresponding to normal diffusion is drawn comparable. The case of ensemble 4, taken rather close to the border of the inner Laplace resonant region, exhibits a somewhat higher diffusion (as already discussed in MCB16). Judging from the results of the Figs. 11 and 12, it appears that the e 2 set of ensembles present faster diffusion. Although this might not be conclusive, mainly due to the lack of a sufficiently large sample of crossing points, the resonant structure, which seems to act as a protection mechanism for the system, appears to be more rapidly altered by the variation of e 2 than by that of the other parameters, namely e 1 or m 3 .
In order to reinforce this conclusion, we computed for the ensemble 9, which appears to be very sensitive to slight variations in e 2 , the variance over the ensemble instead of the cumulative variance. Figure 13 shows that the evolution of the ensemble variances corresponding to a 3 and e 3 are though subdiffusive, closer to the normal diffusion case, revealing that indeed this particular ensemble becomes highly unstable under a small change in the eccentricity of the second planet.
We can thus ensure that a small variation of the initial value of e 2 from the orbital best fit can lead to a significant change in the resonant structure of GJ-876, and may even incur on a completely unstable domain, which could eventually lead to a disruption of the system.

Conclusions
In the present work, we have reviewed the diffusion process in near-integrable Hamiltonian systems with more than two degrees of freedom. We showed that in a simple model like the Arnold Hamiltonian (with parameters not too small), diffusion is inhomogeneous and anisotropic since the system is quite far from the steady state. Therefore, the Brownian a 3 Variance normal diffusion motion approximation is not suitable to modelate the diffusion in our set of experiments. Further, the estimation of a local diffusion coefficient from the variance evolution of the ensemble fails, because it significantly departs from a linear trend. Our results show that large instabilities are not well characterized by a standard diffusion coefficient, since the latter only provides local information about the dynamics around the initial point. In fact, a theory dealing with the correlations produced by stickiness that lead to a time dependent, non-homogeneous and anisotropic diffusion coefficient is still lacking.
Anyway, we showed that diffusion experiments could be very useful to understand the dynamics on chaotic domains. The application of this kind of studies to the well-known planetary system GJ-876 allows us to argue that different dynamical scenarios are possible for the Laplace resonance, in particular when considering allowed values of eccentricities of the planet c. The system could exhibit a very stable inner resonant domain detached from a highly unstable outer one when considering the best fit value of this eccentricity, but after a slight variation of this parameter the inner region would become also unstable and connected with the outer one. In sum, it seems that new (precise) numerical simulations in the parameter space are definitively needed to get a clear understanding of the dynamics of this multiresonant planetary system.