$N$-body chaos, phase-space transport and relaxation in numerical simulations

TL;DR

This paper investigates how chaos in self-gravitating systems depends on the number of particles, revealing a scale-dependent chaos measure that lies between collisional and violent relaxation times, indicating a collective chaos timescale.

Contribution

It demonstrates that the largest Lyapunov exponent decreases with system size, uncovering a new collective timescale in N-body gravitational dynamics.

Findings

Lyapunov exponent decreases as N increases

Lyapunov time is between collisional and violent relaxation times

Identifies a collective chaos timescale in many-body systems

Abstract

Using direct $N$-body simulations of self-gravitating systems we study the dependence of dynamical chaos on the system size $N$. We find that the $N$-body chaos quantified in terms of the largest Lyapunov exponent $\Lambda_{\rm max}$ decreases with $N$. The values of its inverse (the so-called Lyapunov time $t_\lambda$) are found to be smaller than the two-body collisional relaxation time but larger than the typical violent relaxation time, thus suggesting the existence of another collective time scale connected to many-body chaos.

Equations

$${\partial }_{t}f+\mathbf{v}\cdot {\nabla }_{\mathbf{r}}f+\nabla \Phi \cdot {\nabla }_{\mathbf{v}}f=0,$$

$$\rho (r)=\frac{3}{4\pi }\frac{M{r}_c^2}{(r_c^2+r^2{)}^{5/2}},$$

$${\Lambda }_{\max }(t)=\frac{1}{k\Delta t}\sum _k\ln \left|\frac{W(k\Delta t)}{d_0}\right|.$$

$$ϵ=({ϵ}_x{ϵ}_y{ϵ}_z{)}^{1/3},{ϵ}_i=\sqrt{⟨r_i^2⟩⟨v_i^2⟩-⟨r_i{v}_i{⟩}^2},$$

Full text

††volume: 351††journal: Star Clusters: From the Milky Way to the Early Universe††editors: A. Bragaglia, M.B. Davies, A. Sills & E. Vesperini, eds.

1 Introduction

The dynamics of $N$-body gravitational systems, due to the long-range nature of the $1/r^{2}$ force, is dominated by mean-field effects rather than by inter-particle collisions for large $N$. Due to the extremely large number of particles it is often natural to adopt a description in the continuum ($N\to\infty$, $m\to 0$) collisionless limit in terms of the Collisionless Boltzmann Equation (CBE, see e.g. [Binney & Tremaine (2008)]) for the single-particle phase-space distribution function $f(\mathbf{r},\mathbf{v},t)$

[TABLE]

coupled to the Poisson equation $\Delta\Phi(\mathbf{r})=4\pi G\rho(\mathbf{r})$. In many real self-gravitating systems, such as elliptical galaxies, $N$ is large ($N\approx 10^{11}$), and the two-body relaxation time $t_{2b}\propto Nt_{*}/\log N$ associated to collisional processes [Chandrasekhar (1943)] is much larger than the age of the universe. The typical states of such systems are therefore idealized as collisionless equilibria of Eq. (1), approached via the mechanism of violent relaxation first suggested by [Lynden-Bell (1969)], on a time scale of the order of the average crossing time $t_{*}$.

The question whether the continuum limit is effectively meaningful, and how the intrinsic discrete nature (see e.g. [Kandrup (1980)]) of collisionless systems affects their relaxation to (meta-)equilibrium is in principle still open.

Using arguments of differential geometry [Gurzadyan & Savvidy (1986)], conjectured the existence of another relaxation time $t_{*}<\tau<t_{2b}$, linked to the inverse of the exponentiation rate of phase-space volumes, and scaling as $N^{1/3}$ with the system size (see also [Vesperini (1992)]). [Kandrup & Sideris (2001), Sideris & Kandrup (2002), Kandrup & Sideris (2003), Kandrup, Sideris & Bohn (2004)] studied the ensemble properties of tracer orbits in frozen $N$-body potentials (i.e., the potential generated by $N$ fixed particles distributed according to different densities $\rho(r)$). They found that even at large $N$ the discreteness effects are not negligible in models with both integrable and non-integral continuum limit counterparts. Moreover they find that the exponentiation rate of the phase-space volume of initially localized clumps of particles is compatible with the mean Lyapunov exponent associated to the particles. The latter was found to be a slightly increasing function of $N$. In parallel [Hemsendorf & Merritt (2002)] studied the $N$-body chaos in self-consistent Plummer models as a function of $N$, quantified by means of the exponentiation rate of a subset of the $6N$ phase-space volume considering only the position part and finding again a slightly increasing degree of chaoticity.

More recently, [Beraldo e Silva et al. (2019)], using information entropy arguments, suggested the existence of another relaxation scale, associated to discreteness effects in $N$-body models, that scales as $N^{-1/6}$.

Here we revisit this matter further by studying the dynamics of individual tracer particle in frozen and self-consistent equilibrium models as well as the Lyapunov exponents of the full $N$-body problem.

2 Methods

In both frozen and active models we consider the spherically symmetric Plummer density profile

[TABLE]

with total mass $M$ and core radius $r_{c}$. In order to generate the velocities for the active simulations, we use the standard rejection technique to sample the anisotropic equilibrium phase-space distribution function $f$ obtained from $\rho$ with the standard [Eddington (1916)] integral inversion. Throughout this work we assume units such that $G=M=r_{c}=1$, so that the dynamical time $t_{*}=\sqrt{r_{c}^{3}/GM}$ and the scale velocity $v_{*}=r_{c}/t_{*}$ are also equal to unity. Individual particle masses are then $m=1/N$. For the numerical simulations we use a standard fourth-order symplectic integrator with fixed time-step $\Delta t=5\times 10^{-3}$ to solve the particles’ equations of motion and the associated tangent dynamics used to evaluate the maximal Lyapunov exponent with the standard [Benettin, Galgani & Strelcyn (1976)] method, as the large $k$ limit of

[TABLE]

In the above equation $W$ is the norm of the $6N$- or 6-dimensional tangent vector for a self-consistent $N$-body simulation and an individual particle orbit, respectively. In both cases $d_{0}=W(t=0)$.

3 Results and discussion

[Di Cintio & Casetti (2019a)] have computed the largest Lyapunov exponents $\Lambda_{\rm max}$ for the full $N$-body problem and for several tracer orbits for different system sizes $N$, finding that in the first case $\Lambda_{\rm max}$ is a decreasing function of $N$ with a slope between $-1/3$ and $-1/2$ (upward triangles in Fig. 1), contrary to previous numerical estimates (see e.g. [Hemsendorf & Merritt (2002)]). In the second case, when following a single particle in the potential of the others, $\Lambda_{\rm max}$ is almost constant in both active and frozen potentials for tightly bound particles, while decreases for less bound particles with different power-law trends as function of the binding energy.

Curiously, when measuring the average largest Lyapunov exponent of a system of non-interacting particles randomly sampled from a self-consistent isotropic system and propagated in a frozen realization of the latter, $\langle\Lambda_{\rm max}\rangle$ has the same scaling with $N$ (and remarkably close values), as shown by the diamonds in Fig. 1.

Here we compute the so-called emittance, a quantity defined in the context of charged particles beams (see e.g. [Kandrup, Sideris & Bohn (2004)] and references therein) associated to the “diffusion” in phase-space volume, defined by

[TABLE]

where $\langle...\rangle$ denotes ensemble averages. We have evaluated $\epsilon$ for the full $N$-body problem, finding that the global system’s emittance is almost conserved for large $N$, while it appears to grow with time (although the model remains in virial equilibrium) for $N\lesssim 2000$, as shown in Fig. 2 (left pannels).

When we consider instead the emittance of an initially localized cluster of tracers (see Fig. 2, right pannel) in an active $N-$body system, we observe that the latter has (independently of $N$) an increasing beahaviour with time, roughly proportional to $\epsilon_{0}\exp(\langle\Lambda_{\rm max}\rangle t)$ (dashed line), where $\epsilon_{0}$ is its value at $t=0$ and $\langle\Lambda_{\rm max}\rangle$ is the mean maximal Lyapunov exponent of the cluster. We have repeated the same numerical experiment for frozen $N$=body potentials and for time dependent potentials generated by non-interacting particle in a smooth Plummer model finding the same behaviour (See again Fig. 2). Moreover, in direct $N$-body simulations, when we compare the evolution of $\epsilon$ for a subset of active particles with that of the same cluster of non-interacting particles, we find that for large $N$ $\epsilon$ grows in the same fashion for both tracers and active particles. In turn, $\epsilon$ retains generally lower values for initially localized clusters in the frozen case, as a consequence of the conservation of energy and smaller fluctuations of the angular momentum in the latter case (see [Di Cintio & Casetti (2019a)]).

All these results lead us to conclude that the discreteness effects and the intrinsic chaoticity of the $N$-body problem may be linked to another effective relaxation scale shorter than $t_{2b}$. Moreover, it also appears that the conclusions arising from the study of frozen $N$-body models can not be easily extended to real self-consistent systems. The results on the effect of the discreteness chaos on collective instabilities such as the radial-orbit instability (ROI) will be published elsewere [[Di Cintio & Casetti (2019b)]].