Computation of Lyapunov Exponents in General Relativity

Computation of Lyapunov Exponents in General Relativity

arXiv:gr-qc/0302118v2 21 Jun 2003

Xin Wu, Tian-yi Huang
Department of Astronomy, Nanjing University, Nanjing 210093, China

Abstract Lyapunov exponents (LEs) are key indicators of chaos in dynamical systems. In general relativity the classical de?nition of LE meets di?culty because it is not coordinate invariant and spacetime coordinates lose their physical meaning as in Newtonian dynamics. We propose a new de?nition of relativistic LE and give its algorithm in any coordinate system, which represents the observed changing law of the space separation between two neighboring particles (an “observer” and a “neighbor”), and is truly coordinate invariant in a curved spacetime. Key words: chaotic dynamics, relativity and gravitation PACS: 95.10.Fh, 95.30.Sf

Chaos is a popular phenomenon in dynamical systems. One of its main features is the exponential sensitivity on small variations of initial conditions. The exhibition of chaos in the motion of Pluto makes it particularly attractive for scientists to investigate the dynamical behavior of the solar system[1]. In Newtonian mechanics, Lyapunov exponents (LEs), as a key index for meaEmail address: ( Tian-yi Huang ).

Preprint submitted to Elsevier Science

7 February 2008

suring chaos in a dynamical system, can be calculated numerically with either the variational method[2] or the two-particle method[3], and sometimes a mixture of the two[4]. The former is more rigorous but one has to derive the variational equations of the dynamical equations of the system and integrate them numerically with the dynamical equations together. The latter is less cumbersome especially when one wants to compute the maximum LE only, which is the main index for chaos. We will discuss the extension of the two particle method to relativistic models in this letter. ˙ Let q(t) and q(t) be the position and velocity vector of a dynamical system. As a set of initial conditions is selected randomly and the corresponding trajectories are restricted in a compact region in the phase space, the classical de?nition of the maximum LE is λN = lim 1 d(t) ln , t→∞ t d(0) (1)

where the separation d between two neighboring trajectories is required to be ˙ su?ciently small so that the deviation vector (?q, ?q) can be regarded as a good approximation of a tangent vector, and the distance d(t) at time t is of the form d(t) = ˙ ˙ ?q(t) · ?q(t) + ?q(t) · ?q(t). (2)

˙ It is the Euclidian distance in the phase space. q and q are in di?erent dimension, so one has to carefully choose their units to assure that both terms in the expression of d(t) are important. We suggest computing the LE in the con?guration space rather than in the phase space, that is, d′ (t) = ?q(t) · ?q(t) 2 (3)

λ′N = lim

1 d′ (t) ln ′ , t→∞ t d (0)


We argue that both λN and λ′N are e?ective in detecting chaotic behavior of ˙ orbits since ?q(t) = d/dt(?q(t)) and they should have the same Lyapunov exponents. In general relativity the above de?nition and algorithm are questionable. First, there is no uni?ed time for all the reference systems. Secondly, the separation of time and space of the 4-dimensional spacetime is di?erent for di?erent observers. Furthermore, time and space coordinates usually play book-keeping only for events and are not necessary with a physical meaning. Consequently, one would get di?erent values of λN and λ′N in di?erent coordinate systems. Up to now the references for exploring chaos in relativistic models have been mainly interested in studying the dynamical behavior of black hole systems[58] and the mixmaster cosmology[9-13]. Most of them adopt the classical definition of LE. Here we will abandon the variational method in the case of general relativity for it is rather cumbersome to derive variational equations. Actually there exist general expressions for geodesic deviation equations[5] but one has to derive the complicated curvature tensor. Furthermore, in many cases a particle does not follow a geodesic. For instance, a spinning particle in Schwarzschild spacetime doesn’t move along a geodesic[5]. Therefore, we will concentrate on constructing a revised edition of the two-particle method. The classical algorithm of LE lacks invariance in general relativity, which depends on a coordinate gauge and even could bring spurious chaos in some coordinate systems. The chaotic behavior in the mixmaster cosmology[9-13] has been debated for decade or so. Using di?erent time parameterization [14,15], one would get distinct values of λN . Especially, in a logarithmic time, chaos 3

becomes hidden[12]. Chaos, as an intrinsic nature of a given system, should not be a?ected by the choice of a coordinate gauge, and a chaos indicator, LE, should be de?ned as invariant under spacetime transformations. It means that the LE in general relativity should be expressed as a physical or so called proper quantity but not a coordinate quantity[16-19]. Now we consider a particle, called “observer”, moving along an orbit (not necessary to be a geodesic) in the spacetime with a metric ds2 = gαβ dxα dxβ . In this letter Greek letters run from 0 to 3 and Latin letters from 1 to 3. The observer can determine if his motion would be chaotic by observing whether the proper distances from his neighbors are increasing exponentially or not with his proper time. Both the distances and time are observables and should not depend on the choice of a coordinate system. Then we can apply the theory of observation in general relativity[16] to explore the dynamical behavior of the observer. Fig. 1 illustrates the trajectories of the observer and its neighboring particle called “neighbor”. The initial condition of the neighbor is arbitrarily chosen as long as its 4-dimensional distance from the observer is small enough to assure that the neighbor is approximately in the tangent space of the observer. In an arbitrary spacetime coordinate system xα , one can derive the equations of motion of the observer and its neighbor. Here we have to notice that the independent variable of the equations has to be the coordinate time t rather than their proper times τ because the two particles have di?erent proper times but one unique independent variable has to be adopted when integrating numerically their equations of motion together. At the coordinate time t the observer arrives at the point O with the coordinate xα and 4-velocity U α , and the corresponding proper time of the observer is τ . At the same coordinate 4

time t, its neighbor reaches the point P with the coordinate y α along another orbit. A displacement vector δxα = y α ? xα from O to P should be projected into the local space of the observer. The space projection operator of the observer is constructed as hαβ = g αβ + c?2 U α U β (c represents the velocity of ?→ ? light)[16]. OP ′ represents the projected vector δxα = hα δxβ , and its length ⊥ β ?→ ? OP ′ is ?L(τ ) = gαβ δxα δxβ = ⊥ ⊥ hαβ δxα δxβ . (5)

Here gαβ is calculated at xα . ?L(τ ) is the proper distance to the neighbor observed by the observer at his time τ and it is a scalar. Hence the maximum LE in general relativity is de?ned as λR = lim 1 ?L(τ ) ln . τ ?L(0) (6)

τ →∞

Let Στ be the 3-dimensional subspace of the tangent space of the observer at O, which is orthogonal to the 4-velocity U α . It is the point P ′ but not P in Στ as long as P ′ is close enough to the observer O. This tells us that λR is coordinate invariant. The next is an argument for this point. Let us carry out a time transformation t → η. For convenience, we express the projection ?→ ? ? → operation as OP ′ = H OP , where H is the space projection operator and ? → ? → H U O = 0, where a subscript is added for U to describe the 4-velocity. Assume that the observer is located at the point O at η that corresponds to the old coordinate time t, then the neighbor at η would situate at the point Q but not necessary at P . Certainly we have to keep both P and Q inside an ? neighborhood of O (see Fig.1) and ? is considered as a very small quantity. In ? → ? → ? → ? → ? → this case we have U P = U O + O(?). Furthermore, we have OQ = OP ? QP ? → ? → ? → 2 and P Q = ?τP U P + O(?τP ) = k? U P + O(?2 ), where k is a constant and 5

neighbor P P' Q' S
? L(0)



Fig. 1. The trajectories of the observer and its neighbor in spacetime.

the proper time interval of the neighbor between P and Q, ?τP , is in the ? → ? → magnitude of ?. Hence, we get the relation H OQ = H OP + O(?2). This shows that λR is invariant with time transformations. As far as the numerical implementation of the computation of λR is concerned, the following notes are worth noticing. (1) The observer and the neighbor have di?erent proper times, so a coordinate time t should be adopted as the independent variable. The equations of the motion of the particles should be transformed to use t as the independent variable. The variables to be computed step by step are the space coordinates and velocities (with respect to t) of the observer and its neighbor, and the proper time τ and dτ /dt of the observer. In total there are 14 variables. (2) As is well known, renormalization after a certain time interval is essential in this procedure to keep the distance be6

tween the observer and its neighbor small enough. The renormalization must be proceeded in the phase space though our LE is calculated in the con?guration space. (3) One important di?erence during renormalization between the relativistic and classical cases must be noticed. Let Σt be the local 3-space in which all the points have the same coordinate time t. It is evident that the renormalization should be done in Σt otherwise the computation will commit ?→ ? ? → an error. In Fig.1 P is pulled back to M and OM = (?L(0)/?L(τ ))OP . The next integration step for the neighbor will start from the point M. It is obvious ? → that OP is located in Σt because both points O and P are in Σt . On the other ?→ ? hand OP ′ is in the 3-space Στ and a pull back from P ′ to S as renormalization is not correct. Our practice has proven this argument. In order to check the validity of our scheme for calculating the relativistic LE, we are going to reexamine the core-shell system studied by Vieira & Letelier[6]. In Schwarzschild coordinates (t, r, θ, φ), the 4-metric for this system is of the form 2 2 ds2 = ?(1 ? )eA dt2 + eB?A [(1 ? )?1 dr 2 + r 2 dθ2 ] r r ?A 2 2 2 +e r sin θdφ ,


where A and B are functions of r and θ only. Here we adopt nondimensional variables and take c = 1. The metric does not explicitly depend on t and φ and has an energy constant E and an angular momentum constant L, so test particles in free fall are actually in a system with two degrees of freedom with an integral U α Uα = ?1. When A = B ≡ 0, Eq.(7) represents the Schwarzschild spacetime in which test particles in free fall move in regular orbits due to integrability of the system. For simplicity, we discuss the dynamical behavior of geodesics in a black hole plus a dipolar shell, and let 7

A = 2σ?υ, B = γ0 + 4συ ? σ 2 [?2 (1 ? υ 2 ) + υ 2 ],


where ? = r ? 1 and υ = cos θ. The Newtonian limit of the model resembles the Stark problem[20,21], which is fully integrable. However, as far as the relativistic model is concerned, Vieira & Letelier[6] and Saa & Venegeroles[7] have demonstrated strong chaos in Weyl coordinates using the Poincar? surface of section, and Gu?ron & Letelier[8] estimated its maximum e e LE λN = (3.2 ± 0.4) × 10?4 with the classical de?nition of LE (see Eq.(1)) in prolate spheroidal coordinates.

0.40 0.38 0.36

Lyapunov exponent λN (10 )


0.34 0.32 0.30 0.28 0.26 0.24 0.22 0.20 0 2 4 6



coordinate time η (10 )

Fig. 2. The maximum Newtonian Lyapunov exponent, λN , with the time variable η. It becomes one order of magnitude smaller after the time transformation (see Eq.(9)).

We choose parameters as E = 0.975, L = 3.8, σ = 2.5 × 10?4 , and γ0 = σ 2 , 8

4.4 4.2 4.0

Lyapunov exponent λR (10 )


3.8 3.6 3.4 3.2 3.0 2.8 2.6 2.4 0 2 4 6



proper time τ (10 )

Fig. 3. The maximum relativistic Lyapunov exponent, λR , with the proper time τ of the observer. It remains invariant after the time transformation (see Eq.(9)). τ reaches about 9.168 × 105 when η runs 107 .

and the initial conditions of the observer as r = 32, θ =

π , 2

φ = 0, r = 0, ˙

˙ and θ from U α Uα = ?1. As to its neighbor, an initial separation ?r = ?10?8 is adopted, regarded as the best choice[3], and the others remain the same ˙ as the observer’s except θ. The values of E and L are carefully chosen to assure that the trajectories of the observer and its neighbors are bounded in a compact region. We integrate the geodesic equations using two integrators for comparison, Runge-Kutta -Fehlberg 7(8) and the 12th-order Adams-Cowell method, with a coordinate time step 0.01. When t reaches 106 , we ?nd λN = λ′N = (2.2 ± 0.2) × 10?4, which is close to the result of [8] as its integration time amounts to 105 . Meanwhile, we get λR = (2.8 ± 0.2) × 10?4 by Eq.(6). One may notice that λN and λR are in the same magnitude. This is because 9

the Euclidian distance d′ (t) and the invariant proper distance ?L di?er not very much and so do the coordinate time t and the proper time τ . In fact, τ runs about 9.174 × 105 when t passes through 106 . A stronger gravitation by decreasing the angular momentum L will increase the di?erence between λN and λR , but the smallest stable circular orbit in Schwarzschild spacetime is r = 6. To verify the invariance of λR , we do a time transformation t → η = 10t + r 2 /2. (9)

We obtain λN = (2.6 ± 0.2) × 10?5 in the coordinate time η, while λR retains the original value (see Fig.2 and Fig.3). As a further experiment, we slightly change the model to put B/2 ≡ A = 2σ?υ. This system is still chaotic in Schwarzschild coordinates (t, r, θ, φ). To remove the singularity of this metric at the horizon, we go to Lema? cooritre dinates (T, R, θ, φ)[22] as √ √ r? 2 √ |), T = κ(t + 2 2r + 2 ln | √ r+ 2 2κr r + T, R= 3 2 √


where the positive constant κ may be chosen arbitrarily. Our numerical tests display that the larger κ is, the smaller λN becomes. Particularly, the Lema? itre coordinates hide chaos for su?cient large κ if λN is adopted as a chaotic index. However, our λR does not vary with the parameter κ. We would like to emphasize again that the computation of λR can be applied, whether the particles move along geodesics or not. The theory of observa10

tion in general relativity is not relevant to the 4-acceleration of the observer. For example, this method can be used to identify chaos in compact binary systems[23,24]. In this letter we concentrate on calculating an invariant maximum LE in general relativity, but this discussion can be easily extended to ?nd all the relativistic LEs numerically with the technique proposed by Benettin et al.[25] (also see [26]). They suggest choosing an initial set of orthonormal tangent vectors as a base of the tangent space of the observer, then compute the evolution of the volume determined by these vectors. To avoid two vectors getting close to each other under evolution they use Gram-Schmidt procedure to orthonormalize the base after each time step. The only change in a relativistic model is in computing the evolution of a tangent vector, which can be realized by the technique in this letter. As to the inner product in the orthonormalizing procedure and the norm computation, one has to use the Riemanian inner product in place of the Euclidean one. We are very grateful to Dr. Xin-lian Luo of Nanjing University for his helpful discussion. This research is supported by the Natural Science Foundation of China under contract Nos. 10233020 and 10173007.


[1] G. Sussman and J. Wisdom, Science 241 (1988) 433. [2] O. Edward, Chaos in dynamical systems, Cambridge University Press, 1993. [3] G. Tancredi and A. S?nchez, Astron. J. 121 (2001) 1171. a [4] M. Sano and Y. Sawada, Phys. Rev. Lett. 55 (1985) 1082.


[5] S. Suzuki and K.I. Maeda, Phys. Rev. D 55 (1997) 4848. [6] W.M. Vieira and P.S. Letelier, Astrophys. J. 513 (1999) 383. [7] A. Saa and R. Venegeroles, Phys. Lett. A 259 (1999) 201. [8] E. Gu?ron and P.S. Letelier, Astron. and Astrophys. 368 (2001) 716. e [9] G. Contopoulos, N. Voglis and C. Efthymiopoulos, Celest. Mech. Dyn. Astron. 73 (1999) 1. [10] D. Hobill, A. Burd and A. Coley (eds), Deterministic Chaos in General Relativity, Plenum Press, N. York, 1994. [11] L.D. Landau and E.M. Lifshitz, The Classical Theory of Fields, Pergamon Press, Oxford, 1971. [12] M. Szydlowski, Gen. Relativ. and Gravit. 29 (1997) 185. [13] M. Szydlowski and J. Szczesny, Phys. Rev. D 50 (1994) 819. [14] D. Hobill, D. Bernstein, D. Simpkmis and M. Welge, Class. Quantum Grav. 8 (1992) 1155. [15] R. Cushman and J. Sniatycki, Rep. Math. Phys. 36 (1995) 75. [16] R.K. Sachs and H. Wu, General Relativity for Mathematicians, Springer-Verlag, N. York, 1977. [17] T.-Y. Huang, C.-H. Han, Z.-H. Yi and B.-X. Xu, Astron. and Astrophys. 298 (1995) 629. [18] C.-H. Han, T.-Y. Huang and B.-X. Xu, IAUS 141 (1990) 99. [19] J.-H. Tao and T.-Y. Huang, Astron. and Astrophys. 333 (1998) 374. [20] V.I. Arnold, Mathematical Methods of Classical Mechanics, Springer-Verlag, 1989.


[21] K.P. Rauch and M. Holman, Astron. J. 117 (1999) 1087. [22] G. Lema? itre, Ann. Soc. Sci. Bruxelles I A53 (1933) 51. [23] J. Levin, Phys. Rev. Lett. 84 (2000) 3515. [24] J.D. Schnittman and F. A. Rasio, Phys. Rev. Lett. 87 (2001) 121101. [25] G. Benettin, L. Galgani, A. Giorgilli and J.-M. Strelcyn, Meccanica, March 15 (1980) 9. [26] A.J. Lichtenberg and M.A. Lieberman, Regular and Chaotic Dynamics, Springer-Verlag, 1983.




Statistics of finite-time Lyapunov exponents in a random time-dependent potential
Reformulation of QCD in the language of general relativity
Illusions of general relativity in Brans-Dicke gravity
Computation of a few Lyapunov exponents for continuous and discrete dynamical systems
A Real Polynomial Formulation of General Relativity in terms of Connections
Electromagnetic Mass Models in General Theory of Relativity
Form Invariance of Differential Equations in General Relativity
Tests of general relativity in the solar system
Collisions of rigidly rotating disks of dust in General Relativity
Measuring chaos -Lyapunov__ exponents
Geometry of chaos in the two-center problem in General Relativity
Statistical Estimation of Local Lyapunov Exponents Toward Characterizing Predictability in
Lyapunov exponents and transport in the Zhang model of Self-Organized Criticality