{'fraction_non_alphanumeric': 0.04848585968132143, 'fraction_numerical': 0.01529990823392008, 'mean_word_length': 4.035792303814485, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
McGehee regularization of general SO(3)-invariant potentials and applications to stationary and spherically symmetric spacetimes

The McGehee regularization is a method to study the singularity at the origin of the dynamical system describing a point particle in a plane moving under the action of a power-law potential. It was used by Belbruno and Pretorius to perform a dynamical system regularization of the singularity at the center of the motion of massless test particles in the Schwarzschild spacetime. In this paper, we generalize the McGehee transformation so that we can regularize the singularity at the origin of the dynamical system describing the motion of causal geodesics (timelike or null) in any stationary and spherically symmetric spacetime of Kerr-Schild form. We first show that the geodesics for both massive and massless particles can be described globally in the Kerr-Schild spacetime as the motion of a Newtonian point particle in a suitable radial potential and study the conditions under which the central singularity can be regularized using an extension of the McGehee method. As an example, we apply these results to causal geodesics in the Schwarzschild and Reissner-Nordstr\"om spacetimes. Interestingly, the geodesic trajectories in the whole maximal extension of both spacetimes can be described by a single two-dimensional phase space with non-trivial topology. This topology arises from the presence of excluded regions in the phase space determined by the condition that the tangent vector of the geodesic be causal and future directed.


Introduction
Kerr-Schild metrics [Kerr and Schild, 1965] are a well-known Ansatz to solve the Einstein field equations and leads to many physically important exact solutions of the four-dimensional case, such as the Schwarzschild black hole, the Reissner-Nordström, the Kerr black hole, the charged Kerr-Newman black hole, the Vaidya radiating star, Kinnersley photon rocket, pp-waves and also some of their higher dimensional analogues [Málek, 2012]. Kerr-Schild metrics have played a crucial role in the discovery of rotating black holes in higher dimensions [Myers andPerry, 1986, Gibbons et al., 2005] as well as in the recent work on so-called higher order gravities [Anabalón et al., 2009, Anabalón et al., 2011. Also, most static and spherically symmetric spacetimes can be displayed in Kerr-Schild form and their analysis conforms a field that continues giving interesting results nowadays [Parry, 2012, Hackmann et al., 2008. Two of the best known static and spherically symmetric metrics which have been studied extensively and remain an area of current research are the Schwarzschild metric and the Reissner-Nordström metric. Although the behavior of the geodesics in both metrics is well-known, the geodesic equations have a large number of dynamic properties that are still providing new results, such as the characterization of the circular motion in the Reissner-Nordström spacetime for neutral and charged particles [Pugliese et al., 2011a, Pugliese et al., 2011b or the dynamics of the chaotic motion in the Schwarzschild black hole surrounded by an external halo [de Moura and Letelier, 2000]. The dynamical system approach to the analysis of the geodesic flow in these spacetimes and their rotating Kerr generalizations is a novel approach which, besides providing many new and interesting results, also describes known results from a different perspective. Examples are the homoclinic orbits that asymptotically approach the unstable branch of circular orbits [Levin and Perez-Giz, 2008, Misra and Levin, 2010 or the fact that perturbation of the geodesic flow possesses a chaotic invariant set [Moeckel, 1992, Levin, 2000, Suzuki and Maeda, 1999. One of the advantages of this method is that a great amount of information can be obtained without integrating the geodesic equations. Also, by the use of "blow-up" techniques in the dynamical system we can describe the behavior of the geodesic equations near the singularity, of which little is known. One of the most recent works along this line [Belbruno and Pretorius, 2011] has analyzed the null case of the geodesic flow by the use of the McGehee regularization [McGehee, 1981], which is a method designed to deal with the singularities at the center in the motion of Newtonian particles subject to a central power-law potential. In the context of geodesics in Schwarzschild, the limitation of the MacGehee method is the restriction to power-law potentials, which prevents its application to timelike geodesics. This is one of the reasons why null geodesics only where treated in [Belbruno and Pretorius, 2011]. Also, the standard McGehee method involves a somewhat complicated phase space which obscures the analysis. Other approaches to this problem [Stoica and Mioc, 1997] have studied timelike geodesic in Schwarzschild by the use of a variation of the McGehee method. However, the approach is such that one deals with a one-parameter family of energy-dependent dynamical systems in which only one curve in each phase space is relevant. This obviously obscures and complicates unnecessarily the results (in fact, this drawback was not explicitly noticed in [Stoica and Mioc, 1997]).
In this paper we generalize the McGehee regularization so that we can deal with central potentials of a very general form. With this method we can treat not only general causal geodesics in Schwarzschild but also geodesics in the Reissner-Nordström spacetime. In fact, for a substantial fraction of the paper we work in full generality in stationary and spherically symmetric spacetimes of Kerr-Schild form, of which the previous are just particular cases. There are several possible approaches to derive the geodesics equations in such spacetimes. Explicit computation of the Christoffel symbols is tedious and not particularly enlightening. It is particularly cumbersome to incorporate the conserved quantities associated to Killing vectors into the equations. More straightforward and convenient is the use of Hamiltonian methods which, in particular, allows for the incorporation of conserved quantities into the system in a straightforward way. Once we have the geodesic equations for such spacetimes, we can apply the generalized McGehee regularization and subsequently analyze the phase space defined by the geodesic equations, with particular emphasis at the vicinity of the singularity, where new and interesting dynamics appears. A remarkable fact is that the dynamics in the entire maximal extension of the spacetime can be described in a single two-dimensional phase space, which has particular importance in the Reissner-Nordstrom case. The key for this lies in the presence of excluded regions in the phase space arising from the condition that the trajectories correspond to future directed causal geodesics.
The paper is organized as follows: In section 2 we follow a simple way to obtain the geodesic equations for a general stationary Kerr-Schild metric and obtain a simplified Hamiltonian with the Killing conserved quantity already incorporated. In section 3 we particularize to the case of stationary and spherically symmetric Kerr-Schild spacetimes. In particular, we find that the geodesics can be described by a classical Hamiltonian of the form H = T + V with V a spherically symmetric potential. It is not at all clear a priori that this should be possible in the whole Kerr-Schild domain. The Hamilton equations already incorporate all the constants of motion associated to the symmetries. We analyze under which conditions a Hamiltonian trajectory corresponds to a causal, future directed geodesic of the spacetime. These conditions will be translated into excluded regions in the corresponding phase spaces. In Section 4 we generalize the McGehee transformation to radial potentials of very general form and provide a method to choose the appropriate parameter to perform the regularization. This discussion also helps clarifying the original regularization procedure proposed by McGehee. The physical meaning of the generalized McGehee variables is also discussed. In Section 5 we particularize the previous general results to the Schwarzschild spacetime paying particular attention to the collision manifold and to the excluded region for future-oriented geodesics. We recover the known results on null geodesics near the singularity obtained in [Belbruno and Pretorius, 2011] and extend them to timelike geodesics (in fact, all causal geodesics are treated simultaneously). As already mentioned, the understanding of the excluded regions is crucial to have a phase space of physical trajectories with a non-trivial topology capable of dealing with all Kerr-Schild patches of the Kruskal spacetime. Finally, in Section 6 we perform a similar analysis for the maximal extension of the Reissner-Nordström spacetime.
2 Geodesic equations for a general stationary Kerr-Schild metric Throughout this paper, we will consider spacetimes {M = R × (R 3 \ C), g} where C ⊂ R 3 is a closed subset such that M is connected and g is a Lorentzian metric of Kerr-Schild form [Kerr and Schild, 1965]. More specifically, let {x α } = {T, x i } (α, β, · · · = 0, 1, 2, 4 and i, j, · · · = 1, 2, 3) be Cartesian coordinates on R × R 3 and endow M with the Minkowski metric η = −dT 2 + δ ij dx i dx j . Let K be a smooth one-form on M which is null with respect to the metric η and h : M −→ R a smooth function. The metric g being of Kerr-Schild form means that it takes the form It is well-known (and immediate to check) that the inverse metric g −1 is where all Greek indices are raised and lowered with the Minkowski metric η. This expression shows, in particular, that the one-form K is also null in the metric g. We will assume from now on that neither K nor h vanish on a non-empty open set on M.
Our aim in this section is to study the geodesic equations for a Kerr-Schild metric assuming the spacetime to be stationary with Killing vector ξ = ∂ T . It is clear from (1) that ξ is a Killing vector of g if and only where £ denotes Lie derivative. At any point p ∈ M where K| p = 0, let V p ∈ T p M be a vector subspace such that T p M = K| p ⊕ V p and use this direct sum to decompose £ ξ K| p = C| p K + U | p . Inserting this into (2) yields which is equivalent to hU | p = 0 and (£ ξ h + 2Ch)| p = 0. Using the fact that neither h nor K vanish on a non-trivial open set, it follows that ξ = ∂ T is a Killing vector of g if and only if there exists a smooth function C : M −→ R such that £ ξ K = CK and £ ξ h = −2Ch. Let f 0 : R 3 \ C −→ R be any smooth positive function and let f : while the metric g takes the form Dropping the primes, it follows that ξ is a Killing vector for g if and only h and K can be selected to be Lie constant along ξ. We assume this from now on. In the Minkowskian coordinates {x α } let us write K α = ( K, K) where K satisfies K 2 = K 2 := K i K i and Latin indices are raised and lowered with the Euclidean metric δ ij . The Killing vector ξ is timelike on the set {p ∈ M; h K 2 | p < 1}, null on the set {p ∈ M; h K 2 | p = 1} and spacelike on the set {p ∈ M; h K 2 | p > 1}. Note also that we are not assuming K to be future directed or past directed everywhere, so that a priori K may change sign.
In any spacetime (M, g), affinely parametrized geodesics are the solutions of the Hamilton equations of the Hamiltonian defined on the cotangent bundle of M. The Hamilton equations fix p = g(u, ·) where u is the tangent vector to the geodesic. Using the explicit expression (1) for the metric, this Hamiltonian takes the form Given that ξ is a Killing vector, the quantity E := −p(ξ) is conserved along geodesics. Note also that, with this definition, where we have written p = { p, p } and dot means scalar product with δ ij . The Hamiltonian itself is a conserved quantity with the value of H = − 1 2 µ where µ = 0, ±1 depending on whether the geodesic is timelike (µ = 1), spacelike (µ = −1) or null (µ = 0). Inserting (6) and the conserved quantity E into (5) the following Hamiltonian arises naturally which is now defined on the cotangent bundle of R 3 \ C.
The interest of this Hamiltonian lies in the fact (easy to check) that if a curve (T (s), x(s)} is a geodesic in (M, g) with tangent vector u satisfying g(u, u) = −µ and conserved quantity p(u) = −E, then { x(s)} is the projection to the base space R 3 \ C of a solution of the Hamilton equations of (7) satisfying along the curve and T (s) satisfies the ODE which is simply the explicit form for g(u, ξ) = −E in the Cartesian coordinates {T, x}. The converse to this statement will be addressed in the following section in the spherically symmetric case.
3 Geodesic equations for a stationary and spherically symmetric Kerr-Schild metric In this section we want to particularize the problem to the spherically symmetric setting. So, we assume the group of rotations SO(3) acting on M as to be an isometry of g. Note that, for this definition to make sense, the set C must be invariant under the SO(3) action, which we assume from now on. The isometry condition requires £ ζ (hK ⊗K) = 0, for any generator ζ of the group SO(3). Pulling back this relation to the orbits of the isometry group and using the fact that the only symmetric 2-covariant tensor on the sphere which is invariant under SO(3) is a constant times the standard metric on the sphere, it follows that hK ⊗ K pulls back to zero on the SO(3) orbits. Since h does not vanish on open sets, we conclude that K itself pulls back to zero on these surfaces. Given the stationarity condition, this is equivalent to the existence of a smooth function f : The following lemma gives the most general form of g under a mild additional restriction.
Lemma 1. Assume that h is non-zero on a dense set, that the null vector K α does not have any flat zero (i.e. a point where K α and all its derivatives vanish) and that (M, g) is stationary and spherically symmetric. Then h and K α can be chosen so that h( x) is spherically symmetric and Proof. The condition that f has no flat zeros implies that f (and hence K) cannot vanish on a non-empty open set. So, as discussed in Section 2, we can assume £ ξ h = 0 and £ ξ K = 0 where ξ = ∂ T , and that (10) holds. Let S I be the collection of arc-connected components of {f = 0} ⊂ R 3 \ C. On each one of these open sets we haveK = σ I f , where σ I = ±1 is constant on S I . Let S + be the union of components S I with σ I = +1 and S − be the union of components S I with σ I = −1 and assume that both are non-empty. Since S + ∪ S − is dense in R 3 \ C and the latter is connected it follows that there exists a point p ∈ R 3 \ C that can can be approached by a sequence {p + i ∈ S + } and by a sequence {p − i ∈ S − }. SinceK is smooth everywhere, in particular at p, it follows that necessarily f and all its derivatives vanish at p, against assumption. Thus, either S − = ∅ (and we can writeK = f everywhere) or S + = ∅ (and we can writeK = −f everywhere). Consequently, the Kerr-Schild metric takes the form g = η + hf 2 K ⊗ K with K α = (σ, x | x| ). Defining h = hf 2 , and given the spherically symmetric invariance of K α , it follows immediately that g is spherically symmetric if and only if h is spherically symmetric. Dropping the primes in K α and h the lemma follows.
Remark 2. As a consequence of this lemma, the Hamiltonian H in eq. (7) takes the form An important question is to what extent the field equations of this Hamiltonian reproduce the information concerning the geodesics of g. Note first that the Hamilton equation˙ x = ∂H ∂ p reads explicitly˙ which can be written in matrix form as (we denote by x T the transpose of the vector column x and by Id the identity matrix)˙ Now, the relationship between the four-velocity u α of a geodesic and the corresponding fourmomentum p α = g αβ u β is obviously invertible. For a geodesic {T (s), x(s)}, the four-velocity is u =Ṫ (s)∂ T +˙ x(s)∂ x . Lowering indices and using Since E = −p(∂ T ), we also have p = −EdT + p · d x or, equivalently, p The first equation is exactly equation (9) for the case under consideration and must be added to the Hamiltonian system (11) in order to describe the geodesics. Concerning the second equation, its relationship to equation (12) is as follows. First of all, it is immediate to check that any trajectory satisfying (14)-(15) also satisfies the pair of equations (14)-(12). To analyze the converse, observe that the matrix in parenthesis in (12) is invertible for all h = 1. So, given x(s), this equation can be solved uniquely to obtain p(s) and hence, assuming that (14) holds, this solution must be necessarily (15). This shows the equivalence between (14)- (15) and (14) (15) parallel to x which is, for any value of h, Note finally that, at points where h = 1 the dependence ofṪ drops completely from (14). Given a solution { x(s), p(s)} of the Hamiltonian equations of H , it is precisely (16) that allows one to solve forṪ (s) at points satisfying h = 1, and hence must be added to the system. The next lemma shows that the trajectories of the Hamiltonian (11) can be also obtained by solving a much simpler Hamiltonian.
Let { x(s), p(s)} = γ( x 0 , p 0 )(s) and {x(s),p(s)} =γ(x 0 ,p 0 )(s). Both Hamiltonians H andĤ are spherically symmetric and time independent, so there exist constants J,Ĵ, H 0 andĤ 0 such that The respective Hamilton equations implẏ and hencex We next write down explicitly H ( x(s), p(s)) − H 0 = 0. For any vector a, we can compute its square norm as From (21) we have Decomposing p(s) 2 and˙ x(s) 2 according to (22) and inserting (23), a straightforward calculation where the second equality defines V (| x|). For the trajectoryx(s), the form of the HamiltonianĤ immediately implieŝ where the second equality definesV (|x|). Comparing (24) and (25) we conclude that the two trajectories x(s) andx(s) agree if and only if they have initial position, initial velocity and the respective potential functions V (| x|) andV (| x|) agree up to an additive constant c. The condition Since by hypothesis h( x ) is not of the form h( x ) = α| x| 2 β + γ| x| 2 −1 in any neighborhood of x 0 , this equation has as only solution c = 0, J 2 = L 2 and H 0 = 1 2 (E 2 − µ). We conclude that the trajectories x(s) andx(s) agree if and only if x 0 =x 0 ,˙ x| s=0 =ẋ| s=0 , | J| = |L| and H 0 = 1 2 (E 2 − µ). Given the relation (21) between velocity and momentum, the lemma follows.
Remark 4. It is interesting that the HamiltonianĤ is independent of σ, so that we will be able to describe the geodesics in (M, g) both for the case when K is future directed (plus sign) or past directed (negative sign). Moreover, the HamiltonianĤ is a standard Hamiltonian in Newtonian mechanics for a point particle in a central potential. This a substantial simplification over the original problem of solving the geodesic equations in a stationary and spherically symmetric spacetime of Kerr-Schild form, because we can exploit all the information known for trajectories of point particles in Newtonian mechanics under the influence of a radial potential of the form The main consequence of Lemma 3 is, thus, that the spatial part of all geodesics in a stationary and spherically symmetric spacetimes of Kerr-Schild form turns out to be equivalent to the (much simpler) problem of solving the trajectory of a Newtonian point particle in the potential (26).
Once the spatial part of the geodesics is solved, the temporal part is dealt with by solving equation (14) (at points where h = 1) and equation (16) (at points where h = 1). Since we are interested in causal and future directed geodesics we need to find the restrictions on the initial data which guarantee this. The following Proposition summarizes the results above and addresses the issue of future directed initial data for both choices of σ.
Proposition 5. Let M = R × (R 3 \ C) be connected with C ⊂ R 3 closed. The most general stationary and spherically symmetric metric of Kerr-Schild form g = η + hK ⊗ K such that h and K are smooth and with no flat zeros can be written in the form in any domain. Then, the g-geodesic trajectories (T (s), x(s)) with normalized tangent vector correspond exactly to the solutions of where L is an arbitrary constant vector, µ takes the values µ = +1 for timelike geodesics, µ = 0 for null geodesics and µ = −1 for spacelike geodesics anḋ where E is a constant. Moreover, the tangent vector u(s) := (Ṫ (s),˙ x(s)) satisfies In addition, if the time orientation of (M, g) is chosen so that the null vector Proof. The first part of the Proposition is a consequence of Lemma 3 in combination with Remark 3. Note, in particular, that (16) (at points where h( x) = 1) must be used to reconstruct the spacetime trajectory (T (s), x(s)) from the solutions of equations (28)-(29).
Remark 6. Note that when L = µ = 0 we have a 0 ≡ 0 and this Proposition admits the initial dataṙ 0 = 0, E = 0 irrespectively of the value of h 0 . When h 0 = 1, this boundary case corresponds to the situation when the initial tangent four-vector vanishes, and hence the geodesic is a trivial curve. This is consistent with the fact that the zero vector is null and future directed. Admitting trivial curves as null future directed geodesics has the advantage that allows one to treat at once the cases µ = 0 and µ = 1.
Corollary 7. The variation ranges for are independently of the sign of σ and of the function h( x) in the Kerr-Schild metric.
Proof. Immediate from the ranges of variation of E in Proposition 5 and the relation = 1 2 (E 2 −µ).
4 Blow-up of the singularity for radial potentials In his original paper [McGehee, 1981], McGehee proposes a method of blowing-up the singularity by introducing a coordinate system that regularizes the origin for power-law radial potentials where dot is derivative with respect to τ and grad = ∂ ∂x i is the gradient operator. Since the trajectories lie in a plane, this system can be restricted to R 2 without loss of generality. Introducing, as usual, an auxiliary vector variable y this system can be rewritten as a first order system on R 4 as˙ At this point McGehee proposes identifying R 2 with the complex plane C. Writing {x, y} for { x, y} after this identification, the change of coordinates where χ = 2 2+σ has the two properties of (i) regularizing the system at r = 0 and (ii) decoupling the system in the pair of variables {u, v}. Thus the original system transforms into an autonomous dynamical system on the plane {u, v} (with no singularities) together with a pair of first order ODEs in {r, θ} (also free of singularities) which can be integrated afterwards. The new variables take values u, v ∈ R, r ∈ [0, ∞) and θ ∈ [0, 2π).
It is natural to ask whether such a regularization and decoupling procedure also occurs for more general potentials V (| x|). We prove in appendix A that only the power-law and the logarithmic potentials V (| x|) ∝ ln | x| decouple in the variables {u, v}, even after introducing arbitrary functions of r in the transformation (39). Despite this impossibility, the system can still be simplified substantially by a suitable choice of generalized McGehee transformation.
Theorem 8. Let N be an open annulus in C and V : N → R be a radially symmetric function on N × C is equivalent to the system where β is an arbitrary constant. The coordinates {r, θ, u, v} take values in r ∈ (a, b) ⊂ R + , θ ∈ S 1 and u, v ∈ R. The coordinate change is defined by where τ is the flow parameter in (40) and s is the flow parameter in (41).
Remark 9. In the case when the potential V (|x|) is a power-law V (|x|) = |x| −σ the transformation (42) does not agree with the original McGehee transformation (39) even after making the choice β = − σ 2 . The reason lies in the specific choice χ = 2/(2 + σ) made by McGehee. In fact, any choice of non-zero constant χ in the transformation (39) preserves the same properties for the transformed system. We prefer to make the choice χ = 1 because then r measures directly the distance of the particle to the origin. In order to recover the specific form used by McGehee, it is necessary to apply an additional coordinate change r → r χ to the system (41). However, this has no benefits for the dynamical system and has the drawback of obscuring the clear geometric interpretation of r.
Proof. The coordinate change (42) is a particular case of the coordinate change (67) introduced in appendix A with ξ 1 (r) = r, ξ 2 (r) = r β and ξ 3 (r) = r 1−β . In particular, equations (70) and (71) hold with α = 1 and c = 1. Thus, the dynamical system in the new coordinates takes the form (72) which is exactly (41) after setting α = 1, c = 1. N being an open annulus, it must be of the form N = {x ∈ C; x = re iθ with r ∈ (a, b), θ ∈ S 1 } for some 0 < a < b. This proves the claim on the domain of the new coordinates.
Concerning the properties of the new dynamical system we have Theorem 10. With the same conditions as in Theorem 8, the transformed dynamical system admits the following two constants of motion Moreover, if 0 ∈ N and we assume that there is γ > 0 such that |x| γ V (|x|) admits a C 1 extension (as a function of |x|) to |x| = 0, then for any β ≤ − γ 2 the system (41) admits a C 0 extension to Proof. The dynamical system (40) describes the motion of a point particle under the influence of a radial potential V . Thus, the angular momentum L = x ×˙ x and the energy = 1 2˙ x 2 + V ( x) are conserved. In terms of the complex variables {x, y} they take the form [McGehee, 1981]: where is the imaginary part andx is the complex conjugate of x. Applying the coordinate change (42) one finds as claimed. Assume now that ∃γ > 0 such that f (|x|) := |x| γ V (|x|) can be C 1 extended to |x| = 0. The function Λ(|x|) (see expression (40)) is defined to be Inserting this into (41) we see that a sufficient condition for the dynamical system to admit a C 0 extension to r = 0 is that 2β + γ ≤ 0.
The existence of the first integral L can be used to remove v from the equations and reduce the dimensionality of the system, as well as to decouple a two-dimensional subsystem.
Lemma 11. The subset of trajectories of Theorem 8 with constant of motion L are equivalent to the dynamical system defined on (a, b) × S 1 × R. This system is decoupled in the {r, u} variables and admits the first integral Proof. Solve for v in the constant of motion (43) and substitute in the dynamical system (41) and in the expression for (44).

Interpretation of the coordinates
The physical meaning of the coordinates {r, θ, u, v} follows easily from their definition (42): 1. The coordinates r, θ are the standard polar coordinates on the plane.
2. The coordinates u, v are proportional to the radial and the angular components of the velocity. This follows from the first equation in (40) because Thus, u carries all the radial information of the velocity whereas v encodes the angular part of the velocity, The decoupling of the system in the {r, u} variables is adequate since it corresponds to the usual decoupling of the radial motion of a point particle under the influence of a radial potential. Once this motion is solved, the angular motion θ(s) follows by simple integration of θ (s) = Lr(s) −(β+1) .

The regularized reduced dynamical system
We have already discussed in Theorem 10 the range of values for β which regularize the dynamical system (41) at r = 0. The reduced dynamical system (46)-(47) incorporates extra powers of r which are potentially divergent at r = 0. The following lemma determines the range of β which regularizes the reduced system.
Corollary 12. Under the same assumptions as in Theorem 10, the reduced system (46)-(47) for Proof. We already know that β ≤ − γ 2 regularizes the term in Λ of equation (47) at r = 0. For L = 0, equation (48) admits a C 0 extension at r = 0 if and only if 1 + β ≤ 0. This condition also regularizes the first term in equation (47). Thus, β ≤ min{−1, − γ 2 } is a sufficient condition for the existence of a continuous extension to r = 0.
Remark 13. The optimal choice of β for a detailed study of the dynamical system (46)-(47) at r = 0 is β = min{−1, − γ 2 } with γ selected in such a way that |x| γ V (|x|) admits a C 1 extension to |x| = 0 and lim |x|→0 |x| γ V (|x|) = 0. Indeed, a larger value of β is not capable of regularizing the system at r = 0. On the other hand, a smaller value of β overkills the singularity. This has the effect that the invariant submanifold {r = 0} (which is called the collision manifold) has u = 0 as the unique fixed point, and this is always non-hyperbolic. Thus, all details of the phase space structure of the dynamical system at {r = 0} are lost by this choice of β. We will see below an example of this behavior when considering the Schwarzschild limit of the dynamical system describing causal geodesics in the Reissner-Nordström spacetime. As is well-known, the Kruskal spacetime of mass M > 0 outside its bifurcation surface can be covered by four patches, two of them isometric to the advanced Eddington-Finkelstein spacetime and the remaining two to the retarded Eddignton-Finkelstein. These spacetimes [Eddington, 1924, Finkelstein, 1958 consist of the manifold R × R + × S 2 with respective metrics where σ = −1 for the advanced case and σ = 1 for the retarded one. The coordinates take values in V ∈ R, r ∈ R + and dΩ 2 is the round unit metric on the sphere. Furthermore the time orientation is chosen so that V increases along any timelike curve. The coordinate change V = T − σr transforms the metric (51) into ds 2 = −dT 2 + dr 2 + r 2 dΩ 2 + 2M r (dT − σdr) 2 with range of variation T ∈ R, r ∈ R + . Transforming the flat metric dr 2 + r 2 dΩ 2 to Cartesian coordinates brings the metric to Kerr-Schild form where r = √ x · x and the manifold is R × (R 3 \ {0}). The choice of time orientation in (51) implies that the null vector field ∂ T + σ x | x| ∂ x is future directed. This form of the metric was obtained in [Kerr and Schild, 1965]. The case σ = −1 covers the upper and right quadrants of the Penrose diagram of the Kruskal spacetime depicted in Fig. 1 and hence approaches the black hole singularity at r = 0. The case σ = +1 covers the lower and right quadrants of the diagram and approaches the white hole singularity at r = 0. Similarly, the spacetime (52) with σ = −1 also covers the left and upper quadrants of the Kruskal diagram and the spacetime (52) with σ = 1 covers the left-lower quadrants. The only set of points the Kruskal spacetime not covered by these patches is the bifurcation surface at r = 2M .
As noticed in Remark 3, the spatial part of the geodesic equations do not depend of the choice of sign in σ and therefore the dynamical system will also be independent of this choice. This has the following interesting consequence. Consider for instance a future directed causal geodesic stating in the region r < 2M in the lower part of the diagram (for simplicity we call this the white hole region and by black hole region we refer to the domain r < 2M in the upper part of the diagram). This geodesic can be described in the Kerr-Schild metric (52) with σ = 1. After a finite value of its affine parameter, the geodesic will approach r = 2M . This can happen with either v → +∞ or v → v 0 finite. In either case, since the dynamical system only involves the spatial coordinates, this portion of geodesic will have a limit point p in the phase space satisfying r(p) = 2M . Irrespectively of which spacetime point q is approached (even it the point lies on the bifurcation surface), the geodesic can be continued further as a portion of a geodesic in the spacetime (52) with σ = −1 having past endpoint at q. The fact that the dynamical system for the spatial coordinates is independent of σ implies that the trajectory will be described in a single phase space, i.e. the change of spacetime chart will pass fully unnoticed in the phase space of the dynamical system. Thus, we will be able to describe the full geodesic as a single trajectory in the phase space, without having to determine in which spacetime coordinate chart we are working at each portion. As we will see below, this is only possible due to the presence of a excluded region in the phase space. In turn, this excluded region arises as a consequence of the restrictions in the initial data imposed by the condition that the trajectory describes a future directed causal geodesic.

Regularized dynamical system
Lemma 14. The McGehee regularization for the dynamical system that describes the spatial part of the set of geodesic trajectories with angular momentum L in the Kruskal spacetime is where µ = 1, 0, −1 for timelike, null and spacelike geodesics, respectively. The system admits the energy first integral Proof. We can apply Proposition 5 with h = 2M r , which corresponds to the spacetime (52). Equations (28), (29) and (30) becomë where dot means derivative with respect to proper time, affine parameter or arc length depending on the value of µ and the energy is E =Ṫ 2M r − 1 + 2M rṙ , from (31). Given a geodesic, we can rotate the Cartesian coordinates so that the geodesic lies in the plane {x 1 , x 2 } and define x as the complex coordinate x = x 1 + ix 2 ,. The equations of motion (57) becomeẋ = y, This flow is singular at r = 0 and we can apply the McGehee regularization described above. From Corollary 12 and the fact that |x| 3 V (|x|) admits a C 1 extension to x = 0 with non-zero value at this point, we find that the optimal value for regularization is β = min{−1, − γ 2 } = − 3 2 . Applying Lemma 11 the following regularized system is obtained

Excluded regions
The dynamical system (53)-(55) describes all future directed causal geodesics in the Kruskal spacetime. However, not all trajectories in this dynamical system correspond to future directed causal geodesics in this spacetime. The reason is that the set of initial data {r 0 ,ṙ 0 } for future directed causal geodesics is constrained by Proposition 5. Given the relation between u andṙ, this implies that not all points {r, u, θ} in the phase space describe future causal geodesics in the spacetime. We will call the allowed region the set of points in the phase space corresponding to future directed causal geodesics and the excluded region its complement. Let us determine these sets. As discussed in Proposition 5, at points where h < 1, i.e. r > 2M , there is no restriction on the possible values ofṙ 0 , and consequently no restrictions on u arise. On the other hand, when h > 1, i.e. r < 2M , then For h = 1 (r = 2M ),ṙ 0 is restricted to satisfy withṙ 0 = 0 only if L = 0 and the geodesic is null (µ = 0). Given that u = r −βṙ = r 3 2ṙ the allowed region for geodesics in the advanced Eddington-Finkelstein spacetime (i.e. σ = −1) is So, the excluded region is Similarly, for geodesics in the retarded Eddington-Finkelstein spacetime (σ = 1) the excluded region is with (u = 0, r = 2M ) excluded if L = 0 or µ = 0.
As discussed above, the bifurcation surface at r = 2M is not included in either the advanced or retarded Eddington-Finkelstein spacetime. This is the reason why the the point (u = 0, r = 2M ) when either L = 0 or µ = 0 is excluded. Since future causal geodesics with L = 0 in the Kruskal spacetime do cross the bifurcation surface, we must incorporate this set of points of the phase space into the allowed region in order to describe all causal geodesics of the Kruskal spacetime. We conclude that the set of phase space points not describing future directed causal geodesics in the Kruskal spacetime is the intersection of both excluded regions after adding (u = 0, r = 2M ) to the allowed regions, namely The existence of these three types of excluded regions is crucial for describing all future directed causal geodesics in the Kruskal spacetime in a single phase space diagram. Indeed, any trajectory passing through an allowed point in the region r < 2M with u > 0 belongs to the r < 2M region of a retarded Eddington-Finkelstein chart of the Kruskal spacetime, i.e. to the white hole region of the spacetime. A trajectory passing through an allowed point in the region r < 2M with u > 0, belongs to an advanced Eddignton-Finkelstein chart of the Kruskal spacetime, i.e. to the black hole region. A point in the phase space with r > 2M belongs to both charts. The boundary of the allowed region is given by the set of points satisfying which can be rewritten as the set of points with = − µ 2 . In fact, the excluded region can be equivalently defined as the set of points for which < − µ 2 , r ≤ 2M , cf. Corollary 7. The graphical representation of the excluded regions for null and spacelike geodesics is displayed in Fig. 2. Notice that, in the null case, there is no excluded region in the limit L = 0. However, in this case the set of points u = 0 correspond to trivial null geodesics consisting of single points with vanishing tangent vector. For non-null geodesics, the excluded region is always non-empty irrespectively of the value of the angular momentum L.

The collision manifold
The submanifold r = 0 is clearly invariant under the flow. Since r = 0 corresponds to the spacetime singularity, this submanifold is called collision manifold. It can be described globally by the coordinates {(θ, u)} so its topology is R × S 1 . The dynamical system (54)-(55) restricted to the collision manifold reads This system has two lines of critical points: one line of stable points at (θ, u) = (θ 0 , − √ 2M L 2 ) and one line of unstable nodes at (θ, u) = (θ 0 , √ 2M L 2 ), where θ 0 ∈ S 1 is an arbitrary value. The phase space portrait in the collision manifold is shown in Fig. 3. For each value of θ 0 , there is a trajectory extending from u = −∞ and approaching u = − √ 2M L 2 as its future limit point, a trajectory from u = − √ 2M L 2 to u = √ 2M L 2 and a trajectory having u = √ 2M L 2 as its past limit point and extending to u = +∞, all of them with θ = θ 0 . One may wonder how these trajectories relate to causal geodesics in the Kruskal spacetime. To see this, note that any such geodesics must have a finite value of . On the other hand diverges at r = 0 (see (56)). In fact, it does so in the following way The limit when |u| approaches the critical points on the collision manifold depends on the details of how this limit is taken. Since, the trajectories joining the stable critical points to the unstable ones within the collision manifold are interior to the excluded region of the phase phase, it follows that such trajectories are completely unrelated to causal geodesics in the spacetime. This is consistent with the fact that → −∞ on these trajectories, while future directed causal geodesics in the region r < 2M must have ≥ − µ 2 . On the other hand, the trajectories on the collision manifold leaving u = + √ 2M L 2 and those approaching u = − √ 2M L 2 correspond to the limit of trajectories of causal geodesics leaving the white hole singularity at r = 0 and approaching the black hole singularity at r = 0 when their energy diverges to +∞. Thus, this set of trajectories on the collision manifold carries interesting information on the causal geodesics in the spacetime.
To analyze the behavior of the particles near the collision manifold, we can linearize the system at its first order as where u 0 (s) its a solution that corresponds to an orbit in the collision manifold and therefore satisfies The general solution for this differential equation The branch s ∈ (−∞, 0) corresponds to the solution leaving the unstable fixed point at u = + √ 2M L 2 and approaching u −→ +∞, while the branch s ∈ (0, +∞) corresponds to the solution extending from u = −∞ and approaching the fixed point u = − √ 2M L 2 . The first-order linearized system in the variables δr and δu is δr = (δr)u 0 , δu = L 2 (δr) + 3(δu)u 0 .
We can easily solve for δr(s): where δr 0 > 0 is an arbitrary integration constant. Thus, the fixed points are approached exponentially fast in the variable s. To analyze the behavior in the τ variable we recall the relation (42), namely τ (s) = (δr(s)) 5 2 . This integration can be performed explicitly, but it is not particularly enlightening. Instead, we will determine a series expansion of δr(τ ) near τ = 0 where τ is chosen so that the particle reaches the singularity r = 0 at τ = 0 (note that τ < 0 for particles approaching the black hole singularity while τ > 0 for particles leaving the white hole singularity). Define x(s) as so that (δr)(s) = (δr 0 )x(s) and hence dx ds = u 0 x. In terms of x, where the sign depends on the branch we are considering (upper sign for the approach to the black hole and lower sign for the white hole). Then Expanding the right-hand side as a series en x, integrating and inverting we find . A plot of the function δr(τ ) in each of the two branches is given in Fig. 4.
These plots describe the approach to the singularity of very energetic particles in the Kruskal spacetime for different values of the angular momentum. Note that the result is the same for massive and massless particles, as one could expect given the very high energy involved.

Massless particles
The reduced dynamical system with µ = 0 takes the form Its phase space is displayed for different values of L in Fig. 5. The fixed points are (assuming The linearization of the system around these points has the following eigenvalues  Thus, all critical points are hyperbolic. The point (r = 0, u = + √ 2M L 2 ) is an unstable node (repulsor), (r = 0, u = − √ 2M L 2 ) is a stable node (attractor) and (r = 3M, u = 0) is saddle point. This saddle point obviously corresponds to the unstable circular orbit for massless particles.

Massive particles
When µ = 1 the reduced dynamical system is with phase spaces displayed for different values of L in Fig. 6. The fixed points are the second pair under the additional condition L 2 ≥ 12M 2 . For L 2 > 12M 2 all fixed points are hyperbolic, with (r = r + (M, |L|), u = 0) being a center (purely imaginary eigenvalues) and (r = r − (M, |L|), u = 0) being a saddle. When L 2 = 12M 2 , there is a bifurcation in the phase space, which can be visualized in the transition between the third and fourth plots in Fig. 6. Given that r + (M, |L|) is an increasing function of |L| with values ranging from (6M, ∞) we recover the well-known fact that the innermost stable circular orbit (ISCO) in Schwarzschild is r = 6M . The fixed point r − (M, |L|) is a decreasing function of |L| with values ranging from 6M (when |L| → 2 √ 3M ) to 3M when |L| → +∞). We thus recover easily all well-known results for geodesics in Schwarzschild outside the horizon. The approach here, however, is perfectly regular both across the horizon at r = 2M and even at the singularity r = 0. Moreover, it allows us to treat all points in the Kruskal spacetime with a single dynamical system.

The Reissner-Nordström dynamical system
The Reissner-Nordström spacetime [Reissner, 1916, Nordström, 1918 corresponds to the most general solution for Einstein equations with electromagnetic field when spherical symmetry is assumed. The maximal extension of this spacetime is covered by a infinite amount of copies of four basic patches. As in the Kruskal spacetime two of the patches are isometric to the Reissner-Nordström advanced Eddington-Finkelstein spacetime and the remaining two to the Reissner-Nordström retarded Eddignton-Finkelstein. Each one of these spacetimes consist of the manifold R × R + × S 2 with respective metrics where σ = −1 for the advanced case and σ = 1 for the retarded one. The coordinates takes values in V ∈ R, r ∈ R + and dΩ 2 is the round unit metric on the sphere. Performing the same coordinate change as in the Kruskal case, the metric is transformed into its Kerr-Schild form [Kerr and Schild, 1965]: where r = √ x · x and the manifold is R × (R 3 \ {0}). As before we choose the orientation so that the nowhere-zero, null vector ∂ T + σ x | x| ∂ x is future directed. As is well-known, the global structure of the maximal extension of the Reissner-Nordström spacetime depends strongly on whether |Q| > M , |Q| = M or |Q| < M . For the discussion below we concentrate on the latter case because it corresponds to a non-extremal black hole. We emphasize, however, that all the dynamical systems in this Section are valid for any value of Q and M , so they can be used to study geodesics in the extremal black hole or naked singularity cases as well.
When |Q| < M , the basic structure of the maximal extension of the Reissner-Nordström spacetime has now two bifurcation surfaces located, respectively, at the intersection of the two smooth null hypersurfaces with r = r + and the two smooth null hypersurfaces with r = r − , where r ± = M ± M 2 − Q 2 . We call each one of these four hypersurfaces a horizon. The patch with σ = −1 covers the structure unit that goes in ascending way from r = ∞ → r = r + → r = r − → r = 0 (starting from either the right or the left) while the patch with σ = 1 covers the structure unit that goes in an ascending way from r = 0 → r = r − → r = r + → r = ∞ (starting from the right or the left). Notice that there is no causal geodesic starting from a left-most quadrant which, after crossing the null hypersurface r = r + then goes across the portion of the null hypersurface at r = r − lying at the left of the diagram (and similarly for geodesics starting at a right-most quadrant). The reason is that the Killing 1-form ξ := g(ξ, ·) where ξ = ∂ T (defined on a single Eddington-Finkelstein patch) is integrable with orthogonal hypersurfaces foliating the region r − < r < r + with timelike leaves. Let us label these leaves by t. As a consequence we have ξ = Gdt on r − < r < r + where G is a smooth function. Consider the conserved energy for the geodesic, i.e. ∂ T , u = −E where u is the tangent vector. In the region between r − and r + , in order for the geodesic to enter from the left portion of the null hypersurface r = r + and leave across the left portion of the hypersurface r = r − , the geodesic must become somewhere tangent to a hypersurface t = const. At this point we have −E = ξ(u) = Gdt(u) = 0. But E = 0 is impossible for a geodesic starting on the left-most region where ξ is timelike. A similar argument applies to geodesics starting at the right-most quadrant.
As discussed in Section 5, the fact that the dynamical system for the spatial coordinates is independent of σ implies that the trajectory will be described in a single phase space. To understand the geodesic flow along the maximal extension of the Reissner-Nordström we need to analyze the excluded regions. We first write down the regularized dynamical system.

Regularized dynamical system
Lemma 15. The McGehee regularization for the dynamical system that describes the spatial part of the set of geodesic trajectories with angular momentum L in the Kruskal spacetime is where µ = 1, 0, −1 for timelike, null and spacelike geodesics, respectively. The system admits the energy first integral Proof. We can apply Proposition 5 with h = 2M r−Q 2 r 2 , from (60). Equations (28), (29) and (30) Adapting coordinates so that the geodesic lies in the {x 1 , x 2 } plane and defining the complex variable x = x 1 + ix 2 , the equations arė We now apply the McGehee regularization. Since |x| 4 V (|x|) admits a C 1 extension to x = 0, Corollary 12 fixes the optimal value for regularization as β = min{−1, − γ 2 } = −2 and Lemma 11 gives the regularized system (61)-(63) as well as the constant of motion .

Excluded regions
As in the Schwarzschild case, not all trajectories of the dynamical system (61)-(63) correspond to future directed causal geodesics in the Reissner-Nordström spacetime. Proceeding as in Schwarzschild and exploiting Proposition 5, it is straightforward to find that the excluded region of the phase space diagram is In addition, for geodesics in the σ = −1 Eddington-Finkelstein patch u must be lie below the forbidden region (i.e. u ≤ 0 in the strip r − ≤ r ≤ r + ), while geodesics in the σ = +1 Eddington-Finkelstein patch must lie above the forbidden region (i.e. u ≥ 0 in the strip r − ≤ r ≤ r + ). Note that, as in Schwarzschild, the boundary of the allowed region is defined by which correspond to the set of points with = − µ 2 and the excluded region is precisely the set < − µ 2 , r − ≤ r ≤ r + , cf. Corollary 7. The excluded regions for timelike geodesics are displayed in Fig. 8 with representative values of Q and L. The excluded regions for timelike/null geodesics show the same behavior on L as in the Schwarzschild case, i.e. in the null case there is no excluded region in the limit L = 0 but then the line u = 0 consists of a family of trivial geodesics and in the timelike case the excluded region is always non-empty for all values of L.
We can now discuss how the behavior of geodesics across the horizons can be described in the phase space diagram. The crucial information is the restriction for u (which, recall, is proportional toṙ) in the domain r − < r < r + . Let us see how this implies that any geodesic traveling from r = r 0 > r + → r + → r − → r 1 and back to r − → r + must have changed the Kerr-Schild patch along the way (by changing the value of σ). Assume for definiteness that the particle starts in a left-most quadrant. After the crossing of r + the particle must necessarily cross r = r − . This well-known fact can be directly deduced also from the dynamical system because the allowed region in the domain r − < r < r + when σ = −1 has u < 0 everywhere. The crossing of r = r − happens in the right part of the Kerr-Schild patch, as discussed before. Thus the trajectory is still contained in the original Kerr-Schild patch and, at the same time, it has also entered a new patch with σ = −1. Since the collision manifold cannot be attained (see below) the geodesic must necessarily reach a point where u = 0 and cross to positive values of u. From then on, the curve crosses r = r − towards larger values of r. At this crossing, the trajectory necessarily leaves the original Kerr-Schild patch. It does so with u > 0 which is compatible with the fact that we now have σ = −1 and hence the allowed region lies above the forbidden bubble. The geodesic necessarily crosses r = r + and enters a different asymptotic region from which it started. This behavior is of course well-known. What is important here is that a single phase space allows for a complete description of the complicated spacetime trajectory by simply noticing that each time that the forbidden region is encircled, we have moved one step in the ladder in Kerr-Schild patches in the maximal extension of Reissner-Nordström. This is possible only because (i) the phase space has an excluded region and (ii) the crossing at r = r − is not ambiguous, forcing the particle to stay in the same Kerr-Schild patch it started.

The collision manifold
The Reissner-Nordström collision manifold has again topology R × S 1 and global coordinates {(θ, u)}. The dynamical system (62)-(63) restricted to the collision manifold reads which has no fixed points. This is a manifestation of the fact that, in Reissner-Nordström, the singularity has a repulsive behavior, unlike the Schwarzschild case. This is already indicated by the fact that, for |LQ| = 0, lim r→0 = +∞ irrespective of the value of u. So, no physical particles with finite energy can access the collision manifold and, the larger the value of they have, the closer to the collision manifold r = 0 they can get.
To analyze the repulsive nature of the collision manifold in more detail let us linearize the dynamical system to its first order around the collision manifold. Thus, we write r(s) = δr(s), where u 0 (s) is a trajectory on the collision manifold, so that it satisfies The general solution for this differential equation satisfying u 0 (0) = −∞ is u 0 (s) = −|LQ| cot(2|LQ|s).
We can easily solve for δr(s): where δr( π 4|LQ| ) = δr 0 . If follows that, no matter how close to the collision manifold we get (at s = π 4|LQ| where δr(s) is minimum), the trajectory never touches the collision manifold and in fact, it separates from it very quickly. To understand how fast this happens we need to change to the τ variable. Given the relation (42), namely τ (s) = (δr(s)) 3 .
As in the Schwarzschild case we will determine a series expansion of δr(τ ) near τ = 0, where τ is chosen to fulfill τ = 0 when s = π 4|LQ| . Define x(s) as x(s) = 1 sin(2|LQ|s) so that (δr)(s) = (δr 0 )x(s) and hence dx ds = u 0 x. In terms of x Expanding the right-hand side as a series en x, integrating and inverting we find where a = |LQ| (δr 0 ) 3 . It is interesting to note that inserting Q = 0 in (65) we do not recover the Schwarzschild case. This is because the value of the parameter β adapted to Reissner-Nordström is different to that of Schwarzschild. Thus, in the Schwarzschild subcase of Reissner-Nordström we have overkilled the singularity and the fixed points that previously existed at r = 0, u = u ± have both collapsed to u = 0. This collapse can be detected directly on the Reissner-Nordström phase space because the fixed point {u = 0, θ = θ 0 } is no longer hyperbolic when Q = 0. Another way of seeing this is by comparing the shape of the excluded regions of phase space {r, u} for Q = 0 with the shape of the excluded regions in the Schwarzschild phase space. While in the latter case the excluded regions covered a non-empty interval of {r = 0}, the Reissner-Nordström regularization is such that the bubble displayed in Fig. 8 , which stays separated from the collision manifold when Q = 0, just touches the line {r = 0} in the limit Q = 0. Thus the whole line of excluded points in the Schwarzschild regularization has collapsed to a point in the Reissner-Nordström regularization of the Schwarzschild spacetime.  The reduced dynamical system when µ = 0 is The phase portrait can be viewed for different values of L and Q in Fig. 9. The fixed points of this system (assuming LQ = 0) are In the strict inequality case, the two fixed points are hyperbolic, with the point (u, r) = (0, R + (M, Q)) being a saddle point (with real eigenvalues of opposite sign) and the point (u, r) = (0, R − (M, Q)) being a center (purely imaginary eigenvalues). It is straightforward to check that the point (u = 0, r = R + (M, Q)) always lies in the allowed region. The point (u = 0, r = R − (M, Q)) lies in the excluded region as soon as this region is non-empty, i.e. for |Q| < M .
As discussed in Section 6.1.1, curves encircling the excluded region when Q = 0 have periodic properties in the phase space but they are really moving upwards in the Penrose-Carter diagram changing form the white hole patch to the black hole patch as many times as needed. Note that with Q = 0 we recover the Schwarzschild fixed points but, as previously noticed in corollary 12, the phase portrait is nevertheless different because of the different choice of β. This is also the case when µ = 1. When µ = 1 the dynamical system takes the form r = ru, u = r r L 2 − M r + Q 2 − 3L 2 M + 2 L 2 Q 2 + u 2 , with phase spaces displayed for different values of L and Q in Fig. 10. The portraits show very clearly the repulsive nature of the singularity discussed above. The fixed points lie on the line u = 0 and are given by the roots of

Massive particles
The root structure of this polynomial is not uniform in the parameters {M, Q, L}. Let us concentrate for definiteness in the most interesting case L = 0 and 0 < |Q| < M . It turns out that this equation always has one real solution which lies inside the excluded region and corresponds to a hyperbolic critical point that happens to be a center (purely imaginary eigenvalues). Moreover, there exists L 0 (M, Q) > 0 such that, for 0 < |L| ≤ L 0 this is the only root. For |L| = L 0 , there is second root which is double (and hence a non-hyperbolic fixed point for the dynamical system). For |L| > L 0 there are two additional hyperbolic points, both lying outside the excluded region to its right. The one closer to the excluded region is a saddle and the one with largest value of r is a center. The function L 0 (M, Q) is defined as the only positive and real solution of L 6 8Q 2 − 9M 2 + 6L 4 18M 4 − 21M 2 Q 2 + 4Q 4 + 3L 2 8Q 6 − 3M 2 Q 4 + 8Q 8 = 0.
Existence of a unique positive solution of this equation is guaranteed for 0 < |Q| < M . The function L 0 (Q, M ) is displayed in Fig. 11, note that L 0 (Q = 0) = 2 √ 3M . These points are analogous to the well-known critical points in the Schwarzschild spacetime and of course they coincide in the limit Q = 0.

Conclusions
The main result of this work is a dynamical systems method of analyzing causal geodesics in stationary and spherically symmetric Kerr-Schild spacetimes. A remarkable result is that the geodesics can be globally described as the motion of a Newtonian particle in the presence of a radial potential. For spacetimes with singularities at the center we have developed a generalization of the MacGehee transformation that allows for a regularization of the origin and hence for a description of the approach to the singularity in terms of regular variables. In particular, the dynamics at the collision manifold can be analyzed, which gives us useful information for the physical trajectories. We have applied this method to the Schwarzschild and Reissner-Nordström spacetimes. Besides the regularized analysis of the singularity in these spacetimes, we have emphasized the importance of the presence of excluded regions, which in effect makes the phase space diagram acquire a non-trivial topology. This topology and the property that the phase space portrait is independent of whether we deal with an advanced or a retarded Kerr-Schild patch allows one to study in the geodesic motion is spacetimes with complicated global behavior (e.g. the Reissner-Nordström spacetime) in terms of a single two-dimensional phase space portrait.
A Absence of Mcgehee transformation decoupling the system in (u, v) In this appendix we prove a no-go Theorem for McGehee-type transformations capable of decoupling the dynamic equations of a point particle on a plane under the influence of a general radial potential V (|x|). Specifically, we intend to analyze the dynamical systeṁ where {x, y} are coordinates on an open subset of N of C 2 and ∇ = ∂ x 1 + i∂ x 2 . In view of the transformation proposed by McGehee for the power-law potential (39), we define the generalized MacGehee transformation x = e iθ ξ 1 (r), y = e iθ ξ 2 (r)(u + iv), dτ = ξ 3 (r(s))ds, where ξ 1 , ξ 2 , ξ 3 : R + → R are smooth and non-zero functions to be determined. Note that ξ 1 (r) must be invertible for this transformation to be well-defined. The new coordinates {u, v, r, θ} take values on R (for u, v), in R + (for r) and on S 1 (for θ). We note that making ξ 1 complex does not define a more general transformation since it can be reduced to the above one by redefining the variable θ. Given that we are replacing the power-law potential V (|x|) = |x| −σ by a general radial potential, it is reasonable to keep the general structure of the original McGehee transformation and introduce general functions of r in the transformation. In this sense, we can consider (67)  Proof. In the new parameter s, the dynamical system is where prime is derivative with respect to s. Inserting the transformation (67) yields iθ ξ 1 + dξ 1 dr r = ξ 2 ξ 3 (u + iv), iθ ξ 2 + dξ 2 dr r (u + iv) + ξ 2 (u + iv ) = Λ(ξ 2 )ξ 1 ξ 3 .
From the expression for u , the system is uncoupled if and only if ξ 1 (r) 2(1−β) Λ(ξ 1 (r)) = α 3 for some constant α 3 , i.e. Λ(r) is a power-law. Since the potential V is related to Λ by dV dr = −Λ(r)r it follows that the only case for which the generalized McGehee transformation decouples the system is when the potential itself is a power-law or logarithmic (recall that an additive constant is completely irrelevant in the potential V ).