Contribution of individual degrees of freedom to Lyapunov vectors in many-body systems

TL;DR

This paper introduces a diagnostic based on Lyapunov vectors to identify regions in phase space that contribute most to chaos, demonstrated on a mean field model where it highlights the separatrix vicinity.

Contribution

It presents a novel method using the weight δI from Lyapunov vectors to characterize chaotic regions in many-body systems.

Findings

The diagnostic effectively highlights the vicinity of the separatrix.

It works even when the separatrix is barely present.

The method provides insight into the structure of chaos in the model.

Abstract

We use the weight $\delta$I, deduced from the estimation of Lyapunov vectors, in order to characterise regions in the kinetic (x, v) space with particles that most contribute to chaoticity. For the paradigmatic model, the cosine Hamiltonian mean field model, we show that this diagnostic highlights the vicinity of the separatrix, even when the latter hardly exists.

Equations

$$\mathbf{x}\equiv (x_1,x_2,....x_n)$$

$$\frac{\mathrm{d}\mathbf{x}(t)}{\mathrm{d}t}=\mathbf{F}(\mathbf{x}(t)),$$

$$\mathbf{w}\equiv (\delta x_1,\delta x_2,\dots ,\delta x_n).$$

$$\frac{\mathrm{d}\mathbf{w}(t)}{\mathrm{d}t}=\mathbf{J}(\mathbf{x}(t))\mathbf{w}(t),$$

$$\frac{\mathrm{d}\mathcal{J}({\mathbf{x}}_0;t)}{\mathrm{d}t}=\mathbf{J}(\mathbf{x}(t))\mathcal{J}({\mathbf{x}}_0;t).$$

$$\lambda ({\mathbf{x}}_0,{\mathbf{w}}_0)=\underset{t\to \infty }{\lim }\frac{1}{t}\ln \frac{\Vert \mathbf{w}(t)\Vert }{\Vert {\mathbf{w}}_0\Vert }$$

$${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_n.$$

$${\lambda }_i=-{\lambda }_{n-i+1}.$$

$$H=\sum _{\ell =1}^N\frac{{p_{\ell }}^2}{2m}+\frac{1}{N}\sum _{\ell =1}^{N-1}\sum _{j=\ell +1}^{N}V(r_{\ell }-r_j),$$

$$V({\theta }_{\ell }-{\theta }_j)=1-\cos ({\theta }_{\ell }-{\theta }_j),$$

$$\frac{\mathrm{d}{\theta }_{\ell }}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{p}_{\ell }}{\mathrm{d}t}$$

$$\vec{M}=(M_x,M_y)=(M\cos \phi ,M\sin \phi )=\frac{1}{N}\sum _{\ell =1}^N(\cos {\theta }_{\ell },\sin {\theta }_{\ell }),$$

$$h_1({\theta }^{\prime },p^{\prime };\vec{M},\dot{\phi })=\frac{1}{2}(p^{\prime }-\dot{\phi }{)}^2+M(1-\cos ({\theta }^{\prime }-\phi ))$$

$$p-\dot{\phi }=\pm 2\sqrt{M}\cos \frac{\theta -\phi }{2}.$$

$${\dot{M}}_x$$

$${\dot{M}}_y$$

$$\displaystyle\begin{array}[]{ll}{\mathbf{x}}\equiv(x_{1},x_{2},\dots,x_{n}),&\dot{{\mathbf{x}}}(t)={\mathbf{F}}({\mathbf{x}}(t)),\\ {\mathbf{w}}_{i}\equiv(\delta x_{i,1},\delta x_{i,2},\dots,\delta x_{i,n}),&\dot{{\mathbf{w}}}_{i}(t)\,={\mathbf{J}}({\mathbf{x}}(t)){\mathbf{w}}_{i}(t)\,,\hskip 28.45274pt1\leq i\leq n\par\end{array}$$

$$\frac{\mathrm{d}\delta {\theta }_{\ell }}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\delta p_{\ell }}{\mathrm{d}t}$$

$$\delta I_{\ell }=\frac{(\delta {\theta }_{\ell }{)}^2+(\delta p_{\ell }{)}^2}{|\mathbf{w}{|}^2}$$

Full text

Keywords : chaos, Lyapunov exponents, Hamiltonian mean field, kinetic space, molecular dynamics.

I Introduction

Lyapunov exponents are part of the fundamental characterisation of dynamical systems ott ; pikovsky . They measure the rate of divergence between nearby trajectories due to small perturbation in the initial conditions, viz.  they determine the chaoticity of dynamical systems. However, for many-body systems, they relate to the dynamics in the high-dimensional phase space of the dynamics, typically ${\mathbb{T}}^{3N}\times{\mathbb{R}}^{3N}$ for $N$ particles in a 3-dimensional cube with periodic boundary conditions.

In addition, for identical particles, it is usual to describe the system’s evolution not by a single point $({\mathbf{r}}_{1},{\mathbf{r}}_{2},\ldots,{\mathbf{r}}_{N},{\mathbf{p}}_{1},{\mathbf{p}}_{2},\ldots,{\mathbf{p}}_{N})$, but by an empirical measure $N^{-1}\sum_{j=1}^{N}\delta({\mathbf{r}}-{\mathbf{r}}_{j})\delta({\mathbf{v}}-{\mathbf{v}}_{j})$ on the 6-dimensional kinetic space ${\mathbb{T}}^{3}\times{\mathbb{R}}^{3}$, viz. Boltzmann’s $\mu$-space ehrenfest . For any finite $N$, the empirical measure and the single-point in Gibbs’ $6N$-dimensional phase space provide equivalent information on the system microscopic state, and their evolutions under the equations of motion are equivalent.

This is the cornerstone of the mean-field derivation of the Vlasov equation in the $N\to\infty$ limit for smooth inter-particle interactions spohn ; EEDvlasovia13 . Considering that the kinetic description is physically illuminating about the behaviour of particles, we are tempted to search for a particle-related description of the chaoticity of the dynamics. In this paper, we introduce a quantity characterising to what extent each particle, which moves in the kinetic space, contributes to the overall chaoticity in full phase space.

To be specific, we consider a well-known mean-field system, the cosine Hamiltonian mean field model, in which particles move simply on a circle. The one-particle configuration space is the unit circle ${\mathbb{T}}$, the $N$-particle configuration space is ${\mathbb{T}}^{N}$, the one-particle kinetic space is ${\mathbb{T}}\times{\mathbb{R}}$, and the phase space is $({\mathbb{T}}\times{\mathbb{R}})^{N}$.

In Sec. II, we provide a brief reminder on Lyapunov exponents and vectors for continuous-time systems. Sec. III introduces the cosine Hamiltonian mean field model. Sec. IV presents the numerical method and our proposal to investigate sensitive regions in the kinetic space. Sec. V is devoted to our results, and Sec. VI to our conclusion.

II Lyapunov Exponents

Lyapunov exponents are an essential tool for discussing dynamical systems. They provide the distinction between regular and chaotic behaviour by measuring how small perturbations in the system evolve in phase space abraham ; alligood ; bousko ; dorfman ; eckmann-ruelle ; mane . Let

[TABLE]

be the coordinates in the $n$-dimensional phase space (for monatomic gas or liquid models one has $n=6N$ and ${\mathbf{x}}=({\mathbf{r}}_{i},{\mathbf{p}}_{i}),i=1,\ldots,N$ ), in which the system obeys autonomous first-order differential equations

[TABLE]

generating a flow $\Phi({\mathbf{x}}_{0};t)={\mathbf{x}}(t)$ from initial data ${\mathbf{x}}_{0}$ in this space. The vector field ${\mathbf{F}}({\mathbf{x}}(t))$ is the velocity field of the flow. To measure contraction or stretching by ${\mathbf{F}}$ (over short times) and $\Phi$ (over long times) in the neighbourhood of the trajectory ${\mathbf{x}}(t)$ in phase space, consider the difference vector of two trajectories in this space, namely the deviation vector

[TABLE]

In the infinitesimal regime, the evolution equations for this deviation vector are the linear, first variation equations

[TABLE]

with ${\mathbf{J}}:=\partial{\mathbf{F}}/\partial{\mathbf{x}}$ being the $n\times n$ Jacobian matrix of the vector field. These equations also generate the evolution of the Jacobian matrix of the flow ${\mathcal{J}}({\mathbf{x}}_{0};t)=\left.\frac{\partial\Phi({\mathbf{x}}_{0};t)}{\partial{\mathbf{x}}_{0}}\right|_{t}$,

[TABLE]

The evolution equations (4) are integrated with an initial condition ${\mathbf{w}}_{0}=\delta{\mathbf{x}}(0)$. When the elements of ${\mathbf{J}}$ are continuous bounded functions of $t$, the solutions of (4) grow no faster than $\exp(\gamma t)$, for some finite $\gamma$. The Lyapunov exponents are defined by

[TABLE]

when the limit exists. A priori, this limit depends on both the initial data ${\mathbf{x}}_{0}$ and the initial deviation ${\mathbf{w}}_{0}$ since ${\mathbf{J}}({\mathbf{x}}(t))$ depends on the trajectory and ${\mathbf{w}}_{0}$ is the initial data to solve (4). Thus, in an $n$-dimensional system, one has $n$ Lyapunov exponents and each of them refers to the divergence degree of the characteristic directions of the system. All these exponents, with their multiplicity, form the Lyapunov spectrum, which is usually ordered as

[TABLE]

The Largest Lyapunov Exponent (LLE) is denoted by $\lambda_{1}$. The existence of divergent trajectories and therefore chaotic regime may be uncovered by positive values of the exponents.

When the trajectory of interest ${\mathbf{x}}(t)$ visits (densely) an extended domain ${\mathcal{D}}$ in phase space, the Lyapunov spectrum does not depend specifically on ${\mathbf{x}}_{0}$ and simply characterises this domain, which is foliated by manifolds associated with the characteristic directions singled out by the long-time behaviour of the ${\mathbf{w}}_{i}$’s. The asymptotic ${\mathbf{w}}_{i}$’s (up to their normalisation) are the associated Lyapunov vectors posch .

Hamiltonian dynamics are described in terms of pairs of conjugate variables and the form of Hamilton’s equations ensures to the Lyapunov spectrum a special symmetry property, i. e., the Lyapunov exponents come in pairs pikovsky

[TABLE]

In the following sections, we introduce the model, the cosine Hamiltonian mean field model, and the numerical approach to calculate the Lyapunov exponents. Along with these, we also discuss the numerical instruments used to investigate the sensitivity of specific degrees of freedom to the characteristic exponential behaviour associated with the largest Lyapunov exponent. This is built into the definition of the characteristic vectors, which may have arbitrary directions in phase space : our diagnostic attempts at identifying those regions of the kinetic space where particles contribute most to chaoticity. Particles in those regions are also those which will be most important when one attempts to reduce the many-body dynamics to a smaller number of degrees of freedom antoniazzi .

III The cosine Hamiltonian mean field model

We now briefly introduce the $N$-body model used in this paper. The model has a Hamiltonian form

[TABLE]

where $m$, $p_{\ell}$ and $r_{\ell}$ are the mass, momentum and position of particle $\ell$, respectively. The $1/N$ factor is Kac’ factor which makes the total energy an extensive quantity kac-meanfield and ensures that the dynamics of individual particles is well defined in the large $N$ limit under mild conditions spohn ; trocheris ; EEDvlasovia13 ; kie14 .

The cosine Hamiltonian mean field (cosHMF) is a toy model consisting of $N$ classical particles moving on a unit circle and interacting via long range force AntoniRuffo . It is a widely studied system in the literature as it is solvable at equilibrium and allows for fast molecular dynamics simulations as they scale with the number $N$ of particles instead of $N^{2}$ for more generic systems. Moreover, analytical results for the scaling behaviour of the largest Lyapunov exponent with $N$ were obtained by Firpo firpo1998 and recently approached numerically by miranda2018 .

The particle masses $m$ and coupling constant of the interaction can be eliminated by a simple rescaling, so that $p$, $\theta$, $H$, time and Lyapunov exponents are dimensionless. The interaction energy is given by

[TABLE]

with the angles $\theta_{\ell}$ representing the coordinates $r_{\ell}$ of the particles. In this paper, we only consider the attractive case. The equations of motion generated by (9),

[TABLE]

are such that each particle feels the effect of all other ones only through the value of the magnetisation (mean field)

[TABLE]

which is the first Fourier component of the particles spatial density.

This implies that each particle behaves like a charged particle in the field of a wave with amplitude $M=|\vec{M}|$ and phase $\varphi={\mathrm{Arg}}\ \vec{M}$. When the magnetisation is almost constant, particles are almost independent, though they are a priori coupled through their individual contributions to $\vec{M}$ via (13). Then, it makes sense to analyse their motion in the kinetic-theoretical $(\theta,p)$ space where each particle may appear either trapped in the potential well defined by $\vec{M}$, or circulating with respect to this well, viz. in terms of a reduced one-particle effective Hamiltonian

[TABLE]

which generates (11)-(12) for $(\theta_{\ell},p_{\ell})=(\theta^{\prime},p^{\prime})$ after a Galileo transformation to the wave frame when $\dot{\varphi}$ is constant.

As the wave moves at its own velocity $\dot{\varphi}$, it is sensible to understand this trapping in the reference frame of its well EEbook . The border between trapping and circulating behaviour occurs at $h_{1}=e_{\mathrm{s}}$ with the separating energy $e_{\mathrm{s}}=2M$ in the wave frame, and it is drawn in the one-particle $(\theta,p)$ space by the two branches of the separatrix, with equation

[TABLE]

It is well-known that the pendulum has regular dynamics, but trajectories on either side of the separatrix diverge exponentially in time (hence the name separatrix) EEbook .

Note that the dynamics (11)-(12) preserves total energy $H=\sum_{\ell}p_{\ell}^{2}/2+(1-M^{2})N/2$ and total momentum $P=\sum_{\ell}p_{\ell}=\sum_{\ell}\dot{\theta}_{\ell}$. It does not conserve the magnetisation, as

[TABLE]

Actually, for finite values of $N$, the magnetisation cannot be exactly conserved by the dynamics unless it vanishes exactly elskensESAIM .

Note also that the total energy $H$ bears no simple relation with the reduced one-particle energies for individual particles, for two reasons. First, $h_{1}$ is more conveniently expressed in the magnetisation comoving frame. Second, even when $\dot{\varphi}=0$, the potential term in $h_{1}$ sums essentially to twice the total potential energy in $H$, while the kinetic term in $h_{1}$ sums to the total kinetic energy in $H$. In other words, the average of $h_{1}$ over all particles is not the average energy per particle $e:=H/N$.

The equilibrium statistical mechanics of system (9)-(10) can be solved analytically in both canonical and microcanonical ensembles AntoniRuffo . Its caloric curve shows a second order phase transition at the critical energy per particle $e^{\ast}=E^{\ast}/N=3/4$. In terms of the order parameter $|\vec{M}|$, the critical energy separates two types of solutions, those with $|\vec{M}|\neq 0$ ($e^{\ast}<3/4$) and $|\vec{M}|=0$ ($e^{\ast}>3/4$). We refer to these regimes as subcritical and supercritical, respectively.

IV Numerical method

By definition, the calculation of the Lyapunov exponents requires the knowledge of how the flow $\Phi({\mathbf{x}}_{0};t)$, generated by (2), and its local deviations $\{\delta x_{i}\}$, evolving by (4), behave. Starting from a set of orthogonal deviation vectors, they are periodically re-orthonormalised at the same time that the expansion of volumes of different dimensions is computed.

This normalisation rescaling prevents the divergence of the deviation vectors $\delta{\mathbf{x}}(t)$, and the orthogonalisation ensures the calculation of rate of divergence along linearly independent directions in space. The standard numerical approach is available in the references benettin1980 ; shimada1979 ; wolf1985

The perturbations evolve in the tangent space, obeying the linearised equations of motion benettin_1976 ; parker . The appropriate equations are

[TABLE]

where one may set $\|{\mathbf{w}}_{i}(0)\|=1$ for convenience in the Lyapunov exponents calculation. The Jacobian (evaluated along the reference trajectory) and ${\mathbf{w}}_{i}(t)$ refer to the linearised versions of the actual equations of the system. Taking an initial condition ${\mathbf{x}}_{0}$ and an orthonormal basis $(\boldsymbol{\delta}^{(0)}_{1},\boldsymbol{\delta}^{(0)}_{2},\dots,\boldsymbol{\delta}^{(0)}_{n})$ for the linear equations, we carry out the numerical integration of equations (20). In this process, one has a single nonlinear equation and a set of $n$ linearised equations. Every new point ${\mathbf{x}}$ in phase space provides a different Jacobian matrix, so that each time step of integration involves a new linear operator ${\mathbf{J}}$. Since $\|{\mathbf{w}}_{i}(t)\|$ diverges exponentially (with time), we apply at every period $T$ the Gram-Schmidt orthogonalisation procedure to the basis $({\mathbf{w}}_{1}(t),...{\mathbf{w}}_{n}(t))$. This process is iterated $K$ times to estimate the Lyapunov spectrum and vectors.

The steps for obtaining the exponents are thus :

1. i)

for the nonlinear system, choose an initial condition ${\mathbf{x}}_{0}$ ; for the $n$ linearised equations, define an orthonormal set of initial conditions $(\boldsymbol{\delta}^{(0)}_{1},\boldsymbol{\delta}^{(0)}_{2},\dots,\boldsymbol{\delta}^{(0)}_{n})$ ; 2. ii)

the whole set of $1+n$ differential equations is integrated simultaneously for a lapse $T$: ${\mathbf{x}}_{0}\mapsto{\mathbf{x}}(T)$ and $(\boldsymbol{\delta}^{(0)}_{1},\boldsymbol{\delta}^{(0)}_{2},\dots,\boldsymbol{\delta}^{(0)}_{n})\mapsto(\boldsymbol{\delta}^{(1)}_{1},\boldsymbol{\delta}^{(1)}_{2},\dots,\boldsymbol{\delta}^{(1)}_{n})$ ; 3. iii)

intervene on the solutions of linearised equations with the Gram-Schmidt procedure ; the orthogonalisation $(\boldsymbol{\delta}^{(1)}_{1},\boldsymbol{\delta}^{(1)}_{2},\dots,\boldsymbol{\delta}^{(1)}_{n})\mapsto(\boldsymbol{v}^{(1)}_{1},\boldsymbol{v}^{(1)}_{2},\dots,\boldsymbol{v}^{(1)}_{n})$ uses a triangular matrix with unit entries on its diagonal, while the normalisation generates rescalings $(\boldsymbol{v}^{(1)}_{1},\boldsymbol{v}^{(1)}_{2},\dots,\boldsymbol{v}^{(1)}_{n})\mapsto(\boldsymbol{u}^{(1)}_{1},\boldsymbol{u}^{(1)}_{2},\dots,\boldsymbol{u}^{(1)}_{n})$ ; 4. iv)

take the updated basis $(\boldsymbol{u}^{(1)}_{1},\boldsymbol{u}^{(1)}_{2},\dots,\boldsymbol{u}^{(1)}_{n})$ and ${\mathbf{x}}(T)$ as new initial conditions, for the linearised and nonlinear equations respectively, and repeat the process $K$ times ; 5. v)

the Lyapunov exponents are estimated as time averages $\lambda_{i}=\displaystyle\frac{1}{KT}\sum_{k=1}^{K}\ln||\boldsymbol{v}^{(k)}_{i}||$. Stop iterating when $K$ is large enough to get convergence in the values of $\lambda_{i}$’s.

In particular, the Lyapunov exponents calculation for cosHMF model leads to the deviation vector ${\mathbf{w}}\equiv(\delta\theta_{1},\delta\theta_{2},\dots,\delta\theta_{N},\dots,\delta p_{1},\delta p_{2},\dots,\delta p_{N})$, associated to the linearised equations

[TABLE]

that describe the fluctuations evolution ${\mathbf{w}}$ around ${\mathbf{x}}\equiv(\theta_{1},\theta_{2},\dots,\theta_{N},p_{1},p_{2},\dots,p_{N})$.

With unit probability, an arbitrary initial perturbation ${\mathbf{w}}_{0}$ will generate a vector converging to the vector $\boldsymbol{u}_{1}$ in phase space, associated with the largest Lyapunov exponent. We attribute to each particle the weight

[TABLE]

to estimate its individual participation to the divergence rate between nearby trajectories in $2N$-dimensional phase space. This weight (normalised to unity) can be interpreted as the sensitivity level of the region where particle $k$ is in the kinetic space, which leads to the operative notion of Lyapunov modes Borzsak96 ; mila98 ; taniguchi2003B ; bosetti10 ; ginelli11 ; bosetti14 .

V Results

The evolution equations (21)-(22) are solved numerically using the centered leapfrog method with a time step $\Delta t=0.\,05$. We need a symplectic integrator to ensure that the numerical solution preserves the Hamiltonian nature of the model. Moreover, the method is time-reversal invariant, second-order and robust hairer . Figure 1 shows the fluctuations in the energy (9) and the total momentum $P$.

Numerical techniques described in previous section are applied to the cosHMF model in order to obtain the Lyapunov exponents. Figures 2 and 3 show in details the convergence of the numerical procedure and the respective initial conditions for which the equations of motion (11)-(12) are solved.

The initial condition displayed in Figure 2 corresponds to an out of equilibrium configuration called waterbag. This distribution is defined as $f(\theta,p)=1/(\Delta\theta\Delta p)$ if $0<\theta<\Delta\theta$ and $|p|<\Delta p/2$, and $f(\theta,p)=0$ otherwise. This random uniform distribution over the rectangle of area $\Delta\theta\times\Delta p$ provides a state in which the magnetisation is $M_{0}=[(1-\cos{\Delta\theta})^{2}+(\sin{\Delta\theta})^{2}]^{1/2}/\Delta\theta$ and the energy is $e=\Delta p^{2}/24+(1-M_{0}^{2})/2$ rocha-amato-1 .

Figure 3 is obtained from equilibrium initial data. In this case, the initial magnetisation and momentum are determined according to solutions of equilibrium statistical mechanics and generated randomly by Monte Carlo method davidp .

As mentioned previously, the growth of the solutions to (4) is bounded, and, for the same energy and number of particles, the Lyapunov exponents should not depend on a particular trajectory around which the perturbation is defined. To this end, one runs $10$ realisations for both types of initial conditions and show the convergence of $\lambda_{1}$. From here on, our results are generated from the equilibrium initial data as in Figure 3.

Figure 4 displays a full Lyapunov spectrum for the cosHMF model with $N=50$ particles (that means 100 exponent values). We check also whether our numerical method preserves the symplectic nature of the dynamics by plotting the sum $\lambda_{\imath}+\lambda_{2N-\imath+1}$ which must vanish according to (8).

Figure 5 displays a snapshot of the kinetic space for $N=50$ particles in the subcritical region ($e\approx 0.\,55$), using a colour scale and size scale depending on the value of our indicator $\delta I_{\ell}$. The continuous black line is the instantaneous separatrix, with two branches in the typical form of a cat’s eye given by (15). These two branches meet at the instantaneous “X point” associated with the unstable equilibrium of the one-particle Hamiltonian (14).

The colours provide a relative scale of $\delta I_{\ell}$, in which the lowest value is indicated by red and the highest by blue. The mark size is used as absolute scale, proportional to $\delta I$. This figure shows that most of the $\delta I_{\ell}$ values are practically null (less than $10^{-5}$), with few dominant weights located near the separatrix. In this equilibrium state, while particles move due to the dynamics, they essentially remain in the same region, so that the same particles come recurrently close to the cat’s eye, where chaotic motion is prevalent.

A similar pattern is also observed for other values of energy. Figure 6 displays representations of the cosHMF kinetic space for energies $e_{\rm a}\approx 0.\,25$, $e_{\rm b}\approx 0.\,55$, $e_{\rm c}\approx 0.\,70$ and $e_{\rm d}\approx 1.\,00$. The particles highlighted in blue are the particles for which $\delta I_{\ell}\geq 10^{-5}$. We clearly see the correlation of the higher values of $\delta I$ with delimitant regions of the kinetic space. In the subcritical cases (b)(c), the largest values of $\delta I$ are located around the separatrix, between bounded and free particles. Case (a) corresponds to a confined regime where the particles do not have enough energy to escape from their self-consistent potential well. Nevertheless, the largest values of $\delta I$ are situated near the separatrix. With the energy $e_{\rm d}>e^{\ast}$, the system evolved to a disordered phase in which the magnetisation vanishes in the thermodynamic limit and therefore the separatrix is not well defined ; in this case, our numerical results indicate the largest $\delta I_{\ell}$ values distributed near $p=0$, separating positive and negative velocity states.

Actually, even when the thermodynamic limit predicts the order parameter $M=0$ at equilibrium, the instantaneous value of the numerical output $\vec{M}$ for an individual realisation is of the order of $N^{-1/2}$ by the central limit theorem, what defines an instantaneous cat’s eye with a width of the order of $N^{-1/4}$. Particles with a velocity of the order of unity will not be sensitive to this cat’s eye (though they would slowly diffuse in velocity ribeiro , irrespectively of the cat’s eye), but a particle with an O($N^{-1/4}$) velocity will come close enough to the instantaneous X point, and such a particle will have the opportunity to pass the X point (if its energy is high enough) or to reverse its velocity (if its energy is too low to pass). This will generate a dichotomy on the slow particle behaviour, making it most sensitive to chaos.

Moreover, when this particle returns close to the X point (after a time of order O($N^{1/4}\ln N$) as the period of the pendulum diverges near the separatrix), the X point and the separatrix will have moved, because $\vec{M}$ itself fluctuates. Hence the particle may have crossed the fluctuating separatrix between its two passages near the X point, making the new (passing / reversal) dichotomy somewhat decorrelated from the previous one. The velocity of slow particles then approaches a random pattern, similar to the one observed in the slow relaxation of a quasi-stationary state to thermal equilibrium for the wave-particle model closely analogous to the cosHMF firpo-poleni , and in the Chirikov-Taylor standard map misguich98 .

Finally, we discuss the distribution of the weights $\delta I_{\ell}$. A rough understanding of this distribution in a quasi-stationary regime with a mildly varying separatrix may start from the idea that the one-particle dynamics (21)-(22) in tangent space $(\delta\theta_{\ell},\delta p_{\ell})$ results from multiplying this vector with Jacobian matrices generated by the corresponding particle trajectory $(\theta_{\ell}(t),p_{\ell}(t))$ and by the successive magnetisations $\vec{M}(t)$. To the extent that the dynamics exhibits chaos, one may expect the successive stretchings to rescale the vector $(\delta\theta_{\ell},\delta p_{\ell})$ more or less independently, so that $\ln\|(\delta\theta_{\ell}(t),\delta p_{\ell}(t))\|/\|(\delta\theta_{\ell}(0),\delta p_{\ell}(0))\|$ may obey a central limit theorem. This would suggest that the observed weights $\delta I_{\ell}$ could obey a log-normal distribution. Although the argument is oversimplistic, it is compatible with the cumulative distributions observed on Figure 7.

We also use the weight $\delta I$ to identify the particles carrying most of the contribution to the Lyapunov vectors associated with other exponents. In Figure 8, for the same energy as in figures 5, 6b and 7, we highlight the vectors associated with index $kN/4$ ($1\leq k\leq 4$), so that panel 8d corresponds to a null Lyapunov exponent. In agreement with intuition, the smaller exponents give more weight to particles deeper in the potential well.

The distribution of weights was also considered for the first Lyapunov vector by Ginelli et al. ginelli11 , who found a scaling $\sqrt{\delta I_{\ell}}\sim\ell^{-1}$. As their work focused on the largest Lyapunov exponent and its vector, they do not consider higher exponents and associated vectors. The localisation of Lyapunov vectors is also discussed by Bosetti and Posch bosetti14 and Taniguchi and Morris taniguchi2003B , using a global measure of localisation like $N^{-1}\exp(-\sum_{\ell}\delta I_{\ell}\ln\delta I_{\ell})$. Here we consider the detailed distribution of the weights for several Lyapunov vectors. If these weights had a power-law distribution, say $\delta I_{\ell}\approx C\ell^{-k}$ (where $k=2$ would correspond to Figure 6 of Ref. ginelli11 ), then the cumulative distribution function of the weights would scale like $(\delta I)^{1-1/k}$.

Figure 10 presents the values of $\sqrt{\delta I_{\ell}}$ in function of its index $\ell$. As noted in Ref. ginelli11 for the Lyapunov vectors associated with $\lambda_{1}$, we have $\sqrt{\delta I_{\ell}}\sim\ell^{-1}$. The observed distributions for other Lyapunov vectors of the spectrum, Figure 10(b), do not completely rule out the power-law scaling, but do not fully support it either, as well as the cumulative distributions displayed in Figure 11.

VI Concluding remarks

We associate with each particle in an $N$-body system a weight $\delta I_{\ell}$ measuring their respective contribution to the chaoticity of the dynamics, in terms of the components of the first Lyapunov vector, in a macroscopically stationary regime. While our numerical determination of $\delta I_{\ell}$ derives from the Gram-Schmidt procedure for evolving the tangent dynamics (4), this weight refers only to the largest exponent, and therefore corresponds to the covariant Lyapunov vector as well as to the “orthogonal” Lyapunov vector posch . It is thus genuinely intrinsic to the dynamics.

We applied this procedure to molecular dynamics simulations of the paradigmatic cosHMF model. We checked our calculations accuracy by computing the whole Lyapunov spectrum. According to the weight $\delta I$, the most chaotic regions in the kinetic space turn out to be, for the cosHMF model at subcritical values of energy ($e<e^{\star}$), the vicinity of the separatrix, in agreement with previous results ginelli11 . High energy ($e>e^{\star}$) implies a gaslike regime for the particle motion, in which the largest $\delta I$ values are distributed around the particles for which $p=0$ : this domain is also the one where instantaneous one-particle dynamics shows a small, fluctuating separatrix ettoumi-firpo ; ribeiro .

We also discuss the distribution of weights among particles for Lyapunov vectors associated with higher-index exponents.

Thanks to its simplicity and its robustness, our approach may help in identifying chaotic features of more complex interactions antoniES and in providing a way to interpret the behaviour of these systems with a reduced number of degrees of freedom.

Acknowledgements

LHMF thanks the PIIM, CNRS-AMU for their hospitality, Bruno V. Ribeiro and Daniel Santos for productive discussions, the remote access computing laboratory of the DF-UFRPE and CIF-UnB for efficient numerical help, and CAPES for financial support. TMRF was partially financed by CNPq (Brazil). YE enjoyed the hospitality and support from CIFMC/UnB while completing this work. It is a pleasure to gratefully acknowledge fruitful discussions with Dominique F. Escande and members of the teams “turbulence plasma” and “dynamique des systèmes complexes” in Marseille.