Critical escape velocity for a charged particle moving around a weakly magnetized Schwarzschild black hole
Abstract
We discuss charged particles motion in a spacetime of a weakly magnetized static nonrotating black hole. We study under which conditions a charged particle originally revolving around the black hole at a circular orbit after being kicked by another particle or photon can escape to infinity. We determine the escape velocity for particles at the innermost stable circular orbits and discuss the properties of particles moving with nearcritical velocity. We show that in a general case such a motion is chaotic.
pacs:
04.70.Bw, 04.25.g, 04.70.s, 97.60.Lf AlbertaThy0113I Introduction
The mechanism of jets formation in black holes is one of the most intriguing problems of modern astrophysics. The matter in the accretion disk and the black hole rotation might provide sufficient energy to form and support the jets. It is plausible that the key role in the transfer mechanism of this energy to the jets is played by magnetic fields KSKM (); DKBS (); KN (). Both widely discussed models, BlandfordZnajek mechanism BZ (); Par () and Penrose effect for magnetic fields koide:03 (); KSKM (), assume that in the vicinity of the black hole horizon there exists a sufficiently strong magnetic field. This field does not change the black hole geometry, but its interaction with charged particles and plasma is important. Such black holes are called weakly magnetized. Study of such problems in all their complexity requires 3D numerical simulations of the magnetohydrodynamics (MHD) in a strong gravitational field. (For discussion and references, see e.g. Pun ()).
Quite often when dealing with such a complicated problem it might be instructive to consider first its different simplifications, which can be treated either analytically, or by integrating ordinary differential equations. Motion of a charged particle in a weakly magnetized black hole is an important example. In this paper we discuss some aspects of this problem. Namely, we assume that a nonrotating black hole of mass is surrounded by a static axisymmetric test magnetic field, which is homogeneous at infinity, where it has the strength . In the absence of the magnetic field the geodesic equations for the particle are completely integrable. This property is also valid for the motion of the charged particle in a weakly magnetized black hole, provided such a particle moves in the equatorial plane orthogonal to the direction of the magnetic field. This case was studied in details in GP (); AG (); FS (). A similar problem for weakly charged rotating black holes was discussed in AO (). The main result of this study is that in the presence of the magnetic field the innermost stable circular orbit (ISCO) of a charged particle is located close to the black hole horizon. In this sense, the action of the magnetic field on a charged particle is similar to the action of the black hole rotation on a neutral particle. In particular, weakly magnetized black holes may play a role of particle accelerators FJ (); J (), similar to the fast rotating black holes (see e.g. Banados:2009pr (); Banados:2010kn (); Harada:2011pg (); Zaslavskii:2010pw (); Jacobson:2009zg (); Berti:2009bk () and the references therein).
In this paper we consider motion of a charged particle in a weakly magnetized black hole out of the equatorial plane. We focus on the following problem: Suppose the particle revolving around the black hole in the equatorial plane is kicked out of it by another particle or photon. Under which conditions can such a particle leave the black hole vicinity and escape to infinity? In our model the magnetic field far from the black hole is homogeneous. Hence, at a far distance from the black hole, where its gravitational field is weak and can be neglected, the charged particle moves in the homogeneous magnetic field. In this case, the corresponding equations are completely integrable. However, before the particle reaches the spatial infinity, it passes through the region where both of the fields, the gravitational and magnetic, affect its motion, so that the dynamical system looses its complete integrability and the motion of the charged particle may become chaotic.
Several examples demonstrating a similar chaotic behaviour are known. For instance, even in the absence of the magnetic field chaos arises for the motion of a particle when a spherically symmetric metric of a black hole is perturbed. For axially symmetric deformation of a black hole this was demonstrated in Maeda (). There exist also several papers discussing a chaotic particle motion in the MajumdarPapapertou metric, describing the spacetime with two or more extremely charged black holes in equilibrium Yurt (); DFC (). In the vicinity of each of these black holes the spherically symmetric ReissnerNordstrom metric is deformed by the common action of the other black holes. Numerical analysis of a charged particle trajectories in the presence of a toroidal magnetic field in the Schwarzschild spacetime was analyzed in PS (). The results presented there illustrate that under special conditions charged particles moving in the vicinity of the black hole under the action of the toroidal magnetic field can be ejected to infinity along the axis of symmetry. Even in the absence of a black hole, when a charged particle moves in a nonuniform magnetic field it often has either trapped chaotic motion BZ1 (), or chaotic scattering behavior BZ2 (); ABZ (). Chaotic motion of a charged particle in the Ernst spacetime representing a magnetized black hole was analysed in KV ().
In this paper we demonstrate that the charged particle motion near a weakly magnetized black hole is generically chaotic. We find the critical escape velocity for such a particle required to escape to infinity and discuss some properties of the nearcritical motion. The paper is organized as follows: In Sec. II we discuss our model and present an expression for the escape velocity for a neutral particle. In Sec. III we present the equations of motion of a charged particle moving around a weakly magnetized Schwarzschild black hole. Section IV discusses scattering data for our problem. Dimensionless form of the equations and the initial conditions for the particle escape from ISCO orbits are given in Sec. V. In Sec. VI we give several examples of qualitatively different orbits of a charged particle in the weakly magnetized black hole. Basinboundary analysis of the trajectories is applied for an analysis of the charged particle motion in Sec. VII. There we demonstrate the chaotic properties of the trajectories and determine the fractal dimensions in the proper domains. General discussion of the results is given in Sec. VIII. In this paper we use the sign conventions adopted in MTW () and units where .
Ii Escape velocity for a neutral particle
Before considering the escape velocity problem for a charged particle in a weakly magnetized black hole, let us remind the wellknown results for a similar problem in a simpler case when a particle is neutral and the magnetic field is absent. The background Schwarzschild metric is
(1)  
(2) 
Here is the gravitational radius of the black hole. There exist three commuting integrals of motion. Two of them are generated by the Killing vectors
(3) 
reflecting invariance with respect to time translations and rotations around the symmetry axis. The corresponding conserved quantities are the specific energy and the specific azimuthal angular momentum ,
(4)  
(5) 
Here is the mass of the particle, and are its 4velocity and 4momentum, respectively. Here and in what follows, the overdot denotes the derivative with respect to the proper time. The third integral of motion is the square of the specific total angular momentum
(6) 
Here we denoted by the following quantity:
(7) 
Using the normalization condition one obtains
(8) 
The motion of the particle is planar. Let this plane coincides with the equatorial plane. Then and the effective potential for the radial motion takes the form
(9) 
Consider a particle at the circular orbit , where is the local minimum of the effective potential . This orbit exists for . The corresponding specific energy and azimuthal angular momentum are
(10) 
The ISCO is defined by , which corresponds to a convolution point of the effective potential. For ISCO we have and .
Suppose now that the particle at a circular orbit collides with another particle or photon, so that after the collision it will move within a new plane tilted with respect to the original equatorial plane. In a general case, all three types of its motion are possible: (i) bounded motion, (ii) escape to infinity, and (iii) capture by the black hole. The result certainly depends on the details of the collision mechanism. For small values of the transferred energy and momentum the orbit will be only slightly perturbed. However, for larger values of the particle can go away from the initial plane and finally can be captured by the black hole or escape to infinity.
In a general case, as a result of the collision, the particle will have new integrals of motion: , and . For the case of a neutral particle in the Schwarzschild black hole one can easily obtain the conditions of escape in an analytical form. To be able to obtain results which allow rather simple analysis and presentation we simplify the problem and reduce the space of initial data to a oneparameter set by imposing the following restrictions: (i) the azimuthal angular momentum is not changed, and (ii) the initial radial velocity after the collision remains the same, . Under these restrictions there exists only one parameter which determines the motion of the particle, namely the new value of its energy. Under these conditions, as a result of the collision, the particle acquires a velocity in the direction orthogonal to the equatorial plane [see (7)].
After the collision, the total angular momentum and the energy of the particle are
(11)  
(12) 
During the particle motion its polar coordinate changes within the interval .
As a result of the collision, the total angular momentum of the particle changes from its original value to given by Eq. (11). The effective potential defined by the new value of is greater than that before the collision and its extrema correspond to , where
(13) 
and () is a monotonically decreasing (increasing) function of . Thus, immediately after the collision the particle is still at the turning point (), which is located between the maximum and minimum of the effective potential . Therefore, the particle will escape to infinity if , or
(14) 
In particular, for ISCO we have the escape condition , where the last equality corresponds to .
The black hole metric (1) has the evident discrete symmetries
(15) 
These symmetries imply the symmetry of the problem with respect to the following transformation:
(16) 
Iii Charged particles in a magnetized black hole
We consider now the case of weakly magnetized black holes. We assume that a particle has the electric charge and that its motion is affected by the magnetic field in the black hole exterior. Namely, we assume that there exists a magnetic field in the black hole vicinity which is static, axisymmetric and homogeneous at the spatial infinity where it has the strength . According to the procedure given in Wald (); AG (), such a magnetic field can be constructed as follows: Because the metric (1) is Ricci flat, the Killing vectors (3) obey the equation
(17) 
This equation coincides with the Maxwell equation for a 4potential in the Lorenz gauge . The special choice
(18) 
corresponds to a test magnetic field, which is homogeneous at the spatial infinity where it has the strength . The electric 4potential (18) is invariant with respect to the isometries corresponding to the Killing vectors, i.e.,
(19) 
A magnetic field is defined with respect to an observer whose 4velocity is as follows:
(20) 
where
(21) 
and
(22) 
For a local observer at rest in the Schwarzschild spacetime (1) we have , and Eqs. (18) and (20) give
(23) 
At the spatial infinity the magnetic field is directed along the vertical axis. In what follows, we assume that the field is directed upward, therefore, we shall take .
In a curved spacetime a charged particle of mass moving in an external electromagnetic field obeys the equation
(24) 
where is the covariant derivative defined with respect to the metric (1), is the particle 4velocity, , and is its charge.
Denote by the generalized 4momentum of the particle. Then the conserved quantities corresponding to the symmetries of the problem, the specific energy and azimuthal angular momentum, are defined as follows:
(25)  
(26) 
Let us denote
(27) 
Using these conserved quantities and the normalization of the 4velocity vector one gets
(28) 
The and components of Eq. (24) give, respectively
(29)  
(30) 
The normalization condition gives the following first order equation
(31)  
(32) 
This equation is a constraint. If it is satisfied at the initial moment of time, then it is always valid, provided the dynamics of and is controlled by Eqs. (29) and (30).
Let us discuss the symmetry properties of Eqs. (28)–(32). First of all, these equations are invariant under the transformations
(33) 
Therefore, without the loss of generality, we can consider the particle of the positive electric charge (and hence ). To consider a particle of a negative charge one should apply the transformation (33). In other words, a trajectory of a negative charge is related to the positive charge trajectory by the transformation , . However, after making the choice , one needs to study both cases when is positive and negative. They are physically different: The change of the sign of corresponds to the change of the direction of the Lorentz force acting on the particle.
Iv Flat spacetime limit and the scattering data
Suppose a kicked particle escapes to infinity. Let us discuss first the asymptotic properties of such escaped particle. At the asymptotic infinity the gravitational field of the black hole vanishes. Hence, such a charged particle is moving in a practically flat spacetime with homogeneous magnetic field. This type of motion is wellknown (see, e.g., Lan ()). The corresponding equations can be obtained from (28)–(32) by introducing the cylindrical coordinates
(34) 
and taking the limit while keeping finite. The asymptotic form of the equations is
(35)  
(36)  
(37) 
A solution to these equations is wellknown, it represents a helix whose axis is directed along the magnetic field , i.e. parallel to the axis. If the component of the particle velocity vanishes, its trajectory becomes a circle located in a plane. Let be the radius of the circle, then the particle velocity is given by
(38) 
Here is a vector product of the 3vectors and defined in a Euclidean space in the standard way. For the uniform magnetic field defined by (18) we have
(39) 
Therefore, the generalized 3momentum vector of the particle reads
(40) 
and the corresponding angular 3momentum vector about the point of intersection of the axis and the plane is
(41) 
It is directed along the axis, i.e. , where is a unit vector which defines positive direction of . If the center of the circle is located on the axis we have and the component of the angular 3momentum vector reads
(42) 
In our notations (27) we have
(43) 
Thus, for the solution to Eqs. (35)–(37) is
(44) 
where is a constant corresponding to . The general solution to Eqs. (35)–(37) can be obtained by coordinate transformation of (44) by moving the center of the circle to another point on the plane and boosting the solution in the direction. Accordingly, we derive
(45)  
where is a constant corresponding to and is constant velocity along the direction. As a result of this transformation, the azimuthal angular 3momentum reads
(47) 
where is the distance from the axis to the axis of the helix. Thus, one can see that for () the axis is located outside (inside) the circle, while for it passes through the circle.
The energy of the particle does not depend on the location of the circle and can be expressed as follows:
(48) 
If the particle is at rest, we have . This corresponds to and . According to Eq. (47), the last equality implies
(49) 
Here the equality sign corresponds to , i.e., the particle is located on the axis. If and , we have , and
(50) 
Returning back to our main problem of the charged particle escape, we can formulate the corresponding scattering data as the following set: .
V Dimensionless form of the equations
After these remarks we return to our problem. We shall integrate numerically the dynamical equations. For this purpose we introduce the following dimensionless quantities:
(51) 
The and components of the dynamical equation (24) together with the expression for the energy take the form
(52)  
(53)  
(54)  
(55) 
The energy of a particle revolving around the black hole in a circular orbit of radius at the equatorial plane is
(56) 
As we did before, in the case of a neutral particle, we assume that a kick does not change the particle’s azimuthal angular momentum , but only gives the particle the transverse velocity . As a result, the particle energy changes from to
(57) 
We would like to know whether the particle will escape to the asymptotic infinity. To simplify the problem, we shall consider a particle initially moving in an ISCO. In this case, the parameters and are defined by the radius of the orbit as follows (see FS ()):
Here for we have , and for we have . For we have and . For such parametrization the magnetic field, , as well as the specific angular momentum are uniquely specified by the radius of ISCO , while the only one left parameter serves to specify the kick.
Vi Types of the trajectories
Given the orbit radius and the initial energy of the particle after the kick, we integrate the dynamical equations (52)–(53) numerically. As a result of the integration, we can find a trajectory corresponding to the given initial conditions. The dynamical equations were solved using the builtin Mathematica 8.0 function NDSolve
. The integral of motion of the system [see Eqs. (54) and (57)] was used to estimate the accuracy of the numerical solver. For our calculations the energy error was found to be less than .
Results of the numerical integration show that there are three different types of the final particle motion:

The particle is captured by the black hole.

The particle escapes to .

The particle escapes to .
The outcome of the motion is considered a capture when reaches 1. It is considered an escape if reaches . The maximum computation time was chosen . In escape cases, it was found that the cumulative error can reach . The accuracy of the numerical solver can be increased to achieve much better accuracy. While increasing the accuracy is not a problem when few trajectories are to be generated, it can increase the computation time greatly when the equations of motions are intergraded hundreds of thousands of times as in the case of producing basinboundary plots (see below). However, at least in the cases we have studied, increasing the accuracy of the numerical solver does not change the final state of the motion significantly.
To distinguish the final states of the particle we introduce an integer number . It takes the value 0 for the capture case and for escape to , respectively.
If one focuses on a single trajectory, one can see that there exists a variety of its types with qualitatively different behavior of the particle within the domain close to the black hole. For each of these types of the motion the particle may pass a number of times through the equatorial plane. To characterize the dynamical behaviour of the particle in this domain we introduce the “winding” number counting how many times the particle crosses the equatorial plane before it gets captured or escapes to the spatial infinity. The integer number is, evidently, a topological invariant. For our choice of the initial conditions, the particle after the kick starts its motion from the equatorial plane in the positive direction. For this reason, it is evident, that for the winding number is even, while for it is odd. In particular, when a particle goes to without further crossing the equatorial plane, .
In rare cases of the computations the particle stays in the vicinity of the equatorial plane during the numerical lifetime . It crosses the plane many times and forms a compact cloud in the corresponding phase space. However, we expect that such a particle eventually falls into the black hole or escapes to the infinity.
Figures 1–4 illustrate possible “capture” and “escape” trajectories of the charged particle near the magnetized black hole for both at and at cases. They collect several examples of such trajectories for different values of the invariants and .
The fate of a kicked particle was found to be extremely sensitive to the initial conditions: even a very tiny change of these conditions may drastically modify its global behavior. Such an extreme sensitivity is an indication of nonintegrability of the system and its chaotic nature. The dynamical system (24) is Hamiltonian and conservative. Chaos in Hamiltonian systems was extensively studied (see, e.g., ott () and the references therein). However, most of the known tools and results, such as construction of the Poincaré sections and the KAM theorem, are related to (quasi) periodic bounded type of motion, while in our case the particle’s motion can be unbounded. For an analysis of our dynamical system we shall use the methods developed for the study of chaotic dynamical systems which posses several attractors. In our case, we have three different types of asymptotic trajectories. For such a motion a phenomenon of fractal basin boundaries takes place. We shall discuss it in the next section.
Vii Basinboundary analysis
vii.1 Basinboundary plots
The chaotic nature of a dynamical system is exhibited in its sensitive dependence of trajectories on initial conditions. A special case is, what is called, the final state sensitivity. Such a sensitivity takes place if a dynamical system has several coexisting attractors. As a result, for given initial conditions, a dynamical trajectory will typically converge to one of the attractors. Therefore, there must be boundaries of basins of the attractors separating the sets of trajectories with different final states. Such boundaries are often fractals. The larger is the fractal dimension of the basin boundaries the less predictable is the behavior of the dynamical system. To measure fractal dimension one can use the boxcounting dimension. Details and additional information can be found in, e.g. chaos (); ott (). This method, in particular, was used in the analysis of another dynamical system appearing in General Relativity (see, e.g. FrLa ()).
Here we use the basinboundary method for the analysis of the asymptotic behavior of the particle trajectories. The corresponding plots are presented in Fig. 5. As we mentioned, the particle is initially at the ISCO. Plot (a) corresponds to the case and plot (b) to the case . The horizontal axis on these plots shows the dimensional parameter of the ISCO, while the vertical axis shows the dimensionless energy of the particle after the kick. The white domains in these plots correspond to the region forbidden for the particle’s ISCO. The shadowed regions in these plots consist of square pixels of the side size and equal to . We used three different colors for these pixels. The color of a pixel determines the final outcome of the particle motion. These colors are chosen so that dark grey corresponds to the particle capture, grey color corresponds to escape to (), while light grey corresponds to escape to ().
vii.2 Critical escape energy and velocity
The general structure of the plots shown in Fig. 5 can be described as follows: The dark grey region which adjoints to the white region at the bottom of the plots corresponds to the particle capture. The uniformly grey region in the upperright part of the plots corresponds to the particle escape to . The winding number in this region is . This region is restricted from below by a diffuse domain. The uniform region is separated from the diffuse domain by the line which we call the critical escape energy line. Using our numerical results we can estimate the critical escape energy of the particle describing this line. The approximate analytical expression for the case is
(59) 
The maximal relative error of the escape energy is . Plot 6(a) illustrates the escape energy . We can obtain a similar approximate analytical relation for the escape energy for the case. It reads
(60) 
where the maximal relative error is also . This expression is illustrated in plot 6(b). The critical escape velocity as a function of can be derived using these expressions together with Eqs. (56)–(LABEL:lb). Plots (a) and (b) in Fig. 7 illustrate the critical escape velocity for and , respectively.
vii.3 Nearcritical behavior
For a given and energies close but less than the critical one the final state of the particle cannot be strictly predicted. The corresponding nearcritical domain contains the final states of all the three different types. To illustrate it more clearly let us consider a small stripe in these domains. Magnifications of the stripes are shown in Fig. 8. The horizontal axis on these plots shows the dimensional parameter of the ISCO, while the vertical axis shows the dimensionless energy of the particle after the kick which chosen close to for and close to for . These magnified plots demonstrate a linear structure of different regions corresponding to capture and both the types of the escape corresponding to . Similar magnified stripes can be constructed for different values of the magnification factor. The remarkable fact is that each of such plots has similar structure which does not depend on the value of the magnification. In other words, the nearcritical diffuse domain has fractal structure. This fractal structure is a very complicated Cantorset like structure, such that a magnification of any portion of the fractal region reveals similar pattern of the escape and capture regions on a smaller scale and it continues ad infinitum. In the fractal regions the winding number corresponding to either escape or capture can take different values and generally increases with the increasing repetition of the patterns.
vii.4 Fractal dimension of the nearcritical domains
To get a qualitative measure of the complexity of the fractal regions we calculate the boxcounting fractal dimension ,
(61) 
where is the number of squares of the sidelength needed to cover a basin boundary. Such squares are counted only if they contain at least two different colors. The boxcounting fractal dimension gives us a quantitative measure of uncertainty in our numerical computations (see, e.g., ott (); chaos ()). Namely, if our current uncertainty is, say , and we want to reduce it by a factor of by improving the precision of our computations, then the necessary precision of our computations, say , is given by
(62) 
where
(63) 
is the uncertainty exponent. Thus, we need additional digits to achieve the desired precision.
Figure 9 contains plots of vs. for different values of for the fractal structures shown in Fig. 5. The plots illustrate a linear relation for sufficiently small . The fractal dimensions of the two basinsboundaries are
(64)  
(65) 
The fractal dimension is closer to 2 for the case. Thus, the uncertainty exponent is smaller than for the case of , and the corresponding nearcritical domain has more complex fractal structure. In this case an increase in precision requires more digits in computations.
vii.5 Additional details of the basinboundary plots.
We described the main features of the structure and domains in the basinboundary plots. However, these plots contain additional structure which we briefly describe now. First of all, let us mention that for and the values of in the vicinity of and there is an escape lagoon illustrated by the lightgrey color. For the initial data corresponding to this lagoon the charged particle also escapes to infinity but in the direction opposite to the initial kick, that is with . The winding number in this region is . Besides this, there are also smaller size light grey regions which correspond to different values of the winding number. These regions form a wellvisible set of the light grey stripes located to the right of the lagoon. Similar light grey stripes corresponding to the backscattering to are present for but they are much less profound.
Viii Summary
We studied a charged particle motion in the spacetime of a weakly magnetized black hole. We demonstrated that the space of its trajectories has rather a rich structure. There exist three different types of asymptotic behavior: capture and escape to the asymptotic spatial infinity, . There also exists a class of escape trajectories when the charged particle spends a considerable time moving in the vicinity of the black hole close to its equatorial plane, crossing it again and again. For such a particle its escape to infinity has features similar to a diffusion process.
Certainly, our model is rather simplified. We made several assumptions that simplify the problem. These simplifications are of two different types. First, we chose a special type of the magnetic field. In a ‘realistic’ case the magnetic field decreases at far distances. The other simplification was the choice of orbits and parameters of the ‘kicks’. It should be emphasized that the basinboundary plots were constructed for a very special case of the kicking mechanism. For other more general types of kicks the corresponding basinboundary plots might be modified. However, we expect that the following three main features are common: Namely, the motion of the kicked particle is mainly chaotic. There is a critical escape energy line (surface) and the nearcritical domains which have fractal structure. It is interesting to confirm this by direct numerical calculations. Another interesting generalization of the problem is an analysis of the critical escape phenomenon and structure of nearcritical domains for rotating black holes in the presence of a magnetic field.
Acknowledgements.
The authors (V.F and A.S) are grateful to the Natural Sciences and Engineering Research Council of Canada for the financial support. The author (V.F.) thanks also the Killam Trust for its support.References
 J. C. McKinney and R. Narayan, Mon. Not. Roy. Astron. Soc. 375, 523 (2007).
 S. Koide, K. Shibata, T. Kudoh, and D. L. Meier, Science 295, 1688 (2002).
 P. B. Dobbie, Z. Kuncic, G. V. Bicknell, and R. Salmeron, Proceedings of IAU Symposium 259: Cosmic Magnetic Fields: From Planets, To Stars and Galaxies (Tenerife, 2008).
 R. D. Blandford, R. L. Znajek, Mon. Not. Roy. Astron. Soc. 179, 433 (1977).
 K. S. Thorne, R. H. Price, and D. A. Macdonald, Black Holes: The Membrane Paradigm (Yale University, 1986).
 S. Kide, Phys. Rev. D 67, 104010 (2003).
 B. Punsly, Black Hole Gravitohydromagnetics (SpringerVerlag, Berlin, 2001).
 V. P. Frolov and A. A. Shoom, Phys. Rev. D 82, 084034 (2010).
 D. V. Gal’tsov and V. I. Petukhov, Sov. Phys. JEPT 47(3), 419 (1978).
 A. N. Aliev and D. V. Gal’tsov, Sov. Phys. Usp. 32(1), 75 (1989).
 A. N. Aliev and N. Özdemir, Mon. Not. Roy. Astron. Soc. 336, 241 (2002).
 V. P. Frolov, Phys. Rev. D 85, 024020 (2012).
 T. Igata, T. Harada, and M. Kimura, Phys. Rev. D 85, 104028 (2012).
 M. Banados, J. Silk, S. M. West, Phys. Rev. Lett. 103, 111102 (2009).
 M. Banados, B. Hassanain, J. Silk, S. M. West, Phys. Rev. D83, 023004 (2011).
 T. Harada, M. Kimura, arXiv:1109.6722.
 O. B. Zaslavskii, Classical Quantum Gravity 28, 105010 (2011).
 T. Jacobson, T. P. Sotiriou, Phys. Rev. Lett. 104, 021101 (2010).
 E. Berti, V. Cardoso, L. Gualtieri, F. Pretorius, U. Sperhake, Phys. Rev. Lett. 103, 239001 (2009).
 Y. Sota, S. Suzuki, and K.I. Maeda, Classical Quantum Gravity 13 1241 (1996).
 U. Yurtsever, Phys. Rev. D 52, 3176 (1995).
 C. P. Dettmann, N. E. Frankel, and N. J. Cornish, Phys. Rev. D 50, R618 (1994).
 A. R. Prasana and S. Sengupta, Phys. Lett. A 193, 25 (1994).
 J. Büchner and L. M. Zelenyi, J. Geophys. Res. 94, NO. A9, 821 (1989).
 J. Büchner and L. M. Zelenyi, Geophys. Res. Lett. 17, NO. 2, 127 (1990).
 M. AshourAbdalla, J. Büchner and L. M. Zelenyi, J. Geophys. Res. 96, NO. A2, 1601 (1991).
 V. Karas and D. Vokrouhlicky, Gen. Rel. Gravit. 24, No. 7, 729 (1992).
 C. W. Misner, K. S. Thorne and J. A. Wheeler, Gravitation (W. H. Freeman and Co., San Francisco, 1973).
 R. M. Wald, Phys. Rev. D 10, 1680 (1974).
 L. D. Landau and E. M. Lifshitz, The Classical Theory of Fields (Pergamon Press, Oxford, England, 1975).
 E. Ott, Chaos in Dynamical Systems (Cambridge University Press, New York, USA, 1993).
 H. Peitgen, H. Jürgens, amd D. Saupe, Chaos and Fractals: New Frontiers of Science (SpringerVerlag New York, Inc., New York, USA, 1992).
 A. V. Frolov and A. L. Larsen, Classical Quantum Gravity 16 3717 (1999).