Orbit classification in a pseudo-Newtonian Copenhagen problem with Schwarzschild-like primaries

TL;DR

This study explores how relativistic effects and energy levels influence orbital dynamics in a binary black hole system, revealing significant impacts on orbit types and system complexity compared to classical models.

Contribution

It introduces a pseudo-Newtonian model for the Copenhagen problem with Schwarzschild-like primaries and analyzes the effects of relativistic parameters on orbital behavior.

Findings

Jacobi constant significantly affects orbit final states

Schwarzschild radius alters orbit character and fractality

Relativistic effects increase system complexity

Abstract

We examine the orbital dynamics of the planar pseudo-Newtonian Copenhagen problem, in the case of a binary system of Schwarzschild-like primaries, such as super-massive black holes. In particular, we investigate how the Jacobi constant (which is directly connected with the energy of the orbits) influences several aspects of the orbital dynamics, such as the final state of the orbits. We also determine how the relativistic effects (i.e., the Schwarzschild radius) affect the character of the orbits, by comparing our results with the classical Newtonian problem. Basin diagrams are deployed for presenting all the different basin types, using multiple types of planes with two dimensions. We demonstrate that both the Jacobi constant as well as the Schwarzschild radius highly influence the character of the orbits, as well as the degree of fractality of the dynamical system.

Equations

$$\Phi (r,\theta ,\varphi )=\frac{Gm}{r}+\left(\frac{2r_g}{r_{S2}}\right)\left[\left(\frac{r_{S1}}{r_{S2}}\right){\dot{r}}^2+\frac{s(\theta ,\dot{\theta },\dot{\varphi })}{2}\right],$$

$$r_{S1}=r-r_g,$$

$$r_{S2}=r-2r_g,$$

$$s(r,\theta ,\dot{\theta },\dot{\varphi })=r^2\left({\dot{\theta }}^2+{\dot{\varphi }}^2{\sin }^2\theta \right).$$

$$\omega =\sqrt{\frac{2Gm{\left(R-r_S\right)}^2}{R^3\left(R^2-R{r}_S-r_S^2\right)}},$$

$$U(x,y,\dot{x},\dot{y})=\frac{1-\mu }{r_1}+\frac{\mu }{r_2}+\frac{{\omega }^2}{2}\left(x^2+y^2\right)$$

$$\displaystyle+\frac{1}{2}r_{S}\bigg{(}\frac{\left(2r_{1}-r_{S}\right)\left[\dot{x}\left(x-x_{1}\right)+y\dot{y}\right]}{r_{1}^{2}\left(r_{1}-r_{S}\right)^{2}}$$

$$+\frac{\left(2r_2-r_S\right)\left[\dot{x}\left(x-x_2\right)+y\dot{y}\right]}{r_2^2{\left(r_2-r_S\right)}^2}$$

$$\displaystyle+\omega^{2}\left(\frac{r_{1}^{2}}{r_{1}-r_{S}}+\frac{r_{2}^{2}}{r_{2}-r_{S}}\right)\bigg{)},$$

$$r_1=\sqrt{{\left(x-x_1\right)}^2+y^2},$$

$$r_2=\sqrt{{\left(x-x_2\right)}^2+y^2}.$$

$$\mathcal{L}=\frac{1}{2}\left({\dot{x}}^2+{\dot{y}}^2\right)+\omega \left(x\dot{y}-\dot{x}y\right)+U(x,y,\dot{x},\dot{y}).$$

$$\ddot{x}-2\omega \dot{y}=\frac{\partial U}{\partial x}-\frac{d}{dt}\left(\frac{\partial U}{\partial \dot{x}}\right),$$

$$\ddot{y}+2\omega \dot{x}=\frac{\partial U}{\partial y}-\frac{d}{dt}\left(\frac{\partial U}{\partial \dot{y}}\right).$$

$$C_1(x,y)\ddot{x}+C_2(x,y)\ddot{y}+f\left(x,y,\dot{x},\dot{y}\right)=0,$$

$$C_3(x,y)\ddot{x}+C_4(x,y)\ddot{y}+g\left(x,y,\dot{x},\dot{y}\right)=0,$$

$$\ddot{x}=\frac{g\left(x,y,\dot{x},\dot{y}\right)C_2(x,y)-f\left(x,y,\dot{x},\dot{y}\right)C_4(x,y)}{C_5(x,y)},$$

$$\ddot{y}=\frac{f\left(x,y,\dot{x},\dot{y}\right)C_3(x,y)-g\left(x,y,\dot{x},\dot{y}\right)C_1(x,y)}{C_5(x,y)},$$

$$\mathcal{J}=2U-\left({\dot{x}}^2+{\dot{y}}^2\right)-2\left(\dot{x}\frac{\partial U}{\partial \dot{x}}+\dot{y}\frac{\partial U}{\partial \dot{y}}\right)=C.$$

$$U\approx \frac{G{L}^2}{m{c}^2}\frac{1}{r_1^3},$$

$$E=-\frac{1}{r_a+r_p}\pm {\left(\frac{2r_p{r}_a}{r_p+r_a}\right)}^{1/2},$$

$$r_a(r_p)=-r_p-\frac{E+r_p^2}{E^2-2r_p}-\frac{{\left(r_p^4+2E{r}_p^2+2r_p\right)}^{1/2}}{E^2-2r_p},$$

$$r_a(r_p)\approx -\frac{1}{E}+\frac{(2r_p{)}^{1/2}}{E^2}.$$

$$A_{\text{ce}}\sim r_S^{\beta },$$

$$\displaystyle C_{1}(x,y)=1+r_{S}\bigg{[}\frac{\left(2r_{1}-r_{S}\right)\left(x+\mu\right)^{2}}{r_{1}^{2}\left(r_{1}-r_{S}\right)^{2}}$$

$$+\frac{\left(2r_2-r_S\right){\left(x+\mu -1\right)}^2}{r_2^2{\left(r_2-r_S\right)}^2},$$

$$\displaystyle C_{2}(x,y)=C_{3}(x,y)=r_{S}y\bigg{[}\frac{\left(2r_{1}-r_{S}\right)\left(x+\mu\right)}{r_{1}^{2}\left(r_{1}-r_{S}\right)^{2}}$$

$$\displaystyle+\frac{\left(2r_{2}-r_{S}\right)\left(x+\mu-1\right)}{r_{2}^{2}\left(r_{2}-r_{S}\right)^{2}}\bigg{]},$$

$$C_4(x,y)=1+r_{S}y\left[\frac{2r_1-r_S}{r_1^2{\left(r_1-r_S\right)}^2}+\frac{2r_2-r_S}{r_2^2{\left(r_2-r_S\right)}^2}\right],$$

$$\displaystyle C_{5}(x,y)=1+r_{S}\bigg{[}\frac{2r_{1}-r_{S}}{\left(r_{1}-r_{S}\right)^{2}}+\frac{2r_{2}-r_{S}}{\left(r_{2}-r_{S}\right)^{2}}$$

$$\displaystyle+\frac{r_{S}y^{2}\left(2r_{1}-r_{S}\right)\left(2r_{2}-r_{S}\right)}{r_{1}^{2}r_{2}^{2}\left(r_{1}-r_{S}\right)^{2}\left(r_{2}-r_{S}\right)^{2}}\bigg{]},$$

$$f(x,y,\dot{x},\dot{y})=\frac{\left(1-\mu \right)\left(x+\mu \right)}{r_1^3}+\frac{\mu \left(x+\mu -1\right)}{r_2^3}-x{\omega }^2-2\omega \dot{y}$$

$$\displaystyle+r_{S}\bigg{[}\frac{1}{2}\omega\left(-\frac{\left(r_{1}-2r_{S}\right)\left(x+\mu\right)}{\left(r_{1}-r_{S}\right)^{2}}-\frac{\left(r_{2}-2r_{S}\right)\left(x+\mu-1\right)}{\left(r_{2}-r_{S}\right)^{2}}\right)$$

$$\displaystyle+\frac{x+\mu}{r_{1}^{4}\left(r_{1}-r_{S}\right)^{3}}\bigg{[}r_{1}^{2}\left(r_{1}-r_{S}\right)\left(2r_{1}-r_{S}\right)\left(\dot{x}^{2}+\dot{y}^{2}\right)$$

$$\displaystyle-\left(3r_{1}^{2}-3r_{1}r_{S}+r_{S}^{2}\right)\left(\dot{x}\left(x+\mu\right)+y\dot{y}\right)^{2}\bigg{]}$$

$$\displaystyle+\frac{x+\mu-1}{r_{2}^{4}\left(r_{2}-r_{S}\right)^{3}}\bigg{[}r_{2}^{2}\left(r_{2}-r_{S}\right)\left(2r_{2}-r_{S}\right)\left(\dot{x}^{2}+\dot{y}^{2}\right)$$

$$\displaystyle-\left(3r_{2}^{2}-3r_{2}r_{S}+r_{S}^{2}\right)\left(\dot{x}\left(x+\mu-1\right)+y\dot{y}\right)^{2}\bigg{]}\bigg{]},$$

$$g(x,y,\dot{x},\dot{y})=y\left(\frac{1-\mu }{r_1^3}+\frac{\mu }{r_2^3}\right)-y{\omega }^2+2\omega \dot{x}$$

$$\displaystyle+r_{S}y\bigg{[}-\frac{1}{2}\omega^{2}\left(\frac{r_{1}-2r_{S}}{\left(r_{1}-r_{S}\right)^{2}}+\frac{r_{2}-2r_{S}}{\left(r_{2}-r_{S}\right)^{2}}\right)$$

Full text

††pubyear: 2019††volume: 487††pagerange: 2340–2353

1 Introduction

Black hole binary systems have become the subject of intense research since direct confirmation of their existence was provided by the recent discovery of the gravitational wave signals emitted from about a dozen black hole binary merger events by the LIGO collaboration (Abbott et al., 2016a, b). The black holes involved in these discoveries span a mass range going from $~{}10M_{\odot}$ to $~{}100M_{\odot}$, and are all consistent to have initially formed from the death of massive stars.

There is also strong observational support of the existence of a different population of so-called super-massive black holes (SMBHs), with masses spanning from $10^{5}$ to $10^{10}M_{\odot}$, residing at the center of virtually all galaxies (Beckmann & Shrader, 2010). It is expected that some of these SMBHs will pair up as binaries following the merger of their host galaxies (Begelman et al., 2010). There is indeed evidence of several active galaxies with double nucleus (Komossa et al., 2003; Müller et al., 2015), while the eventual inspiral and merger of some of these SMBH binaries constitutes a prime gravitational wave source for the planned LISA observatory (Amaro-Seoane et al., 2012).

Under general conditions following a galaxy merger, the recently formed SMBH binary will typically be left with a binary separation of the order of 1 pc. This separation is at least two orders of magnitude larger than that needed for gravitational wave emission to be an efficient mechanism to drive the binary towards coalescence. This theoretical difficulty is known as the final parsec problem (Milosavljević & Merritt, 2002). It is considered that three-body interactions between the binary and individual stars in the circumnuclear cluster might play an important role for bridging this separation gap. Most of the studies that have been done of these systems have been based on non-relativistic analysis (Mukhopadhyay, 2002). It is clear, however, that general relativistic effects should be dominant at least for a restricted region of the parameter space available to the interacting stars.

General relativity is a complex but successful geometric theory of gravitation that can be easily applied for the analytical modelling of some symmetric astrophysical scenarios of single compact objects (Stephani et al., 2009). However, for the simplest case of two massive bodies (i.e., SMBHs), there are no known exact solutions, only approximations or numerical solutions (see e.g., Poisson & Will, 2014; Baumgarte & Shapiro, 2010). To avoid the high complexity of the equations of general relativity, analytical post-Newtonian or pseudo-Newtonian (PN) approaches are available (Uryũ & Eriguchi, 1999; Steklain & Letelier, 2006; Blanchet, 2009). The use of such approaches allows one to obtain approximate solutions and to reproduce some relativistic effects with some degree of accuracy.

Pseudo-Newtonian approaches have been widely used in astronomy and astrophysics (Mukhopadhyay, 2002), motivated by the fact that its general predictions of particle dynamics do accurately agree with general relativistic computations (Chakrabarti & Mondal, 2006). Since the seminal paper by Paczyński & Wiita (Paczyński & Wiita, 1980), various authors have proposed different PN potentials to model accretion disks around black holes (see e.g., Løvås, 1998; Artemova et al., 1996; Semerak & Karas, 1999). In this context, it has been shown that these pioneer potentials have serious limitations when applied to non-circular trajectories, and hence, several alternative potentials have been developed recently that extend the range of validity of this kind of tools (see e.g., Witzany & Lämmerzahl, 2017; Tejeda & Rosswog, 2013; Wegg, 2012).

On the other hand, the study of the motion of test particles around binary black hole systems based on reduced three-body problem approaches (Roche approximations) is a relatively recent line of research (see e.g., Dubeibe et al., 2017). In point of fact, in a recent paper, some of us have given a general theoretical framework based on the Paczyński & Wiita potential for describing the gravitational field of two non-Newtonian primaries (Zotos et al., 2018). The possibility of using PN potentials for studying the dynamics around binary black holes is justified by the fact that all PN approaches are just modifications to the Newtonian physics designed for mimicking some particular general-relativistic behaviours. In the present work, we examine the nature of motion of the planar PN Copenhagen problem, in the case of Schwarzschild-like primaries. The primaries are modelled with the aid of the Tejeda-Rosswog potential (Tejeda & Rosswog, 2013), which can reproduce exactly general trajectories of test particles in presence of a Schwarzschild black-hole. By adopting this PN potential we intend to capture and replicate some general relativistic effects, relevant to the Copenhagen problem. Moreover, we expect these effects to become more prominent and obvious, as we increase the value of the Schwarzschild radius. With the present work we expect to contribute to the understating of the role played by general relativistic effects in the dynamics of test particles around SMBH binaries. In particular, our numerical analysis has tangible astrophysical applications, since it describes the planar motion of a test particle in a binary system of two equally massed black holes, while it also demonstrates the influence of the general relativistic effects (i.e., the Schwarzschild radius) on the character of the orbits of the test particle

It should be stressed that for studying this particular problem, the following assumptions were made: i) The primaries move in circular closed loop orbits about their common mass center; ii) the total potential of the combined gravitational field can be expressed as the sum of the individual potentials of each primary; iii) the motion of the test particle takes place on the $(x,y)$ with $z=0$; and iv) the binary separation is larger than the length scale over which the gravitational wave emission become relevant for the orbital dynamics.

The present article is, to a considerable extent, a continuation of the work initiated in Zotos et al. (2018), in which we determine the character of motion of a test particle, moving in a binary system of two compact objects. In Section 2 the mathematical formulation of the astrophysical system, along with the disaggregation of the critical energy levels, are presented. In Section 3, the basin diagrams are sketched on the configuration plane and then compared to the distribution of close encounter and escape times, spent by each orbit. Next, we discuss the parametric variation of the percentages of all types of orbits, classified in our study. We close this section by showing how the Schwarzschild radius affects the extent of the main types of basins: bounded, close encounter and escape. The most important conclusions of this numerical study are summarised in Section 4.

2 Pseudo-newtonian model

Recently, in Tejeda & Rosswog (2013) a new PN potential (TR herafter) has been introduced, which reproduces many relativistic features with better accuracy than many commonly used simple potentials, such as the Paczyńsky-Wiita (PW) (Paczyński & Wiita, 1980) and the Nowak-Wagoner (NW) (Nowask & Wagoner, 1991) potentials. In particular, the TR potential was derived in a self-consistent manner by taking the low-energy limit of the general relativistic equation for geodesic motion of test particles around a Schwarzschild black hole. This is in contrast to commonly used pseudo-Newtonian potentials (e.g., PW and NW) that are introduced in an ad-hoc basis in order to capture some particular relativistic feature. In particular, the PW potential was conceived with the express purpose of capturing the location of the marginally stable circular orbit (also known as innermost stable circular orbit or ISCO), but other relevant features are not so well reproduced (e.g., epicyclic frequencies and pericentre precession), not to mention that general particle motion (non-circular trajectories) cannot be reproduced with this potential. On the other hand, the TR potential captures several key relativistic effects exactly (e.g., the location of the marginally stable and marginally bound circular orbits, the radial dependence of the energy and angular momentum of circular orbits, the geometry of general orbits), and at the same time gives an accurate description of other relevant relativistic features (e.g., pericentre precession, epicyclic frequencies). More information on the performance of the TR potential as well as a comparison with other pseudo-Newtonian potentials can be found in Tejeda & Rosswog (2013).

It should also be noted that all three types of potentials (PW, NW and TR) coincide far away from the compact object, as they converge to the Newtonian potential, but in the relativistic regime their dynamical properties are totally different.

The TR potential, in spherical coordinates $(r,\theta,\phi)$, can be written as

[TABLE]

where $r_{g}=Gm/c^{2}$ corresponds to the gravitational radius, while

[TABLE]

This potential contains two terms of which the first represents the usual, classical Newtonian potential, while the second term dictates the contributions due to the kinetic energy.

In this context, the two body problem can be formulated in the following way: it is assumed that both primaries have the same mass $(m_{1}=m_{2})$ and they freely move in a circle of radius $R/2$ around the barycenter, then we equate the radial force of gravity, exerted by each primary body on the other, with the centrifugal force, such that the angular velocity must satisfy

[TABLE]

where $r_{S}=2Gm/c^{2}$ denotes the Schwarzschild radius. It should be noted, that by setting $r_{S}=0$ Eq. (3) reduces to the classical form.

Using the standard notation for the normalization of constants in the planar circular restricted three-body problem, i.e., setting $G(m_{1}+m_{2})=1$ for the masses and $x_{1}+x_{2}=1$ for the distance between the primaries, the origin always coincide with the barycenter in the synodic system, such that introducing the parameter $\mu$, we have that $Gm_{1}=1-\mu$, $Gm_{2}=\mu$, $x_{1}=-\mu$ and $x_{2}=1-\mu$. It should be noted, that the value $\mu=1/2$ corresponds to the Copenhagen problem111For more details regarding the dimensionless units, we refer the reader to the reference Szebehely (1967)..

Hence, the PN effective potential that describes a relativistic binary system of two compact objects, in the equatorial plane, can be written in cartesian coordinates as

[TABLE]

with

[TABLE]

Now, we use Eq. (4) to construct the Lagrangian per unit mass in a synodic reference frame for the planar RTBP (see e.g., Szebehely, 1967)

[TABLE]

The respective equations of motion for the third body can be calculated with the aid of the Euler-Lagrange equations, as

[TABLE]

Due to the explicit dependence of the effective potential on the velocities, we get a system of two equations of the form

[TABLE]

whose solutions are given by

[TABLE]

where $C_{5}=C_{1}C_{4}-C_{2}C_{3}$.

The Jacobi integral is the only constant of motion and reads as

[TABLE]

Note, that unlike in the classical RTBP, the last term is not zero due to the velocity dependent terms in the PN effective potential (4).

In Appendix A we provide the explicit formulae for all the terms entering the equations of motion as well as the Jacobi integral.

The binary system of the two compact objects has five equilibrium points, all of them lying on the plane $(x,y)$. The firth three libration points, that is $L_{i}$, with $i=1,...,3$, are collinear equilibrium points on the horizontal $x$-axis, while the other two $L_{4}$ and $L_{5}$ lie on the vertical $y$-axis (see Fig. 1). These equilibrium points divide the plane of motion $(x,y)$ into two areas: (i) the interior region (IR) for $x(L_{3})<x<x(L_{2})$ and the exterior region (ER) for $x<x(L_{3})$ or $x>x(L_{2})$.

The Jacobi values at the points of equilibrium are critical values and they are denoted as $C_{i}$, $i=1,...,4$. The geometry of the Hill’s regions configuration on the plane $(x,y)$ is a function of these critical energy levels. More specifically:

• •

Energy interval I $(C>C_{1})$: The test particle can move either very close or very far from the primaries, while both IR and ER are disconnected.

• •

Energy interval II $(C_{2}=C_{3}<C<C_{1})$: The third body is allowed to circulate around the primary bodies in the IR, because the channel around the equilibrium point $L_{1}$ is open.

• •

Energy interval III $(C_{4}=C_{5}<C<C_{2}=C_{3})$: The test particle can move from IR to ER and vice versa, due to the open necks around the saddle points $L_{2}$ and $L_{3}$.

• •

Energy interval IV $(C<C_{4}=C_{5})$: The third body is allowed to move on the entire configuration space $(x,y)$, due to the complete absence of the energetically not allowed regions of motion.

3 Classification of the initial conditions

All orbits that we are going to classify start on the plane $(x,y)$ with starting velocities $\dot{x_{0}}=0$, while in all cases $\dot{y_{0}}$ is derived through the Jacobi integral. The character of the orbits will be displayed by colouring the corresponding starting conditions $(x_{0},y_{0})$, according to the final state of the orbits. The computational methodology and all the numerical criteria we used for the taxonomy of the orbits are explained in detail in Section 4 of Zotos et al. (2018).

The following subsections correspond to the different Hill’s regions configurations and of course to different Jacobi values. We chose to present the nature of motion of the configuration space $(x,y)$ for those Jacobi values, where the orbital content is interesting. It should be noted, that we do not present results for the case where $C>C_{1}$ (energy interval I), mainly for saving space, simply because the orbital content for this case is not interesting.

3.1 Energy interval II $(C_{2}=C_{3}<C<C_{1})$

The first case corresponds to the second energy interval $C_{2}=C_{3}<C<C_{1}$, that is when the throat (channel) near the equilibrium point $L_{1}$ is open, which means that the test particle can move around both primary bodies, inside the IR. In Fig. 2 we present the orbit classification for $C=3.6$, for the classical Copenhagen problem (panel (a)) and for the PN Copenhagen problem (panel (c)). In all the following similar figures the diagrams on the first row correspond to $r_{S}=0$, while those on the second row correspond to $r_{S}=0.01$. We chose to present the results for the extreme case with $r_{S}=0.01$ so as to observe the highest impact of the PN dynamics on the nature of the orbits.

When $r_{S}=0$ it is seen that inside the IR there exist two stability islands of regular bounded motion, which are formed by starting conditions of orbits moving around only one of the primary bodies. Just outside these bounded basins there are thin layers of chaotic trapped orbits, while the rest of the IR is filled with a highly chaotic mixture of close encounter orbits. One can argue, that relatively close to the primaries we can identify local basins of close encounter orbits. The ER is dominated by a unified area corresponding to escaping motion, while at the left side there is a small stability basin, composed of starting conditions of orbits circulating around $P_{1}$ and $P_{2}$.

For the PN problem, with $r_{S}\neq 0$, the orbital structure changes as follows:

• •

The area of the bounded motion in both IR and ER is smaller.

• •

The chaotic layer, around the stability islands, disappears.

• •

The close encounter basins in the IR are more prominent.

Parts (b) and (d) of Fig. 2 show the way in which the times for close encounter and escape are distributed on the $(x,y)$ plane, using tones of red and blue, respectively. By inspecting the range of values at both colour bars it becomes evident that the time of escape of the orbits is almost unperturbed by the increase of $r_{S}$. On the other hand, the close encounter time of the orbits in the case of the PN problem is about ten times shorter, with respect to the classical problem. By the term time of escape $(t_{\rm esc})$ we refer to the time needed, so as the test particle to cross the escape radius $R_{\rm esc}=10$ (with origin at $(0,0)$). In the same vein, the close encounter time $(t_{\rm ce})$ is the time required for a test particle to enter the disk with radius $R_{c}=r_{S}+r_{\rm ce}$, around each primary, with velocity pointing inwards. As in Zotos et al. (2018), $r_{\rm ce}=10^{-5}$.

3.2 Energy interval III $(C_{4}=C_{5}<C<C_{2}=C_{3})$

The third energy interval $C_{4}=C_{5}<C<C_{2}=C_{3}$ corresponds to the case where the two channels, in the vicinity of the points of equilibrium $L_{2}$ and $L_{3}$, are open and the test particle, coming from the IR, can enter the ER and vice versa. The character of motion on the $(x,y)$, for $C=3.3$, is revealed in Fig. 3. In the classical RTBP we see that the stability islands inside the IR are still there, while there no stability islands of regular motion around both primaries. With a closer look at panel (a) of Fig. 3 one can identify a small portion of trapped chaotic starting conditions. Furthermore, starting conditions of close encounter orbits exist also in the ER, forming thin trails, leading to the primaries.

When $r_{S}=0.01$ we observe the following changes:

• •

Regular bounded motion is only possible in the IR, however, the corresponding area of the bounded basins is smaller.

• •

There aren’t signs of chaotic behaviour.

• •

Starting conditions which lead to close encounter to one of the primaries form well-defined basins, which extend beyond the IR.

Once more, according to parts (b) and (d) of Fig. 3, far from the primaries, the diagrams of escape times are identical in both cases (classical and PN problems). On the opposite, the close encounter time spent by the orbits in the PN case is about 2.4 times shorter, compared to the close encounter time of the same starting conditions for the classical version of the problem.

3.3 Energy interval IV $(C<C_{4}=C_{5})$

The last energy interval corresponds to the case where the test particle can freely move to the plane $(x,y)$, without any restrictions of energetically not allowed regions of motion. From a dynamical point of view, this energy range is the most interesting region because the orbital content changes drastically, with respect to the energy level.

Fig. 4 shows the nature of orbits when $C=2$. In both the classical and the PN cases the areas of regular bounded motion seem to be about the same. The only apparent difference is about the basins of the close encounter motion which in the case of the PN dynamics are much more defined. Also for this energy level the distribution of the time of escape of the orbits is exactly the same in both cases, while in the case of the PN dynamics the orbits need about 2.4 times less time to perform a close encounter to one of the primaries, with respect to the classical version of the problem.

The colour diagrams in Fig. 5 display the classification of the orbits when $C=-0.5$. For this energy range escaping motion is the dominant type of motion, while bounded motion around the primaries is no longer possible. Close encounters are allowed in both cases, however, in the PN problem the corresponding basins are larger, implying that in this case, the phenomenon of close encounters is more probable. As for the distribution of the close encounter and escape time of the orbits the diagrams shown in parts (b) and (d) of Fig. 5 indicate that in both cases the rates are quite similar.

When $C=-2$ we observe in Fig. 6 that the geometry of the plane $(x,y)$ changes again drastically. In both, the classical and the PN cases, stability islands corresponding to bounded regular motion around both primaries reappear, while at the boundaries of the bounded basins we observe the presence of starting conditions corresponding to chaotic orbits. It deserves mentioning, that when $r_{S}=0.01$ the amount of trapped chaotic starting conditions is lower, with respect to the classical RTBP. Moreover, the close encounter basins are more prominent in the PN case. The diagrams in parts (b) and (b) of Fig. 6 with the distributions of the close encounter and escape time of orbits suggest that in both cases the rates are fairly similar, however in the case of the PN problem the orbits seem to escape a little bit faster with respect to the classical problem.

At this point, we would like to explain an issue regarding the geometry of the different basin types on the configuration space $(x,y)$. In Zotos et al. (2018) the colour-coded diagrams had several symmetries (for example the area of the close encounter basins to primary $P_{1}$ was equal to the area of the close encounter basins to primary $P_{2}$), which were manifested by the use of polar coordinates. However, in the present study, the respective basins diagrams are completely asymmetric because in this case, we cannot use again polar coordinates. This is because the effective potential (4) is also a function of the velocities $\dot{x}$ and $\dot{y}$.

3.4 An overview analysis

Previously, we presented the nature of orbits on the $(x,y)$ plane for specific values of $C$. Now we want to obtain more information on how the Jacobi constant affects the character of motion. For this purpose, we present in Fig. 7(a-d) numerical outcomes on the $(x,C)$ space, for four characteristic values of $r_{S}$, when all orbits have $y_{0}=\dot{x_{0}}=0$, while $\dot{y_{0}}$ is always calculated through the Jacobi integral (10). The corresponding diagrams showing the distributions of the close encounter and escape time of the orbits are given in Fig. 8(a-d).

As the value of $r_{S}$ increases the most notable changes on the $(x,C)$ plane are:

• •

The areas of bounded regular and chaotic motion in both IR and ER are reduced.

• •

The extent of the close encounter basins increases.

• •

The fractal-like regions on the $(x,C)$ plane are reduced. More precisely, by computing the uncertainty dimension $D_{0}$ (see e.g., Aguirre et al., 2001, 2009), we found that the degree of fractality of the $(x,C)$ plane is indeed reduced. Our computations suggest that when $r_{S}=0$ $D_{0}=1.4123$, when $r_{S}=10^{-4}$ $D_{0}=1.4082$, when $r_{S}=10^{-3}$ $D_{0}=1.3871$, while when $r_{S}=10^{-2}$ $D_{0}=1.3107$.

• •

The time spent by the orbits in a close encounter seems to be reduced, while the escape time of the orbits seems to be unaffected by the shift of $r_{S}$.

From the diagrams shown in Fig. 7 we can extract additional information, regarding the influence of the relativistic effects (i.e., the Schwarzschild radius) on the character of the orbits. In Fig. 9 we present, as a function of $C$, the parametric evolution of the rates of all types of orbits. It is seen, that the most interesting energy range is for $-2<C<4$. In particular, we observe that in this range with an increasing value of $r_{S}$ the percentage of orbits corresponding to close encounters increases, whereas the number of escaping orbits is reduced. Moreover, we see that for extremely low $(C<-5)$ or high $(C>5)$ escaping motion completely dominates, despite of the particular value of $r_{S}$.

In the same vein, we present, as a function of $C$, in Fig. 10 the evolution of the average escape and close encounter times of the orbits, with starting conditions on the $(x,C)$ plane. Note that in both diagrams the scale of the times is logarithmic. One can see, that there are specific energy intervals where there is a clear hierarchy, regarding the influence of $r_{S}$. More specifically, for $-0.5<C<3.5$ $<\log_{10}\left(t_{\rm esc}\right)>$ is reduced, with increasing values of $r_{S}$. Similarly, for $1<C<4$ $<\log_{10}\left(t_{\rm ce}\right)>$ is also reduced with increasing value of $r_{S}$, while for $C>4.2$ the trend is reversed. For the rest of the energy ranges the influence of $r_{S}$ is not very clear and we cannot deduce safe conclusions.

Our numerical analysis ends with Fig. 11 in which we present how $r_{S}$ influences the areas of the main basin types (bounded, escape and close encounter). One can see, that in all cases the most extended basin corresponds to escaping motion, while the area of this type of basin is practically unperturbed by the shift of $r_{S}$. Moreover, the area of the bounded ordered orbits decreases for increasing value of $r_{S}$. On the other hand, for the close encounter motion, the corresponding numerical value of the area on the $(x,C)$ planes displays a linear increase, in a log-log scale. This means that we have a power law between the value of the area of close encounter orbits and the Schwarzschild radius. A similar behaviour (a power law) has also been observed in Nagler (2004) and Nagler (2005), for the classical RTBP, and also in Zotos et al. (2018), for the Paczyński-Wiita potential.

In Fig. 11 a power law $A_{\text{ce}}=r_{S}^{\gamma}$ (where $A_{\text{ce}}$ is the area corresponding to close encounter orbits), with $\gamma=0.296$ is found. The power law behaviour can be explained by the following scaling arguments.

Using Kepler’s problem we can roughly approximate the RTBP in the case where the third body is sufficiently close to one of the primaries, just before a close encounter occurs. In addition, all rotation terms can be neglected. Assuming $r_{S}\ll r_{1}\ll r_{2}$, the problem reduces to Einstein’s computation of Mercury’s orbit, with an additional $1/r_{1}^{3}$ term that determines the precession rate of the Keplerian ellipse in the radial potential. For $r_{1}\ll 1$, we have $1/r_{1}^{3}\gg 1/r_{1}$, such that the relativistic correction term dominates and according to Einstein (1915) one can approximate the effective potential as

[TABLE]

with $L^{2}$ the square of the angular momentum, which is constant.

The energy $E$ of the (precessing) Keplerian ellipse is given by

[TABLE]

where $r_{p}$ denotes the pericenter and $r_{a}$ the apocenter distance.

Following the arguments used in Nagler (2004, 2005), we solve Eq. (12) for $r_{a}$ that yields

[TABLE]

where $E$ is the total energy.

In general, a close encounter occurs when the third body enters the circle of radius $R:=r_{1}$ around the singularity, $r_{p}\leq R$. Hence, for $R\ll 1$, $A_{\text{ce}}\approx 2\pi r_{a}(0)\left(r_{a}(0)-r_{a}(R)\right)$. Assuming that $r_{p}\ll 1$, then Eq. (13) displays, as a function of $r_{p}$, an approximated square root behaviour

[TABLE]

Since we define close encounter through the condition $R=r_{S}$ the close encounter area $A_{\text{ce}}$, as a function of $r_{s}$, is given by

[TABLE]

with the exponent $\beta=1/2$. This approximation explains qualitatively the power law behavior, observed in Fig. 11 but deviates to some extent from the numerical value of the exponent $\gamma=0.296$.

4 Concluding remarks

The scope of the article was a numerical investigation of the type of the orbits of the planar PN Copenhagen problem with Schwarzschild-like primaries, such as super-massive black holes. For describing the dynamical properties of the primaries we used the Tejeda-Rosswog (TR) potential, which reproduces many relativistic effects (e.g., the innermost stable circular orbit, non-circular trajectories, epicyclic frequencies and pericentre precession) with better accuracy, than many other simple potentials, such as the Paczyński & Wiita (PW) and Nowak-Wagoner (NW). By integrating large sets of orbits, we obtained the several basin types on different versions of two-dimensional (2D) planes, by means of colour-coded basin diagrams. Moreover, we explored how the Jacobi constant $C$, as well as the relativistic effects (i.e., the Schwarzschild radius $r_{S}$) influence the nature of the orbits. It was demonstrated, that the both these quantities strongly affect the final states of the orbits, as well as the degree of fractality of the dynamical system.

The following list summarizes the main outcomes of our study:

Regardless of the $r_{S}$ value, escaping orbits predominate for both extremely low and high values of $C$. 2. 2.

In the PN Copenhagen problem, the basins corresponding to starting conditions of close encounter orbits are larger and better defined, while the amount of close encounter orbits grows, with increasing values of $r_{S}$. 3. 3.

The relativistic effects (i.e., the Schwarzschild radius) also affects the bounded regular motion of the system. In particular, the area of the ordered basins corresponding to ordered orbits circulating around one or both primaries seems to reduce, as the value of $r_{S}$ increases. 4. 4.

For $-0.5<C<3.5$, it was found that the time of escape of the orbits is reduced (at least by 15 time units) with increasing value of $r_{S}$. On the contrary, the close encounter time of the orbits either increases (for $C>4.5$, at least by 25 time units) or decreases (for $1<C<4$, at leat by 30 time units) with increasing value of $r_{S}$, depending on the particular energy range. 5. 5.

The fractal degree of the system is reduced with increasing value of $r_{S}$, since all the basin boundaries become smoother and sharper, while at the same time all the noisy (chaotic) regions are confined. The drop of the degree of fractality was also confirmed by the computation of the uncertainty dimension.

In Zotos et al. (2018) we conducted a similar numerical analysis classifying initial conditions of orbits, using the PW potential. On this basis, it would be very interesting to discuss the main similarities and differences of the PW and TR potentials. Here it should be noted, that the PW effective potential is only a function of the coordinates $x$ and $y$, so we were able to use a specialized choice of initial conditions, which express the symmetries of the system. The TR effective potential on the other hand (see Eq. (4)), is also a function of the velocities $\dot{x}$ and $\dot{y}$, which means that we cannot use again the same type of initial conditions. Therefore, we cannot directly relate the numerical outcomes of the basin diagrams on the $(x,y)$ plane. Nevertheless, a comparison can be made by using the numerical information of the classification of orbits on the $(x,C)$ plane.

For the PW and TR potentials we have the following main similarities:

• •

The stability islands of non-escaping regular motion appear at the same locations, i.e., between the two primaries and outside them.

• •

The area of the close encounter basins grow with increasing value of the Schwarzschild radius. Additionally, the increase of the close encounter area follow, for both potentials, a linear pattern (on a log-log scale).

• •

Both the average escape and the close encounter time of the orbits evolve very similarly (displaying similar patterns) for both types of the potentials, as a function of the Jacobi constant.

For the PW and TR potentials we have the following main differences:

• •

In the case of the PW potential the area of the stability islands seems to be unaffected by the shift of the Schwarzschild radius. On the contrary, using the TR potential we observed that the area of the bounded basins decreases, with increasing value of the Schwarzschild radius.

• •

At the left hand side of the $(x,C)$ diagrams and between the stability basins the area of the close encounter basins is larger in the PW potential, with respect to the corresponding diagrams of the TR potential for all values of the Schwarzschild radius.

• •

For the PW potential, we found that for relatively high values of the Schwarzschild radius all three main basin types (bounded, escaping and close encounter) seem to converge to the same rate. In the case of the TR potential, the area of the bounded and close encounter basins do converge, while that of the escaping orbits remains almost unperturbed.

For numerically integrating the grids of starting conditions we used a Bulirsch-Stoer routine in the standard version of FORTRAN 77 (e.g., Press, 1992), with double precision. In all calculations the error, regarding the conservation of the Jacobi constant, was of the order of $10^{-13}$, or even smaller. Using an Intel Quad-Core i7 vPro 4.0 GHz processor, we needed between 3 and 11 hours of CPU time, for the taxonomy of orbits, per grid. All the graphics of the article have been constructed by using the 11.3 version of the Mathematica*®* software (Wolfram, 2003).

Acknowledgments

The authors would like to express their warmest thanks to the anonymous referee for all the apt suggestions and comments which improved both the quality as well as the clarity of the paper.

Appendix A Detailed presentation of the involved equations

Below we present the analytic expressions of all the coefficients, entering the equations of motion:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

while the velocity dependent functions $f$ and $g$ are given by

[TABLE]

and

[TABLE]

The explicit form of the Jacobi integral of motion reads

[TABLE]

As can be noted from the previous expressions, in the limit $r_{S}\to 0$, the equations of motion (7) and the Jacobian (10) both reduce to their Newtonian counterparts.