Chaos in Complex Astrophysical Systems

Ana Nunes¹ & Nuno Pereira²

¹Departamento de Física, Faculdade de Ciências da Unversidade de Lisboa,

Campo Grande, Edifício C1, 1749-016 Lisboa, Portugal.

² Área Departamental de Matemática, Escola Superior de Tecnologia e Gestão,

Instituto Politécnico de Beja,

Rua Afonso III, nº 1, 7800-050 Beja, Portugal.

Phone: +351 284311540, Fax: +351 284311542

Email: nuno.pereira@estig.ipbeja.pt

Abstract

The dynamics of some astrophysical systems, such as star clusters (N  10²-10⁶) or planetary systems (N  10), can be modelled by a set of differential equations known as the N-body gravitational problem. The impossibility of solving analytically the general problem when N  3 [1] requires the use of numerical integration when studying these systems.

The sensitivity of the trajectories of stars to small changes (perturbations) on the initial conditions was first reported by Miller [2]. The improved computational power now available allowed a systematic study of the typical time scale of this instability in terms of the characteristic crossing time of the system [3].

In this work we present results obtained from numerical integration of the N-body gravitational problem and the associated variational equations of motion, using the NNEWTON package [4]. We use the “Lyapunov Characteristic Indicator” [5] to estimate the time scale of the instability associated with the exponential growth of perturbations (variations) of the initial conditions.

Our preliminary results are in good agreement with those of [6] and show a simple relation between the time scale of the instability, the number of particles and the crossing time of the system. 

The N-Body Problem

Definition

 N point masses with known position and velocities for

time t_o.

 Interaction between masses according to the Newton Law of Gravity.

 The dynamics of each particle is determined by the interaction with the N-1 particles of the system.

Problem: what are the positions and velocities of the N particles at given instant of time t  t_o?



The N-Body Problem

Application to the Study of

Real Astrophysical Systems

 Celestial Mechanics (N  10): the study of planetary systems and satellite trajectories.

 Stellar Dynamics (N  10² – 10⁶): the study of star clusters.

 Galactic Dynamics (N  10¹⁰ – 10¹²): formation and evolution of structures with a large number of stars.

 Cosmology (N  “#{Universe}”): large scale structure of the Universe; clusters of galaxies.



The N-Body Problem

Mathematical Formulation

 System of 3N 2^nd order Differential Equations.

 System of 6N 1^st order Differential Equations.

 Initial Conditions: 6N coordinates – 3D vectors of position and of velocity for each particle.

 12 Integrals of Motion: the order of the system is reduced to 6N-12.

Problem: The impossibility of solving analytically the general problem when N  3 [1].



The N-Body Problem

Numerical Resolution

Direct Methods: the force in one particle is calculated taking into account the contribution of all the N-1 particles of the system.

 Advantages: spatial resolution, exact force calculation.

 Problems: Computational effort proportional to N².

Software: The NNEWTON Package [4]

 N-Body Integrators: NNEWTON Programs.

 Initial Conditions Generator: NN-VIRIAL (systems in virial equilibrium).

 Analysis Tools:

SOLEXACT2: exact solution of the 2-body problem.
NN-ELT: energy and angular momentum calculations.



The N-Body Problem

Numerical Resolution

NNEWTON Package/ Parallel Algorithm

The use of supercomputers (or clusters of PC’s) and parallel algorithms allows us study larger systems and to follow their evolution for longer (physical) times.

 Data obtained on the CRATY-T3D located at the Edingburgh Parallel Computer Center (TRACS Program, 1997). The program is a parallel version of the NNEWTON integrator [7].

 The N² dependence of the computational time is well established in this example. This is the “signature” of the direct methods.



The Exponential Instability

Definition

 Reference system s and perturbed system S=s+

 The trajectories of the systems on the 6N-dim phase space diverge exponencially.

 First observed by Miller (1964) [1] on systems with

N  32.

Problem: What is the dependence between the time scale t_e, the number of particles N, and the crossing time of the system, t_cr(=linear dimension of the system/v_rms)?



The Exponential Instability

Mesuring the evolution of perturbations:

Variational Equations

 First order variational equations of motion for particle i [4,5]:

 Variations of each particle i defined by

 Average variation (for each integration step):



The Exponential Instability

Lyapunov Characteristic Indicator

 We consider a finite dynamical system with state vector q, and the solution of the corresponding variational equation, q. An exponential growth is caracterized by

 The parameter  - Lyapunov Characteristic Exponent - is defined by

 As t   the N-body problem is asymptotically integrable and so =0. We therefore define a“Lyapunov Characteristic Indicator” as in [5], by

where t is the simulated time.

 The time scale of the exponential instability is defined by



Case Study

Simulation of N-Body Systems N=4/8/16

 Initial Conditions:

unitary masses
particles distributed in a box in the configurational space and in the velocity space.
virial equilibrium (2T-U=0)
total energy E = -1
perturbation in the x coordinate of particle 1: x₁=10^-6.
time of simultation t=100

 Programs: NN-VIRIAL, NNEWTON.



Case Study

Data Analysis (I)

 The numerical data produced by the NNEWTON program is processed by the PROVAR program [4].

 For each time iteration de values of r_a and v_a are calculated and plotted as a function of the time of simulation.

 The time evolution of r_ashows an exponential growth.

 The observed spikes are associated with the occurrence of close encounters between particles of the system.



Case Study

Data Analysis (II)

 Although the initial perturbation was introduced on a spatial coordinate of only one particle, this perturbation will “infect” the spatial and the velocity coordinates of every particle.

 The characteristic exponential growth of the perturbations is independent of the “metric”. In the above figure we plot the time evolution the quantity v_a.

 The observed spikes are associated with the occurrence of close encounters in the system (in our case, binary encounters).



Case Study

Data Analysis (III)

 Time scale estimation: linear regression on the numerical data.

 Crossing time calculation: t_cr= R/v_rms, (R is the linear dimension of the system, and v_rms the root mean squared velocity).

 Assuming that we have

 After adjusting the data we obtain



Conclusions (I)

The exponential instability and its implications to numerical simulations

“The exponential separation of trajectories was found as one of those cases in which something unexpected happens in a calculation, and a discovery was made as a result of tracking it down. (…) I foolishly thought that a good way to check a new gravitational N-body program would be to test it for microscopic reversibility. That produced a disaster”.

Miller (1993)[8]

“It may be that numerical simulations of adequate accuracy give consistent results only because they are all equally inaccurate (by the standards which would be necessary to compute reliable positions and velocities for individual stars.”

Heggie (1991)[6]



Conclusions (II)

 The positions and velocities of individual particles are quite unreliable after a few crossing times: this is a period much shorter than that over which we would like to study the dynamics of real N-body systems, such as star clusters.

 Statistical approach to the results of N-body simulations: for instance the determination of the average rate at which stars escape.

Work in Progress

 Simulation of larger systems with the use of parallel algorithms and parallel computers.

 Statistics of collisions.



References

[1] Poincaré, H. (1892): Les Methodes Nouvelles de la Mecanique Celeste. Gauthier-Villars, Paris, Vol. (I)

[2] Miller, R. H. (1964): Irreversibility in Small Stellar Dynamics Systems. Astrophysical Journal, 140, 250.

[3] Kandrup, H. E. (1994): Chaos, Regularity, and Noise in Self-Gravitating Systems. An invited plenary talk at: The Seventh Marcel Grossmann Meeting.

[4] Pereira, N. S. A. (2001): Master’s Thesis. Faculty of Sciences, University of Lisbon.

[5] Heggie, D. C. (1991): Chaos in the N-Body Problem of Stellar Dynamics. Predictability, Stability, and Chaos in N-Body Dynamical Systems, Edited by A. E. Roy, Plenum Press, New York.

[6] Goodman, J., Heggie, D. C., Hut, P. (1993): On the Exponential Instability of N-Body Systems. The Astrophysical Journal, 415:715-733.

[7] Pereira, N. S. A. (1999): A Parallel N-Body Integrator Using MPI. In: Vector and Parallel Processing – VECPAR’98, pp 626-635, José M. L. M. Palma, Jack Dongarra, Vincent Hernández (Eds.), Lecture Notes in Computer Science 1573, Springer-Verlag.

[8] Miller, R. H. (1993): Core Motions and Global Chaotic Oscillations. In: Ergodic Concepts in Stellar Dynamics, pp 137-150, V. G. Gurzadyan, D. Pfenniger (Eds.), Lecture Notes in Physics 430, Springer-Verlag.



Simulation of N-Body Systems N=4/8/16

Data Analysis (I)

Data Analysis (II)

 Although the initial perturbation was introduced on a spatial coordinate of only one particle, this perturbation will “infect” the spatial and the velocity coordinates of every particle.

 The characteristic exponential growth of the perturbations is independent of the “metric”. In the above figure we plot the time evolution the quantity va.

Data Analysis (III)

Work in Progress

References

 The characteristic exponential growth of the perturbations is independent of the “metric”. In the above figure we plot the time evolution the quantity v_a.