Chaotic motion in pp-wave spacetimes

J. Podolský^{∗}, K. Veselý

Department of Theoretical Physics, Faculty of Mathematics and Physics, Charles University, V Holešovičkách 2, 180 00 Prague 8, Czech Republic.

September 23, 1998

Abstract

We investigate geodesics in non-homogeneous vacuum pp-wave solutions and demonstrate their chaotic behavior by rigorous analytic and numerical methods. For the particular class of solutions considered, distinct “outcomes” (channels to infinity) are identified, and it is shown that the boundary between different outcomes has a fractal structure. This seems to be the first example of chaos in exact radiative spacetimes.

PACS: 04.20.Jb; 04.30.-w; 05.45.+b; 95.10.Fh

Keywords: pp-waves, chaotic motion

∗E–mail: podolsky@mbox.troja.mff.cuni.cz

1 Introduction

Over the past few decades the realization that the behavior of some nonlinear dynamical systems is extremely sensitive to initial conditions has changed our view on time evolution in physics, biology, and even economics. In the context of general relativity, which is a nonlinear dynamical theory par excellence, the first systems for which chaotic behavior of solutions to the Einstein equations has been recognized and studied were spatially homogeneous but anisotropic cosmological models. The “chaotic cosmology” programme initiated in 1968 by Misner directed attention to Bianchi IX (Mixmaster) models. The oscillatory dependence of their cosmological scale factors in different spatial directions on time has been studied both analytically and numerically in many works (see relevant contributions in [1], [2], and references therein). For a long time, however, it has been discussed how conclusively and rigorously the existence of chaos in these relativistic systems has been demonstrated as gauge invariant measures of chaotic behavior are required. Complicated nonlinear effects also occur in systems with coupled gravitational and scalar or Yang-Mills fields where the number of degrees of freedom is increased compared to purely gravitational systems (see for example [3]-[10]).

Another important type of problems providing a nonlinear dynamical system in the context of general relativity is the study of geodesic motion in a given spacetime. It is well known that Newtonian many-body systems are chaotic and it is of primary interest to investigate similar situations in Einstein’s theory. In particular, the chaotic behavior of geodesics in the relativistic analogue of the two fixed-centres problem was examined by Contopoulos [11], Dettmann et al [12] and Yurtsever [13], who investigated null and timelike geodesics in a spacetime consisting of two (fixed) extreme Reissner-Nordström black holes. These results were generalized by Cornish and Gibbons [14] to the Einstein-Maxwell-dilaton two-centre spacetimes. Bombelli and Calzetta [15] and Letelier and Vieira [16] studied chaotic geodesic motion in perturbed Schwarzschild spacetimes. The chaotic behavior of a spinning test particle in Schwarzschild spacetime was demonstrated by Suzuki and Maeda [17]. Karas and Vokrouhlický [18], Sota et al [19] and Vieira and Letelier [20], [21] studied the behavior of test particles in some static axisymmetric spacetimes. Chaotic motion has also been demonstrated in (topologically nontrivial) Robertson-Walker spacetimes by Lockhart et al [22] and Tomaschitz [23].

In our work [24] we announced that geodesic motion in the well-known class of pp-waves, which are vacuum type-N solutions representing the simplest exact gravitational wave spacetimes, is also chaotic. Here we present a detailed analysis. In the next section we investigate the geodesic equations for a non-homogeneous case. The chaotic behavior of timelike, null and spacelike geodesics in these spacetimes is then established by analytic and numerical fractal methods in sections 3 and 4, respectively. Some concluding remarks are given in section 5.

2 Geodesics in pp-waves

The metric of the widely known class of vacuum pp-waves [25], plane-fronted gravitational waves with parallel rays, can be written in the form

ds^{2} = 2 dζdζ̄ − 2 dudv − (f + f̄) du^{2} ,  (1)

where f(u, ζ) is an arbitrary function of the retarded time u and the complex coordinate ζ spanning the plane wave surfaces u = const. The non-trivial curvature tensor components are proportional to f_{,ζζ} so that (1) represents the Minkowski universe when f is linear in ζ. The case f = d(u)ζ^{2} describes the well-known plane gravitational wave (the “homogeneous” pp-wave) which has been thoroughly investigated, see [25] for Refs. This textbook example of an exact radiative solution has also been used for the construction of sandwich and impulsive waves [26]-[38].

Here we study motion in non-homogeneous vacuum pp-waves. The geodesic equations for the metric (1) are

ζ̈ + \frac{1}{2} f̄_{,ζ̄} U^{2} = 0 ,  (2)
u̇ = U = const ,  (3)
v̈ + (f_{,ζ} ζ̇+ f̄_{,ζ̄} ˙ζ̄) U + \frac{1}{2}(f + f̄)_{,u} U^{2} = 0 ,  (4)

where dot denotes d/dτ with τ being an affine parameter along the geodesic. In addition, we also assume a condition normalizing the tangent to the geodesic such that U_{μ}U^{μ} = ϵ where ϵ = −1, 0, +1 for timelike, null or spacelike geodesics, respectively. In particular, τ is a proper time for timelike geodesic observers. This condition can explicitly be written as

v̇ = \frac{1}{U} [ζ̇˙ζ̄ − \frac{1}{2}(f + f̄) U^{2} − \frac{ϵ}{2}] .  (5)

Here we assume U ≠ 0 (for U = 0 the geodesic equations can simply be integrated yielding only trivial null geodesics ζ = ζ_{0}, u = u_{0}, v = v_{1}τ +v_{0} and spacelike geodesics ζ = \frac{1}{\sqrt{2}} exp(iϕ_{0})τ +ζ_{0}, u = u_{0}, v = v_{1}τ + v_{0}, where ζ_{0}, u_{0}, v_{1}, v_{0} and ϕ_{0} are constants). By differentiating Eq. (5) with respect to τ and using (2) we immediately obtain (4) which can thus be omitted. It suffices to find solutions of (2) since v(τ) can then be obtained by integrating Eq. (5) while u = Uτ + u_{0}.

Now we concentrate on the complex equation (2) which has the same form for timelike, null and spacelike geodesics. The first integral of (2) for f independent of u is

ζ̇˙ζ̄ + U^{2}Re f = 2E ,  (6)

where E is a real constant. Using (6) we can further simplify (5) into v(τ) = U^{−1}(2E −ϵ/2) τ − 2U ∫ Re f(ζ(τ)) dτ. Moreover, it indicates that the Hamiltonian for the system of equations (2) written in two real coordinates x and y introduced by ζ = x + iy is

H = \frac{1}{2} (p_{x}^{2} + p_{y}^{2}) + V (x, y) ,  (7)

with a potential V (x, y) = \frac{1}{2}U^{2} Re f. For non-homogeneous pp-wave spacetimes given by f ∼ ζ^{n}, n = 3, 4, ⋯, the corresponding polynomial potential is called an “n-saddle”. Its shape for n = 3 and n = 5 is shown in Fig. 1. The metric (1) of these spacetimes is invariant under rotation ζ → ζ̃ = ζ exp(i2π/n) so that for any geodesic there exist other n − 1 geodesics differing only by rotations. There is also a symmetry ζ → ζ̄ corresponding to y → −y.

For simplicity, from now on we shall assume the simplest case n = 3. Since any constant multiplicative factor of V can be removed by a suitable rescaling of the parameter τ, we can without loss of generality consider a potential of the form

V (x, y) = \frac{1}{3}x^{3} − xy^{2} ,  (8)

called a “monkey saddle”. Surprisingly, the Hamiltonian (7), (8) is a special case of famous Hénon-Heiles Hamiltonian [39] which is known to be a “textbook” example of a chaotic system. Here, however, quadratic terms in V are missing. This particular case of the Hénon-Heiles Hamiltonian has rigorously been investigated by Rod [40].

3 Analytic demonstration of chaos in pp-waves

Rod described in detail the chaotic ^{1} behavior of orbits in the Hamiltonian system given by (7), (8). He concentrated on bounded orbits. These only appear in the positive energy manifolds H(x, y, p_{x}, p_{y}) = E > 0 for which the level surfaces of the potential (8) have the form of a “monkey saddle” (see Fig. 2). The homogeneity of V guarantees that the orbit structure for any two positive values of E is isomorphic modulo a constant scale factor and adjustment of time: x → x̃ = λx, y → ỹ = λy and τ → τ̃ = τ/\sqrt{λ} results in E → Ẽ = λ^{3}E. Therefore, without loss of generality we can restrict attention to any arbitrary value of E, say E = \frac{1}{3}. Geodesics for all other values of E can simply be obtained by rescaling space coordinates and the affine parameter.

There are three simple geodesics L_{j}(τ), j = 1, 2, 3, for a given E which will be described in detail in the next section. In projection to the (x, y)-plane they follow the axes y = 0 and y = ±\sqrt{3}x of the three symmetric channels of the connected region bounded by the level curves V (x, y) = E > 0 (for notation see Fig. 2). The boundary consists of three disjoint branches given by y = ±\sqrt{\frac{1}{3}(x^{2} − 1/x)} with asymptotes x = 0, x = ±\sqrt{3}y: V_{1} intersecting the x-axis, V_{2} above and V_{3} below the x-axis, respectively.

1He called it “pathological” since his paper [40] was written long before the term “chaos” came into vogue.

In order to obtain a description of the structure of bounded orbits Rod first constructed three basic periodic orbits Π_{j}(τ) in each of the three channels (Fig. 2) having perpendicular intersections with L_{j}(τ). Rod proved that these three periodic orbits are unstable and in fact are isolated invariant sets for the flow generated by the equations of motion. He also showed the existence of two additional periodic orbits Π_{4} and Π_{5} which follow the same trajectory in the (x, y)-plane, only Π_{5}(τ) = Π_{4}(−τ).

The analysis of chaotic behavior was then translated in [40] from the (x, y)-plane projection into a topological description of the asymptotic sets for the periodic orbits Π_{1}, Π_{2}, Π_{3} in the three-dimensional energy manifold H = E. The region in which these bounded orbits occur can be decomposed into three cells R_{j} such that, e.g., R̄_{1} is the compact, connected subset of the (x, y)-plane bounded by the lines D̄_{1}, D̄_{2} and Σ̄_{1} , see Fig. 2 (here the bar denotes a projection of the corresponding set from the energy manifold H(x, y, p_{x}, p_{y}) = E to the (x, y)-plane). Formally, D̄_{1} = {(x, y) ∈ L_{2} with y ≥ 0}, D̄_{2} and D̄_{3} are defined analogously using the symmetry and Σ̄_{1} denotes the minimum distance line segment between V_{2} and V_{3} (this intersects the branches V_{2} and V_{3} in points given by y = ±x); Σ_{2}, Σ_{3} are again defined analogously. Thus, Σ_{j} and the two D_{k} are boundaries of the isolating block R_{j} for the flow containing the unstable periodic solution Π_{j}. Each of the three cells R_{j} contains one such orbit and no other bounded orbits; Π_{j} is the only invariant set in R_{j}, i.e., it is isolated. Note also that any orbit that leaves the central region R_{0} = R_{1} ∪ R_{2} ∪ R_{3} never enters it again due to the structure of the acceleration field in the channels beyond Σ_{j}. Such orbits “go to infinity” through the j-channel.

Subsequently, Rod presented a description of orbits asymptotic to the basic periodic orbits Π_{j} as τ → ±∞. It is appropriate to localize these sets of asymptotic orbits by their intersection with the boundaries D_{k} which are closed topological two-discs (Fig. 3). Any point P in the disc represents a point in the phase space section whose spatial coordinates are all points (x, y) ∈ D̄_{k} while impulses are restricted by the energy condition implying p_{x}^{2} + p_{y}^{2} = r^{2} ≡ 2[E − V (x, y)]. The radius r monotonically increases from 0 to \sqrt{2E} as the points (x, y) ∈ D̄_{k} are taken between the intersection of D̄_{k} with V and the origin (0, 0). Therefore, we can identify any point P = (p_{x}, p_{y}) of the disc with an orbit passing through the point (x(r), y(r)) ∈ D̄_{k} with a given velocity (ẋ, ẏ) = (p_{x}, p_{y}). Using this we can indicate in D_{k} important points, areas and their boundaries relevant to the chaotic behavior of orbits, see Fig. 3 for D_{3} (figures for D_{1} and D_{2} can simply be obtained by symmetry). In the figure, any arc a^{+} (j) denotes the boundary of points representing orbits that go to Σ_{j} (and then to infinity) as τ increases; orbits given by a^{+} (j) approach Π_{j} as τ → +∞. Similarly a^{−}(j) describes orbits asymptotic to Π_{j} as τ → −∞. Points l_{j} denote the intersections that L_{j}(τ) has with the circle of velocities at the origin (orbits pointed towards l_{j} continue directly to the j-channel as τ → +∞).

The analysis in [40] showed that the asymptotic sets which are bordered by a^{±}(i) intersect transversally on D_{k} thus defining open sets W (m, n) called “windows” and that the intersection of the closure of any window with its bounding sets is uncountable. This gives the existence of orbits that “connect” the three unstable periodic orbits Π_{j}: they are homoclinic (asymptotic to the same periodic orbit in both time directions) or heteroclinic (asymptotic to two different periodic orbits, one in each time direction). It is the existence of bounded orbits which temporarily enter the region R_{j}, winding past Π_{j} and then emerging for recycling through the next region, and the existence of homoclinic and heteroclinic orbits that illustrates complicated chaotic structure of the flow in the central region R_{0}.

In order to describe the topology of all possible orbits in phase space it is convenient to introduce a set of sequences (finite, infinite, or biinfinite), s ≡ ⋯, s_{k}, s_{k+1} , s_{k+2} , ⋯, where s_{k} ∈ {1, 2, 3}, s_{k} ≠ s_{k+1} . Any sequence then corresponds to a sequence ⋯, R_{s_{k}}, R_{s_{k+1}} , R_{s_{k+2}} , ⋯ of blocks through which the orbit is running in the prescribed order as τ goes from −∞ to +∞. Using this symbolic dynamics, Rod defined thirteen disjoint orbit classes intersecting R_{0}:

1. bounded orbits with s = {s_{k}}_{−∞}^{−∞} which do not leave R_{0}

2. homoclinic orbits with s = s_{1}, s_{2}, ⋯, s_{1} that come from Π_{s_{1}} and go through the finite sequence to the same Π_{s_{1}}

3. heteroclinic orbits with s = s_{1}, s_{2}, ⋯, s_{f} that come from Π_{s_{1}} and go to different Π_{s_{f}}

4. orbits given by s = s_{1}, s_{2}, ⋯ that come from Π_{s_{1}}

5. orbits given by s = ⋯, s_{−2}, s_{−1} that go to Π_{s_{−1}}

6. orbits that come from Π_{s_{1}} and go through a finite sequence s = s_{1}, s_{2}, ⋯, s_{f} to Σ_{s_{f}}

7. orbits that come from Σ_{s_{b}} and go through a finite sequence s = s_{b}, ⋯, s_{−2}, s_{−1} to Π_{s_{−1}}

8. orbits that come from Π_{j} and go to Σ_{j} intersecting only the region R_{j} for j = 1, 2, 3

9. orbits that come from Σ_{j} and go to Π_{j} intersecting only the region R_{j} for j = 1, 2, 3

10. orbits that come from Σ_{s_{1}} and go through the sequence s = s_{1}, s_{2}, ⋯

11. orbits that go through the sequence s = ⋯, s_{−2}, s_{−1} to Σ_{s_{−1}}

12. unbounded orbits that come from Σ_{s_{1}} through a finite sequence s = s_{1}, s_{2}, ⋯, s_{f} to Σ_{s_{f}}

13. orbits that come from Σ_{j} to the same Σ_{j} intersecting only the region R_{j} for j = 1, 2, 3

Clearly, considering the reversibility of time, the classes 5., 7., 9. and 11. are equivalent to 4., 6., 8. and 10, respectively. Also, the union of orbit classes 1. through 5. together with the periodic orbits Π_{1}, Π_{2}, Π_{3} form the invariant set of bounded orbits. Orbits that are unbounded in only one time direction are in classes 6. through 11., and those that are unbounded in both time directions are in classes 12. and 13.

Finally, the chaotic structure of bounded orbits was analysed in [40] by (a topological version of) the Smale horseshoe map. With this technique Rod showed that to any biinfinite sequence of symbols {1, 2, 3} there exists an uncountable number of bounded orbits running through the blocks R_{j} in a given sequence. (Similarly, all the orbit classes 4. through 13. each contains an uncountable number of orbits.) Also, the flow admits at least a countable number of homoclinic and heteroclinic orbits.

Rod remarked that the results could be refined if the basic orbits Π_{j} were known to be hyperbolic so that they would admit stable and unstable asymptotic manifolds. Consequently, to each periodic symbol sequence there would correspond a countable collection of periodic orbits. The difficulties in proving the hyperbolicity of Π_{j} were subsequently overcome (in an even more general context) in [41]: the “monkey saddle” potential studied above is a special case of “Example A” of [41] given by V (x, y) = \frac{1}{3} ∏_{i=1}^{3} (x − ρ_{i}y) with ρ_{1} = 0, ρ_{2} = \sqrt{3}, ρ_{3} = −\sqrt{3}.

In [42], summarizing and generalizing some previous results [43], the Hamiltonian (7), (8) (as a particular case of the Hénon-Heiles Hamiltonian) was presented as an example of a system for which the Smale horseshoe mapping can be explicitly embedded as a subsystem into the flow along the nondegenerate homoclinic and heteroclinic orbits to hyperbolic unstable periodic orbits. The complex behavior of nearby orbits then implied the nonexistence of global second analytic integral. This completed the proof of the chaotic nature of the studied system in the sense of a rigorous definition of chaos, cf. [44].

4 Numerical demonstration of chaos in pp-waves

The equations of motion resulting from (7), (8) have a very simple form

ẍ = y^{2} − x^{2},  ÿ = 2xy .  (9)

However, their explicit analytic solutions can only be found in very special cases. Of course, there are (unstable) geodesics with E = 0 given by x = 0 = y. Less trivially, letting y = 0 for all τ, one gets radial geodesics L_{1}(τ) through the y = 0 channel (analogous geodesics L_{2}(τ) and L_{3}(τ) through remaining two channels can simply be obtained by rotation, see Fig. 2). Solutions of this type starting at τ = 0 from the level curve V_{1} = E (i.e., x(0) = ^{3} \sqrt{3E}, ẋ(0) = 0) are given by

x(τ) = ^{3} \sqrt{3E} [1 − \sqrt{3} \frac{1 − cn(βτ)}{1 + cn(βτ)}] ,  (10)

where β^{2} = 3^{−1/6} 2 ^{3} \sqrt{E} and cn is the Jacobian elliptic function with modulus k^{2} = (2 + \sqrt{3})/4, i.e., k = sin(\frac{5}{12}π). The function x(τ) monotonically decreases from x(0) = ^{3} \sqrt{3E} across x(τ_{b}) = 0 to x(τ_{s}) = −∞ (where the singularity is located). Interestingly, the proper times τ_{b} and τ_{s} are finite. In fact, one can easily calculate that βτ_{s} = 2K(k) ≈ 5.53612629 where K(k) is the complete elliptic integral of the first kind, and βτ_{b} = F (k, φ) ≈ 1.84537543 where F (k, φ) is the elliptic integral of the first kind with sinφ = 3^{−1/4}2/(1 + \sqrt{3}). Clearly, τ_{s} = 3τ_{b}.

In particular, for E = 0 the solution with y = 0 is simply

x(τ) = −\frac{6}{(τ − τ_{s})^{2}} ,  (11)

which describes (for τ ≥ τ_{s}) a geodesic emerging from the singularity x = −∞ at τ = τ_{s} and approaching the origin x = 0 asymptotically as τ → ∞, or (for τ ≤ τ_{s}) a geodesic that falls from the origin to the singularity. Again, if we set the initial condition x(0) = x_{0} < 0 then τ_{s} = \sqrt{−6/x_{0}} so that the time needed to fall from x_{0} to the singularity is finite.

We also considered to localize the trajectories of the basic periodic orbits Π_{j}. Due to the symmetry one can concentrate on Π_{1} only which can be described as a function x(y). Assuming E = \frac{1}{3}, it must be a solution of a non-linear equation

\frac{1 + 3xy^{2} − x^{3}}{1 + x^{′2}} x′′ + 3xyx′ = \frac{3}{2}(y^{2} − x^{2}) ,  (12)

where x′ = dx/dy, such that x(−y) = x(y) and x′(0) = 0. Then Π_{1} is given by

x = a + by^{2} + cy^{4} + dy^{6} + ⋯,  (13)

where b = \frac{3}{4}a^{2}/(a^{3}− 1), c = \frac{1}{16}(2 + 13a^{3})/(a^{3}− 1)^{2}, d = \frac{1}{320}a(505a^{6}− 437a^{3}− 68)/(a^{3}− 1)^{4}. We found the value of a numerically and then calculated the remaining constants using the above relations: a = −0.5152, b = −0.1751, c = 0.0107, d = −0.0012.

The explicit geodesics presented above are very special. In order to get the global picture of a motion one has to perform the integration of Eqs. (9) numerically. Typical geodesics in the studied spacetime are shown in Figs. 4 and 5.

In Fig. 4 we present geodesics starting from the branch V_{2}. The curves are orthogonal to V_{2}, proceed first downwards and “fan out” from both sides of Π_{1} (and also Π_{3}). Such behavior illustrates some of the analytic results presented in [40] and indicates that all Π_{j} are unstable.

Other geodesics that we obtained by numerical integrations are shown in Fig. 5. Their initial conditions are chosen such that the geodesics start at τ = 0 from a circle x^{2} + y^{2} = ρ^{2} in the (x, y)-plane. Due to a scaling property of the “monkey saddle” potential (see section 3) we can, without loss of generality, assume ρ = 1 (all other geodesics except those intersecting the origin can simply be obtained by rescaling). In Fig. 5a we present geodesics with ẋ = 0 = ẏ at τ = 0 and in Fig. 5b geodesics starting with non-vanishing (but same) velocities. (Note that these conditions are different from the approach adopted in the previous section since the geodesics are not on the same “energy manifolds” E = const. However, E is not the energy of particles and photons and there is no physical reason to sort all the geodesics according to E no matter how useful it proved to be from the mathematical point of view).

We observe that each unbounded geodesic escapes to infinity (i.e., falls to the singularity) through only one of the three channels in the potential. Choosing non-zero initial velocities makes more geodesics prefer one of the channels (cf. Fig. 5a and Fig. 5b) but does not significantly change the character of the motion.

In fact, all unbounded geodesics through the three channels oscillate around the corresponding “basic” radial geodesics L_{j}(τ) discussed above. Let us assume geodesics through the first channel centered by L_{1}(τ) given by (10) (similar results for the second and the third channel follow by symmetry). As they approach the singularity at x = −∞, we may assume |y| ≪ |x| (this is also justified by our numerical simulations). Then the asymptotic solution to Eq. (9) is x(τ) ≈ −6(τ_{s} − τ)^{−2} and therefore

y(τ) ≈ \sqrt{τ_{s} − τ} (A cos[\frac{\sqrt{47}}{2} ln(τ_{s} − τ)] + B sin[\frac{\sqrt{47}}{2} ln(τ_{s} − τ)]) ,  (14)

where A and B are arbitrary constants. As the geodesics approach the singularity, their frequency of oscillations around L_{1}(τ) grows to infinity while the amplitude of oscillations tend to zero. We call this effect a “focusation”.

The main objective of this section, however, is to establish the chaotic behavior of geodesics in pp-waves by numerical means. This may be more illustrative than the formal analysis presented in the previous section. Chaos is usually indicated by a sensitive dependence of possible outcomes on the choice of initial conditions. The standard approach, called a fractal method, was advanced in the papers by Cornish, Dettmann, Frankel, Levin and others (see for example [2], [11]-[14]). It starts with a definition of several different outcomes, i.e., types of ends of all possible trajectories. Subsequently, a set of initial conditions is evolved numerically until one of the outcome states is reached. Chaos is established if the basin boundaries which separate initial conditions leading to different outcomes are fractal. As we shall now demonstrate, we observe exactly these structures in the studied system.

It is natural to parametrize the unit circle from which the geodesics in Fig 5a start at τ = 0 (with vanishing velocities) by x(0) = cosϕ, y(0) = sinϕ, ϕ ∈ [−π, π). We have already pointed out that all unbounded geodesics approach the singularity at infinite values of x and y through only three distinct channels. These represent possible outcomes of our system and we assign them symbol j which takes one of the corresponding values, j ∈ {1, 2, 3} (thus, for example, j = 1 means that the geodesic approaches the infinity at x = −∞, y = 0 through the first channel centered by L_{1}(τ) as τ → τ_{s} > 0). From Fig 5a we observe that the behavior of the function j(ϕ) depends very sensitively on ϕ in certain regions. We calculated j(ϕ) numerically and we display the results in Fig 6 ^{2}. Also, we plot in the same diagram the function τ_{s}(ϕ) which takes the value of the parameter τ when the singularity is reached by a given geodesic.

2When the values of j were colour coded we obtained nice fractal pictures. Unfortunately, here we could present only their black-and-white versions which are not so impressive. Therefore, we chose simply to plot the function j(ϕ).

The boundaries between the outcomes appear to be fractal. This is confirmed on the enlarged detail of the image and the enlarged detail of the detail etc. up to the sixth level. In Fig. 7 we show such zooming in of the fractal interval localized around the value ϕ ≈ 0 (there are two similar fractal intervals around ϕ ≈ \frac{2}{3}π and ϕ ≈ −\frac{2}{3}π corresponding to the other outcome channels). On each level the structure has the same property, namely that between two larger connected sets of geodesics with outcome channels j_{1} and j_{2} ≠ j_{1} there is always a smaller connected set of geodesics with outcome channel j_{3} such that j_{3} ≠ j_{1} and j_{3} ≠ j_{2}. Similarly as in [11], [14], the structure of initial conditions of these three types of orbits resembles three mixed Cantor sets, and this fact is again a manifestation of chaos.

The above structure of j(ϕ) has its counterpart in the fractal structure of τ_{s}(ϕ), see Fig. 6 and Fig. 7. We observe that the value of this function increases considerably on each discontinuity of j(ϕ), i.e., on any fractal basin boundary between the different outcomes (in fact, τ_{s} is infinite there). Thus, there is an infinite number of peaks, each of which corresponds to an orbit asymptotic to one of the periodic orbits Π_{j} as τ_{s} → ∞ (these orbits never “decide” on a particular outcome, and so never escape to infinity). Also, the value of τ_{s} increases in the non-chaotic regions of ϕ as one zooms into the higher levels of the fractal. This can be understood physically. Most geodesics fall into the singularity directly. In Fig. 6 they form the three largest connected sets of length Δϕ_{0} ≈ 1.83 on which j(ϕ) is constant. Three smaller connected sets of length Δϕ_{1} ≈ 0.248 in the first level of the fractal structure correspond to geodesics that approach the singularity after one “bounce” on the potential wall V_{j} in the central region. Their values of τ_{s} must naturally be larger. (Note that they also contain geodesics L_{j}(τ) given by ϕ = 0, ϕ = \frac{2}{3}π and ϕ = −\frac{2}{3}π for j = 1, 2 and 3, respectively; we calculated explicitly their times τ_{s} = 2K(k) ≈ 5.5361/β ≈ 5.1519 which agrees with the values in Fig. 6). In the second level of the fractal there are sets of length Δϕ_{2} ≈ 8.6 × 10^{−3} with geodesics that reach the singularity after two “bounces” having even larger values of τ_{s} etc. We found numerically that Δϕ_{3} ≈ 3.23 × 10^{−4}, Δϕ_{4} ≈ 1.18 × 10^{−5}, Δϕ_{5} ≈ 4.28 × 10^{−7}, Δϕ_{6} ≈ 1.56 × 10^{−8}, etc.

Using this we can estimate the fractal dimension D of the above structure. We observe that Δϕ_{i}/Δϕ_{i+1} → r = 27.4 ± 0.2 for large values of i. Since the “fractal pattern” is doubled in each step we get D = ln2/ lnr = 0.209 ± 0.001 so that the dimension is clearly non-integer.

The motion can also be visualized by the time evolution of a ring of free test particles in the (x, y)-plane which are in rest initially. In Fig. 8 we observe that the circle is deformed in a complicated way. In fact, the circle forms loops that escape to different outcome channels as τ grows. They approach the singularity and one can easily check that as the particles forming the original circle move in different channels, their relative proper distance grows. Of course, there are also isolated particles of the initial ring that will forever remain near the origin approaching asymptotically the basic periodic orbits Π_{j}.

Finally, in order to confirm once more the chaotic behavior of geodesics in pp-waves we return back to the Rod analysis. We checked his main result concerning the motion in the energy manifold E = const. by numerical simulation of the disc D_{3}. This is shown in Fig. 9 which significantly improves Rod’s schematic sketch (cf. Fig. 3). There are two narrow chaotic bands with a fractal structure. Each boundary in the band between two different outcomes represents some bounded geodesic as τ → ∞. Clearly, there is an uncountable number of such geodesics. A similar set for τ → −∞ can simply be obtained by a reflection with respect to the y-axis. An interesting feature of the disc is that the fractal bands bend as they approach the boundary of the disc which represents the geodesics crossing the origin x = 0 and y = 0. In Fig. 10 we present some geodesics of this type. Again, the outcome is extremely sensitive on initial conditions. It also demonstrates that the higher the level of the fractal, the greater is the number of “bounces” that a geodesic undergoes before falling into the singularity.

5 Final remarks

In this paper we have demonstrated by invariant analytic and numerical methods the chaotic behavior of geodesic motion in non-homogeneous vacuum pp-wave spacetimes. It was established for all types of geodesics: timelike, null and spacelike. This seems to be the first explicit demonstration of chaos in exact radiative spacetimes (note that a chaotic interaction of particles with particular classes of linearized gravitational waves on given backgrounds has been studied in [15], [16], [45]-[48]). In fact, pp-wave solutions represent the simplest models of exact gravitational waves in general relativity. However, so far most work has concentrated on homogeneous pp-waves for which the motions are integrable (Eq. (2) is linear) and hence non-chaotic.

Although we investigated in detail only the spacetime (1) with the function f(ζ) of the form f ∼ ζ^{3} corresponding to the “monkey saddle” potential V given by (8), as we indicated in [24] the results might be carried out to f ∼ ζ^{n} with n ≥ 4, i.e., to the general n-saddle potential V = \frac{1}{n} Re ζ^{n} (cf. Fig. 1 for n = 5). In particular, the decomposition of the central region into topological isolating blocks and the existence of basic unstable periodic solutions Π_{j} in each channel j = 1, 2, ⋯, n is immediate [40]. These periodic solutions are hyperbolic [41] and the existence of nondegenerate homoclinic and heteroclinic orbits is established [42], [43]. As an example we present our numerical results for n = 5: Fig. 11a shows unbounded geodesics starting from a unit circle from rest and Fig. 11b displays corresponding functions j(ϕ) and τ_{s}(ϕ). There are now five outcome channels and again, the structure is fractal.

Our work demonstrates that geodesic motion in spacetimes even as simple and well-known as pp-waves can be complex. Hopefully, it will initiate investigation of chaos in other exact radiative space-times.

Acknowledgments

We acknowledge the support of grants No. GACR-202/96/0206 and No. GAUK-230/1996 from the Czech Republic and Charles University. We thank the developers of the software system FAMULUS (Tomáš Ledvinka in particular) which we used for computation and drawing of all the pictures. Also, we thank Jerry Griffiths for his help with the manuscript.

References

[1] Hobill D, Burd A and Coley A (ed) 1994 Deterministic Chaos in General Relativity (New York: Plenum Press)

[2] Cornish N J and Levin J J 1997 Phys. Rev. Lett. 78 998; Phys. Rev. D 55 7489

[3] Choptuik M W 1993 Phys. Rev. Lett. 70 9

[4] Calzetta E and El Hasi C 1993 Class. Quantum Grav. 10 1825

[5] Larsen A L 1994 Class. Quantum Grav. 11 1201

[6] Blanco S, Domenech G, El Hasi C and Rosso O A 1994 Gen. Rel. Grav. 26 1131

[7] Blanco S, Costa A and Rosso O A 1995 Gen. Rel. Grav. 27 1295

[8] Calzetta E and Gonzales J J 1995 Phys. Rev. D 51 6821

[9] Cornish N J and Levin J J 1996 Phys. Rev. D 53 3022

[10] Barrow J D and Levin J J 1998 Phys. Rev. Lett. 80 656

[11] Contopoulos G 1990 Proc. R. Soc. Lond. A 431 183; 1991 435 551

[12] Dettmann C P, Frankel N E and Cornish N J 1994 Phys. Rev. D 50 R618

[13] Yurtsever U 1995 Phys. Rev. D 52 3176

[14] Cornish N J and Gibbons G W 1997 Class. Quantum Grav. 14 1865

[15] Bombelli L and Calzetta E 1992 Class. Quantum Grav. 9 2573

[16] Letelier P S and Vieira W M 1997 Class. Quantum Grav. 14 1249

[17] Suzuki S and Maeda K 1997 Phys. Rev. D 55 4848

[18] Karas V and Vokrouhlický D 1992 Gen. Rel. Grav. 24 729

[19] Sota Y, Suzuki S and Maeda K 1996 Class. Quantum Grav. 13 1241

[20] Vieira W M and Letelier P S 1996 Phys. Rev. Lett. 76 1409

[21] Vieira W M and Letelier P S 1996 Class. Quantum Grav. 13 3115

[22] Lockhart C M, Misra B and Prigogine I 1982 Phys. Rev. D 25 921

[23] Tomaschitz R 1991 J. Math. Phys. 32 2571; 1993 34 1022; 1993 34 3131

[24] Podolský J and Veselý K 1998 Phys. Rev. D 58 081501 [25] Kramer D, Stephani H, MacCallum M A H and Herlt E 1980 Exact Solutions of Einstein’s

Field Equations (Cambridge: Cambridge University Press)

[26] Bondi H, Pirani F and Robinson I 1959 Proc. R. Soc. Lond. A 251 519

[27] Bondi H and Pirani F 1989 Proc. R. Soc. Lond. A 421 395

[28] Podolský J and Veselý K 1998 Czech. J. Phys. 48 871

[29] Aichelburg P C and Sexl R U 1971 Gen. Rel. Grav. 2 303

[30] Penrose R 1972 General Relativity ed L. O’Raifeartaigh (Clarendon Press: Oxford)

[31] Ferrari V and Pendenza P 1989 Gen. Rel. Grav. 22 1105

[32] D’Eath P D and Payne P N 1993 Phys. Rev. D 46 658

[33] Balasin H and Nachbagauer H 1996 Class. Quantum Grav. 13 731

[34] Aichelburg P C and Balasin H 1997 Class. Quantum Grav. 14 A31

[35] Griffiths J B and Podolský J 1997 Phys. Lett. A 236 8

[36] Balasin H 1997 Class. Quantum Grav. 14 455

[37] Podolský J and Veselý K 1998 Phys. Lett. A 241 145

[38] Steinbauer R 1998 J. Math. Phys. 39 2201

[39] Hénon M and Heiles C 1964 Astronom. J. 69 73.

[40] Rod D L 1973 J. Differential Equations 14 129

[41] Rod D L, Pecelli G and Churchill R C 1977 J. Differential Equations 24 329

[42] Churchill R C and Rod D L 1980 J. Differential Equations 37 23

[43] Churchill R C and Rod D L 1976 J. Differential Equations 21 39; 21 66

[44] Churchill R C 1994 Deterministic Chaos in General Relativity ed D Hobill, A Burd and A Coley (New York: Plenum Press)

[45] Varvoglis H and Papadopoulos D 1992 Astron. Astrophys. 261 664

[46] Cornish N J and Frankel N E 1997 Phys. Rev. D 56 1903

[47] Chicone C, Mashoon B and Retzloff D G 1997 Class. Quantum Grav. 14 699

[48] Kokubun F 1998 Phys. Rev. D 57 2610

Figure Captions

Fig. 1 The shape of the potential V (x, y) = \frac{1}{n} Re ζ^{n}, where ζ = x+iy, describing geodesics in corresponding non-homogeneous vacuum pp-wave spacetimes (here for n = 3 and n = 5).

Fig. 2 For n = 3, the “monkey saddle” potential V (x, y) = \frac{1}{3}x^{3} − xy^{2} plays a crucial role in famous Hénon-Heiles chaotic system. The level surface V = E consists of disjoint branches V_{j}, j = 1, 2, 3, defining three channels to infinity, each centered by a special geodesic L_{j}(τ). The central region can be decomposed into three cells R_{j} bounded by D_{k} and Σ_{j} such that the only bounded orbit in R_{j} is an unstable periodic orbit Π_{j}(τ). Also, the trajectory of two additional periodic orbits Π_{5}(τ) = Π_{4}(−τ) which we found by numerical simulations is drawn by a dashed line. See the text for more details.

Fig. 3 The boundary D_{3} between the isolating blocks R_{2} and R_{3} is a closed topological twodisc. Any point (p_{x}, p_{y}) in this schematic figure represents an orbit passing with a velocity (ẋ, ẏ) = (p_{x}, p_{y}) through the point (x(r), 0) ∈ D̄_{3} such that (for E = \frac{1}{3}) x = ^{3} \sqrt{1 − \frac{3}{2}r^{2}} where r = \sqrt{p_{x}^{2} + p_{y}^{2}} . Arc boundaries of orbits that go to infinity through the j-channel as τ → ±∞ are denoted by a^{±}(j). Their intersections form “windows” W (m, n) containing homoclinic and heteroclinic orbits asymptotic to and from the basic periodic orbits Π.

Fig. 4 Geodesics starting from the branch V_{2} “fan out” from both sides of Π_{1} and Π_{3}.

Fig. 5 Geodesics starting from a unit circle in the (x, y)-plane a) from rest, b) with initial velocity ẋ = 0.4, ẏ = 0.8. They escape to infinity (where the singularity is localized) through only three channels in the potential and a sensitive dependence of these three possible outcomes j ∈ {1, 2, 3} on the choice of initial conditions is observed.

Fig. 6 Plot of functions j(ϕ) and τ_{s}(ϕ) which labels the three possible outcomes and the value of τ when the singularity is reached by a given geodesic, respectively, on ϕ ∈ [−π, π) parametrizing the initial position on a unit circle, x(0) = cosϕ, y(0) = sinϕ; ẋ(0) = 0 = ẏ(0). Boundaries separating different outcomes are fractal establishing chaos. Each peak in τ_{s}(ϕ) which coincides with discontinuity in j(ϕ) corresponds to some bounded orbit.

Fig. 7 The fractal structure of j(ϕ) and τ_{s}(ϕ) was confirmed by zooming in the interval around the value ϕ ≈ 0 up to the sixth level. Between any two connected sets representing geodesics with outcome channels j_{1} and j_{2} ≠ j_{1} there is always a smaller connected set with j_{3} ≠ j_{1} and j_{3} ≠ j_{2}.

Fig. 8 The time evolution of a ring of free test particles in the (x, y)-plane, initially in rest. The circle is deformed in a complicated fractal way with different segments moving to different outcome channels. Notice, for example, a similarity of the patterns at τ = 4.0 and τ = 7.0. However, at τ = 7.0 the lines are in fact doubled and consist of particles coming from different parts of the original circle.

Fig. 9 Exact form of the disc D_{3} which we obtained by numerical simulation significantly improving Rod’s schematic sketch presented in Fig. 3. The chaotic bands in the dics of radius \sqrt{\frac{2}{3}} ≈ 0.8165 are very narrow and also have fractal structure which is clear from the enlarged details. The coding is such that white colour denotes the outcome channel j = 1, black corresponds to j = 2 and grey corresponds to j = 3. Each boundary in the fractal basin represents a geodesic asymptotic to some Π_{j} as τ → +∞.

Fig. 10 Typical geodesics starting from the origin x(0) = 0 = y(0) with velocities ẋ(0) = \sqrt{\frac{2}{3}} cosψ, ẏ(0) = \sqrt{\frac{2}{3}} sinψ for a) ψ ∈ (−0.28, +0.28), b) ψ ∈ (0.024, 0.028), c) ψ ∈ (0.0254, 0.0256). Again, the outcome is extremely sensitive on initial conditions.

Fig. 11 Motion in pp-waves given by the structural function f ∼ ζ^{5}: a) geodesics starting from a unit circle from rest approach infinity through five outcome channels, b) corresponding functions j(ϕ) and τ_{s}(ϕ) also have fractal structure.