© The Authors 2026

1 Introduction

Gravitational systems, even for relatively low numbers of degree of freedom, are known to be characterised by the coexistence of regular and chaotic orbits, the latter typically defined as exhibiting strong sensitivity to initial conditions (Contopoulos 2002). The main open issues surrounding chaos in the gravitational N-body problem can be summarised as follows: (1) the relation between Hamiltonian chaos and the onset of dynamical instabilities in the few body systems (Milani & Nobili 1992; Batygin & Laughlin 2008; Mogavero et al. 2023); (2) the dependence of the extent of chaos on the specific orbital structure of a given energy landscape in galactic potentials (Carpintero & Aguilar 1998; Karanis & Caranicolas 2001; Kandrup & Siopis 2003); and (3) the scaling of chaos with the number of degrees of freedom of the systems at hand (Miller 1971; Gurzadyan & Savvidy 1986; Kandrup & Smith 1991; Goodman et al. 1993; Kandrup & Sideris 2001; El-Zant et al. 2019; Di Cintio & Casetti 2019; Portegies Zwart et al. 2023).

Usually, the degree of chaos in an autonomous dynamical system that can be expressed in Hamiltonian from ℋ(x₁, p₁, ..., x_N, p_N) = const is estimated by evaluating numerically its largest Lyapunov exponent (e.g. see Lichtenberg & Lieberman 1992), formally defined by $${\lambda }_{max}=\underset{t\to \mathrm{\infty }}{lim}\underset{\parallel {\mathbf{W}}_0\parallel \to 0}{lim}\frac{1}{t}\ln \frac{\parallel \mathbf{W}(t)\parallel }{\parallel {\mathbf{W}}_0\parallel }.$$(1)

In the equation above, W = (δx₁, δp₁, ..., δx_N, δp_N) denotes the state vector in the tangent space, 𝒮_T, to the phase space, 𝒮, of the system. The dynamics of the tangent vectors is given by $$\delta \dot{\mathbf{x}}=\delta \mathbf{p};\delta \dot{\mathbf{p}}=-{\mathbf{D}}_V^2(\mathbf{x})\delta \mathbf{x},$$(2)

where x = (x₁, ..., x_N), with analogous definitions for δx, p, and δp. We also have $${\mathbf{D}}_V^2(\mathbf{x}{)}_{jk}={\frac{{\mathrm{\partial }}^{2}V(\mathbf{x})}{\mathrm{\partial }{x}_j\mathrm{\partial }{x}_k}|}_{\mathbf{x}};j,k=1,2,\dots ,N,$$(3)

where V(x) is the potential part of the Hamiltonian ℋ = p²/2 + V(x).

If computing the Hessian matrix D²V is particularly cumbersome or if the system is non-autonomous (i.e. ℋ depends explicitly on time), we would typically replace the norm of the tangent vector, ∥W(t)∥, in Eq. (1) with the phase-space distance between two initially nearby trajectories, ξ(t | x₀) and ξ′(t | x₀ + δx₀) (Benettin et al. 1980). In this approach, the numerically estimated λ_max may depend on the initial normalisation of W₀ unless periodic renormalisation is applied.

This procedure consists of evolving the two trajectories for a finite time interval and then redefining the perturbed trajectory so that its separation from the reference trajectory is reset to the initial magnitude, δx₀, while preserving its direction. Explicitly, at each renormalisation step, we have $${\xi }^{\mathrm{\prime }}(t)\to \xi (t)+\delta x_0\frac{{\xi }^{\mathrm{\prime }}(t)-\xi (t)}{\parallel {\xi }^{\mathrm{\prime }}(t)-\xi (t)\parallel }.$$(4)

This renormalisation ensures that the evolution remains within the linear (tangent) regime; namely, that the perturbation stays infinitesimal. Without such a procedure, the separation between trajectories eventually saturates at a characteristic scale of the system (e.g. the system size or the typical impact parameter between particles), thereby invalidating the measurement of exponential divergence, as noted by several authors (Goodman et al. 1993; Valluri & Merritt 2000; El-Zant et al. 2019; Portegies Zwart et al. 2023).

More generally, the full Lyapunov spectrum can be computed by propagating the variational equations up to a time, τ, after which the tangent vectors are orthonormalised (e.g. via a Gram–Schmidt procedure) to prevent their alignment with the most unstable direction (Wiesel 1993; Christiansen & Rugh 1997; Quarles et al. 2011). The finite-time Lyapunov exponents, λ_i(τ), are then obtained by accumulating the logarithmic growth factors associated with the orthonormalisation via $${\lambda }_i(\tau )=\frac{1}{\tau }({\lambda }_{i-1}(\tau )+\ln \parallel {\mathbf{z}}_i-\sum _{j=1}^i⟨{\mathbf{z}}_i,{\mathbf{z}}_{j-1}^{\mathrm{\prime }}⟩{\mathbf{z}}_{j-1}^{\mathrm{\prime }}\parallel ),$$(5)

up to the desired precision dictated by the chosen convergence criterion. In Eq. (5), z_i = (δx_i, δp_i) denotes the i-th deviation (tangent) vector, initially normalised to unity, and λ₀ = 0.

The largest Lyapunov exponent is only the leading member of the Lyapunov spectrum. While λ_max controls the fastest infinitesimal separation between nearby trajectories, the subleading positive exponents quantify additional unstable directions in tangent space (Ishida et al. 2006; Zhou & Wang 2020). Equivalently, the sum of the first k exponents gives the exponential growth rate of an infinitesimal k-dimensional volume element. Therefore, two systems with comparable λ_max may still differ substantially in the dimensionality and isotropy of their instability. This distinction is particularly relevant for Hamiltonian systems with mixed phase space, where local stretching, transport across phase space, and finite-time trapping near regular structures need not be controlled by a single exponent.

Using the largest Lyapunov exponent as a metric of chaos to investigate its scaling in the continuum limit of the N–body problem (i.e. N → ∞ with asymptotically vanishing individual particle mass), leads to puzzling findings. For a given model, determined by a potential density pair supported by an equilibrium phase-space distribution, f, single particle trajectories, largely independent of their specific energies, will increasingly resemble their parent realisations in the smooth potential of the continuum model, while their largest Lyapunov exponents depend weakly on the number of simulation particles, N.

Kandrup & Sideris (2001) carried out frozen N-body experiments (i.e. test particles propagated through the potential generated by a distribution of N fixed particles sampling a given density profile) and concluded that chaos associated with discreteness effects in the N-body problem should be viewed very differently from the chaos associated with a bulk, possibly non-integrable, potential. More recently, Di Cintio & Casetti (2019) studying test orbits in live direct N-body simulations with N up to ≈10⁵ reached similar conclusions. In particular (see their Figs. 7 and 12) the weak scaling of orbital chaoticity with the number of degrees of freedom (i.e. number of particles at a fixed normalisation) has a different trend for different energies.

With respect to the largest Lyapunov exponent, λ_max, of the full N-body problem, whereas some studies suggest that it increases for increasing N see (Goodman et al. 1993; Hemsendorf & Merritt 2002; Portegies Zwart et al. 2022; Portegies Zwart & Boekholt 2023; Asano & Portegies Zwart 2026, some others (Beraldo e Silva et al. 2019; Di Cintio & Casetti 2019, 2020) indicate instead that λ_max is inversely proportional to N, at least for N ≳ 10⁴. This was previously suggested in the semi-analytical estimations of Gurzadyan & Savvidy (1986) and Gurzadyan & Kocharyan (2009). These manifestly contradicting findings are not reconciled at present. The use of a regularised version of the Coulombian 1/r potential, such as the usual $1/\sqrt{r^2+{ϵ}^2}$, can only partially explain why models with softened interactions have systematically lower Lyapunov exponents due to the absence of hard scattering events¹ among particles. Moreover, convergence studies at fixed system parameters with decreasing values of the softening length ϵ, yield similar values of the Lyapunov exponents below some critical ϵ (see Fig. 2 of Kandrup & Sideris 2001 and Fig. 6 of Di Cintio & Casetti 2019).

These apparently conflicting results indicate that the use of the Lyapunov exponent alone might not be enough to characterise chaos in gravitational N-body systems. Of course, up to this point we have discussed the largest λ. The full spectrum may differ substantially between realisations with different N and similar λ_max. For example Di Cintio & Trani (2025) showed that in the three-body problem with relativistic corrections included in the force calculations, the largest Lyapunov exponents can be comparable to (or even smaller than) their classical counterparts, while the full spectrum may differ significantly between the two cases.

It is therefore worth exploring other metrics of chaos in addition to the Lyapunov exponents and related quantities such as the small alignment indexes (SALI, Skokos 2001), generalised alignment indexes (GALI, Skokos et al. 2007), or the mean exponential growth factor of nearby orbits (MEGNO, Cincotta & Simó 2000; Mestre et al. 2011). For example, several studies in the context of simple dynamical systems (e.g. see Chernov 1997; Papenbrock 2000; Shiozawa & Tokuda 2024) used the Kolmogorov-Sinai (KS, Kolmogorov 1958, 1959; Sinai 1959) entropy 𝒮_KS (see Lichtenberg & Lieberman 1992). Once a coarse graining of time and phase-space is defined, the KS entropy quantifies the rate of information production along a trajectory, as inferred from the conditional probabilities that the trajectory occupies a given phase-space cell, j, at time, t_n+1, given its history up to time t_n. Remarkably, Pesin (1977) showed that for smooth dynamical systems with an invariant measure, such as Hamiltonian systems obeying Liouville’s theorem, the KS entropy is bounded from above by the sum of the positive Lyapunov exponents, $${\mathcal{S}}_{\mathrm{K}\mathrm{S}}\le {\mathcal{S}}_{\mathrm{K}\mathrm{S}}^{+}=\sum _{{\lambda }_i>0}{\lambda }_i.$$(6)

When the dynamics is sufficiently chaotic, this bound becomes an equality, known as Pesin’s identity. In systems with mixed phase space, where regular and chaotic regions coexist, this relation should be understood as applying to the chaotic component in the asymptotic limit.

Another commonly used entropy measure is the Shannon information entropy (Shannon & Weaver 1949). Given a partition of phase space into cells k, with occupation probabilities, p_k(t), for instance, estimated from an ensemble of trajectories at a fixed time t), the Shannon entropy is $${\mathcal{S}}_{\mathrm{S}\mathrm{h}}(t)=-\sum _k{p}_k(t)\ln p_k(t),$$(7)

which quantifies how broadly the ensemble is distributed over the chosen coarse-grained representation. The Shannon and KS entropy are conceptually similar in that both depend on coarse graining and capture information content. However, 𝒮_Sh(t) is an instantaneous property of a probability distribution (and can be tracked in time if the distribution evolves), whereas 𝒮_KS is an asymptotic entropy rate associated with the time evolution of a single trajectory, defined from the statistics of sequences of its visits in phase space cells.

The use of entropy-based quantities as indicators of relaxation and chaos has been rather successful in dynamical astronomy (e.g. Romeo 1990). In a series of studies (Núñez et al. 1996; Cincotta et al. 1999; Giordano & Cincotta 2018; Cincotta et al. 2021), the information entropy was employed to quantify the dynamical stability of ensembles of trajectories in the restricted three-body problem and to identify periodic and chaotic patterns in astronomical data.

These diagnostics are complementary to frequency map analysis, which characterises orbital regularity through the Fourier decomposition of time series and the (in)stability of the associated fundamental frequencies. Originally introduced in planetary dynamics (Laskar 1990), frequency map analysis has since been widely adopted in accelerator beam dynamics (Shatilov et al. 2011) and galactic dynamics (Valluri et al. 2012; Vasiliev 2013; Bajkova et al. 2023; Woudenberg & Helmi 2025). More recently, Hyman et al. (2025) introduced a framework for orbital-complexity analysis based on permutation entropy and statistical complexity measures, while Canbaz (2025) proposed a novel indicator called a permutation entropy of the power spectrum, combining permutation entropy with Fourier-based methods. Similarly, Bajkova et al. (2025) computed the information entropy of the discrete Fourier transform of time series of radial Galactic distances to classify the chaoticity of globular cluster orbits perturbed by the Galactic bar.

In this work, we investigate the relationship between different chaos indicators by comparing the largest Lyapunov exponent and the Shannon entropy in two idealised systems: a test particle in the Hénon-Heiles potential and in the potential generated by N particles distributed according to a Plummer profile. We assess the degree to which these indicators yield consistent classifications of regular and chaotic dynamics, and we identify regimes where their interpretations diverge. We also examine the dependence of these indicators on the particle number, N, with the aim of clarifying the apparently contradictory results on the relation between chaoticity and the number of degrees of freedom in the N-body problem.

The rest of the paper is structured as follows. In Sect. 2, we introduce the governing equations, discuss the numerical integration and describe the procedure to evaluate the Lyapunov exponents and the information entropy. In Sect. 3, we present our numerical simulations and discuss the results. Finally, in Sect. 4 we draw our conclusions and interpret our findings in light of previous work.

2 Methods

2.1 Models

In this work, we consider two models, specifically, the Hénon & Heiles (1964, hereafter HH) system and a spherical N-body system in which the motion of a single particle is studied within the time-dependent gravitational potential generated by the remaining N − 1 particles.

2.1.1 Hénon-Heiles system

The HH system is defined by the Hamiltonian, $${\mathcal{H}}_{\mathrm{H}\mathrm{H}}=\frac{p_x^2+p_y^2}{2}+\frac{x^2+y^2}{2}+x^{2}y-\frac{y^3}{3},$$(8)

and it was originally conceived as a toy model for the dynamics in the meridional plane of an axisymmetric galactic potential. It has since become a paradigmatic example in nonlinear dynamics, exhibiting a mixed phase space; namely, a nontrivial coexistence of regular and chaotic orbits on the same energy surface (Weiss et al. 2003). In this work, we adopted the conventional rescaling in which the particle mass m and the coefficients of all terms in the 2D potential are set to unity. With this normalisation, the critical energy above which all trajectories escape is E_crit = 1/6. The equations of motion derived by Eq. (8) and the associated variational equations, defined by Eq. (3), is expressed as $$\ddot{x}=-x-2xy,$$(9) $$\ddot{y}=-x^2+y^2-y,$$(10)

and $$\ddot{\delta x}=-(1+2y)\delta x-2x\delta y,$$(11) $$\ddot{\delta y}=-2x\delta x-(1-2y)\delta y,$$(12)

respectively. Figure 1 illustrates a typical regular and chaotic orbit in the HH system, corresponding to E = 0.04 and E = 0.16.

2.1.2 N-body system

We probe the dynamics of individual orbits in spherical equilibrium self gravitating N-body models evolving test particles in the time-dependent gravitational potential exerted by the remaining N − 1 bodies. The equation of motion for the i-th particle of mass, m, is $${\ddot{\mathbf{r}}}_i=-Gm\sum _{i\ne j=1}^N\frac{{\mathbf{r}}_i-{\mathbf{r}}_j}{r_{ij}^3},$$(13)

where G is the gravitational constant and r_ij = ∥r_i − r_j∥. The associated variational equations in tangent space, required for the computation of Lyapunov exponents (see, e.g. Rein & Tamayo 2016), is expressed as $$\ddot{\delta }{\mathbf{r}}_i=-Gm\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N[\frac{\delta {\mathbf{r}}_i-\delta {\mathbf{r}}_j}{r_{ij}^3}-3({\mathbf{r}}_i-{\mathbf{r}}_j)\frac{⟨\delta {\mathbf{r}}_i-\delta {\mathbf{r}}_j,{\mathbf{r}}_i-{\mathbf{r}}_j⟩}{r_{ij}^5}].$$(14)

In the N-body experiments, the initial positions and velocities of the N particles are drawn from the widely adopted Plummer (1911) model, with total mass, M, and scale radius, a, and defined by the following spherical density-potential pair: $$\rho (r)=\frac{3M}{4\pi a^3}{(1+\frac{r^2}{a^2})}^{-5/2};\mathrm{\Phi }(r)=-\frac{GM}{\sqrt{r^2+a^2}},$$(15)

supported by the ergodic phase-space distribution function, $$f(E)=\frac{\sqrt{2}}{378{\pi }^{2}G{a}^2{\sigma }_0}{\left(\frac{-E}{{\sigma }_0^2}\right)}^{7/2}.$$(16)

Here, E denotes the specific energy of a particle and σ₀ = $\sqrt{GM/6a}$ is the central velocity dispersion. All quantities are expressed in N-body units, such that G = M = a = 1 and each particle has equal mass, m = M/N. Consequently, the typical crossing time $t_{\mathrm{c}}\equiv \sqrt{a^3/GM}=1$. We simulate different sets of initial conditions with N between 10⁴ and 10⁶. As N increases, the total mass is held fixed, consequently, the individual particle mass decreases accordingly.

We approximate the dynamics by evolving N − 1 independent “field” particles in the smooth potential (15), while computing the trajectory and associated tangent dynamics of a “test” particle according to Eqs. (13)–(14). In this framework, the variational vectors of the field particles, $\ddot{\delta }{\mathbf{r}}_j$, evolve in the Plummer potential according to $$\ddot{\delta }{\mathbf{r}}_j=-Gm[\frac{\delta {\mathbf{r}}_j}{{(r_j^2+a^2)}^{3/2}}-3{\mathbf{r}}_j\frac{⟨\delta {\mathbf{r}}_j,{\mathbf{r}}_j⟩}{{(r_j^2+a^2)}^{5/2}}].$$(17)

Fig. 1 Illustration of two trajectories in the HH system. Top panel: low-energy, regular trajectory with a nearly zero largest Lyapunov exponent, λmax. Bottom panel: high-energy, chaotic trajectory with finite λmax.

2.2 Numerical integration

For both the HH system and the N-body system, we integrated the equations of motion, Eqs. (9) and (13), together with their associated variational equations (Skokos & Gerlach 2010), using the fourth-order symplectic integration scheme discussed in Laskar & Robutel (2001, see also Yoshida 1990, 1993; Kinoshita et al. 1991 and references therein). In all runs, we adopted a fixed time step, Δ, t that depends on the initial orbital energy. The time step was determined through preliminary trial-and-error integrations, requiring that over an integration time of 5× 10³ time units the relative energy error remains below 10⁻⁹ in double precision for the HH system. In the N-body runs, for the dimensionless equations of motion considered here, the resulting optimal timesteps typically lie in the range 2 × 10⁻⁵ ≤ Δt ≤ 10⁻². As we are following individual test trajectories in time dependent N-body potentials, we did not employ any softening or regularisation of the gravitational interactions.

2.3 Evaluation of the Lyapunov exponents

Evaluating the largest Lyapunov exponent of a given system by simply time-dependently computing (1) up to a desired convergence is typically unfeasible as it is in principle plagued by round-off errors resulting in a rapid blow-up of the numerical solution. Following Benettin et al. (1976) in the numerical calculations we obtain λ_max as the limit of $${\lambda }_{max}(t)=\frac{1}{L\mathrm{\Delta }t}\sum _{k=1}^L\ln \frac{\parallel \mathbf{W}(k\mathrm{\Delta }t)\parallel }{\parallel {\mathbf{W}}_0\parallel },$$(18)

for a (large) time t = LΔt with the additional precaution of periodically re-normalising W to its initial magnitude (see e.g. Skokos 2010 and references therein). When computing λ_max for an individual particle in N–body experiments, we also average over different experiments where the other N − 1 background particles are sampled from the same model in a similar fashion to what done in Sartorello et al. (2025) for particles in time-dependent gravitational potential subjected to noise.

2.4 Information entropy

We computed the information entropy in two distinct ways. First, for individual orbits we evaluate Eq. (7) from the empirical phase-space density in the full state vector (x, y, p_x, p_y)(t), obtained by coarse-graining the trajectory samples into a multi-dimensional histogram. Using ${\hat{p}}_k$ to denote the fraction of stored samples that fall in bin, k, our estimator is $${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}\equiv -\sum _k{\hat{p}}_k\ln {\hat{p}}_k,$$(19)

where the sum runs over all bins (with the convention that empty bins contribute zero). To mitigate the dependence on an arbitrary bin choice, we estimated a reference bin width in each coordinate using the Freedman-Diaconis rule and then scanned a range of nearby binnings to verify the stability of the inferred entropy. We enforced a minimum sampling quality by requiring an average occupancy of at least five samples per non-empty bin (i.e. N_samp/N_occ ≥ 5), where N_samp is the number of stored phase-space samples and N_occ is the number of occupied bins. For consistent entropy comparisons, a common binning is adopted for all trajectories within a given comparison set. We also explored nearby binning configurations and verified that the qualitative trends, as well as the quantitative values within uncertainties, remain unchanged. In this definition, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ quantifies how broadly the trajectory populates the accessible phase space (at the chosen coarse-graining). This is therefore a measure of phase-space exploration, rather than of instantaneous information-production rate.

Second, we evaluated the mutual information entropy introduced by Núñez et al. (1996, see also Fraser & Swinney 1986), which we denote with ℐ_N. Although its formal definition resembles that of a mutual information Shannon entropy, the quantity ℐ_N is not a mutual information measure in the probabilistic sense. Instead, it is constructed by calculating the Shannon entropy of normalised arc-length weights along two trajectories in phase space and by combining the resulting entropies in a way that is sensitive to the dynamical de-correlation between nearby realisations. Given a reference trajectory ξ(t) and a nearby realisation, ξ′(t), initialised at a separation of ∥ξ′(0) − ξ(0)∥ = d₀, Núñez et al. (1996) define $${\mathcal{I}}_{\mathrm{N}}(t)={\mathcal{S}}_{\mathrm{N}}(\xi ,{\xi }^{\mathrm{\prime }};t)-\frac{1}{2}[{\mathcal{S}}_{\mathrm{N}}(\xi ;t)+{\mathcal{S}}_{\mathrm{N}}({\xi }^{\mathrm{\prime }};t)],$$(20)

where the entropies are computed from normalised path-length rates along the curves. For a single trajectory, ξ(t), we define the instantaneous speed in phase space, $${\dot{\ell }}_{\xi }(t)=\parallel \dot{\xi }(t)\parallel ,$$(21)

the accumulated path length up to time, t, $$L_{\xi }(t)={\int }_0^t{\dot{\ell }}_{\xi }\left(t^{\mathrm{\prime }}\right)\mathrm{d}{t}^{\mathrm{\prime }},$$(22)

and the corresponding normalised length density, $${\epsilon }_{\xi }(t)=\frac{{\dot{\ell }}_{\xi }(t)}{L_{\xi }(t)}.$$(23)

The path-length entropy for an individual trajectory is then the Shannon entropy of this time-density, $${\mathcal{S}}_{\mathrm{N}}(\xi ;t)=-{\int }_0^t{\epsilon }_{\xi }(t^{\mathrm{\prime }})\ln {\epsilon }_{\xi }(t^{\mathrm{\prime }})\mathrm{d}{t}^{\mathrm{\prime }}.$$(24)

For the joint entropy estimator 𝒮_N(ξ, ξ′; t), we define $${\dot{\ell }}_{\xi {\xi }^{\mathrm{\prime }}}(t)=\sqrt{{\dot{\ell }}_{\xi }^2(t)+{\dot{\ell }}_{\xi }^2(t)},$$(25) $$L_{\xi {\xi }^{\mathrm{\prime }}}(t)={\int }_0^t{\dot{\ell }}_{\xi {\xi }^{\mathrm{\prime }}}(t^{\mathrm{\prime }})\mathrm{d}{t}^{\mathrm{\prime }},$$(26) $${\epsilon }_{\xi {\xi }^{\mathrm{\prime }}}(t)=\frac{{\dot{\ell }}_{\xi {\xi }^{\mathrm{\prime }}}(t)}{L_{\xi {\xi }^{\mathrm{\prime }}}(t)},$$(27)

and substitutes ε_ξξ′ into Eq. (24) to obtain 𝒮_N(ξ, ξ′; t). In contrast to ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$, ℐ_N(t) is sensitive to dynamical divergence: as the perturbed trajectory deviates from the reference one, the joint and marginal length densities evolve differently, producing a characteristic growth and eventual saturation of ℐ_N(t) in chaotic regimes (Núñez et al. 1996). In this sense, ℐ_N plays a role analogous to that of a Lyapunov exponent as a diagnostic of instability, while remaining a finite-time quantity. A practical distinction is that the mutual information entropy as defined by Núñez et al. (1996) cannot, in general, be computed from tangent-vector dynamics alone: its definition depends on the accumulated path length of a genuinely distinct trajectory, which is not preserved under the renormalisation steps required by the Lyapunov-exponent computation algorithm.

Fig. 2 Largest Lyapunov exponent, λmax (top panel), and Shannon entropy, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ (bottom panel), versus energy for 2000 realisations of the HH system. In the top panel, the blue and red dashed lines mark the trends λmax ∝ $\sqrt{E}$ and λmax ∝ E3.4 reported by Shevchenko & Mel’nikov (2003). In the bottom panel, the curves ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ ∝ ln E and ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ ∝ 4 ln E are shown for comparison. Both λmax and ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ increase with energy, with a change in slope around E ≃ 0.08, marking the transition from predominantly regular to chaotic trajectories.

Fig. 3 Average mutual information entropy, ℐN (top panel), and Shannon entropy (bottom panel) versus the largest Lyapunov exponent, λmax, for different values of the energy, E, in the range 0.04 ≤ E ≤ 0.16 in the HH system. Error bars indicate the 1σ scatter over the ensemble, shown by the coloured points. In both panels, dashed lines represent the best-fitting relation between the entropy and λmax.

Fig. 4 Poincaré sections (y, py) of the HH Hamiltonian for different values of the energy, E, defined by crossings of the plane x = 0. The colour coding represents the largest Lyapunov exponent, λmax, of the orbits.

Fig. 5 Poincaré sections (y, py) of the HH Hamiltonian for different values of the energy, E, defined by crossings of the plane x = 0. The colour coding represents the mutual information entropy, ℐN (Eq. (20)), of the orbits.

3 Numerical experiments and results

3.1 HH system

It is well known that in the HH potential the largest Lyapunov exponent λ_max, averaged over an energy surface, increases with E. Early computations by Benettin et al. (1976) suggested an exp E trend for both λ_max and the KS entropy, S_KS. Subsequent studies, for example Shevchenko & Mel’nikov (2003) reported a broken power-law trend instead, with a sharp change in slope at E ≈ 0.08 for the averaged Lyapunov exponent.

In the energy range 10⁻³ ≤ E < 1/6, we evaluated both the largest Lyapunov exponent and the Shannon entropy on a grid of 10³ realisations. The top panel of Fig. 2 shows λ_max as a function of orbital energy E. Orange and blue points highlight the two regimes, indicated by the corresponding dashed lines, where λ ∝ $\sqrt{E}$ and λ ∝ E^3.4, respectively. The Shannon entropy, $\hat{\mathcal{S}}\mathrm{S}\mathrm{h}$, for the same orbits is shown in the bottom panel, and like-wise exhibits a broken trend with E, consistent with ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ ∝ ln E for E ≲ 0.08 and ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ ∝ ln E⁴ for E ≳ 0.08. We note that, in principle, for a given value of E, due to the coexistence of wildly chaotic and nearly regular orbits, the values of the Lyapunov exponents for different trajectories can differ significantly, resulting in highly non-trivial distributions P(λ) (e.g. Vallejo et al. 2008).

In Fig. 3, we present the mutual information entropy, ℐ_N (top panel), in addition to the Shannon entropy ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ (bottom panel) against the mean Lyapunov exponent. We find that ℐ_N scales linearly with λ_max, while the Shannon entropy is compatible with a ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ ∝ ln λ_max relation, in agreement with the picture emerging from Fig. 2.

The correspondence between the mutual information entropy and the Lyapunov exponent can be visualised by colour-coding these indicators in the Poincaré sections of the trajectories at fixed values of E. In Figs. 4 and 5, we show the sections in (y, p_y) defined by x = 0, colour-coded by λ_max and ℐ_N, respectively, for different energies. Both indicators, despite the greater scatter of ℐ_N at nearly equal λ_max, reveal how for E ≳ 0.8 the fraction of accessible phase-space is indeed occupied by a rapidly growing fraction of chaotic orbits though with a persistent measure of regular islands, even for the E = 0.16, closer to the E = 1/6 threshold case.

Fig. 6 Tracer particle orbit with initial energy E ≈ −0.66 in different N-body realisations of an isotropic Plummer model with increasing N, and in the smooth-potential (continuum) limit. The colour coding represents time in units of tc. As N increases, the orbit approaches its continuum counterpart.

Fig. 7 Evolution of the largest Lyapunov exponent λmax for different values of N in the N-body system. Solid lines show the mean over 10 realisations and the shaded regions indicate the 1σ dispersion.

3.2 N-body models

We performed a series of N-dependent numerical experiments on test orbits for different values of the initial particle energy. As mentioned above, numerical studies on frozen N-body realisations of simple galactic potentials (see e.g. Kandrup & Sideris 2001; Sideris & Kandrup 2002; Kandrup et al. 2004) seem to suggest that individual particle trajectories, for increasing N approach more and more closely their counterparts in the parent smooth system, while their largest Lyapunov exponent does not show a significant decreasing trend with N, as we would naively expect.

In the range 10⁴ ≤ N ≤ 10⁶, we observe a similar tendency even in our time-dependent N–body potentials for the orbits of test particles initialised with exactly the same phase space position approaches in realisations of the potential with increasing values of N. As an example, in Fig. 6, we show the trajectory for a particle of initial energy per unit mass E ≈ −0.66 in Plummer models sampled by N = 10⁴, 3 × 10⁴, 10⁵, 3 × 10⁵, and 10⁶. Notably, as N increases, the fluctuations of the orbital inclination become less and less important, being already negligible for N = 10⁶. A similar behaviour is observed for other values of the orbital eccentricity and initial energy (not shown here), with highly eccentric orbits at fixed E (i.e. having lower values of angular momentum) being generally more subjected to wild fluctuations in the inclination, even at N ≈ 10⁵.

We note that λ_max evaluated as function of a given value of the energy E depends weakly on N (see also Fig. 2 in Kandrup & Sideris 2001, Fig. 12 in Di Cintio & Casetti 2019 and Fig. 10 in Sartorello et al. 2025). This is evidenced in Fig. 7 where we tracked the evolution of the mean finite time largest Lyapunov exponent for the same orbital initial condition (i.e. test particle with equal combinations of r and v), propagated in different realisations of the Plummer model at different values of N. Using the same number of independent samplings of the underlying model (shaded areas) the mean values of λ_max over the ensemble relax to comparable values over 10³ dynamical times, regardless of N.

We emphasise that unlike in the frozen-potential experiments of Kandrup and collaborators, the energy of the test particle is not conserved for any N, owing to the time-dependent potential generated by the N moving field particles. It is therefore necessary to quantify the range of energies explored along each trajectory.

In the top panel of Fig. 8, we show λ_max as a function of the median energy over the full integration. Orbits at lower energies (i.e. particles deep in the potential well) have larger values of their largest Lyapunov exponent, whereas the orbits exploring higher energies (i.e. weakly bound particles) systematically show smaller λ_max. This trend is consistent with the relation between the averaged λ_max and the initial energy found in noisy-potential experiments mimicking an increasing number of field particles (up to 10¹³) by Sartorello et al. (2025), as well as in low-resolution direct N-body simulations by Di Cintio & Casetti (2019). Furthermore, as N increases, the scatter in λ_max about its relation with the median energy decreases.

Because the particle energy is not conserved, we compared λ_max against both the initial and median energy at fixed N. As shown in the lower panels of Fig. 8, both λ_max and ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ decreases monotonically with both quantities. Notably, for weakly bound particles (E ≳ −0.3), the initial and median energies (green crosses) are nearly identical, as expected for orbits confined to the lower-density outskirts of the cluster.

With this established, the top and middle panels of Fig. 9 show the violin plots of λ_max and the information entropy, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$, respectively, including their statistical spread, for orbits with initial energy E = −0.45 (corresponding to a circular orbit of radius a in the smooth Plummer potential). The orbits are evolved across ensembles of discrete realisations with increasing N.

The distribution of the largest Lyapunov exponent shows no significant dependence on N over this range, while its spread narrows markedly with increasing N. At first sight, this appears inconsistent with the ln(ln N) scaling proposed by Goodman et al. (1993). However, their estimate was derived assuming the propagation and amplification of independent perturbations through strong two-body encounters in a fully self-gravitating N-body system. In contrast, our setup follows a test particle evolving in the live potential generated by N particles and, therefore, it does not capture the same self-consistent amplification mechanism. In contrast, the entropy decreases with N and exhibits a broader spread, particularly at N ~ 10⁶. This trend is also evident in the bottom panel, which shows ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ as a function of λ_max for individual realisations. At fixed N, larger Lyapunov exponents are broadly associated with higher entropy, although the correlation is weaker than it is for the HH system.

Fig. 8 Top panel: largest Lyapunov exponent, λmax, versus the orbital energy, E, of the test particle in the N-body system, where E is taken as the median along the trajectory. Colors indicate the number of particles, N. Lower panels: λmax (upper sub-panel) and Shannon entropy, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ (lower sub-panel), versus the initial (blue crosses) and median (green crosses) orbital energy in a realisation of an isotropic Plummer model at fixed N = 105. Both λmax and ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ decrease with increasing orbital energy and approach a saturation level for the most tightly bound orbits.

Fig. 9 Top panels: distributions of the largest Lyapunov exponent, λmax (upper sub-panel), and Shannon entropy, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ (lower sub-panel), for different values of N in the N-body system. Error bars indicate the median and the associated 1σ scatter. In the upper sub-panel, the dashed line denotes the scaling proposed by Goodman et al. (1993). Bottom panel: ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ versus λmax for individual realisations of the N-body system. The largest Lyapunov exponent shows no significant dependence on N, whereas the Shannon entropy decreases with increasing N.

4 Discussion and conclusions

We investigated the relation between the largest Lyapunov exponent and different formulations of information entropy in the paradigmatic Hénon & Heiles (1964) 2D potential and for individual particle orbits in time-dependent N-body potentials. Our main results can be summarised as follows.

For the HH model, the Shannon entropy averaged over trajectories on the same energy surface scales logarithmically with energy, with different pre-factors in the ≲0.08 and ≳0.08 regimes (Fig. 2). This change in the pre-factor reflects the transition from weakly to strongly chaotic motion across the well known onset of widespread chaos in the HH system. This behaviour is consistent with that of the largest Lyapunov exponent and of the mutual information entropy. The latter scales approximately linearly with λ_max, whereas ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ exhibits a logarithmic dependence (Fig. 3).

For single-particle orbits in N-body realisations of a given gravitational potential (here, a spherical Plummer model), λ_max is larger for more tightly bound orbits (i.e. lower energies), while it decreases for less bound trajectories; the same qualitative trend is observed for the Shannon entropy (Fig. 8). At fixed initial conditions, increasing the resolution of the N-body system (i.e. larger N at fixed total mass) produces little variation in λ_max, even as the dynamics approaches the continuum limit (Fig. 9). In contrast, the Shannon entropy decreases monotonically with N, largely independently of E, in qualitative agreement with frequency-map analyses of frozen N-body systems (Kandrup & Sideris 2001; Sideris & Kandrup 2002). This discrepancy suggests that λ_max alone is not a sufficient indicator of chaoticity in gravitational N-body systems, at least if chaoticity is understood in terms of finite-time phase-space exploration. The maximum Lyapunov exponent measures the fastest local rate of infinitesimal stretching in tangent space. In contrast, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ measures the finite-resolution volume effectively sampled by an orbit over the integration time. These two quantities are expected to correlate when local stretching efficiently produces transport across the accessible region of phase space and when the dominant stretching direction is representative of the overall instability. However, this connection need not hold in general. In this sense, the entropy indicator is not a replacement for the Lyapunov exponent, but a complementary diagnostic sensitive to the macroscopic region explored by the trajectory. Although ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ is not an estimator of the KS entropy per se, it is motivated by the same connection between instability, information production, and phase-space exploration. Its behavior could therefore reflect dynamical information that is not contained in λ_max alone, but is instead associated with the collective contribution of the positive Lyapunov exponents.

We conjecture that as N increases, λ₁ ≡ λ_max remains approximately constant (at least up to N = 10⁶), while the second and third exponents λ₂ and λ₃ decrease significantly beyond N ~ 10⁵. This interpretation is supported by Fig. 6 (N = 10⁶), where perturbations in the precession frequency persist, but fluctuations of the orbital plane are strongly suppressed, indicating a near conservation of the angular momentum. We note that in a related context, the inclusion of post-Newtonian corrections can decrease the leading exponent, while increasing subdominant ones in certain configurations of the three-body problem (Di Cintio & Trani 2025).

If the decrease in ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ with N is driven by the suppression of sub-leading unstable directions, or more generally, by the progressive smoothing of the finite-N potential, then we should expect ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ to approach an asymptotic value at a sufficiently large N. This limiting value need not be connected to the saturated value of λ_max. Thus λ_max may remain approximately constant, while ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ decreases and eventually saturates to the value associated with the corresponding smooth-potential orbit, for the adopted integration time and coarse graining. In a smooth spherical Plummer potential, this baseline is expected to reflect regular motion constrained by the integrals of the mean-field problem, rather than continued chaotic exploration.

A direct test would require computing the full Lyapunov spectrum (Eq. (5)) for the test particle. However, this entails orthonormalisation in a 6N-dimensional phase space, which requires evolving and orthonormalising a 6N × 6N tangent basis and becomes computationally prohibitive already at N ≈ 10⁴. An alternative is to track the deformation of an initially hyperspherical ensemble of trajectories, as proposed by Sartorello et al. (2025).

In light of these results, we speculate that the N-body chaos driven by close encounters (sometimes referred to as punctuated chaos, see Portegies Zwart et al. 2023; Boekholt et al. 2023), which appears to increase with N (Portegies Zwart & Boekholt 2023), primarily reflects the behaviour of the largest Lyapunov exponent, which probes local instability driven by close encounters, rather than global phase-space transport. In contrast, a dynamical entropy measure, being sensitive to global mixing through the full Lyapunov spectrum, would decrease with N, owing to the suppression of subdominant exponents. We emphasise that trajectory-based entropy measures provide a viable alternative to λ_max only in low-dimensional systems (e.g. driven 1D models or 2D non-integrable potentials such as the HH system), where the sum of the positive part of the Lyapunov spectrum closely traces λ_max. In higher dimensional phase spaces, quantities such as ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ can instead serve as diagnostics indicating whether a full spectral analysis is required.

Finally, alternative estimates of dynamical entropy can be obtained from the compressibility of time series (e.g. from trajectory data), for instance via Lempel-Ziv-type complexity measures (Lempel & Ziv 1976; Puglisi et al. 2003). A natural extension of this work is to compare, for fixed system parameters (e.g. energy), trajectory-based entropy estimates with those derived from compression schemes. This approach may be particularly useful when variational equations are unavailable or computationally expensive or when dealing with observational data (e.g. minor-body ephemerides or Galactic-centre S-stars), where only finite time series are accessible.

Acknowledgements

We thank the anonymous referee for their insightful comments, which have improved the manuscript. PFDC acknowledges the support from the MUR PRIN2022 project “Breakdown of ergodicity in classical and quantum many-body systems” (BECQuMB) Grant No. 20222BHC9Z. AAT and PFDC are grateful to Anna Lisa Varri and the University of Edinburgh for their hospitality during the 2nd “Chaotic rendezvous” meeting, and to the IFPU for hosting the workshop “Chaos and nonlinearity in dynamical astronomy”. AAT dedicates this work to the memory of his father, Vincenzo Trani, a mathematics teacher who inspired his interest in science. This study was completed during a period of personal recovery, and he gratefully acknowledges the support received during this time.

References

All Figures

	Fig. 1 Illustration of two trajectories in the HH system. Top panel: low-energy, regular trajectory with a nearly zero largest Lyapunov exponent, λmax. Bottom panel: high-energy, chaotic trajectory with finite λmax.
In the text

Fig. 2 Largest Lyapunov exponent, λmax (top panel), and Shannon entropy, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ (bottom panel), versus energy for 2000 realisations of the HH system. In the top panel, the blue and red dashed lines mark the trends λmax ∝ $\sqrt{E}$ and λmax ∝ E3.4 reported by Shevchenko & Mel’nikov (2003). In the bottom panel, the curves ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ ∝ ln E and ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ ∝ 4 ln E are shown for comparison. Both λmax and ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ increase with energy, with a change in slope around E ≃ 0.08, marking the transition from predominantly regular to chaotic trajectories.

In the text

Fig. 3 Average mutual information entropy, ℐN (top panel), and Shannon entropy (bottom panel) versus the largest Lyapunov exponent, λmax, for different values of the energy, E, in the range 0.04 ≤ E ≤ 0.16 in the HH system. Error bars indicate the 1σ scatter over the ensemble, shown by the coloured points. In both panels, dashed lines represent the best-fitting relation between the entropy and λmax.

In the text

	Fig. 4 Poincaré sections (y, py) of the HH Hamiltonian for different values of the energy, E, defined by crossings of the plane x = 0. The colour coding represents the largest Lyapunov exponent, λmax, of the orbits.
In the text

	Fig. 5 Poincaré sections (y, py) of the HH Hamiltonian for different values of the energy, E, defined by crossings of the plane x = 0. The colour coding represents the mutual information entropy, ℐN (Eq. (20)), of the orbits.
In the text

	Fig. 6 Tracer particle orbit with initial energy E ≈ −0.66 in different N-body realisations of an isotropic Plummer model with increasing N, and in the smooth-potential (continuum) limit. The colour coding represents time in units of tc. As N increases, the orbit approaches its continuum counterpart.
In the text

	Fig. 7 Evolution of the largest Lyapunov exponent λmax for different values of N in the N-body system. Solid lines show the mean over 10 realisations and the shaded regions indicate the 1σ dispersion.
In the text

Fig. 8 Top panel: largest Lyapunov exponent, λmax, versus the orbital energy, E, of the test particle in the N-body system, where E is taken as the median along the trajectory. Colors indicate the number of particles, N. Lower panels: λmax (upper sub-panel) and Shannon entropy, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ (lower sub-panel), versus the initial (blue crosses) and median (green crosses) orbital energy in a realisation of an isotropic Plummer model at fixed N = 105. Both λmax and ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ decrease with increasing orbital energy and approach a saturation level for the most tightly bound orbits.

In the text

Fig. 9 Top panels: distributions of the largest Lyapunov exponent, λmax (upper sub-panel), and Shannon entropy, ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ (lower sub-panel), for different values of N in the N-body system. Error bars indicate the median and the associated 1σ scatter. In the upper sub-panel, the dashed line denotes the scaling proposed by Goodman et al. (1993). Bottom panel: ${\hat{\mathcal{S}}}_{\mathrm{S}\mathrm{h}}$ versus λmax for individual realisations of the N-body system. The largest Lyapunov exponent shows no significant dependence on N, whereas the Shannon entropy decreases with increasing N.

In the text