Shannon entropy diffusion estimates: sensitivity on the parameters of the method

In the present effort, we revisit the Shannon entropy approach for the study of both the extent and the rate of diffusion in a given dynamical system. In particular, we provide a theoretical and numerical study of the dependence of the formulation on the parameters of the method. We succeed in deriving not only a diffusion coefficient, DS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{S}$$\end{document}, but also an estimate of the macroscopical instability time for the system under study. Dealing with a toy model, namely a 4D symplectic application that represents the dynamics around a junction of resonances of different order, and an a particular case of the planar three-body problem, the HD20003 planetary system, we obtain numerical evidence that DS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{S}$$\end{document} is a robust measure of the diffusion rate, no significant dependence on the free parameter of the entropy formulation (the size of the elements of the partition) being observed. Moreover, successful results concerning the estimation of macroscopical instability times obtained from DS\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_{S}$$\end{document} are presented in both cases.


Introduction
The Shannon entropy approach to characterize the chaotic diffusion in multidimensional systems was already introduced in Giordano and Cincotta (2018), Cincotta and Shevchenko (2019). Therein, both the entropy S and the diffusion rate D S (or S in the very first formulation) defined through the time-derivative of S were shown to provide a proper measure of unstable dynamics despite the nature of the transport process induced by the dynamics. In a recent work Beaugé and Cincotta (2019), successful estimates of the escape or disruption times in the planar restricted three-body problem were derived by means of the inverse of D S . In the latter approach, some variations of the original formulation were considered in order to allow computationally feasible numerical experiments.
Chaotic diffusion is a typical global instability observed in different fields of astronomy, physics and chemistry, among others. In dynamical astronomy, diffusion plays a relevant role in planetary and galactic dynamics, as reviewed and discussed in Brunetti et al. (2011), Contopoulos and Harasoula (2013), Martí et al. (2016), , Maffione et al. (2015Maffione et al. ( , 2018. As shown in several investigations concerning weak chaos, the diffusion proceeding along the homoclinic tangle of a single resonance (or resonance web) over very large motion times could be approximated by a nearly uncorrelated process (see, for example, Lega et al. 2003;Froeschlé et al. 2005Froeschlé et al. , 2006Guzzo et al. 2005;Lega and Froeschlé 2008;Efthymiopoulos and Harsoula 2013;Cincotta et al. 2014), leading then to a slow drift of the actions in such a way that the action variance scales nearly linear with time.
In the chaotic regime and for relatively short or physical motion times, diffusion could deviate from the normal regime as discussed in . Indeed, in this scenario the variance evolution is highly sensitive to stickiness and this phenomena seriously affects its time rate, as we will show. As a consequence, any estimation of a diffusion coefficient computed over a finite time interval by means of a normal diffusion assumption would, in general, fail.
In the present work, we extend and reinforce the Shannon entropy formulation, encouraged by the results presented in Giordano and Cincotta (2018), Beaugé and Cincotta (2019). In particular, we will focus on the dependence of both S and D S on the parameters involved in the dynamical system as well as on the size of the partition, which in fact is a free parameter in the computation of the entropy. We provide theoretical and numerical evidence that D S is almost independent of the partition size provided that it is chosen in some optimal way. Moreover, a macroscopical instability time τ inst can be derived by means of the inverse of D S . The estimated τ inst can later be compared with the times obtained by direct numerical simulations.
To this end, we first consider the simple 4D symplectic map introduced in Gelfreich et al. (2013) to represent the dynamics in the neighborhood of a junction of resonances of different order. This model becomes an interesting and adequate scenario to perform many series of numerical experiments at low computational cost.
With the help of the obtained results for the symplectic application, a similar approach is given in the context of the planar three-body problem, namely the HD 20003 exoplanetary system (Mayor et al. 2011;Gillon et al. 2017), where two close-in Neptune-mass planets are near the 3/1 commensurability of mean motions.

The Shannon entropy
Herein, we briefly summarize the formulation already given in Giordano and Cincotta (2018), , Cincotta and Shevchenko (2019) and recently applied to the restricted planar three-body problem in Beaugé and Cincotta (2019), where the Shannon entropy is proposed as an alternative way to estimate a diffusion coefficient. Further theoretical background on the Shannon entropy can be found, for instance, in Shannon and Weaver (1949), Katz (1967), Wehrl (1978), Lesne (2014), Arnold and Avez (1989).
Let us consider a discrete dynamical system defined through actions I 1 , I 2 ∈ R and phases ϑ 1 , ϑ 2 ∈ S 1 , such that both actions can also be defined in the torus. Indeed, if I 1 , I 2 are periodic with periods I m 1 , I m 2 , the action plane can be reduced to the unit square with opposite sides identified T . Let S γ = {(I 1 , I 2 ) i , i = 1, . . . , N s } ⊂ T be the section defined by the intersections of any trajectory γ of the full 4D map with |ϑ 1 | + |ϑ 2 | < δ 1. Let α be a partition of T , i.e., a collection of r bi-dimensional cells that cover T . The elements of the partition are assumed to be measurable and disjoint. If n k is the number of intersecting points in S γ falling in the cell a k and if r 0 ≤ r denotes the number of non-empty cells, then the entropy can be recast as Notice that for a given partition and any γ , the entropy is always bounded. In particular, for (1) it is 0 ≤ S(γ , α) ≤ ln r 0 . The minimum of S is reached when γ is confined to a single element of α, which corresponds to the case of stability and, in such a case, it is S = 0. On the other hand S reaches its maximum, S = ln r 0 , when n k = N s /r 0 for all k, that is, when γ fills uniformly T , i.e., in case of ergodicity.
Following Cincotta and Shevchenko (2019) if γ r is a random trajectory, r 0 = r and (2) reduces to which follows after assuming a Poissonian distribution for n k , with mean value N s /r 1. Thus, the last relation is a good approximation to the entropy in case of random motion.
On the other hand, for a chaotic orbit γ where in general the n k distribution is not Poissonian and r 0 ≤ r , if the restriction N s /r 0 1 holds, then the entropy reduces to where β measures the deviation of the distribution density from a Poissonian one (see Cincotta and Shevchenko 2019;Cincotta et al. 2021 for details). Again, if βr 0 /N s is small, it is S ≈ S 0 = ln r 0 . Accordingly to Giordano and Cincotta (2018), Cincotta and Shevchenko (2019), Beaugé and Cincotta (2019), recalling thatṠ and assuming that δr where δ I 2 j (t) denotes either δ I 2 1 , δ I 2 2 or (δ I 2 1 + δ I 2 2 )/2 and · is the average over the interval (t, t + δt), on introducing S ≡ r 0 (t)Ṡ, then a local diffusion coefficient can be estimated as where σ is a normalization constant with dimensions of [I 1 I 2 ] that depends on the full domain where the diffusion occurs (not restricted to T ). If |I j (t)| ≤ I m j , j = 1, 2 for all times, then this domain is (−I m 1 , I m 1 ) × (−I m 2 , I m 2 ), while if |I j (t)| > I m j then the diffusion spreads over a region defined by I jmax , I jmin , j = 1, 2, being the latter the maximum and minimum values of the actions over the full motion time. Let ( I ) M be the diameter of the region where the diffusion spreads, then provides the area in action space of the unit cell on the full domain where the motion takes place. Therefore, a global diffusion coefficient can be estimated as Further details about this approach can be found in Giordano and Cincotta (2018), where 2 denotes a mean square displacement (for instance the squared distance between the initial and given boundary values of the actions), K being a numerical constant of the order of unity that takes into account the diffusion path in action plane and any eventual mild dependence of D S on r as we will show.

The choice of the partition
The entropy formulation includes a free parameter, the size of the partition r . Clearly, its selection depends on the total number of intersecting points with the section S, N s , and on its distribution density, which is determined by the dynamics. At first sight, the numerical computation of the entropy, for a finite number of iterates N , requires that N s /r 1. However, this condition should apply in the case of random dynamics, so that the number of intersections per cell, n k , follows a Poissonian distribution with mean value (and variance) λ = N s /r , as already discussed in , Cincotta and Shevchenko (2019). This means that r 0 (t) = r for t = N large enough or r small enough.
This would not be the case in the scenarios we are interested in, namely, systems with divided phase space where the instability would be confined to some domain of action plane. In such a case, if r 0 r and the density of intersecting points over the r 0 cells is nearly uniform, with n k ≈ λ 0 = N s /r 0 , the distribution of the n k approaches δ(n k − λ 0 ). Indeed, as discussed in , when the motion is dense with respect to some domain in T , |S − ln r 0 | ∼ O(N −2 ) instead of O(N −1 ) that applies for a Poissonian distribution.
For the computation of the entropy, the restriction λ 1 allows small values of r . However, a necessary condition for the approximation (6) is that r 0 (t) grows with time, so r 0 (t) r in such a way that δr 0 (t) could take into account increasing values of the actions. A relative small value of r would lead to r 0 (t * ) = r after a motion time t * ≤ N , the entropy values become S(t ≥ t * ) = ln r and thenṠ(t ≥ t * ) = 0. Therefore, for the computation of both S and D S , r should be large enough in such a way thatṠ is positive for unstable motion.
Assume full ergodicity enclosed in some domain T d ⊆ T . If d ≡ diam(T d ) diam(T ) = 1, then for r large enough, most of the cells would be empty, the N s intersections are confined just to r 0 (t) r cells that should cover T d . Thus, we could replace the restriction λ 1 by λ 0 1, with r 0 r but allowing λ ∼ 1, implying S ≈ ln r 0 < ln r . In other words, the condition λ 0 1 states that, in mean, after a given motion time each occupied cell reaches an invariant measure λ 0 .
A suitable selection of r could be, in case of a nearly dense distribution in T d , r N s , so r 0 (t) N s r . Therefore, the final value of the normalized entropy should satisfŷ and we can empirically adopt a threshold bound,Ŝ <Ŝ c . On the other hand, the minimum size of the cells, i.e., an upper bound for r , depends on the smallest variation of the actions that would be regarded as diffusion. Therefore, if δ I denotes the minimum change in the actions that we are interested to measure, then r ≤ δ I −2 .
Notice that D S should be invariant under a change of partition, r →r . Indeed, if the diffusion is localized on T d such that the distribution of the n k approaches a δ-function, then the density ρ(I 1 , I 2 ) ≈ ρ 0 for (I 1 , I 2 ) ∈ T d then from (8) where ( I ) 2 M denotes the maximum variation of the actions on phase space and thusD S = D S .
As defined, the entropy depends on the partition, so under this changeS becomes that reduces toS Then, if pq > 1,S > S, whileS < S if pq < 1, but in general it will beS ≈ S for a wide range of values of pq and large r .
Recall that all these analytical estimates on the optimal choice of the partition rest on the assumption that the density of iterates approaches a δ-function; otherwise, numerical simulations are required.

A simple model
where ε, η are real parameters and a 2 , a 3 ∈ Q. As discussed in Gelfreich et al. (2013), this map represents the motion around the junction of the resonances I 1 + a 2 I 2 = 2πk/η, a 3 I 2 + a 2 I 1 = 2πk /η, k, k ∈ Z whenever η 2, ε 1, and it is 4D generalization of the Chirikov standard map family (see, for instance, Chirikov 1969Chirikov , 1979Kaneko and Bagley 1985).
The map is invariant under the transformation with q, p the minimum integer numbers such that q/a 2 ∈ Z and pa 2 = ra 3 , r being also an integer number. Thus, we restrict the action space and therefore the opposite sides of D are identified. In what follows, we simply adopt a 2 = 1/2, a 3 = 5/4. A more detailed study of this map is given in Cincotta et al. (2021).
The global structure of the map and its resonance web for given values of the parameters is presented in Fig. 1 where a dynamical map by means of the MEGNO (Cincotta and Simó 2000;Cincotta et al. 2003;Cincotta and Giordano 2016) is given after adopting the section ϑ 1 = ϑ 2 = 0. Light blue indicates stable motion, while dark blue chaotic one. The resonance pattern is clearly observed from the plot, in particular the primary resonances I 1 + a 2 I 2 = 2πk/η, a 3 I 2 + a 2 I 1 = 2πk /η. The Mean Exponential Growth factor of Nearby Orbits (MEGNO) provides a fast estimate of the maximum Lyapunov exponent; it approaches a constant value, 2, for stable motion, while it increases linearly with time for chaotic motion, being the maximum Lyapunov exponent its time rate.
We observe that for the ensembles located very close to the resonance crossings (those depicted in violet and yellow), the diffusion is bounded around the junction and its extent is comparable in both directions, while, for the remaining ensembles, it spreads over the resonances a 2 I 1 + a 3 I 2 = 2πkη −1 and I 1 + a 2 I 2 = 2πk η −1 .
It is customary to characterize and measure diffusion through the evolution of the variance of a fast action I f ; thus, here we compute the variance over the section (i.e., the action plane) (see  and references therein for alternative ways of computing the variance). To this aim, we adopt an integer time interval N 0 N , define k l = l N 0 , l ∈ N, and consider discrete motion times k such that k l−1 < k ≤ k l . If the n p trajectories of the ensemble present m l 1 intersections with S at different times κ l ∈ (k l−1 , k l ], then the variance over the section is computed as where is the mean value of the ensemble. For an initial condition lying on a single resonance, I f ( j) denotes the values of the actions of the ensemble taken along the direction of the resonance line while, on an overlapped domain, I 2 f = I 2 1 + I 2 2 . For the computation of the variance (15), I 1 , I 2 are neither normalized nor restricted to any given domain in action plane as, for instance, we did in Fig. 1.
In Fig. 2 (left panel), we plot the evolution of y 2 taking N 0 = 10 4 and considering For the adopted values of the parameters, the variances of the actions for the four ensembles evolve, in mean, in different ways. Besides oscillations, they increase over the full time span as t b ; in general with b = 1, this is particularly evident for the ensemble located far away from the origin (plotted in blue), for which b < 1, while for the ensemble taken close to the origin b > 1 (plotted in violet) at least for the time-span considered. Figure 2 (right panel) reveals this behavior where y 2 /k is computed, and it provides a crude estimate of the diffusion rate in case b ≈ 1. Let us mention that in every case, for all the ensembles, it is m l ≈ 200. Since the diffusion is anisotropic and inhomogeneous, the evolution of the variance for each ensemble is expected to behave in different ways.
A fit of the evolution of the variance by a power law like y 2 = Dt b leads to values of b below and above 1; however, the discussion concerning the nature of the diffusion, i.e., whether it is normal (b ∼ 1) or not, and the estimation of the numerical values of b and D is postponed to a forthcoming work Cincotta et al. (2021). Indeed, the question regarding the nature of diffusion in 4D maps for finite motion times is in fact an open problem. The above examples represent a typical behavior of the variance in the general case (see  and references therein), the diffusion, for relatively small values of the perturbation parameters and comparatively short motion times, looks like anomalous and thus it seems to be cumbersome to get a diffusion coefficient by a linear fit, i.e., assuming a normal transport process. As mentioned, all these topics will be discussed in detail in a future effort.
In this direction, in Giordano and Cincotta (2018) the Shannon entropy was introduced to provide a measure of both the extent of the diffusion and its time rate, despite of its nature, that is, whether normal or anomalous.

Numerical computation of the Shannon entropy
Let us consider the same set of parameters and initial conditions that led to the action excursions shown in Fig. 1. In those simulations, the approximation of a nearly constant density of iterates over the section could be a assumed. We introduce then in T two different partitions, r = 1000 × 1000 and r = 3000 × 3000 elements, and compute both S and D S for all the experiments, considering again N 0 ≡ δt = 10 4 , n p = 10 3 and N = 10 7 iterates of the map (14), the total number of iterates being n p × N and N s = m l N /N 0 the intersections with S. In any case, experimentally we get m l ≈ 200, so N s ≈ 2 × 10 5 , and accordingly to (12) witĥ S c = 0.75, we should take 2 × 10 5 < r 11.7 × 10 6 (500 × 500 < r 3500 × 3500). Notice then that both adopted partitions fulfill the criterion (12). From now on, S ≡ S/ ln r such that 0 ≤ S ≤ 1. Figure 3 presents the evolution of the entropy with the integer time k, N 0 ≤ k ≤ N , for the four ensembles taken above. Note the mild dependence of the entropy on the partition size, for the entropy exhibits a very small shift upward when considering smaller cells, but in any case the derivativeṠ seems to be non-sensitive to the partition size as it was already discussed.
The entropy for the ensemble presenting the more extended diffusion is always higher than the rest, while for the ensemble that evidence a bounded diffusion around the origin, the entropy takes smaller values. In any case, S displays a close to logarithmic evolution with k, revealing a quasi-linear increase of r 0 with time so that its time derivative measures the rate of occupancy of the partition's cells.   Fig. 1 for the two adopted partitions. We observe that, in particular, for r = 3000 × 3000, D S evolves in a nearly constant fashion and it takes final values similar to those of y 2 /k (compare with Fig. 2). In any case, a diffusion rate can be estimated despite the size of the partition, namely ∼ 10 −7 . Let us mention that, as already discussed in Giordano and Cincotta (2018), equivalent overlapping diffusion extents (corresponding to different initial ensembles) would lead to a quite different section variance evolution, i.e., to a different dynamical characterization of the diffusion. Instead, both S and D S present similar temporal evolutions as well as final values after N iterates as Figs. 3 and 4 show, which is consistent with the comparable diffusion areas observed in Fig. 1 for the four initial ensembles.
Let us consider the ensemble in Fig. 1 whose intersections with S are depicted in blue. Figure 5 presents the final values of S and D S at N = 10 7 and the same values of the remaining parameters for different partitions r = m × m, where m ranges from m = 1000 to 4900. In the plot for S, we also include the approximation given in (13) with S ≈ 0.69 for r = 1000 × 1000 and 1000 pq = m × m. The results show that both S and D S are nearly invariant with respect to the size of the partition whenever r > N s . Indeed, within the full range of m, the maximum variation of the entropy is S ≈ 0.014, it increases from m = 1000 up to m ≈ 2000 but for m > 2000 it slightly decreases, reaching similar values at both extreme values of m. This reduction of the entropy for large m values is due to the fact that for resolutions higher than ∼ 2000 −2 the density of the N s intersections with S starts to reveal its non-uniform character, i.e., when S deviates from the expected theoretical value as m increases. Indeed, the approximation ρ(I 1 , I 2 ) ≈ ρ 0 requires that the surface density of points n ≈ N s /d 2 1 (d being the diameter of the diffusion region), and thus, the average distance between intersecting points n −1/2 is small with respect to 1/m.
On the other hand, D S is bounded in the interval 1.2 × 10 −7 < D S < 3.8 × 10 −7 and it decreases with increasing r in a mild way. Let us notice that the computed diffusion coefficient differs at most in a factor 3 when r does so in a factor 25. However, since the optimal partitions are those with r 3500, we get the estimate D S ≈ 2.5 ± 1.5 × 10 −7 . Therefore, the order of D S is well determined despite the size of the partition and the estimate given by (12) provides a natural bound to the choice of r whenever the distribution density is nearly constant in T d .
Let us proceed with further experiments now adopting slightly higher values of the parameters, ε = 0.65, η = 0.75, such that a higher diffusion rate would be expected. We pick up 150 initial conditions on the segment I 2 (0) = 0.07I m 2 , 0 ≤ I 1 (0) ≤ 0.35I m 1 on the action plane, which is depicted as a red line onto the MEGNO contour plot for a small region of T in the left panel of Fig. 6. We take small ensembles of 100 nearby initial conditions centered at any of the above defined values of I 1 (0), I 2 (0) and iterate M up to N = 5 × 10 8 to finally determine a mean escape time, t esc defined as the instant when either |I 1 (t)| ≥ I m 1 or |I 2 (t)| ≥ I m 2 for any of the initial conditions in the ensemble. In other words, we define t esc as the time when the trajectories leave the region D = (−I m 1 , I m 1 ) × (−I m 2 , I m 2 ). Similarly, we compute D S for larger ensembles of n p = 1000 initial conditions in the neighborhood of the same initial values on the red segment, taking a partition r = 2000×2000 and iterating the map up to N = 5 × 10 6 . Using the definition given in (11), we derive an instability time τ inst introducing a mean square displacement 2 = (I m 1 − I 1 (0)) 2 + (I m 2 − I 2 (0)) 2 . As we have already seen from Fig. 1, the diffusion mostly spreads over the homoclinic tangles of the primary resonances and thus the above expression for 2 would overestimate the true mean displacement on T , so we set K = 1/4 (the proper value of K could be adjusted once for all by making a couple of values of t esc and τ inst to coincide). Figure 6 right panel shows the results for both t esc and τ inst revealing the efficiency of this approach to estimate macroscopic instability times. Indeed, the plot shows a very good agreement between t esc and the instability time derived from the diffusion coefficient for values of the former below the total time considered N = 5 × 10 8 . For initial conditions leading to a very small D S no escape takes place, the motion remaining encompassed in D at least for N iterations. In such cases, we get values of τ inst > N . Since the determination of t esc is done up to It is worth mentioning that this simple model, a 4D application, allows us to perform numerical simulations involving for instance quite large values of r and n p that could not be the scenario in more complex systems as the one we discuss in the next section.

Application to the planar 3BP
Herein, we present an application of the Shannon entropy technique in the framework of the planar (non-restricted) 3BP. In particular, we have chosen a specific exoplanetary system, the HD 20003, and solved the trajectories of the bodies through long-term integrations of the exact equations of motion for given sets of initial conditions (ICs). At first, the resulting solutions are analyzed in terms of the Shannon entropy, the dependence of the method on its parameters being discussed. Further, we have focused on observing the concordance between the outcomes from the instability times derived from the Shannon entropy in comparison with the system's lifetime obtained from direct N-body simulations.
The HD 20003 exoplanetary system was reported in Gillon et al. (2017), Mayor et al. (2011) and more recently in Udry et al. (2019). It consists of two close-in Neptune-mass planets whose mean motion ratio n 1 /n 2 is found near the 3/1 commensurability. The central star is a late-G dwarf of mass m 0 = 0.875M . In the literature, the planets are named HD 20003 b and HD 20003 c, but throughout this section we designate them with the subscripts 1 and 2, respectively. In Udry et al. (2019), the values of the nominal solution of the system are given, as well as the associated errors. These nominal values are gathered in Table 1, where we include the respective values for the semi-major axes (a), eccentricities (e), the orbital periods (P), the arguments of pericenter (ω) and the time of passages at periastron (T P ) (actually, with these latter quantities we have computed the initial values of the mean anomalies (M)).
The system can be initially described by a 4-degree of freedom (4-dof) Hamiltonian  the planar case), the system is reduced to a 3-dof Hamiltonian Ferraz-Mello (2005). We considered then a model where the Hamiltonian could be split in two parts, H = H 0 + H 1 , in the framework of the Jacobi canonical coordinate system of positions and momenta (ρ 1 , ρ 2 , p 1 , p 2 ), with the Keplerian (unperturbed) part H 0 written as in Andrade-Ines and Robutel (2018), where the mass factors are μ 1 = G(m 0 + m 1 ) and μ 2 = Gm 2 (m 0 + m 1 ), and the reduced masses m 1 = (m 0 m 1 )/(m 0 + m 1 ) and m 2 = m 2 (m 1 + m 2 )/(m 0 + m 1 + m 2 ). Recall that 2 Ferraz-Mello (2007). In all these definitions, a i indicates the i-th body's Jacobian osculating semimajor axis. The perturbation H 1 can be recast as (see Andrade-Ines and Robutel 2018) where κ 0 = −m 1 /(m 0 + m 1 ) and κ 1 = m 0 /(m 0 + m 1 ). The perturbation H 1 allows trigonometric series expansion which can be obtained by classical analytical methods. Though it is not a concern of the present work, further details on such a development for the disturbing function as well as the derivation of the relations (17) and (18) can be found in Andrade-Ines and Robutel (2018).
In the numerical simulations for the HD 20003 system, we applied an N-body routine (see, e.g., Beaugé and Cincotta 2019) for solving the classical Newtonian 3-body equations of motion. On adopting the Jacobi reference frame, we numerically integrated the equations of motion in rectangular coordinates of the positions and momenta (ρ 1 , ρ 2 , p 1 , p 2 ) and, by applying the classical transformations, we finally recovered the values of the osculating orbital elements Ferraz-Mello (2007).
The left panel in Fig. 7 shows a MEGNO contour plot map (Y ) in the (a 1 , e 1 ) plane for each trajectory for a grid of 400×400 ICs integrated over a 10 4 yrs time-span, the outer planet being fixed at its nominal position given in Table 1. The plot displays the particular structure of phase space associated with the 3/1 mean motion resonance (MMR), which is highlighted by a red line. We observe a central region of regular trajectories (Y ≈ 2) that widens as the inner planet eccentricity grows. The nominal position of the system (hereafter called (C0)) is highlighted by a black dot, and it lies considerably away from the assumed resonant region and yet stable motion. Also, we have indicated by a white dot another IC, (C1), for which the initial value of the inner planet eccentricity is e 1 = 0.4589. Such an IC is seen to be located inside a completely chaotic region of the phase space. The plot on the right shows a zoom of the full MEGNO map around both initial conditions. Let us note an indentation of the phase space through some "discontinuous" structures that we have associated with high-order resonances. (Two of them are identified by vertical lines in the surroundings of (C0) and (C1).) By means of the N-body routine, we performed long-term numerical simulations using both ICs, (C0) and (C1). Figure 8 displays the time evolution of the semi-major axis (upper plot) and the eccentricity (middle plot) of the outer planet for both ICs. The distinct character of these two single solutions is quite evident. Indeed, in the case of (C0), even though the trajectory is chaotic accordingly to its MEGNO value, we observe that such character does not affect the overall dynamics, both a 2 and e 2 evolving in a steady fashion. Instead, (C1) exhibits a strong chaotic behavior that drives the system into eventual close approaches (that will later culminate in a collision between the orbiting bodies at t ≈ 4 × 10 5 yrs). The intermittent behavior in the evolution of a 2 may be justified by incursions of the trajectory into high-order resonances' domains.
With respect to the partition of the actions' phase space, we adopted beforehand a rectangular region in the (a, e)-plane of side dimensions (2 a, 2 e) and centered at the respective initial pair of elements (a 0 , e 0 ) for each trajectory. Though the partition box will actually be defined in the action plane (L 2 , G 2 ), our initially regarding a partition on the (a, e))-plane permits the adoption of the Hill's criterion for global stability. In this direction, we considered a minimum distance h = a 2 − a 1 beyond which the system could be considered macroscopically stable and close encounters would not occur at least over an extended time-span (possibly comparable with the host star's age) (Marzari 2014). Thus, from Marchal and Bozis (1982) it is: We then adopted the boundaries a 1 = a 2 = h and e 1 = e 2 = 0.5, both singular cases e − e < 0 and e + e > 1 being avoided by internal computational conditions. In But in our experiments the actual partition box is taken onto the action plane (L 2 , G 2 ), its size being determined by the maximum and minimum values of both L 2 and G 2 corresponding to the extreme values of a and e. Inside this region, we define the partition with r = r a × r e bi-dimensional cells of individual area (L 2 max − L 2 min )(G 2 max − G 2 min )/r . We emphasize that a partition is defined in the action plane of both planets and consequently, a single diffusion coefficient D S is estimated for each body which is a remarkable feature of the method.
The partition boxes corresponding to the outer planet for (C0) and (C1) are depicted in the bottom panel in Fig. 8, where the black dots indicate the successive intersections of each trajectory with the action plane. The diffusive character of (C1) is clearly revealed, the iterates spreading over regions beyond the predefined boundaries. Moreover, a rather fast departure from the center of the box occurs as compared to (C0), for which a quite concentrated configuration accounts for a lower diffusion. We emphasize that, in our method, eventual excursions beyond the boundaries of the box are not taken as instantaneous escape/collision of the trajectories. In fact, for the computation of the entropy, those orbital points that escape Fig. 9 Time evolution of S/ ln r (upper panels) and the estimated D S,2 (lower panel) obtained with an ensemble of initial conditions around (C0). We present the outcomes for three values of r . All panels show the individual curves obtained from both the original and the "ghost" trajectories (brown) and the resulting evaluation for the whole ensemble (black) from the partition box are reintroduced in the box area following the scheme discussed in Beaugé and Cincotta (2019).
Let us now turn to our main scope, the dependence of the Shannon approach on the partition size. To such an end we have integrated, not only the IC but a set of n e = 5 planetary pairs with a and e with small displacements (δa 2 , δe 2 ) of order 10 −6 around the respective initial values for both (C0) and (C1).
Thus, we considered the actual system plus a set of n e = 5 "ghost" systems and followed their evolution by recourse to the numerical integration of the equations of motion for a total span of T = 10 5 yrs, with a setpsize h = 10 −1 yr. For the n e + 1 systems, we computed both the Shannon entropy S and the diffusion coefficient D S . The results are displayed in Figs. 9 and 10 for (C0) and (C1), respectively.
As discussed above, in the above experiments the density of iterates seems to obey neither a Poissonian nor a delta distribution, so we relay in numerical experiments rather than in any analytical estimates.
For each IC, we present the outcomes obtained with three different values for the partition size, namely r (1) = 100 × 100, r (2) = 800 × 800 and r (3) = 1600 × 1600. Let us note that the integration of the whole bunch of trajectories in the 3BP is so demanding that restrains the choice of a larger partition, r (3) being the maximum value allowed by the CPU's memory in our case.  Fig. 9 obtained for trajectories integrated from IC (C1). The brown curves show the j-th "ghost" system of the ensemble, with 1 ≤ j ≤ 5. The single black curves represent the behavior corresponding to the whole ensemble of systems In Fig. 9, we observe that the computations for both S and D S corresponding to each individual trajectory (brown curves) cannot be discriminated from those obtained for the whole ensemble (indicated with a dotted black line) regardless the size of the partition. As expected, D S → 0 since the motion around this IC is stable as revealed in Fig. 8.
In Fig. 10 instead, the temporal evolution of S and D S for the ensemble of systems (in black) departs from the curves for the different individual trajectories (in brown) which can be interpreted in terms of the chaotic nature of the IC (C1). Indeed, slight variations in the initial values selected in a small neighborhood of (C1) result in trajectories that visit a different set of cells each. As a consequence the S values corresponding to the ensemble as a whole account for a larger number of occupied cells. The plots in the bottom panel confirm that an ensemble of systems needs to be considered in order to obtain a reliable value for the diffusion coefficient, contrary to what was found in a similar application to the restricted three-body problem Beaugé and Cincotta (2019).
Furthermore, let us remark that from our numerous experiments we observe that the outcoming values of D S varies slightly with the size of the partition for r values larger than ∼ 100 × 100.
As in the former application, an instability time, τ ins , can be obtained by (11) from the computed diffusion coefficient. In this case, we adopted first K = 1, with 2 = (L 2 max − L 2 0 ) 2 + (G 2 max − G 2 0 ) 2 , where L 2 0 and G 2 0 are the respective initial values defining the center of the partition box in the action plane. We carried out a series of 882 simulations by integrating the IC (C1) varying the parameters (T , h, r , δa, δe, n e ) in order to analyze how such selections affect the results obtained with the Shannon approach. The sensitivity of the computed instability time for the outer planet to the variation of the total number of points obtained by the numerical integration, T /h, the number of cells in both a and e, r a and r e so that r = r a × r e , the size of the dispersion in the initial values for a and e in the ensemble, δa and δe, and the number of ghost systems considered, n e , can be observed in Fig. 11 in case of IC (C1). We recall that after defining the values δa and δe the concomitant values δL 2 and δG 2 are to be computed. Each one of the total of 882 simulations consisted of picking up at random a time-span of integration T ∈ [3 × 10 4 , 10 5 ] yrs, time-steps h ∈ [10 −2 , 2 × 10 −1 ] yr, taking n e = 5, 14 ghost systems, r a , r e ∈ [500, 1600], and δa, δe ∈ [10 −6 , 10 −3 ]. In the light of our experiments no correlation is observed between τ ins and the free parameters, the average value of τ ins ≈ 2 × 10 6 with a dispersion of one order of magnitude. Therefrom, we concluded that the Shannon formulation seems to provide a quite robust estimate of the instability time, regardless the parameters involved in the procedure, since the observed fluctuations are within the intrinsic dispersions arising in N-body simulations, as we show below.
We compared the Shannon estimates of a macroscopic instability time-scale with the life-time of the system obtained from an N-body simulation. We are concern about the intrinsic numerical dispersion as a consequence of taking a small neighborhood around a  (11) given unstable, chaotic initial condition. To this end, we numerically integrated up to T = 10 9 yrs a set of 841 nearby ICs close to (C1), namely a 29 × 29 grid in the range [−10 −4 , 10 −4 ] in both directions of the (a 1 , e 1 )-phase plane. Thus, the histogram on the left panel in Fig. 12 was generated, t dis being defined as the instant at which the integration process stops, because either an eventual collision or some kind of planetary scattering occurs. We observe an almost symmetric histogram around t dis ∼ 10 6 yr, with a dispersion 1 order of magnitude.
For the sake of comparison, we have included in the right panel of Fig. 12 the histogram corresponding to the estimates of the instability time τ inst by (11), for two different values of K . Such a coefficient accounts for any eventual dependence of D S on the diffusion path in the action plane. We notice that while for K = 1 we obtain a more probable value of τ inst ∼ 2 × 10 6 yrs, when adjusting the K value there results τ inst ∼ 7 × 10 5 yrs which is in quite good agreement with the more probable value of the disruption time. Notice should be taken that for these two histograms the system was integrated at most for a total motion time of 1.5 × 10 5 yrs.
The results from the N-body integrations exhibit a "natural" dispersion in the disruption times of the system which can be interpreted in terms of the intrinsic chaotic nature of the solutions, which affects both the Shannon and N-body approaches. Furthermore, let us emphasize the qualitative and quantitative concordance between the histograms in Fig. 12.

Discussion
In the present effort, we briefly review the Shannon entropy approach to measure chaotic diffusion and investigate the dependence of the estimates on the parameters. In particular, we provide theoretical and empirical arguments for the optimal choice of the partition, a free parameter of the method. The results given along this work suggest that both the entropy and the diffusion rate show a mild dependence on the size of the cells.
Quite remarkable is the agreement between an instability time derived from D S and one resulting from straight long-term numerical simulations. While the definition of τ inst involves the constant K , its value is of the order of unity as shown in this work. Moreover, in Cincotta et al. (2021) similar values of K are derived in other applications; the same 4D map but with different values of the parameters and location of the initial conditions and a distinct 3BP model for the HD 181433 system.
The application of the Shannon entropy method to the HD20003 planetary system shows that this approach is a suitable technique to quantify diffusion processes due to the interaction of MMRs. The results presented in Figs. 11 and 12 show the independence of the derived instability time-scale on the parameters used for its estimation. The observed dispersion in the estimate of the disruption time is just a natural consequence of the intrinsic character of unstable chaotic motion; a quite similar result is obtained with N-body simulations.