Geodesic dynamics in Chazy-Curzon spacetimes

TL;DR

This paper numerically analyzes geodesic motion in Chazy-Curzon spacetimes, finding exclusively regular trajectories that imply the existence of Carter's constant, and observes oscillatory motions similar to the MacMillan problem.

Contribution

It provides the first thorough numerical study of geodesics in Chazy-Curzon metrics, revealing regular motion and potential conserved quantities.

Findings

All geodesics exhibit regular motion.

Evidence suggests the existence of Carter's constant.

Oscillatory motions similar to the MacMillan problem are observed.

Abstract

In the last decades, the dynamical studies around compact objects became a subject of active research, partially motivated by the observed differences in the profiles of the gravitational waves depending on the dynamics of the system. In this work, via the Poincar\'e section method, we conduct a thorough numerical analysis of the dynamical behavior of geodesics around Chazy-Curzon metrics. As the main result, we find only regular motions for the geodesics in all cases, which suggest the existence of the so-called Carter's constant in this kind of exact solutions. Moreover, our simulations indicate that in the two-particle Chazy-Curzon solution, some oscillatory motions take place as in the classical MacMillan problem.

Equations

$$d{s}^2=-e^{2\psi }d{t}^2+e^{-2\psi }[e^{2\gamma }(d{\rho }^2+d{z}^2)+{\rho }^{2}d{\varphi }^2],$$

$$0$$

$$\psi =-\frac{m}{\sqrt{{\rho }^2+z^2}},$$

$$\gamma =-\frac{m^2{\rho }^2}{2({\rho }^2+z^2{)}^2}.$$

$$\psi =-\frac{m_1}{\sqrt{{\rho }^2+(z-a{)}^2}}-\frac{m_2}{\sqrt{{\rho }^2+(z+a{)}^2}},$$

$$\gamma$$

$$\underset{\rho ,z\to \infty }{\lim }\psi =0,\mathrm{and}\underset{\rho ,z\to \infty }{\lim }\gamma =0,$$

$$C=-\frac{m_1{m}_2}{2a^2}.$$

$$\frac{m_1{m}_2}{2a^2}\left[\frac{(z^2-a^2)}{|z^2-a^2|}-1\right]=0,$$

$$P_0=m,P_1=0,P_2=-m^3/3,P_3=0,\mathrm{and}{P}_4=19m^5/105,$$

$$P_0=2m,P_1=0,P_2=2a^{2}m-8m^3/3,P_3=0,\mathrm{and}$$

$$P_4=2a^{4}m-64a^2{m}^3/7+608m^5/105,$$

$$2\mathcal{L}=e^{-2\psi }\left[e^{2\gamma }\left({\dot{\rho }}^2+{\dot{z}}^2\right)+{\rho }^2{\dot{\varphi }}^2\right]-e^{2\psi }{\dot{t}}^2,$$

$$\dot{t}=e^{-2\psi }E,\dot{\varphi }=\frac{e^{2\psi }L}{{\rho }^2},$$

$$\Phi (\rho ,z)=e^{-2\gamma }\left(E^2-\frac{e^{4\psi }{L}^2}{{\rho }^2}-e^{2\psi }\right).$$

$$\ddot{\rho }+({\dot{\rho }}^2-{\dot{z}}^2)({\gamma }_{,\rho }-{\psi }_{,\rho })+2\dot{\rho }\dot{z}({\gamma }_{,z}-{\psi }_{,z})+e^{-2\gamma }\left[E^2{\psi }_{,\rho }+(\rho {\psi }_{,\rho }-1)\frac{L^2{e}^{4\psi }}{{\rho }^3}\right]=0,$$

$$\ddot{z}-({\dot{\rho }}^2-{\dot{z}}^2)({\gamma }_{,z}-{\psi }_{,z})+2\dot{\rho }\dot{z}({\gamma }_{,\rho }-{\psi }_{,\rho })+e^{-2\gamma }{\psi }_{,z}\left[E^2+\frac{L^2{e}^{4\psi }}{{\rho }^2}\right]=0,$$

$$\ddot{\rho }+{\dot{\rho }}^2({\gamma }_{,\rho }-{\psi }_{,\rho })+e^{2(2\psi -\gamma )}(\rho {\psi }_{,\rho }-1)\frac{L^2}{{\rho }^3}+e^{-2\gamma }{E}^2{\psi }_{,\rho }=0.$$

$$\alpha =\frac{16m^2{e}^{2m\left[\frac{m{\rho }^2}{({\rho }^2+z^2{)}^2}-\frac{2}{\sqrt{{\rho }^2+z^2}}\right]}}{({\rho }^2+z^2{)}^6}\hat{A}$$

$$\hat{A}$$

Full text

1 Introduction

Unlike in Newton’s theory of gravitation, there is no exact solution for a double-star system in orbital motion in general relativity. This is due to the fact that in the strong field regime the spacetime describing binary systems is strongly nonlinear and time-dependent Poisson2014 . Nowadays, in literature, this problem is tackled in two different ways: by means of approximation methods (e.g. post-Newtonian approximations) or resorting to numerical relativity Blanchet2001 . Despite these partial solutions, there is always an imperative need to find exact solutions describing astrophysical scenarios, even if the solutions are not completely realistic. Exact solutions can be used as initial conditions in realistic situations or to fine-tune numerical codes. For example, the Wald’s electro-vacuum exact solution in Schwarzschild spacetime Wald1974 it has been extensively used to test general relativistic force-free codes (see e.g. Paschalidis2013 ).

Since the seminal papers of Chazy Chazy1924 and Curzon Curzon1925 , there have been many attempts to find a suitable astrophysical exact solution of Einstein and Einstein-Maxwell equations modeling two or more compact objects (see e.g.. the works of Papapetrou Papapetrou1945 , Majumdar Majumdar1947 , Bonnor Bonnor1966 , Hartle & Hawking Hartle1972 , Kramer & Neugebauer Kramer1980 , Emparan Emparan2000 , and Manko et al. Manko2014 ). In all cases, the equilibrium between the constituents is reached either by electrostatic repulsion or by the so-called Weyl struts. In particular, the two-particle Chazy-Curzon metric comes from the superposition of two Newtonian point sources at different positions on the symmetry axis Griffiths2009 . Such superposition gives place to a conical singularity which holds the two sources apart in a static configuration Einstein1936 . At this point, it is important to note that notwithstanding the simplicity and practicality of the Chazy-Curzon metrics, the true nature of the sources has been evaded for decades, Abdelqader2013 and therefore the dynamics of timelike geodesics around this metrics has never been considered in the literature.

The dynamics of test particles orbiting around two compact objects is of physical interest in at least two contexts. First, it could serve as a first approximation to the problem of the existence (or not) of a third isolating integral of motion in the two-body problem in general relativity. In other words, the numerical evidence of chaotic or regular orbits could shed lights on the existence of the so-called Carter’s constant Carter1968 ; Dubeibe2007 . Second, previous studies on the dynamics of strong gravity sources have shown that the sensitive dependence in the geodesic dynamics shall be displayed in the gravitational waveforms, that is because the gravitational waveforms can be written in phase space coordinates and hence the differences on nearby trajectories are translated into differences of nearby waveforms Cornish2001 . Taking into account the previous motivation, in the present paper, we shall study the dynamics around two compact objects using the Chazy-Curzon metric via the Poincaré sections method. Given the setup, this study can be considered as a first step and a starting point toward understanding the general relativistic three-body problem.

The paper is organized as follows: In Section 2 we present the solutions derived by Chazy and Curzon along with its most relevant properties. The effective potential and the equations of motion are introduced in Section 3. In Section 4, we briefly analyze the results for the geodesic dynamics in both cases, and finally, in Section 5 we present our conclusions and discuss some possible applications of the results obtained.

2 The Chazy-Curzon solution

For the sake of completeness, we shall begin by deriving the Chazy-Curzon solution Chazy1924 ; Curzon1925 . To do so, let us start with the simplest metric for a static axisymmetric spacetime, which was first derived by Weyl Weyl1917 and reads as

[TABLE]

where the metric functions $\psi$ and $\gamma$ are only dependent on Weyl’s canonical coordinates $(\rho,z)$.

With this choice of the metric, Einstein’s field equations can be reduced to

[TABLE]

where equation (4) is just the classical Laplace’s equation in cylindrical coordinates. The simplest non-trivial solution of the overdetermined system of equations (2-4) can be calculated by setting the classical expression for the gravitational potential of a point mass source as a solution of (4), that is

[TABLE]

and hence, substituting into Eq. (2) or (3), we get

[TABLE]

Now, it should be noted that due to the linearity of Laplace’s equation it is possible to construct a new solution corresponding to the superposition of two Newtonian point sources with masses $m_{1}$ and $m_{2}$ located at $z=a$ y $z=-a$, respectively. In this case, the metric functions are explicitly given as

[TABLE]

and

[TABLE]

where $C$ is an integration constant that can be determined from the asymptotic flatness condition

[TABLE]

and is given by

[TABLE]

On the other hand, from the elementary flatness condition, $\lim\limits_{\rho\rightarrow 0}\gamma=0$, we get

[TABLE]

which is satisfied if and only if $z>a$ and $z<-a$, while a singularity takes place in the interval $-a<z<a$. Such singularity is interpreted as a Weyl strut, which holds the two particles apart and does not exert a gravitational field Griffiths2009 . The term strut comes from the fact that, at the lowest order, the stress is approximately equal to $m_{1}m_{2}/4a^{2}$, as expected from Newtonian theory.111At this point, it is important to note that such strut can be removed by adding spin to the sources, see e.g. Dietz1982 and Hernandez1993 , however, these metrics will be considered further in a forthcoming paper.

The Chazy-Curzon solution can be interpreted in terms of its Newtonian limit through the multipole moments. The knowledge of multipole moments let us infer the physical meaning of each parameter in the solution, and hence, get an idea of the kind of source represented by the particular solution. Although there are many methods for finding multipole moments, the Hoenselaers & Perjés method Hoenselaers1990 , along with the corrections introduced by Sotiriou & Apostolatos Sotiriou2004 , allows us to calculate all moments in an efficient and accurate manner.

Following the Hoenselaers & Perjés procedure, we find that the first five multipolar gravitational moments $P_{n}$ for a single particle solution (5,6) are

[TABLE]

whence it follows that the parameter $m$, denotes the total mass, the total angular momentum is zero, and all higher-order mass moments are proportional to increasing powers of $m$, i.e., we may infer that the solution describes a static non-spherical source.

On the other hand, the first five multipolar gravitational moments for the two-particle solution (7,8) with $m_{1}=m_{2}=m$, read as

[TABLE]

such that the total mass of the source is $2m$, the total angular momentum is zero with $P_{i}=0$, for $i$ even, and all higher mass moments are proportional to increasing powers of $m$ and $a$. It means that the two-particle solution represents a static pair of non-spherical sources separated by a fixed distance $2a$.

3 Effective potential and equations of motion

In general, the dynamics of a test particle in General Relativity is determined by the space-time curvature, i.e., from the metric. The test particles can exhibit in general two types of motion: bounded and unbounded orbits and, as in the Newtonian orbital dynamics, each type of motion is solely specified by means of the effective potential. Such potential can be derived as follows (see e.g. Dubeibe2016 ): (i) From (1) and the relation $2{\cal L}=g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}$, we get

[TABLE]

with ${\cal L}$ the relativistic Lagrangian and $\dot{x}^{\mu}$ the four-velocity of the test particle. (ii) From the cyclic coordinates $t$ and $\phi$, the associated four-velocities read as

[TABLE]

where the constants $E$ and $L$ are related to the Killing vectors $\xi_{t}$ and $\xi_{\phi}$, representing the total energy and the angular momentum of the test particle, respectively. (iii) For test particles, the Lagrangian satisfies the relation $2{\cal L}=-\delta$, with $\delta=1$ for massive particles. (iv) Using the previous relations and some straightforward algebra, the effective potential can be defined as

[TABLE]

It should also be noted that, from the definition of $\Phi$, the motion is restricted to regions where $\Phi(\rho,z)\geq 0$.

On the other hand, the equations of motion resulting from the Lagrangian formalism can be written as

[TABLE]

while, in the case of equatorial motion in axisymmetric static spacetimes, the equation of motion reduces to

[TABLE]

4 Results and discussions

4.1 Case 1: One-particle solution

Let us start by considering the motion of a test particle in presence of a single Chazy-Curzon source, such that the metric functions are given by (5) and (6) with $m=1$. A typical contour plot of the respective effective potential (16) is shown in Fig. 1 for $L=4$ and different values of $E$. The main difference with the case of fixed $E$ is that as the value of $L$ increases, the confinement region moves to the $+\rho$ direction. As can be noted from Fig. 1, there are three feasible types of trajectories: bounded, unbounded and plunging orbits.

In Fig. 2, a confined orbit is presented for initial conditions inside a closed contour of the effective potential, with initial conditions and parameters $L=4,E=0.472,\rho_{0}=20,z_{0}=-1,\dot{z}_{0}=0,\dot{\rho}_{0}=0.08$ and $\phi_{0}=\pi/6$. The periodicity of this orbit will be determined later through the study of the Poincaré sections. It should be noted that when setting the value of the initial condition on the $z$-coordinate $z_{0}=0$ with $\dot{z}_{0}=0$, the motion is restricted to the equatorial plane.

An additional zone of possible motion is the one that surrounds the origin (see Fig. 1). For test particles inside this region, we observed that their orbits approach the source in an unusual way: if a test particle moves toward the source along the $z=0$ axis it falls directly to the origin, but if it is not the case, the test particle is repelled and then moves to the origin along the $\rho=0$ axis. A detailed analysis of this behavior was carried out by Gautrean & Anderson who after computing the Kretschmann scalar $\alpha=R_{\mu\nu\sigma\tau}R^{\mu\nu\sigma\tau}$, showed that the singular behavior of this invariant quantity depends on the approaching direction to the origin Gautrean1967 .

For the sake of completeness, we have reproduced the calculation of the Kretschmann scalar, obtaining

[TABLE]

with

[TABLE]

A sketch of the behavior of the effective potential $\Phi(\rho,z)$ and the Kretschmann scalar $\alpha$ as we approach the origin along the $\rho$ or $z$-axis, is presented in Fig. 3. In accordance with the results predicted by the geodesic motion, it can be noted that in the first case (upper quadrant of Fig. 3), the Kretschmann scalar $\alpha$ tends to zero, while the effective potential $\Phi(\rho,z)$ tends to infinity, which can be interpreted as a repulsive gravitational phenomenon (see e.g. Herrera2005 ); in the second case (lower quadrant of Fig. 3), the Kretschmann scalar $\alpha$ tends to infinity, while the effective potential $\Phi(\rho,z)$ tends to zero, which can be interpreted as an intrinsic singularity. From the previous discussion, it is easy to understand the directional character of the singularity in the one-particle Chazy-Curzon solution.

In order to fully characterize the geodesic motion of test particles in presence of the one-particle Chazy-Curzon solution, we study the geodesic dynamics by means of the Poincaré surfaces of section. In Fig. 4 we show the Poincaré section $z=0$ in phase space ($\rho,P_{\rho}$) for the set of parameters of energy and angular momentum used in Fig. 2. It can be noted that the available phase space is filled with regular islands and no chaos is observed in the system. The same behavior exhibited in Fig. 4, occurs for all the different values of $E$ and $L$ used in this study. Our numerical results suggest the existence of only integrable geodesics.

4.2 Case 2: Two-particle solution

Following the same procedure outlined in the previous subsection, now we consider the motion of test particles in presence of a double Chazy-Curzon source with the metric functions given by (7) and (8) and $m_{1}=m_{2}=a=1$. In Fig. 5 for $L=8$ and using different values of $E$, we show four contour plots of the effective potential (16). Like in the previous case, as the value of $L$ increases, the confinement region moves to the $+\rho$ direction. The structure of the zero velocity surfaces is very similar to the one presented in Fig. 1, except for the fact that in the inner contour there appear two open regions, hence, the three types of trajectories of the one-particle Chazy-Curzon solution are still possible, but, with two available sources to fall in.

As an example of bounded orbit, in Fig. 6 we show the trajectory followed by a test particle whose motion is confined to the zero velocity surface presented in the interval $\rho\in(16.8,29.5)$ and $z\in(-4.75,4.75)$. In essence, this is the same type of bounded periodic orbit found in the one-particle Chazy-Curzon solution. Additionally, in Fig. 7 we show two examples of orbits falling toward the sources. From this figure, it can be seen that there is a reflection symmetry about the $\rho$-axis, which is a consequence of the symmetry properties $(\psi(\rho,z)=\psi(\rho,-z)\land\gamma(\rho,z)=\gamma(\rho,-z))$ of the metric functions. Also, the repulsive character of the Weyl strut can be easily recognized from the same figure. By setting $L=0$, the repulsive potential between the primaries disappears, and some oscillatory type of motions take place, similar to the ones in the classical MacMillan problem MacMillan1911 (see Fig. 8).

Concerning the dynamics of the test particles in presence of the two-particle Chazy-Curzon solution, we perform the same analysis using the Poincaré surfaces of section. In Fig. 9 we show the Poincaré section $z=0$ in phase space ($\rho,P_{\rho}$) for the set of parameters of energy and angular momentum used in Fig. 6. From the Kolmogorov-Arnold-Moser (KAM) theory Holmes1982 , the set of ordered points represent periodic orbits, while closed curves should correspond to the quasi-periodic orbits. The same behavior exhibited in Fig. 6, occurs for the different values of $E$ and $L$ swept in this study, i.e. almost ten thousand different combinations. As in the case of a single Chazy-Curzon source, our numerical results suggest the existence of only integrable geodesics in the two-particle Chazy-Curzon solution.

5 Conclusions

In the present paper, we have conducted a thorough numerical analysis of the dynamical behavior of geodesics around Chazy-Curzon sources. We start with the one-particle Chazy-Curzon metric and then we move forward to the two particle solution, finding that the multipolar structure suggests that the one-particle Chazy-Curzon solution describes a static non-spherical source, which possesses a directional singularity. On the other hand, the multipolar structure of the two-particle solution shows us that this solution represents a static pair of non-spherical sources separated by a fixed distance $2a$. The two-particle Chazy-Curzon solution satisfies the asymptotic flatness conditions, but the elementary flatness condition is not satisfied in the interval $-a<z<a$, where a conical singularity (Weyl strut) takes place. The strut behaves as a repulsive field in the two-particle Chazy-Curzon solution for $L\neq 0$, but setting $L=0$, the repulsive potential between the primaries disappears and some oscillatory motions take place, as in the circular Sitnikov problem.

Once we had an idea of the physical meaning of both solutions, the equations of motion for time-like test particles were derived by means of the Lagrangian formalism. In each case, the system of equations was integrated using a Runge-Kutta-Fehlberg algorithm with an adaptive step size. With this method, the relative error in the system’s energy conservation allows us to determine the maximum integration time (in most of the cases $\approx 10^{4}$) in which the error tolerance is below $10^{-10}$. Given the parameters $E$ and $L$ and the initial conditions $\rho_{0},z_{0},\dot{\rho}_{0}$, the third integral of motion $2{\cal L}=-\delta$ determines $\dot{z}_{0}$. Taking into account that the effective phase-space is three-dimensional, for fixed values of energy and angular momentum, Poincaré sections are a convenient tool in dynamical systems theory to analyze the regularity or chaoticity of motion. It is important to note that unlike in the classical, post-Newtonian or Pseudo-Newtonian system (see e.g. Dubeibe2017a ; Dubeibe2017b ; Zotos2018 ), in the general relativistic case the Lyapunov exponents must not be used for the analysis, due to the non-invariance of these exponents in general relativity Motter2003 .

Concerning the geodesic dynamics, in the present study, we used approximately ten thousand different combinations of values of $E$ and $L$ for each solution. We have scanned the parameter space for energy and angular momentum using the interval $[-1,1]$ for $E$, with a step size of $\Delta E=0.02$, and $[-20,20]$ for $L$, with a step size of $\Delta L=0.4$ (or less if necessary), these intervals define closed contours in the effective potential which is a guaranty for bounded orbits. In both cases (single Chazy-Curzon source and the two-particle Chazy-Curzon solution), the number of initial conditions integrated for each Poincaré section is of 50, which is enough to conclude that the phase space is filled by regular orbits. In other words, our findings suggest the existence of the so-called Carter’s constant in both systems, which shall provide the fourth conserved quantity necessary to uniquely determine all orbits in each spacetime. We hope our contribution to be useful for a further understanding of the general relativistic three-body problem, in particular, our results could be used as initial conditions in realistic binary simulations or to fine-tune numerical codes.

Acknowledgements.

We thank the anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions. One of the authors (FLD) gratefully acknowledges the financial support provided by Universidad de los Llanos and Colciencias (Colombia), under Grants No. 8840 and 8863.