# Constraining deviations from spherical symmetry using $\gamma$-metric

Hrishikesh Chakrabarty<EMAIL_ADDRESS>School of Astronomy
and Space Sciences, University of Chinese Academy of Sciences (UCAS), Beijing
China  Yong Tang<EMAIL_ADDRESS>School of Astronomy and Space Sciences,
University of Chinese Academy of Sciences (UCAS), Beijing China  School of
Fundamental Physics and Mathematical Sciences, Hangzhou Institute for Advanced
Study, UCAS, Hangzhou 310024, China  International Center for Theoretical
Physics Asia-Pacific, Beijing/Hangzhou, China  National Astronomical
Observatories, Chinese Academy of Sciences, Beijing, China

###### Abstract

The $\gamma$-spacetime metric is a static and axially symmetric vacuum
solution of the Einstein equation. This spacetime represents a naked
singularity and it has an extra parameter $\gamma$ which signifies deviations
from spherical symmetry. In this work, we study the possibility of
constraining the deformation parameter with astrophysical observations. We
start with gravitational lensing in the weak and strong-field limits and
calculate the respective deflection angles to show that only strong-field
lensing observations will be able to constrain $\gamma$ independently. Later
we study two other classical tests of gravity: Shapiro time delay and
precession of perihelion. We show that, out of these two experiments, the
deformation parameter affects the observables only in perihelion shift.

##  I Introduction

The $\gamma$-metric, also known as Zipoy-Voorhees spacetimegamma1 ; gamma2  is
an exact vacuum solution of Einstein’s field equations. The metric belongs to
the Weyl class of spacetimes and is static, axially symmetric and
asymptotically flat. In Erez-Rosen coordinates erez-rosen  the spacetime is
given by the line element

$\displaystyle{\rm d}s^{2}=$
$\displaystyle-f^{\gamma}dt^{2}+f^{\gamma^{2}-\gamma}g^{1-\gamma^{2}}\left(\frac{dr^{2}}{f}+r^{2}d\theta^{2}\right)$
(1) $\displaystyle+f^{1-\gamma}r^{2}\sin^{2}\theta d\phi^{2},$

where

$\displaystyle f=1-\frac{2M}{r},$ (2) $\displaystyle
g=1-\frac{2M}{r}+\frac{M^{2}\sin^{2}\theta}{r^{2}}.$

There are two parameters that characterize the spacetime: $M>0$ is related to
the mass of source and $\gamma>0$ quantifies deformation from spherical
symmetry as for $\gamma>1$ ($\gamma<1$), the spacetime is oblate (prolate).
The spacetime is spherically symmetric for $\gamma=1$, and it corresponds to
Schwarzschild solution. The total Arnowitt-Deser-Misner (ADM) mass measured by
an observer at infinity is $M_{\rm ADM}=\gamma M$.

The metric contains a genuine curvature singularity at $r=2M$ for $\gamma\neq
1$, meaning the spacetime is a naked singularity Virbhadra:1996cz ;
Papadopoulos:1981wr . However, this surface at $r=2M$ must be regarded as an
infinitely redshifted one which may show features similar to the Schwarzschild
event horizon observationally Abdikamalov:2019ztb . For this reason, the
$\gamma$-metric can be considered a “black hole mimicker.” It was shown in
Chakrabarty:2017ysw  that the singularity at $r=2M$ can be resolved in
theories of conformal gravity leaving only the infinitely redshifted surface.

Geodesic motion of massive and massless particles in this spacetime was
studied in Herrera:1998eq ; Chowdhury:2011aa ; Boshkayev:2015jaa ;
Toshmatov:2019bda .The $\gamma$-metric has also been used to describe the
exterior of a general relativistic disk, for details see Gonzalez:2003fe .
Recently, oscillation of neutrinos and its lensing was studied in
$\gamma$-metric and it was shown that constraints on $\gamma$ can be obtained
from the detection of extra-solar neutrinos. There have been studies on
optical properties and shadows of $\gamma$-metric which verified the “black
hole mimicking” property of the spacetime Abdikamalov:2019ztb ; Shaikh:2021cvl
. However, the $\gamma$-metric is non-integrable in general, as shown in
Lukes-Gerakopoulos:2012qpc  and leads to interesting chaotic behavior for the
motion of test particles. The fact that $\gamma$-metric mimics a black hole
and the particle motion scenario is very much similar to that in Schwarzschild
spacetime enables it as a well-motivated and simple candidate to study toy
models of astrophysical scenarios where the exterior of a compact massive
compact object is not given by a usual black hole line element but by the
deformed metric.

Despite this interesting motive, there have not been many efforts to study
different types of experiments in $\gamma$-spacetime. In this work, we discuss
the possibility of constraining the $\gamma$ parameter with astrophysical
observations in the weak and strong field limits. We start with the
calculation of the deflection angle in $\gamma$-metric for gravitational
lensing in the strong and weak-field limit. The weak lensing deflection angle
can be directly applied to solar-system experiments. However, we show that the
parameter $\gamma$ cannot be independently constrained with weak lensing
observations. On the other hand, the strong lensing observations would be able
to constrain $\gamma$ with future very long baseline interferometry (VLBI)
projects. Finally, we shift our attention to other two solar-system
experiments: radar echo delay and precession of perihelia.

The rest of the paper is organized as follows: In Sec. II, we study the strong
and weak field limit of gravitational lensing in $\gamma$-metric. Then in sec.
III and Sec. IV we investigate the Shapiro time-delay and precession of
perihelia in $\gamma$-metric. Finally, we summarize and discuss our findings
in Sec. V. Throughout the paper we use natural units setting $G=c=1$.

##  II Gravitational lensing

In this section, we would calculate the angle of deflection of light in
$\gamma$-metric and see if any observational constraints can be obtained.
First, we shall work in the weak field to find an analytical expression for
the deflection angle. Later, we extend our calculations to the strong field
limit using a method first described in Bozza:2002zj ; Bozza:2010xqn .

###  II.1 Weak-field limit

The world lines of light rays in a curved spacetime can be described by giving
the coordinates $x^{\alpha}$ as functions of any one of a family of affine
parameters $\lambda$. Now, the null vector $v^{\alpha}$ is a tangent to the
world line given by

$v^{\alpha}\equiv\frac{dx^{\alpha}}{d\lambda}.$ (3)

The $\gamma$-metric is independent of $t$ and $\phi$ coordinates and hence we
shall have two conserved quantities along the light ray trajectories,

$\displaystyle e$
$\displaystyle\equiv-\bm{\varepsilon}.\mathbf{v}=f^{\gamma}\left(\frac{dt}{d\lambda}\right)\Rightarrow\left(\frac{dt}{d\lambda}\right)=\frac{e}{f^{\gamma}},$
(4) $\displaystyle l$
$\displaystyle\equiv\bm{\eta}.\mathbf{v}=r^{2}f^{1-\gamma}\left(\frac{d\phi}{d\lambda}\right)\Rightarrow\left(\frac{d\phi}{d\lambda}\right)=\frac{l}{r^{2}f^{1-\gamma}},$

where $\eta^{\alpha}=(0,0,0,1)$ and $\varepsilon^{\alpha}=(1,0,0,0)$ are
Killing vectors of the field. A third integral can be written considering the
requirement that the tangent vector be null

$\mathbf{v}.\mathbf{v}=g_{\alpha\beta}\frac{dx^{\alpha}}{d\lambda}\frac{dx^{\beta}}{d\lambda}=0.$
(5)

Now, on the equatorial plane ($\theta=\pi/2$), the previous equation can be
written explicitly as

$f^{\gamma}\left(\frac{dt}{d\lambda}\right)^{2}-\frac{g^{1-\gamma^{2}}}{f^{1+\gamma-\gamma^{2}}}\left(\frac{dr}{d\lambda}\right)^{2}-f^{1-\gamma}r^{2}\left(\frac{d\phi}{d\lambda}\right)^{2}=0.$
(6)

From Eq. (4), we replace $(dt/d\lambda)$ and $(d\phi/d\lambda)$ to obtain

$\frac{e^{2}}{f^{\gamma}}-\frac{g^{1-\gamma^{2}}}{f^{1+\gamma-\gamma^{2}}}\left(\frac{dr}{d\lambda}\right)^{2}-f^{\gamma-1}\left(\frac{l^{2}}{r^{2}}\right)=0.$
(7)

Multiplying by $(f^{\gamma}/l^{2})$, it gives

$\frac{1}{u^{2}}=\frac{1}{l^{2}}\left[f^{\gamma^{2}-1}g^{1-\gamma^{2}}\right]\left(\frac{dr}{d\lambda}\right)^{2}+W_{\rm
eff}(r),$ (8)

where $u=l/e$ is the impact parameter and

$W_{\rm eff}(r)=\frac{1}{r^{2}}f^{2\gamma-1}$ (9)

is the effective potential felt by the photons. From Eq. (4), we know

$\left(\frac{d\phi}{d\lambda}\right)^{2}=\left(\frac{l}{r^{2}f^{1-\gamma}}\right)^{2}.$
(10)

Dividing Eq. (10) by Eq. (8), we get

$\left(\frac{d\phi}{dr}\right)^{2}=\frac{f^{\gamma^{2}-1}g^{1-\gamma^{2}}}{r^{4}f^{2(1-\gamma)}\left(\frac{1}{u^{2}}-\frac{f^{2\gamma-1}}{r^{2}}\right)}=\mathcal{Z}.$
(11)

This equation gives us the change in azimuthal angle with respect to the
radial coordinate and it can be integrated to obtain the deflection angle
$\Delta\phi$. The magnitude of the total angel swept out as the light ray
proceeds in from infinity and back out again. $\Delta\phi$ is just twice the
angle swept out from a turning point $r=r_{1}$ to infinity. Therefore,

$\Delta\phi=2\int d\phi=2\int^{\infty}_{r_{1}}\sqrt{\mathcal{Z}}dr,$ (12)

and the turning point $r_{1}$ is the radius where $1/u^{2}=W_{\rm
eff}(r_{1})$. The integral here is very complicated and cannot be solved
analytically. Hence, we would consider a weak-field limit where $M/r<<1$. To
obtain an analytical solution, we introduce a new variable $w$, such that

$r=\frac{u}{w}\Rightarrow dr=-\frac{u}{w^{2}}dw.$ (13)

With the new variable, the functions of the metric coefficients in the weak-
field limit,

$f\simeq g\simeq\left(1-\frac{2M}{u}w\right)$ (14)

Now using Eq. (13) and Eq. (14) in the integral, we obtain

$\delta\phi=2\int^{w_{1}}_{0}\frac{1+\frac{M}{u}w}{\left[1+(2\gamma-1)\frac{2M}{u}w-w^{2}\right]^{1/2}}dw,$
(15)

where the limit $w_{1}$ is the value of $w$ at which the denominator vanishes.
This integral can be solved analytically using a standard algebraic
manipulation software. The result is

$\displaystyle\Delta\phi=$ $\displaystyle
2\Bigg{[}-a_{1}\sqrt{1+a_{1}w-w^{2}}$ (16)
$\displaystyle-\arctan\left(\frac{a_{2}-2w}{2\sqrt{1+a_{1}w-w^{2}}}\right)\Bigg{]}_{0}^{w_{1}},$

where $a_{1}=M/u$ and $a_{2}=(2\gamma-1)(2M/u)$. Taking the limit, we get

$\Delta\phi=\pi+\frac{4\gamma M}{u}$ (17)

The deflection angle $\delta\phi_{\rm def}$ is related to $\Delta\phi$ by the
following relation

$\delta\phi_{\rm def}=\Delta\phi-\pi.$ (18)

Thus, the deflection angle in the weak field limit is

$\delta\phi_{\rm def}^{\gamma}=\frac{4\gamma M}{u}=\frac{4M_{\rm ADM}}{u}.$
(19)

Now, we can see from the deflection angle that even if the expression involves
the parameter $\gamma$, it is unlikely to constrain it by observing
gravitational lensing in the weak-field limit. The gravitational mass of the
system is $M_{\rm ADM}=\gamma M$ and the observations will only provide
information about it, unable to break the degeneracy. In the following
section, we shall see if it is possible to break this degeneracy in the strong
field limit.

###  II.2 Strong field limit: Bozza’s method

Figure 1: Left panel: Coefficients of the deflection angle (Eq. (28)) $\bar{a},\bar{b}$ and $u_{m}$ with respect to $\gamma$ evaluated at $u-u_{m}=0.003$. Right panel: The deflection angle $\alpha(\theta)$ with respect to $\gamma$ evaluated at $u-u_{m}=0.003$.  Parameters |  $\gamma=$ 0.95 |  $\gamma=$ 0.97 |  $\gamma=$ 1.0 (Schwarzschild) |  $\gamma=$ 1.02 |  $\gamma=$ 1.05
---|---|---|---|---|---
$\theta_{\infty}(\mu{\rm arcsec})$ |  16.81 |  16.85 |  16.9 |  16.94 |  16.98
$s(\mu{\rm arcsec})$ |  0.02303 |  0.0221 |  0.0211 |  0.02057 |  0.01984
$p_{m}$ (magnitude) |  6.714 |  6.758 |  6.821 |  6.859 |  6.912
$u_{m}/R_{s}$ |  2.58 |  2.589 |  2.598 |  2.603 |  2.609
$\bar{a}$ |  1.0159 |  1.0091 |  1.0 |  0.9944 |  0.9869
$\bar{b}$ |  -0.4148 |  -0.4097 |  -0.4002 |  -0.3927 |  -0.3801
Table 1: Estimates of the observables and the coefficients $\bar{a},\bar{b}$
and $u_{m}$ of strong lensing for the black hole at the center of our galaxy
considering $\gamma$-metric as the spacetime geometry of the exterior.
$\theta_{\infty},s$ and $p_{m}(p)=2.5\log p$ are defined in the text. Here
$R_{s}=\frac{2GM}{c^{2}}$ is the Schwarzschild radius.

In this subsection, we attempt to obtain an analytical expression for
gravitational lensing in the strong field limit in $\gamma$-metric in the
equatorial plane. This analytical procedure was developed in Bozza:2001xd  for
Schwarzschild spacetime and later extended to general cases in Bozza:2010xqn .
Other works on gravitational lensing of black holes and naked singularities
can be found in Claudel:2000yi -Virbhadra:2008ws . We shall first describe the
method briefly and then apply it to $\gamma$-metric. For convenience, we shall
keep the notation used in Bozza:2010xqn .

Let us first express the metric in the equatorial plane as

$ds^{2}=-A(r)dt^{2}+B(r)dr^{2}+C(r)d\phi^{2},$ (20)

where $A(r)$, $B(r)$ and $C(r)$ are metric coefficients. To check the
existence of photon sphere around a compact object described by
$\gamma$-metric, we require that the equation Virbhadra:2002ju ;
Claudel:2000yi

$\frac{C^{\prime}(r)}{C(r)}=\frac{A^{\prime}(r)}{A(r)}$ (21)

admits at least one positive solution and the largest root of this equation
will be called photon sphere, $r_{m}$. Also note that $A,B,C,A^{\prime}$ and
$C^{\prime}$ must be positive for $r>r_{m}$.

The strong field expansion takes the photon sphere as the starting point. For
$\gamma$-metric, the photon sphere is at $r_{m}=(2\gamma+1)M$
Abdikamalov:2019ztb . Note that, when $\gamma\neq 1$, the coordinates are not
spherical and hence photon sphere will not be a sphere. However, in the
equatorial plane, this issue does not make any difference and the strong
lensing method can be applied to $\gamma$-spacetime.

A photon coming from infinity with an impact parameter $u$ will be deviated
while approaching the compact object, reaching a minimum distance $r_{0}$ and
then emerging in another direction. The approach phase will be symmetrical to
the departure phase with the time reversed. Now, using the conservation of
angular momentum, the closest approach distance can be related to the impact
parameter $u$ by the equation

$u=\sqrt{\frac{C_{0}}{A_{0}}},$ (22)

where the subscript $0$ indicated that the functions are evaluated at $r_{0}$.
Similar to the previous section, the geodesic equation can be easily used to
extract the quantity

$\frac{d\phi}{dr}=\frac{\sqrt{B}}{\sqrt{C}\sqrt{\frac{C}{C_{0}}}\frac{A_{0}}{A}-1},$
(23)

which gives the angular shift of the photon as a function of the radial
coordinate. Now, in the strong field limit, the deflection angle can be
calculated as a function of the closest approach Virbhadra:1998dy

$\alpha(r_{0})=I(r_{0})-\pi,$ (24)

where

$I(r_{0})=\int_{r_{0}}^{\infty}\frac{\sqrt{B}}{\sqrt{C}\sqrt{\frac{C}{C_{0}}\frac{A_{0}}{A}-1}}dr.$
(25)

We can see that for a vanishing gravitational field, $\alpha(r_{0})$ vanishes
identically. In the weak-field limit, the integrand can be expanded to the
first order in the gravitational potential reproducing Eq. (15). Decreasing
the impact parameter increases the deflection angle and at some point, it will
exceed $2\pi$ resulting in a loop around the compact object spacetime.
$r_{0}=r_{m}$ corresponds to the impact parameter $u=u_{m}$, and at this point
the deflection angle diverges meaning the photon gets captured.

The divergence here is logarithmic Bozza:2001xd ; Bozza:2010xqn . The
analytical expansion for the deflection angle close to the divergence has the
form

$\alpha(r_{0})=-a\log\left(\frac{r_{0}}{r_{m}}-1\right)+b+\mathcal{O}(r_{0}-r_{m}),$
(26)

where the coefficients depend on the metric functions evaluated at $r_{m}$.

The deflection angle, $\alpha(r_{0})$ can be expressed as a function of the
angular separation of the image from the lens, $\theta$. Angular separation of
the lens from the image is defined as $\theta={u}/{D_{OL}}$, where $D_{OL}$ is
the distance between the lens and the observer. In terms of this variable the
deflection angle is expressed as

$\alpha(\theta)=-\bar{a}\log\left(\frac{\theta
D_{OL}}{u_{m}}-1\right)+\bar{b}+\mathcal{O}(u-u_{m}).$ (27)

The coefficients $\bar{a}$ and $\bar{b}$ are as follows

$\displaystyle\bar{a}$
$\displaystyle=\frac{a}{2}=\frac{R(0,r_{m})}{2\sqrt{\beta_{m}}},$ (28)
$\displaystyle\bar{b}$
$\displaystyle=-\pi+b_{R}+\bar{a}\log\left(\frac{2\beta_{m}}{y_{m}}\right).$

Here, $R(0,r_{m})$, $\beta_{m}$, $y_{m}$ are again functions of metric
coefficients evaluated at $r_{0}=r_{m}$. The derivation of the deflection
angle along with the forms of the coefficients in the strong field limit is
rather long and is omitted in the main text. For reference we have added in
Appendix A the derivation, which alternatively can also be found in
Bozza:2010xqn . Below, we apply this procedure to $\gamma$-metric.

Our task is to find the coefficients $\bar{a}$, $\bar{b}$ and $u_{m}$ to check
how the deflection angle varies with respect to the parameter $\gamma$. We
rewrite the metric coefficients again as $A,B$ and $C$ for $\gamma$-metric

$\displaystyle A(r)=\left(1-\frac{2M}{r}\right)^{\gamma},$ (29) $\displaystyle
B(r)=\left(1-\frac{2M}{r}\right)^{\gamma^{2}-\gamma-1}\left(1-\frac{2M}{r}+\frac{M^{2}}{r^{2}}\right)^{1-\gamma^{2}},$
$\displaystyle C(r)=\left(1-\frac{2M}{r}\right)^{1-\gamma}r^{2}.$

We can calculate the radius of the photon sphere by solving the equation
$\alpha=0$ (first equation of (67) of Appendix A). This yields

$r_{m}=(2\gamma+1)M.$ (30)

Now, $\beta$ (Eq. (72)) at $r_{0}=r_{m}$ is

$\beta_{m}=\frac{(2\gamma-1)(2\gamma+1)\left[\left(\frac{2\gamma-1}{2\gamma+1}\right)^{\gamma}-1\right]^{2}}{4\gamma^{2}\left(\frac{2\gamma-1}{2\gamma+1}\right)^{\gamma}}.$
(31)

The integral Eq. (77) cannot be solved exactly to find $b_{R}$. Therefore we
expand the integrand in powers of $(\gamma-1)$ and evaluate the single
coefficients. We get

$b_{R}=b_{R,0}+b_{R,1}(\gamma-1)+\mathcal{O}(\gamma-1)^{2}.$ (32)

Here $b_{R,0}$ is the value of the coefficient for Schwarzschild spacetime
Bozza:2010xqn

$b_{R,0}=0.9496.$ (33)

$b_{R,1}$ represents the correction term for $\gamma$-metric. [see Appendix
for expression].

Finally, we can compute the coefficients of lensing in the strong field limit.
For $M_{\rm ADM}=\gamma M=1$,

$\displaystyle\bar{a}$
$\displaystyle=\sqrt{\left(\frac{4\gamma^{2}}{4\gamma^{2}-1}\right)^{1-\gamma^{2}}},$
(34) $\displaystyle\bar{b}$
$\displaystyle=-\pi+b_{R}+\sqrt{\left(\frac{4\gamma^{2}}{4\gamma^{2}-1}\right)^{1-\gamma^{2}}}\log\left(\left(2-\frac{1}{2\gamma^{2}}\right)\left[1-\left(\frac{2\gamma-1}{2\gamma+1}\right)^{-\gamma}\right]^{2}\right),$
$\displaystyle u_{m}$
$\displaystyle=\left(2+\frac{1}{\gamma}\right)\left(\frac{2\gamma-1}{2\gamma+1}\right)^{\frac{1-2\gamma}{2}}.$

These parameters are functions of $\gamma$ only, since we have fixed the
observable mass. Now we can calculate the dependence of the deflection angle
on the parameter $\gamma$. In Schwarzschild spacetime, it was shown in
Bozza:2010xqn  that the most external image appears where $\alpha(\theta)$
falls below $2\pi$ and this happens at the location $u-u_{m}=0.003264$. To see
the dependence of the strong lensing parameters and deflection angle on
$\gamma$, we evaluate them at $u=u_{m}+0.003$ and plot the result in Fig. 1.

We can see from the left panel of Fig. 1 that $u_{m}$ increases while
$\bar{a}$ decreases with increasing $\gamma$. $\bar{b}$ on the other hand
initially decreases with increasing $\gamma$ reaching a minimum value at
$\gamma=0.89$ and later again increases. The panel shows the decrease of
deflection angle $\alpha(\theta)$ with increasing $\gamma$. Both these plots
show an asymptotic behaviour when we approach $\gamma\rightarrow 0.5$. This is
precisely because of the disappearance of the photon sphere at $\gamma=0.5$ as
the photon sphere radius becomes $r_{m}=2M$ and it coincides with the
singularity which is also an infinitely redshifted surface at the equatorial
plane.

Let us now have a look at one realistic numerical example. We have significant
evidence that the center of our galaxy hosts a supermassive black hole with
mass $M_{\rm BH}\simeq 4\times 10^{6}{\rm M_{\odot}}$. Gravitational lensing
by such a black hole was discussed in detail by Virbhadra and Ellis in
Virbhadra:1999nm  for the Schwarzschild case and Bozza in Bozza:2002zj
comparing modified spacetimes. Considering a distance of $D_{OL}=8.5{\rm kpc}$
between the black hole and the sun, they showed that the separation between
each set of relativistic images with respect to the central lens would be
$\theta_{\infty}=17\mu{\rm arcsecs}$. Now we are on the verge of achieving
such a resolution through actual VLBI projects such as Event Horizon Telescope
EventHorizonTelescope:2019uob . Here we follow the same procedure Bozza:2002zj
; Virbhadra:1999nm  for $\gamma$-metric and estimate the quantities required
for a complete strong field limit reconstruction.

The observables for gravitational lensing are discussed in Virbhadra:1999nm ;
Bozza:2002zj ; Virbhadra:2002ju . Basically we are interested in two
quantities

$\displaystyle s=\theta_{1}-\theta_{\infty},$ (35) $\displaystyle
p=\frac{\mu_{1}}{\sum_{n=2}^{\infty}\mu_{n}}.$

Here, $s$ and $p$ represents the separation between the first image and the
others and the ratio between the flux of the first image and that coming from
the others respectively. Clearly, $\theta_{1}$ is the angular separation for
the first image and $\theta_{\infty}$ signify the asymptotic position by the
other set of images. $\mu_{n}$ is the flux, where the subscript represents the
corresponding image. The parameters $s$ and $p$ can be related to the
coefficients of strong lensing as

$\displaystyle
s=\theta_{\infty}\exp^{\frac{\bar{b}}{\bar{a}}-\frac{2\pi}{\bar{a}}}$ (36)
$\displaystyle p=\exp^{\frac{2\pi}{\bar{a}}}.$

So, just by measuring the angular separation and a ratio of flux, one can
reconstruct the full strong field limit expansion of the deflection angle.

In Table 1, we show the estimates of these parameters for different values of
$\gamma$, including $\gamma=1$ which is the Schwarzschild case. Looking at the
table, it is immediately clear that the easiest quantity to evaluate is the
minimum impact parameter $u_{m}$. Once we achieve a microarcsecond resolution
in the next years, it would become possible to distinguish between a
Schwarzschild black hole and a naked singularity represented by
$\gamma$-spacetime metric. On the other hand, to fit all other coefficients
into any compact object model, we need to separate at least the outermost
relativistic image from the others. This can be done only with an increased
optical resolution at least by two orders of magnitude concerning actual
observational projects Bozza:2002zj .

##  III Shapiro time delay

Another interesting relativistic effect in the propagation of light rays is
the apparent delay in propagation time for a light signal passing near the
sun. This is important because radar-ranging techniques can measure this delay
and give constraint on $\gamma$-metric. This effect is called the Shapiro time
delay effect.

The idea is to measure the time required for radar signals to travel to an
inner planet or satellite in two circumstances: a) when the signal passes very
near the Sun and b) when the ray does not go near the Sun. The time required
to travel for light $t_{0}$ between two planets sitting far away from the Sun
is given by

$t_{0}=\int_{-l_{1}}^{l_{2}}dy,$ (37)

where $dy$ is the differential distance in the radial direction in the solar
system and $l_{1}$ and $l_{2}$ are the distances of the the planets from the
Sun. When the radar signal travels close to the sun, the previous equation
should be modified as

$t=\int_{-l_{1}}^{l_{2}}\frac{dy}{v}=\int_{-l_{1}}^{l_{2}}\sqrt{\frac{g_{rr}}{g_{tt}}}dy.$
(38)

Here, $v=\sqrt{g_{rr}/g_{tt}}$ is the speed of light in presence of the
gravitational field, $g_{tt}$ and $g_{rr}$ are the metric components. Now the
time difference is

$\Delta
t=t-t_{0}=\int_{-l_{1}}^{l_{2}}\left(\sqrt{\frac{g_{rr}}{g_{tt}}}-1\right)dy.$
(39)

The radial coordinate can be expressed as $r=\sqrt{R_{\odot}^{2}+y^{2}}$,
where $R_{\odot}$ is the radius of the Sun. Now Eq. (39) becomes an integral
of variable $y$, so we can write

$\Delta
t=\int_{-l_{1}}^{l_{2}}\left(\sqrt{\frac{g_{rr}(\sqrt{R_{\odot}^{2}+y^{2}})}{g_{tt}(\sqrt{R_{\odot}^{2}+y^{2}})}}-1\right)dy.$
(40)

Now $g_{rr}$ and $g_{tt}$ can be replaced in the above integral

$\Delta
t=\int_{-l_{1}}^{l_{2}}\left(\sqrt{f^{\gamma^{2}-2\gamma-1}g^{1-\gamma^{2}}}-1\right)dy$
(41)

Explicitly, in the equatorial plane

$\Delta
t=\bigints_{-l_{1}}^{l_{2}}\left(\sqrt{\left(1-\frac{2M}{\sqrt{R_{\odot}^{2}+y^{2}}}\right)^{\gamma^{2}-2\gamma-1}\left(1-\frac{2M}{\sqrt{R_{\odot}^{2}+y^{2}}}+\frac{M^{2}}{R_{\odot}^{2}+y^{2}}\right)^{1-\gamma^{2}}}-1\right)dy.$
(42)

Now the integration can be performed considering a Earth-Sun-Mars system.
Relevant quantities here are the distance to the Earth from the Sun $R_{\rm
E}=l_{1}=1.525\times 10^{13}{\rm cm}$, the distance to the Mars from the Sun
$R_{\rm M}=l_{2}=2.491\times 10^{13}{\rm cm}$, mass of the Sun
$M_{\odot}=\gamma M=1.989\times 10^{33}{\rm g}$ and the radius of the Sun
$R_{\odot}=6.955\times 10^{10}{\rm cm}$.

For Schwarzschild metric, with the above mentioned values of the parameters
the radar echo delay has the value

$\Delta t_{RD}^{(Sch)}\simeq 4M_{\odot}\ln(4l_{1}l_{2}/R_{\odot}^{2})\simeq
2.4927\times 10^{-4}{\rm s}$ (43)

In $\gamma$-spacetime, the radar echo delay expression can be written as

$\Delta
t_{RD}^{\gamma}=2\bigints_{-l_{1}}^{l_{2}}\left(\sqrt{\left(1-\frac{2(M_{\odot}/\gamma)}{\sqrt{R_{\odot}^{2}+y^{2}}}\right)^{\gamma^{2}-2\gamma-1}\left(1-\frac{2(M_{\odot}/\gamma)}{\sqrt{R_{\odot}^{2}+y^{2}}}+\frac{(M_{\odot}/\gamma)^{2}}{R_{\odot}^{2}+y^{2}}\right)^{1-\gamma^{2}}}-1\right)dy$
(44)

First let us check the weak field limit of the above integral. We expand the
expression under the square root in terms of $M/R_{\odot}$ and keep only the
terms of order $\mathcal{O}(M/R_{\odot})$. The integral becomes

$\displaystyle\Delta t_{RD}^{\gamma}$
$\displaystyle=2\int_{-l_{1}}^{l_{2}}\frac{2\gamma(M_{\odot}/\gamma)}{R_{\odot}}\left(1-\frac{y^{2}}{R_{\odot}^{2}}\right)^{-1}dy$
(45)
$\displaystyle=2\int_{-l_{1}}^{l_{2}}\frac{2M_{\odot}}{R_{\odot}}\left(1-\frac{y^{2}}{R_{\odot}^{2}}\right)^{-1}dy$
$\displaystyle\simeq
4M_{\odot}\ln\left(\frac{4l_{1}l_{2}}{R_{\odot}^{2}}\right)=\Delta
t_{RD}^{Sch},$

which is nothing but the expression for Shapiro time-delay in Schwarzschild
spacetime. Therefore, the light rays in the time delay experiment will not
experience the effect of $\gamma$ on them.

Let us now verify this in the strong-field regime. We numerically integrate
Eq. (44) without any approximations. The result of the integration is shown in
Fig. 2, which clearly shows that the Shapiro time delay in the $\gamma$ metric
is independent of the parameter $\gamma$.

Figure 2:  Numerical integration of Eq. (44). This Plot shows the time delay
with respect to gamma when a radar signal travels back and forth through path
close to the sun where the exterior of the sun assumed to be described by
$\gamma$ metric.

##  IV Perihelion precession

To find the rate precession of perihelion in $\gamma$-metric, let us now
concentrate on the motion of massive particles. The metric is independent of
$t$ and $\theta$ coordinates and hence we shall again have two timelike
Killing vectors

$\displaystyle e$
$\displaystyle\equiv-\bm{\varepsilon}.\mathbf{v}=f^{\gamma}\left(\frac{dt}{d\tau}\right)\Rightarrow\left(\frac{dt}{d\tau}\right)=\frac{e}{f^{\gamma}},$
(46) $\displaystyle l$
$\displaystyle\equiv\bm{\eta}.\mathbf{v}=r^{2}f^{1-\gamma}\left(\frac{d\phi}{d\tau}\right)\Rightarrow\left(\frac{d\phi}{d\tau}\right)=\frac{l}{r^{2}f^{1-\gamma}}.$

A third constant of motion is given by the considering the requirement that
the tangent vector to be timelike for massive particles

$\mathbf{v}.\mathbf{v}=g_{\alpha\beta}u^{\alpha}u^{\beta}=-1.$ (47)

On the equatorial plane, this equation can be written as

$-f^{\gamma}\dot{t}^{2}+\left(\frac{f^{\gamma^{2}-\gamma}g^{1-\gamma^{2}}}{f}\right)\dot{r}^{2}+f^{1-\gamma}r^{2}\dot{\phi}^{2}=-1,$
(48)

where $\dot{t}=dt/d\tau$, $\dot{\phi}=d\phi/d\tau$ and $\tau$ is the proper
time. Now we can replace $\dot{t}$ and $\dot{\phi}$ using Eq. (46) and the
equation becomes

$-f^{-\gamma}e^{2}+\left(\frac{f^{\gamma^{2}-\gamma}g^{1-\gamma^{2}}}{f}\right)\dot{r}^{2}+\frac{l^{2}}{r^{2}f^{1-\gamma}}=-1.$
(49)

With some algebraic manipulations, this equation can be written as

$\dot{r}^{2}+A(r)^{-1}\frac{l^{2}}{r^{2}f^{1-\gamma}}=A(r)^{-1}\left(f^{-\gamma}e^{2}-1\right),$
(50)

where

$A(r)=\frac{f^{\gamma^{2}-\gamma}g^{1-\gamma^{2}}}{f}.$ (51)

Now we employ a change of variable of the form $r=1/w$ and replace $d/d\tau$
with $(lw^{2}/f^{1-\gamma})d/d\phi$ and obtain

$\left(\frac{dw}{d\phi}\right)^{2}+\frac{f^{1-\gamma}}{A(w)}w^{2}=\frac{f^{2(1-\gamma)}}{l^{2}A(w)}\left(f^{-\gamma}e^{2}-1\right).$
(52)

Arranging the terms

$\left(\frac{dw}{d\phi}\right)^{2}+w^{2}=\frac{e^{2}}{l^{2}}X(w)-\frac{1}{l^{2}}Y(w)-w^{2}Z(w),$
(53)

where

$\displaystyle X(w)$ $\displaystyle=\frac{f^{2-3\gamma}}{A(w)},$ (54)
$\displaystyle Y(w)$ $\displaystyle=\frac{f^{2(1-\gamma)}}{A(w)},$
$\displaystyle Z(w)$ $\displaystyle=\frac{f^{1-\gamma}}{A(w)}-1.$

By taking the derivative of the previous equation with respect to $\phi$, we
find

$\frac{d^{2}w}{d\phi^{2}}+w=F(w),$ (55)

which is the orbit equation. Here,

$F(w)=\frac{1}{2}\frac{dG(w)}{dw},$ (56)

and

$G(w)=w^{2}Z(w)+\frac{e^{2}}{l^{2}}X(w)-\frac{1}{l^{2}}Y(w).$ (57)

Using the method as shown in Harko:2009qr ; Capistrano:2019qdv , we work with
the previous orbit equation. The deviation angle $\delta\phi$ can be found
analytically using

$\delta\phi=\pi\frac{dF(w)}{dw}\Bigg{|}_{w_{0}},$ (58)

along with the constraint $F(w_{0})=w_{0}$ for a near-circular orbit. In the
weak field regime ($Mw<<1$), we can expand $X(w)$, $Y(w)$ and $Z(w)$ in powers
of $w$ to write $G(w)$ and $F(w)$

$\displaystyle G(w)$
$\displaystyle=\frac{w^{2}\left(3\left(3\gamma^{2}-4\gamma+1\right)e^{2}M^{2}-3(\gamma-1)^{2}M^{2}\right)}{l^{2}}$
(59) $\displaystyle+\frac{w\left(4(\gamma-1)e^{2}M-2\gamma
M+4M\right)}{l^{2}}$
$\displaystyle+\frac{e^{2}-1}{l^{2}}-(\gamma-1)(\gamma+1)M^{2}w^{4}+2Mw^{3}$
$\displaystyle F(w)$
$\displaystyle=\frac{M\left(-\gamma+2(\gamma-1)e^{2}+2\right)}{l^{2}}$
$\displaystyle+\frac{Mw\left(3(\gamma-1)(3\gamma-1)e^{2}M-3(\gamma-1)^{2}M\right)}{l^{2}}$
$\displaystyle-2(\gamma-1)(\gamma+1)M^{2}w^{3}+3Mw^{2}$

So the circular orbits will be given by the roots of the equation
$F(w_{0})=w_{0}$. Explicitly

$A_{0}+A_{1}w_{0}+A_{2}w_{0}^{2}+A_{3}w_{0}^{3}=w_{0},$ (60)

where

$\displaystyle A_{0}$
$\displaystyle=\frac{M\left(-\gamma+2(\gamma-1)e^{2}+2\right)}{l^{2}},$ (61)
$\displaystyle A_{1}$
$\displaystyle=\frac{M^{2}\left[3(\gamma-1)(3\gamma-1)e^{2}-3(\gamma-1)^{2}\right]}{l^{2}},$
$\displaystyle A_{2}$ $\displaystyle=3M,$ $\displaystyle A_{3}$
$\displaystyle=-2(\gamma-1)(\gamma+1)M^{2}.$

Neglecting higher order corrections, the solution to the equation
$F(w_{0})=w_{0}$ can be written as

$w_{0}=\frac{A_{0}}{1-A_{1}},$ (62)

which reduces to $w_{0}=M/l^{2}$ for $\gamma=1$. Using the above results, the
perihelion precession $\delta\phi$ can be written as

$\displaystyle\delta\phi$ $\displaystyle=\pi\frac{dF(w)}{dw}\Bigg{|}_{w_{0}}$
(63)
$\displaystyle=\pi\Bigg{[}\frac{M^{2}\left[3(\gamma-1)(3\gamma-1)e^{2}-3(\gamma-1)^{2}\right]}{l^{2}}$
$\displaystyle\ \ \ \ \ \ \ \ \ \
-6(\gamma-1)(\gamma+1)M^{2}w_{0}^{2}+6Mw_{0}\Bigg{]}.$

For $\gamma=1$, we recover the classical General Relativistic result
$\delta\phi_{\rm GR}=6\pi\left(M^{2}/l^{2}\right)$.

For the $\gamma$ metric, we report the variation of perihelion precession
angle of Mercury around the sun (left panel) and S2 star around Sgr A* (right
panel) with respect to the deformation parameter $\gamma$ in Figure 3. In
Table 2, we show the numerical estimates of perihelion precession for Mercury
for values of $\gamma=0.95,0.97,1.0,1.02$ and $1.05$. We can put constraints
on the parameter $\gamma$ in the solar system using the measured values of
perihelion precession for Mercury $\delta\phi_{\rm prec}$. The observed value
of the perihelion shift of Mercury is $42.98\pm 0.04$ arcsec/century
Will:2005yc ; Will:2014kxa . Using this data, $\gamma$ can be constrained at
$1.0000\pm 0.0005$ where we clearly recover the Schwarzschild limit. On the
other hand, for S2, the measured value of the orbital precession per orbit is
GRAVITY:2020gka  $\delta\phi_{\rm prec}=12.1\times(1.10\pm 0.19)$
arcminute/orbit. This constrains $\gamma$ at $0.96\pm 0.08$ for Sgr A* on the
plane of the orbit in the weak field limit. We still recover the Schwarzschild
limit. However, note that this analysis is preliminary. For better
constraints, full numerical simulations with all orbital parameters should be
done. With the increase in precision of the measured value in the coming
years, we would be able to obtain a more stringent constraint on the
deformation parameter. Similarly, an analysis in the strong field regime would
also require full numerical simulations of the orbit equations, which then can
be applied to the motion of S stars to constraint deviation from spherical
symmetry.

$\gamma$ |  $\delta\phi$ (arcsec/century)
---|---
0.95 |  47.45
0.97 |  45.62
1.0 |  42.98
1.02 |  41.28
1.05 |  38.83
Table 2: Estimates of precession of perihelion of Mercury for different values
of $\gamma$.

Figure 3:  Perihelion shift as a function of the deformation parameter for
Mercury (left) and S2 star (right) with respect to $\gamma$.

##  V Outlook and discussion

In the present work, we tested the applicability of the $\gamma$-metric for
astrophysical observations. We studied strong and weak gravitational lensing,
radar echo delay and perihelion precession in the context of the considered
metric spacetime. We demonstrated that most of the weak field experiment
observables do not show any dependence on $\gamma$, except perihelion
precession. However, strong lensing shows a significant effect on light
bending in $\gamma$-metric.

We started with gravitational lensing in the weak field limit and calculated
the deflection angle of photons. The expression of the deflection angle
depends on the deformation parameter $\gamma$ and the parameter related to the
mass $M$. This dependence comes in the form $\gamma M$ which is the observable
mass of the system. Therefore, it is unlikely to independently constrain
$\gamma$ with weak lensing experiments. On the other hand, the mechanism in
strong lensing is quite involved and we can see a clear dependence of the
observational quantities in the deformations parameter. This can be seen in
Fig. 1. The deflection angle in strong lensing increases with decreasing
$\gamma$ and blows up at $\gamma=0.5$. This is because of the disappearance of
the photon sphere at the corresponding value of $\gamma$.

We then concentrated on two classic solar system experiments, namely Shapiro
delay and perihelion precession. We showed that light rays in the time delay
experiment are not affected by the deformation parameter. Finally, we studied
the precession of perihelion in this spacetime and showed that perihelion
shift depends on the deformation parameter, i.e. with increasing $\gamma$ the
shift in perihelion decreased as can be seen in Fig. 3. Using the currently
accepted value of the perihelion shift of Mercury, we obtained numerical
constrain on the deformation parameter which is $\gamma=1.0000\pm 0.0005$.
Similarly, data from S2 star contraints $\gamma$ at $0.96\pm 0.08$. Shapiro
time delay is still a weak field effect in the solar system and it is expected
that the deformation parameter of a modified spacetime will affect the travel
time of a signal if the signal passes close to the massive object. However, in
our analysis of $\gamma$-metric, we can see that the integrand of the time
delay integral does not depend on the non-trivial $g_{\phi\phi}$ term of the
metric and the contribution from $\gamma$ in $g_{tt}$ and $g_{rr}$ seems to
contribute only towards the observable mass of the massive object. On the
other hand, in strong lensing and perihelion precession, the deviation from
the Schwarzschild value comes from the $g_{\phi\phi}$ term.

Finally, concerning the possibility of constraining deviations from spherical
symmetry through observation of gravitational lensing, our results show that
only precise measurements of observables in strong field lensing would allow
us to distinguish a black hole from the $\gamma$-spacetime metric. On the
other hand, more precise observation of the perihelion shift of Mercury and S
stars in the future would allow us to put much tighter bounds on $\gamma$ in
the respective systems.

###### Acknowledgements.

HC thank Prof. Bobomurat Ahmedov, Dr. Daniele Malafarina and Dr. Ahmadjon
Abdujabbarov for useful discussions. HC also acknowledges support from
University of Chinese Academy of Sciences through Special Research Assistant
fellowship (SRA). YT is supported by Natural Science Foundation of China
(NSFC) under Grants No. 11851302, the Fundamental Research Funds for the
Central Universities and Key Research Program of the Chinese Academy of
Sciences, Grant No. XDPB15.

##  Appendix A Deflection angle in the strong field limit: Bozza’s method

To find the coefficients of the deflection angle in Eq. (28), first two new
variables are defined as

$\displaystyle y$ $\displaystyle=A(r)$ (64) $\displaystyle z$
$\displaystyle=\frac{y-y_{0}}{1-y_{0}}$

where $y_{0}=A_{0}=A(r_{0})$. In terms of these new variables the integral in
deflection angle (25) becomes

$\displaystyle I(r_{0})=\int_{0}^{1}R(z,r_{0})f(z,r_{0})dz,$ (65)
$\displaystyle
R(z,r_{0})=\frac{2\sqrt{By}}{CA^{\prime}}(1-y_{0})\sqrt{C_{0}},$
$\displaystyle
f(z,r_{0})=\frac{1}{\sqrt{y_{0}-[(1-y_{0})z+y_{0}]\frac{C_{0}}{C}}}$

where all functions without the subscript $0$ are evaluated at
$x=A^{-1}[(1-y_{0})z+y_{0}]$.

Here the function $R(z,r_{0})$ is regular for all values of $z$ and $r_{0}$
but $f(z,r_{0})$ diverges for $z\rightarrow 0$. The argument under the square
root in $f(z,r_{0})$ is expanded to the second order in $z$,

$\displaystyle f(z,r_{0})\sim f_{0}(z,r_{0})=\frac{1}{\sqrt{\alpha z+\beta
z^{2}}},$ (66)

where,

$\displaystyle\alpha$
$\displaystyle=\frac{1-y_{0}}{C_{0}A_{0}^{\prime}}(C_{0}^{\prime}y_{0}-C_{0}A_{0}^{\prime}),$
(67) $\displaystyle\beta$
$\displaystyle=\frac{(1-y_{0})^{2}}{2C_{0}^{2}A_{0}^{\prime
3}}\Big{[}2C_{0}C_{0}^{\prime}A_{0}^{\prime
2}+\left(C_{0}C_{0}^{\prime\prime}-2C_{0}^{\prime
2}\right)y_{0}A_{0}^{\prime}$ $\displaystyle-
C_{0}C_{0}^{\prime}y_{0}A_{0}^{\prime\prime}\Big{]}.$

In case of a non-zero $\alpha$, the leading order divergence is $z^{1/2}$,
which can be integrated to give a finite result. The integral diverges when
$\alpha=0$, since the leading order divergence is $z^{-1}$. Analyzing the form
of $\alpha$, we can see that it vanishes for $r_{0}=r_{m}$, resulting in
capture of each photon having $r_{0}<r_{m}$.

Now to solve the integral, it is split into two parts, $I_{D}$ (contains the
divergence) and $I_{R}$ (regular):

$\displaystyle I(r_{0})$ $\displaystyle=I_{D}(r_{0})+I_{R}(r_{0})$ (68)
$\displaystyle I_{D}(r_{0})$
$\displaystyle=\int_{0}^{1}R(0,r_{m})f_{0}(z,r_{0})dz,$ $\displaystyle
I_{R}(r_{0})$ $\displaystyle=\int_{0}^{1}g(z,r_{0})dz,$

where

$g(z,r_{0})=R(z,r_{0})f(z,r_{0})-R(0,r_{m})f_{0}(z,r_{0}).$ (69)

We can see that $I_{R}(r_{0})$ is the original integral with the divergence
subtracted. Now both the integrals should be solved separately and then summed
up rebuilding the deflection angle.

Let us first handle the integral $I_{D}(r_{0})$. It can be solved exactly. The
results is

$I_{D}(r_{0})=R(0,r_{m})\frac{2}{\sqrt{\beta}}\log\frac{\sqrt{\beta}+\sqrt{\alpha+\beta}}{\sqrt{\alpha}}.$
(70)

We are interested only in the terms up to $\mathcal{O}(r_{0}-r_{m})$, so we
expand $\alpha$ as

$\alpha=\frac{2\beta_{m}A_{m}^{\prime}}{1-y_{m}}(r_{0}-r_{m})+\mathcal{O}(r_{0}-r_{m})^{2},$
(71)

where

$\beta_{m}=\frac{C_{m}(1-y_{m})^{2}(C_{m}^{\prime\prime}y_{m}-C_{m}A_{m}^{\prime\prime})}{2y_{m}^{2}C_{m}^{\prime\prime}}.$
(72)

We then substitute it in $I_{D}(r_{0})$

$I_{D}(r_{0})=-a\log\left(\frac{r_{0}}{r_{m}}-1\right)+b_{D}+\mathcal{O}(r_{0}-r_{m})$
(73)

where

$\displaystyle a=\frac{R(0,r_{m})}{\sqrt{\beta_{m}}},$ (74) $\displaystyle
b=\frac{R(0,r_{m})}{\sqrt{\beta_{m}}}\log\frac{2(1-y_{m})}{A_{m}^{\prime}r_{m}}.$

$I_{D}(r_{0})$ yields logarithmic leading order divergence of the deflection
angle.

Now let us concentrate on the regular term. To find $b$ we need to add an
analogous term coming from the regular part of the integral to $b_{D}$. First
we expand $I_{R}(r_{0}-r_{m})$

$I_{R}(r_{0})=\sum_{n=0}^{\infty}\frac{1}{n!}(r_{0}-r_{m})^{n}\int_{0}^{1}\frac{\partial^{n}g}{\partial
r_{o}^{n}}\Bigg{|}_{r_{0}=r_{m}}dz$ (75)

and evaluate the single coefficients.

Not subtracing the divergent part would have resulted in an infinite
coefficient for $n=0$, but the other coefficients would be finite. However,
$g(z,r_{0})$ is regular for $z=0$ as $r_{0}\rightarrow r_{m}$. We are only
interested to terms up to $\mathcal{O}(r_{0}-r_{m})$, so just retaining the
leading order term

$I_{R}(r_{0})=\int_{0}^{1}g(z,r_{m})dz+\mathcal{O}(r_{0}-r_{m}).$ (76)

Then we have

$b_{R}=I_{R}(r_{m})$ (77)

which is the term needed to be added to $b_{D}$ to get the regular
coefficient. So finally,

$b=-\pi+b_{D}+b_{R}.$ (78)

The term $b_{R}$ can be evaluated numerically for all metric forms since the
integrand does not contain any divergences. However, for Schwarzschild metric,
it is solved exactly. For $\gamma$-metric we expand the integrand over
$(\gamma-1)$ and try to obtain an analytical expression.

Now to go from $\alpha(r_{0})$ to $\alpha(\theta)$ we expand the equation for
impact parameter

$\displaystyle u_{m}=\sqrt{\frac{C_{m}}{y_{m}}},$ (79) $\displaystyle
u-u_{m}=c(r_{0}-r_{m})^{2},$

where

$c=\frac{C_{m}^{\prime\prime}y_{m}-C_{m}A_{m}^{\prime\prime}}{3\sqrt{y_{m}^{3}C_{m}}}=\beta_{m}\sqrt{\frac{y_{m}}{C_{m}^{3}}}\frac{C_{m}^{\prime
2}}{2(1-y_{m})^{2}}.$ (80)

Using this relation, the deflection angle can be expressed as a function of
$\theta$

$\alpha(\theta)=-\bar{a}\log\left(\frac{\theta
D_{OL}}{u_{m}}-1\right)+\bar{b}$ (81)

where

$\displaystyle\bar{a}=\frac{a}{2}=\frac{R(0,r_{m})}{2\sqrt{\beta_{m}}}$ (82)
$\displaystyle\bar{b}=-\pi+b_{R}+\bar{a}\log\frac{2\beta_{m}}{y_{m}}.$

So finally, in order to calculate the deflection angle as a function of
$\theta$, we need to

  * •

solve Eq. (21) to find $r_{m}$.

  * •

write $\beta_{m}$ and $R(0,r_{m})$ from Eq. (72) and Eq. (65) respectively

  * •

compute $b_{R}$ from Eq. (77) by a an expansion of the parameters of the
metric

  * •

compute $u_{m}$, $\bar{a}$ and $\bar{b}$ from Eq. (79) and (82) respectively.

The only integral involved in the process is calculating $b_{R}$.

## References

  * (1) D. M. Zipoy. 1966. J.Math.Phys.,7,1137
  * (2) B. H. Voorhees, “Static axially symmetric gravitational fields,” Phys. Rev. D 2, 2119 (1970).
  * (3) G. Erez, N. Rosen “Bull. Res. Council Israel 8F (1959) 47.
  * (4) K. S. Virbhadra, “Directional naked singularity in general relativity,” [arXiv:gr-qc/9606004 [gr-qc]].
  * (5) D. Papadopoulos, B. Stewart and L. Witten, “Some Properties of a Particular Static, Axially Symmetric Space-time,” Phys. Rev. D 24, 320-326 (1981)
  * (6) A. B. Abdikamalov, A. A. Abdujabbarov, D. Ayzenberg, D. Malafarina, C. Bambi and B. Ahmedov, “Black hole mimicker hiding in the shadow: Optical properties of the $\gamma$ metric,” Phys. Rev. D 100, no.2, 024014 (2019)
  * (7) H. Chakrabarty, C. A. Benavides-Gallego, C. Bambi and L. Modesto, “Unattainable extended spacetime regions in conformal gravity,” JHEP 03, 013 (2018)
  * (8) L. Herrera, F. M. Paiva and N. O. Santos, “The Levi-Civita space-time as a limiting case of the gamma space-time,” J. Math. Phys. 40, 4064-4071 (1999)
  * (9) A. N. Chowdhury, M. Patil, D. Malafarina and P. S. Joshi, “Circular geodesics and accretion disks in Janis-Newman-Winicour and Gamma metric,” Phys. Rev. D 85, 104031 (2012)
  * (10) K. Boshkayev, E. Gasperin, A. C. Gutierrez-Pineres, H. Quevedo and S. Toktarbay, “Motion of test particles in the field of a naked singularity,” Phys. Rev. D 93, no.2, 024024 (2016)
  * (11) B. Toshmatov and D. Malafarina, “Spinning test particles in the $\gamma$ spacetime,” Phys. Rev. D 100, no.10, 104052 (2019)
  * (12) G. A. Gonzalez and P. S. Letelier, “Exact general relativistic thick disks,” Phys. Rev. D 69, 044013 (2004) [arXiv:gr-qc/0311078 [gr-qc]].
  * (13) R. Shaikh, S. Paul, P. Banerjee and T. Sarkar, “Shadows and thin accretion disk images of the $\gamma$-metric,” [arXiv:2105.12057 [gr-qc]].
  * (14) G. Lukes-Gerakopoulos, “The non-integrability of the Zipoy-Voorhees metric,” Phys. Rev. D 86, 044013 (2012)
  * (15) V. Bozza, “Gravitational lensing in the strong field limit,” Phys. Rev. D 66, 103001 (2002)
  * (16) V. Bozza, “Gravitational Lensing by Black Holes,” Gen. Rel. Grav. 42, 2269-2300 (2010)
  * (17) V. Bozza, S. Capozziello, G. Iovane and G. Scarpetta, “Strong field limit of black hole gravitational lensing,” Gen. Rel. Grav. 33, 1535-1548 (2001)
  * (18) K. S. Virbhadra and G. F. R. Ellis, “Schwarzschild black hole lensing,” Phys. Rev. D 62, 084003 (2000)
  * (19) K. S. Virbhadra and G. F. R. Ellis, “Gravitational lensing by naked singularities,” Phys. Rev. D 65, 103004 (2002)
  * (20) K. Akiyama et al. [Event Horizon Telescope], “First M87 Event Horizon Telescope Results. II. Array and Instrumentation,” Astrophys. J. Lett. 875, no.1, L2 (2019)
  * (21) T. Harko, Z. Kovacs and F. S. N. Lobo, “Solar system tests of Hořava-Lifshitz gravity,” Proc. Roy. Soc. Lond. A 467, 1390-1407 (2011)
  * (22) A. J. S. Capistrano, P. T. Z. Seidel and L. A. Cabral, “Effective apsidal precession from a monopole solution in a Zipoy spacetime,” Eur. Phys. J. C 79, no.9, 730 (2019)
  * (23) C. M. Will, “The Confrontation between General Relativity and Experiment,” Living Rev. Rel. 17, 4 (2014)
  * (24) C. M. Will, “Was Einstein right?: testing relativity at the centenary,” Annalen Phys. 15, 19-33 (2005)
  * (25) R. Abuter et al. [GRAVITY], “Detection of the Schwarzschild precession in the orbit of the star S2 near the Galactic centre massive black hole,” Astron. Astrophys. 636, L5 (2020) [arXiv:2004.07187 [astro-ph.GA]].
  * (26) C. M. Claudel, K. S. Virbhadra and G. F. R. Ellis, “The Geometry of photon surfaces,” J. Math. Phys. 42, 818-838 (2001)
  * (27) K. S. Virbhadra, D. Narasimha and S. M. Chitre, “Role of the scalar field in gravitational lensing,” Astron. Astrophys. 337, 1-8 (1998)
  * (28) K. S. Virbhadra and C. R. Keeton, “Time delay and magnification centroid due to gravitational lensing by black holes and naked singularities,” Phys. Rev. D 77, 124014 (2008)
  * (29) K. S. Virbhadra, “Relativistic images of Schwarzschild black hole lensing,” Phys. Rev. D 79, 083004 (2009)


11institutetext:  Department of Physics, School of Science, Aristotle
University of Thessaloniki, 541 24, Thessaloniki, Greece
11email<EMAIL_ADDRESS>22institutetext: Facultad de Ciencias Humanas
y de la Educación, Universidad de los Llanos, Villavicencio, Colombia
22email<EMAIL_ADDRESS>33institutetext: Academic Department of
Mathematics, Universidade Tecnológica Federal do Paraná (UTFPR),
Av. Sete de Setembro, 3165, Curitiba, Brazil
33email<EMAIL_ADDRESS>44institutetext: Nonlinear Analysis and Applied
Mathematics (NAAM)-Research Group, Department of Mathematics,
Faculty of Science, King Abdulaziz University, P.O. Box 80203, Jeddah 21589,
Saudi Arabia
44email<EMAIL_ADDRESS>
# Orbit classification in a disk galaxy model with a pseudo-Newtonian central
black hole

Euaggelos E. Zotos Orbit classification in a disk galaxy model with a pseudo-
Newtonian central black holeOrbit classification in a disk galaxy model with a
pseudo-Newtonian central black hole Fredy L. Dubeibe Orbit classification in a
disk galaxy model with a pseudo-Newtonian central black holeOrbit
classification in a disk galaxy model with a pseudo-Newtonian central black
hole André F. Steklain Orbit classification in a disk galaxy model with a
pseudo-Newtonian central black holeOrbit classification in a disk galaxy model
with a pseudo-Newtonian central black hole Tareq Saeed Orbit classification in
a disk galaxy model with a pseudo-Newtonian central black holeOrbit
classification in a disk galaxy model with a pseudo-Newtonian central black
hole

We numerically investigate the motion of stars on the meridional plane of an
axially symmetric disk galaxy model, containing a central supermassive black
hole, represented by the Paczyński-Wiita potential. By using this pseudo-
Newtonian potential we can replicate important relativistic properties, such
as the existence of the Schwarzschild radius. After classifying extensive
samples of initial conditions of trajectories, we manage to distinguish
between collisional, ordered, and chaotic motion. Besides, all starting
conditions of regular orbits are further classified into families of regular
orbits. Our results are presented through modern color-coded basin diagrams on
several types of two-dimensional planes. Our analysis reveals that both the
mass of the black hole (in direct relation with the Schwarzschild radius) as
well as the angular momentum play an important role in the character of orbits
of stars. More specifically, the trajectories of low angular momentum stars
are highly affected by the mass of the black hole, while high angular momentum
stars seem to be unaffected by the central black hole. Comparison with
previous related outcomes, using Newtonian potentials for the central region
of the galaxy, is also made.

###### Key Words.:

Galaxies: kinematics and dynamics – Galaxies: structure, chaos – Black hole
physics

##  1 Introduction

Galactic dynamics is of paramount importance not only because of the evident
astronomical interest in classifying and analyzing the nature of galaxies
(Binney & Tremaine, 2008; Manos & Athanassoula, 2011; Zotos, 2012), but also
to investigate the nature of dark matter that constitutes a significant part
of these objects (Zotos & Caranicolas, 2013). Exploring the nature of orbits
on different models of galaxies can provide a hint to the nature of dark
matter present in these structures. Moreover, the phase space structure of a
particular model shall depend on its basic components, i.e., disk, halo, and
the nucleus or central bulge (Zotos, 2012; Zotos & Carpintero, 2013; Zotos,
2014), such that its phase space is as unique as a fingerprint. In general
terms, orbits can be classified according to their dynamical behavior into
regular and chaotic (Caranicolas, 1996; Caranicolas & Papadopoulos, 2003;
Zotos, 2012; Zotos & Carpintero, 2013; Zotos, 2014), while regular orbits can
be subclassified into different families depending on its resonances (Martinet
& Mayer, 1975; Manabe, 1979; Greiner, 1987; Lees & Schwarzschild, 1992), which
can be analyzed using different technics (Binney & Spergel, 19828, 1984;
Laskar, 1993; Carpintero & Aguilar, 1998; Zotos & Carpintero, 2013) that in
essence, decomposes a trajectory function into its constituent frequencies.
Also, it is important to note that a typical galaxy has billions of stars,
making simulations computationally expensive, therefore several analytical
models are proposed based on Newtonian theory and General Relativity.

In Newtonian theory, globular and spherical galaxies are represented by the
popular models of Plummer (Plummer, 1911) and King (King, 1966). Highly
flattened axisymmetric galaxies are usually modeled by the so-called Toomre
model (Toomre, 1963), while thickened-disk are modeled by Miyamoto and Nagai
potentials (Miyamoto & Nagai, 1975). The Miyamoto-Nagai model is a stationary
and axially symmetric potential, which despite the deviations in its density
profile in comparison with realistic models, provides a simple and reasonable
model for a disk. Concerning the galactic nuclei, most of the analytical
models are spherical, finite mass, and its isotropic velocity dispersion can
be calculated analytically (see e.g., Jaffe, 1983; de Zeeuw, 1985; Hernquist,
1990; Zhao, 1996). Newtonian models for massive black hole nuclei have been
also proposed aiming to reproduce the final stage for the evolution of active
galactic nuclei Rees (1984) or to mimic realistically the nuclei of galaxies
(Häring & Rix, 2004). The linearity of Newtonian equations allows the
superposition of different models to study, e.g., galactic models composed of
a disk and a bulge (Binney & Tremaine, 2008).

In General Relativity besides the well-known spherical solutions, several
exact disk-like solutions have been proposed to model galaxies, including
static thin disks (Bonnor & Sackfield, 1968; Morgan & Morgan, 1969; Ujevic &
Letelier, 2004) and thick disks (González & Letelier, 2004). In order to
obtain these models is usual to start with the metric and calculate the
energy-momentum from the metric. A different approach uses the energy-momentum
tensor as a source of the Einstein equations but this is highly nontrivial. In
the general relativistic framework, the composition of different models is not
always possible due to the nonlinearity of the Einstein equations.
Nevertheless, there are some special cases where the superposition of exact
solutions is possible (Lemos & Letelier, 1993; Semerak, 2002, 2004). A review
of these and other relativistic models can be found in (Semerak, 2002; Vogt &
Letelier, 2005).

There is compelling evidence of supermassive black holes in the center of
several galaxies, including Milky Way (Liu et al., 2019), where general
relativity is found to provide a natural scenario for its study. On the other
side, Newtonian models are computationally cheaper and can be easily built and
analyzed. A balance between the two extremes can be achieved using Newtonian
equations that include pseudo-potentials which reproduce some relativistic
effects. In particular, some of the most renowned pseudo-Newtonian potentials
have been introduced originally to study accretion disks in Black Holes
(Paczyński & Wiita, 1980; Semerak & Karas, 1999). The simplest and at the same
time practical of these potentials is the one derived by Paczyński and Wiita
(Paczyński & Wiita, 1980), which can correctly reproduce the locations of both
the marginally bounded circular orbit and the last stable circular orbit of
the Schwarzschild metric (Abramowicz, 2009).

Since is a well-known fact that pseudo-Newtonian potentials can have a
remarkable effect on the dynamics of different systems (Steklain & Letelier,
2006; Dubeibe et al., 2017; Zotos & Steklain, 2019; Zotos et al., 2019), in
the present paper we shall use the pseudo-Newtonian Paczyński-Wiita potential
(Paczyński & Wiita, 1980), for modeling the central part of a galaxy hosting a
supermassive black hole, while the flat thin disk is represented by the
Miyamoto-Nagai potential. Aiming to obtain direct outcomes on how the mimicked
relativistic properties of the supermassive black hole affect the character of
the trajectories of the stars, we study the influence of the Schwarzschild
radius and the angular momentum of a test particle on the orbit
classification, i.e., collisional, ordered, and chaotic. Therefore, we would
be able to obtain direct outcomes on how the relativistic properties of the
supermassive black hole (such as the Schwarzschild radius) affect the
character of the trajectories of the stars, by comparing results from previous
works (see e.g., Zotos, 2014; Zotos & Carpintero, 2013), where the central
nucleus was modeled by using a simple Newtonian Plummer potential (Plummer,
1911). Regular and chaotic orbits will be identified by using the SALI method
introduced by (Skokos et al., 2004), while regular orbits are further
classified into families of regular orbits (Carpintero & Aguilar, 1998; Zotos
& Carpintero, 2013).

The article is structured as follows: in Section 2 we present the main
properties of the galaxy model. In the following Section 3 we explain the
computational methods we use for obtaining the classification of the orbits
which is presented in Section 4. Our article ends with Section 5, where we
emphasize the main outcomes of our analysis.

##  2 Description of the galaxy model

Our galaxy model describes the motion of stars on the meridional $(R,z)$ plane
and consists of two components. The first component is a flat thin disk,
represented by the Miyamoto-Nagai potential (Miyamoto & Nagai, 1975)

$\Phi_{\rm d}(R,z)=-\frac{GM_{\rm
d}}{\sqrt{R^{2}+\left(k+\sqrt{z^{2}+s^{2}}\right)^{2}}},$ (1)

where, $M_{\rm d}$, $k$, and $s$ are the mass, the scale length and the scale
height of the disk, respectively.

For the description of the central black hole, we use the Paczyński-Wiita
potential (see e.g., Paczyński & Wiita, 1980; Abramowicz, 2009).

$\Phi_{\rm bh}(R,z)=-\frac{GM_{\rm bh}}{\sqrt{R^{2}+z^{2}}-r_{s}},$ (2)

where $M_{\rm bh}$ is the mass of the black hole. This is a very simple and
practical pseudo-Newtonian potential, which however can realistically
replicate and model general relativistic effects, associated with the motion
of test particles around a non-spinning black hole. Moreover, the Paczyński-
Wiita potential alters the classical Newtonian potential $1/r$ to
$1/(r-r_{s})$, thus introducing the relativistic Schwarzschild radius $r_{s}$.

Due to the fact that the total potential $\Phi_{t}(R,z)=\Phi_{\rm d}+\Phi_{\rm
bh}$ has an axial symmetry, we have the conservation of the $z$-component
$(L_{z})$ of the total angular momentum. This automatically implies that the
motion of the test particle can be described using the effective potential

$\Phi_{\rm eff}(R,z)=\Phi_{t}(R,z)+\frac{L_{z}^{2}}{2R^{2}},$ (3)

and the set of the following equations of motion

$\displaystyle\ddot{R}=-\frac{\partial\Phi_{\rm eff}}{\partial R},$ (4)
$\displaystyle\ddot{z}=-\frac{\partial\Phi_{\rm eff}}{\partial z}.$ (5)

The energy integral of the system is given through the Hamiltonian

$H(R,z,\dot{R},\dot{z})=\Phi_{\rm
eff}(R,z)+\frac{1}{2}\left(\dot{R}^{2}+\dot{z}^{2}\right)=E,$ (6)

where $E$ is the test particle’s conserved value of the total orbital energy.
Here it should be noted that the test particle can access those areas of the
phase space where $E\geq\Phi_{\rm eff}$.

In our computations we use a system of units in which the unit of length is 1
kpc, the unit of time is $10^{6}$ yr (1 Myr), the unit of mass is
$2.22508\times 10^{11}{\rm M}_{\odot}$, the velocity unit is 978.564 km s-1,
the angular momentum unit (per unit mass) is 978.564 km kpc-1 s-1, while the
unit of energy (per unit mass) is $9.576\times 10^{5}$ km2s-2. Then for the
gravitational constant we have that $G=1$. In the above units the values of
the constant parameters are: $M_{\rm d}=0.5$ (corresponding to about $10^{11}$
M${}_{\odot})$, $k=3$, and $s=0.175$. On the other hand, the parameters
$L_{z}$ and $M_{\rm bh}$ are treated as free parameters.

The Schwarzschild radius is related to the mass of the black hole through the
relation

$r_{s}=\frac{2GM_{\rm bh}}{c^{2}}.$ (7)

In the adopted system of galactic units, the value of the speed of light is
roughly equal to 300 velocity units. Moreover, the values of the mass of the
central black hole lie in the interval $[10^{8}{\rm M}_{\odot},10^{10}{\rm
M}_{\odot}]$, which is a typical range of values in disk galaxies (Binney &
Tremaine, 2008). Consequently, from Eq. (7) we derive that the values of the
Schwarzschild radius lie in the interval $[10^{-8},10^{-6}]$. However, such
low values will certainly induce computational malfunctions, since almost all
numerical integrators cannot cope with motion extremely close to the
singularity111Our previous experience indicates that the Bulirsch-Stoer
integrator we use in our study can adequately follow the evolution of a test
particle up to a radius equal to $10^{-5}$ around the singularity, while for
smaller radii it displays several well known malfunctions, such as time step
“overflow”.. On this basis, in order to eliminate all ill behaviors during the
numerical integration, we define a close encounter radius $R_{c}$ around the
singularity as $R_{c}=10^{3}r_{s}$. At this point, it should be noted that
although there are some studies in which the motion of test particles can
reach a closer position to the galactic center (see e.g. Sukova & Semerak
(2013)), in the current study the close encounter radius does not affect the
orbital dynamics of the test particle because the numerical artifact solely
stops the numerical integration, therefore in practice, its direct effect is
the loss of information in the basin diagrams around the origin in a disk of
the order $[10^{-5},10^{-3}]$ units of length for $M_{bh}=[10^{8}{\rm
M}_{\odot},10^{10}{\rm M}_{\odot}]$, respectively.

Figure 1: The geometry of the zero velocity curves (black lines), near the
central region of the galaxy, for $M_{\rm bh}=0.05$ and (a): $L_{z}=0.1$, (b):
$L_{z}=0.5$. The energetically allowed and forbidden regions of motion are
indicated using white and gray color, respectively, while the red line
corresponds to the Schwarzschild radius.

Figure 2: White areas correspond to open ZVCs, where collision to the central
black hole is possible, while gray areas denote closed ZVCs, where the test
particle is not energetically allowed to approach the Schwarzschild radius.

In Zotos et al. (2018) we have shown that as long as $r_{s}>0$ there exists an
additional energetically forbidden circular region around the black hole.
Furthermore, the radius of the circular region is equal to the Schwarzschild
radius. The interesting question is the following: is the test particle always
able to reach the central black hole? The answer is no. Whether the test
particle can approach the Schwarzschild radius or not, strongly depends on the
value of its angular momentum. More precisely, if its angular momentum is low
enough then it can approach the Schwarzschild radius of the black hole and
display a collision event. On the other hand, if its angular momentum is
relatively high then it can never reach the boundaries of the black hole. In
panel (a) of Fig. 1, we present the first case, where the approach to the
central black hole is possible. In this case, the zero-velocity curves (ZVCs),
defined as $\Phi_{\rm eff}(R,z)=E$, are open, thus allowing the test particle
to approach the black hole. In panel (b) of the same figure, we see the case
where the test particle cannot reach the central region of the galaxy. This is
because in this case the ZVCs are closed, thus acting as potential barriers.
In Fig. 2 we show the regions of open (white) and closed (gray) ZVCs, as a
function of $(L_{z},M_{\rm bh})$.

##  3 Computational methodology

Figure 3: A collection showing the most important types of motion, when
$M_{\rm bh}=0.05$. (a): box orbit; (b): 2:1 orbit; (c): 1:1 orbit; (d): 2:3
orbit; (e): 4:3 orbit; (f): chaotic orbit.

Our work aims to determine how the mass $M_{\rm bh}$ of the black hole, along
with the value of the angular momentum $L_{z}$, affects the character of
motion of the test particles. For this task, we perform an orbit
classification, by numerically integrating sets of $1024\times 1024$ initial
conditions of orbits. In all cases, the value of the total orbital energy is
fixed to $E=-0.05$, which corresponds to the maximum possible value of $R$
coordinate (about 10 kpc for a typical disk galaxy (see e.g., Binney &
Tremaine, 2008)). We decided to keep fixed the value of the energy because its
influence on the character of orbits moving in the meridional plane has
already been revealed in Zotos (2016).

All initial conditions of the trajectories were numerically integrated for
$10^{4}$ time units, which correspond to about 1 Hubble time (14 billion
years) while using a constant time step, equal to $10^{-3}$. During the
numerical integration, we also compute the Smaller Alignment Index (SALI)
(Skokos, 2001), so as to be able to distinguish between chaotic and regular
motion. This method allows classifying the orbits according to the numerical
value obtained after evolving two orthonormal deviation vectors $\@vec{w_{1}}$
and $\@vec{w_{2}}$, which are periodically normalized to prevent overflow
errors. More precisely, if SALI $>10^{-4}$ the orbit is categorized as
regular, while if SALI $<10^{-8}$ it is classified as chaotic, or if the
result belongs to the interval $10^{-4}<$ SALI $<10^{-8}$, it is classed as
sticky, and the orbit requires a longer time of integration to be classified
222We refer the reader to (Skokos, 2001) for a detailed explanation of the
method.. The index is formally defined as
$\mathrm{SALI}\equiv\min\left(d_{-},d_{+}\right)$, with

$d_{\mp}\equiv\left\|\frac{\@vec{w}_{1}}{\left\|\@vec{w}_{1}\right\|}\mp\frac{\<EMAIL_ADDRESS>(8)

Once an orbit was classified as a regular one, then we performed an additional
categorization, thus classifying regular initial conditions into regular
families. Our analysis indicates that several types of regular orbits appear
in our disk galaxy model, while the most important ones are: (i) box orbits;
(ii) 1:1 orbits; (iii) 2:1 orbits; (iv) 2:3 orbits; (v) 4:3 orbits; and (vi)
higher resonant orbits. Specifically, in all studied cases, the relative
percentage of higher resonant orbits is always less than 0.1% and therefore we
may argue that these orbits do not have any significant contribution to the
overall orbital dynamics of the galaxy. The notation of orbits (box orbits and
$n:m$ resonant orbits) in the axially symmetric galaxy model is according to
the initial works of Carpintero & Aguilar (1998) and Zotos & Carpintero
(2013). In panels (a)-(e) of Fig. 3 we provide the shapes of the main regular
orbits of the system, while in panel (f) of the same figure we give an example
of a typical chaotic orbit when $M_{\rm bh}=0.05$. In all cases, the ZVC is
the black thick curve circumscribing the orbits. For all regular types of
trajectories, we tried to present regular orbits with starting conditions very
close to the respective parent periodic orbits, to be able to clearly
recognize their shapes.

A double-precision Bulirsch-Stoer integrator, written in `FORTRAN 77` (e.g.,
Press et al., 1992), has been applied for the numerical integration of the
starting conditions of the trajectories. During the numerical integration, the
test particle’s total energy was sufficiently conserved, while the
corresponding observed error was of the order of $10^{-14}$. For the
classification of the initial conditions (per grid) the required CPU time,
using a Quad-Core Intel i7 4.0 GHz processor was varying between 5 and 22
hours, strongly depending on the number of collision orbits. Version 12.1 of
Mathematica® software has been deployed for developing all the graphical
illustration of the article.

##  4 Orbit classification

Figure 4: Basin diagrams on the $(R,z)$ plane for $L_{z}=0.001$, with (a):
$M_{\rm bh}=0.0005$ and (b): $M_{\rm bh}=0.05$.

Figure 5: Distribution of the collision times for $L_{z}=0.001$, with (a):
$M_{\rm bh}=0.0005$ and (b): $M_{\rm bh}=0.05$.

Figure 6: Basin diagrams on the $(R,z)$ plane for $L_{z}=0.1$, with (a):
$M_{\rm bh}=0.0005$ and (b): $M_{\rm bh}=0.05$.

Figure 7: Basin diagrams on the $(R,z)$ plane for $L_{z}=0.5$, with (a):
$M_{\rm bh}=0.0005$ and (b): $M_{\rm bh}=0.05$.

Figure 8: Basin diagrams on the $(R,L_{z})$ plane, with (a): $M_{\rm
bh}=0.0005$ and (b): $M_{\rm bh}=0.05$.

Figure 9: Evolution of the percentages of all main types of orbits on the
$(R,L_{z})$ plane, as a function of the angular momentum, when (a): $M_{\rm
bh}=0.0005$ and (b): $M_{\rm bh}=0.05$.

In this section, we will present the nature of the motion of the test particle
as a function of the mass of the black hole and the angular momentum. For
illustrating the orbital properties of the galactic system, hereafter we shall
follow the graphical approach introduced in Nagler (2004, 2005), thus
presenting color-coded basin diagrams, where each pixel corresponds to a
unique trajectory and it is colored according to the final state of the test
particle.

All trajectories have initial conditions $(R_{0},z_{0})$, while for a given
total orbital energy $E_{0}$ the initial velocities are given by

$\displaystyle\dot{R_{0}}=\frac{z_{0}}{d_{0}}f(R_{0},z_{0};E_{0}),$ (9)
$\displaystyle\dot{z_{0}}=-\frac{R_{0}}{d_{0}}f(R_{0},z_{0};E_{0}),$ (10)

where

$f(R_{0},z_{0};E_{0})=\sqrt{2\left(E_{0}-\Phi_{\rm eff}(R_{0},z_{0})\right)},$
(11)

while $d_{0}=\sqrt{R_{0}^{2}+z_{0}^{2}}$.

###  4.1 Low angular momentum

We begin our analysis, considering the case where the test particle (star)
takes low values of the angular momentum. In Fig. 4(a-b) we present the basin
diagrams on the $(R,z)$ plane, when $L_{z}=0.001$. We consider only the $z>0$
part because due to the axial symmetry of the galaxy the orbital structure of
the $z<0$ part is mirror-symmetrical. Panel (a) corresponds to the case where
the mass of the black hole is $M_{\rm bh}=0.0005$. The first type of motion
which we encounter as we move away from the center is the box orbits. In fact,
there are two main basins occupied by initial conditions of box orbits, while
between these basins we have the presence of the 2:1 resonant orbits. Both
basins of box orbits are “polluted” by a mixture of chaotic and collision
orbits, while tiny stability islands of higher resonant orbits are also
present. Very close to the boundaries of the ZVC we observe the presence of
2:1 and 1:1 resonant orbits, while there are no signs of other resonant
orbits, such as 4:3.

The case with $M_{\rm bh}=0.05$ is displayed in panel (b) of Fig. 4. Here we
see that the pollution of collision orbits we have seen in panel (a) is much
stronger leading to well-defined basins of collision. Inside the collision
basin near the central black hole, we observe the presence of a stability
island corresponding to 1:1 resonant orbits. On the other hand, higher
resonant orbits seem to be completely absent for this value of $M_{\rm bh}$.
Moreover, all initial conditions belonging to chaotic orbits (according to the
value of SALI) lead, sooner or later, to a collision with the central black
hole, and therefore there are no remaining trapped chaotic trajectories.

In order to illustrate more clearly the difference between the two cases
regarding the collision orbits, we present in Fig. 5(a-b) the distributions of
the corresponding collision times of the trajectories. In panel (a),
corresponding to $M_{\rm bh}=0.0005$ one can see that the vast majority of the
starting conditions need a significant amount of time (more than 2500 time
units) for reaching the central black hole. However, we distinguish some very
thin basins of initial conditions for which the corresponding collision time
is very low (less than 300 time units). On the other hand, in panel (b) of
Fig. 5 we see that almost all the initial conditions lead to a collision with
the black hole within the first 1000 time units of the numerical integration.
Therefore, one may reasonably ask: why the orbits collide so quickly in the
second case? Fortunately, the answer is very easy. In both cases, the ZVCs are
open near the center and therefore the test particle can approach the central
black hole. The main difference is that in the second case (where $M_{\rm
bh}=0.05$) the Schwarzschild radius is 100 times larger which implies that the
collision channel is much wider. Consequently, the test particles need
significantly less time to find the open throat and collide with the
singularity.

###  4.2 Intermediate angular momentum

The second case under investigation concerns test particles moving with
intermediate values of angular momentum. The basin diagrams on panels (a) and
(b) of Fig. 6 show the orbital structure of the meridional $(R,z)$ plane, when
$L_{z}=0.1$. Panel (a) corresponds to $M_{\rm bh}=0.0005$. Here we see that
collision motion is not possible. This is also seen in the diagram of Fig. 2,
where one can extract the valuable information that for $(M_{\rm
bh},L_{z})=(0.0005,0.1)$ the corresponding ZVCs are closed. We may argue that
for this set of values of $M_{\rm bh}$ and $L_{z}$ almost the entire $(R,z)$
plane is covered by box and 2:1 orbits, while all other types of orbits are
extremely limited or even absent.

For a much higher value of the mass of the black hole, we see in panel (b) of
Fig. 6 that extended collision basins emerge close and far from the central
singularity, thus leading to a considerable reduction of the area
corresponding to box orbits. Besides, one can see that now the test particle
can perform also 1:1 and 4:3 resonant orbits.

###  4.3 High angular momentum

The last scenario under consideration is the case with high angular momentum
stars. Fig. 7(a-b) shows the basins diagrams on the $(R,z)$ plane, when
$L_{z}=0.5$. Now we observe that the value of the mass of the black hole
affects only resonant orbits. This implies that high angular momentum stars,
moving on box and 2:1 orbits are hardly affected by the central black hole.
Furthermore, this type of stars cannot even approach the singularity, taking
into account that for such high values of the angular momentum the ZVCs are
always closed, thus not allowing collision with the black hole. Another
interesting fact is that the 1:1 and 2:3 resonant types of motion seem to be
more prominent, compared to all other previous cases, while on the other hand,
the 4:3 resonant orbits are barely noticeable. Moreover, it should be noted
that with increasing value of the angular momentum, while keeping constant the
value of the total orbital energy, the energetically allowed region on the
meridional $(R,z)$ plane is considerably reduced.

###  4.4 Overview analysis

It would be very informative to present from another perspective the influence
of the angular momentum $L_{z}$ as well as of the mass of the black hole
$M_{\rm bh}$ on the nature of the motion of the test particles.

In Fig. 8(a-b) we present the character of the trajectories, with initial
conditions on the $(R,L_{z})$ plane, while for all orbits we set $z_{0}=0$. It
is seen, that when the value of the mass of the black hole is low, collision
orbits are very limited and appear only for very low values of the angular
momentum $(L_{z}<0.01)$, that is when stars can approach very close to the
galactic center. On the other hand, when $M_{\rm bh}$ is high enough collision
is still possible, even for relatively large angular momentum values
$(L_{z}<0.1)$. Box orbits occupy most of the $(R,L_{z})$ plane, while the
second most populated type of motion is the 2:1 resonant family, which lies
around 4-5 kpc. Furthermore, one can see, that for $M_{\rm bh}=0.05$, the 4:3
resonant orbits appear at about $L_{z}=0.1$, while the 2:3 resonant orbits
exist mainly for $L_{z}>0.2$. Isolated initial conditions of higher resonant
orbits can be observed scattered inside the area of box orbit.

The parametric evolution of the percentages of all the main types of orbits,
as a function of the angular momentum, is given in panels (a) and (b) of Fig.
9. From this figure, it is seen that the orbital content of the system is
affected by the value of the mass of the black hole mainly for low values of
the angular momentum $(L_{z}<0.2)$, while for higher values the pattern is
almost unaffected by the shift on the value of $M_{\rm bh}$. We also see, that
for large angular momentum values $(L_{z}\to 1)$ box orbits completely
dominate, while at the same time resonant orbits, such as the 2:3 and 4:3
orbital families disappear.

Figure 10: Evolution of the average collision time of the trajectories, with
initial conditions on the $(R,L_{z})$ plane, as a function of the angular
momentum. The horizontal black dashed line indicates the total time of the
numerical integration, while the vertical red dashed lines denote the maximum
value of $L_{z}$, for which the ZVCs are open.

Before closing, we would like to present some additional information regarding
the collision time of the trajectories and the connection between the angular
momentum and the mass of the black hole. In Fig. 10 we display the evolution
of the average collision of orbits as a function of the angular momentum, for
three values of $M_{\rm bh}$. We see, that the average collision time is
reduced, with increasing value of the mass of the black hole. Moreover, in all
three cases, the collision time increases, with increasing value of the
angular momentum, and tends to the total time of the numerical integration
($10^{4}$ time units). This is anticipated because with an increasing value of
$L_{z}$ the width of the collision channel of the ZVCs, near the galactic
center, is reduced. The vertical red dashed lines indicate the maximum values
of the angular momentum for which a collision is possible (evidently, for
higher values of $L_{z}$ the corresponding ZVCs are closed). We would also
like to note that these critical values of the angular momentum, obtained
through the numerical integration, coincide with the respective values
obtained earlier (see Fig. 2), using analytical methods.

##  5 Concluding remarks

In this work, we numerically investigated the motion of stars on the
meridional $(R,z)$ plane of an axially symmetric disk galaxy model, containing
a central supermassive black hole. For modeling the black hole we used the
Paczyński-Wiita potential which can replicate important relativistic
properties, such as the existence of the Schwarzschild radius $r_{s}$. After
classifying extensive samples of initial conditions of trajectories, we
managed to distinguish between collisional, ordered, and chaotic motion.
Furthermore, all starting conditions of regular orbits are further classified
into families of regular orbits. Our results are presented through modern
color-coded basin diagrams on several types of two-dimensional planes, while
we also monitored the evolution of the respective percentages.

The following list summarizes the main outcomes, regarding the motion of
stars:

  1. 1.

We have seen that the parameters of both the angular momentum and the mass of
the black hole determine if the stars can approach or not to the central
singularity. In general terms, open ZVCs (which allow collision with the black
hole) are present only when $L_{z}<0.2$.

  2. 2.

For low values of the angular momentum, only the most basic types of orbits
exist (such as box, 1:1, and 2:1), while higher resonant orbits (such as the
2:3 and the 4:3 orbital families) appear only for higher values of the
$L_{z}$.

  3. 3.

For relatively low values of the mass of the black hole $(M_{\rm bh}\simeq
0.0005)$ we encountered the phenomenon of trapped chaos, where chaotic orbits
remained trapped inside the ZVCs for at least one Hubble time, even though
collision to the central black hole is energetically allowed. On the other
hand, for high values of $M_{\rm bh}$, there is no evidence of trapped chaos
and all chaotic orbits lead fast to a collision.

  4. 4.

Low angular momentum stars (with $L_{z}<0.1$) are highly affected by the mass
of the black hole. More specifically, for low values of $M_{\rm bh}$ they can
either move on a box, 2:1, or chaotic orbit, while a collision is less
possible. On the contrary, for high values of $M_{\rm bh}$ the possibility of
collision is much higher, while regular motion seems to be shared between box
and 2:1 types of orbits.

  5. 5.

High angular momentum stars (with $L_{z}>0.5$) seem to be less affected by the
value of the mass of the black hole since they cannot even approach the
central regions of the galaxy. Our analysis indicates that the vast majority
of high angular momentum stars move either on box or 2:1 trajectories.

In two earlier works (Zotos, 2014; Zotos & Carpintero, 2013) we have also
investigated the nature of motion in a disk galaxy model, using the same
Miyamoto-Nagai model. However, the central nucleus was modeled, using a
typical Plummer potential (Plummer, 1911), while in the present study we use a
pseudo-Newtonian Paczyński-Wiita potential (Paczyński & Wiita, 1980).
Therefore, we can directly compare the differences between the classical
Newtonian and the pseudo-Newtonian models. In both cases, the main types of
regular orbits are the same, while the influence of the angular momentum seems
to be similar. In particular, we see that high angular momentum stars move
mainly on regular orbits, while low angular momentum stars can display chaotic
motion, by passing close to the central region. The main difference between
the two cases is the fact that in the case of the Plummer potential we
encountered additional types of resonant orbits (e.g., 5:4 and 6:5). Thus, we
may argue that in the pseudo-Newtonian case of the Paczyński-Wiita potential
the orbital content of regular motion is not as rich as in the case of the
Plummer potential. Additionally, in the case of the Plummer potential several
types of resonant orbits (such as the 2:3, 4:3, 5:4 and 6:5) were
simultaneously present. On the contrary, in the case of the Paczyński-Wiita
potential, we observe a hierarchy regarding the appearance of resonant orbits.
Specifically, for low values of $L_{z}$ we have the presence of the 4:3
orbits, while the 2:3 orbits appear for much higher values of the angular
momentum. Finally, for the Newtonian case of the Plummer potential we
encountered many more types of higher resonant orbits $n:m$ (with $n,m>5$)
inside the area of the box orbits, while in the case of the Paczyński-Wiita
potential all these higher resonant orbits are completely absent.

Consequently, we believe that our findings can be relevant not only from a
theoretical point of view but also from an astrophysical perspective. In
particular, the identification of resonant orbits could shed lights on the
formation and origin of kinematic moving groups (MGs) (Martinez-Medina et al.,
2016), while the presence of chaotic or collisional orbits could certainly
modify the structure and hence the evolution of the galaxy with the increase
of the central black hole mass. Taking into account the positive and
encouraging results of the present analysis, it is in our plans to expand our
orbit classification into three dimensions. We would aim to determine how the
mass of the black hole in combination with the angular momentum of the stars,
affects their trajectories inside the three-dimensional $(x,y,z)$ space.
Moreover, it would also be interesting to determine the influence of the mass
of the black hole on the network of periodic orbits of the galaxy.

## Acknowledgments

This project was funded by the Deanship of Scientific Research (DSR) at King
Abdulaziz University, Jeddah, Saudi Arabia, under grant no. KEP-17-130-38. The
authors, therefore, acknowledge with thanks DSR technical and financial
support. FLD acknowledges financial support from the Universidad de los
Llanos. Our warmest thanks go to the anonymous referee for the careful reading
of the manuscript as well as for all the apt suggestions and comments, which
allowed us to improve both the quality and the clarity of the paper.

## References

  * Abramowicz (2009) Abramowicz, M.A. 2009. A&A, 500, 213
  * Binney & Tremaine (2008) Binney, J., Tremaine, S. 2008. Galactic Dynamics. Princeton Univ. Press, Princeton, USA
  * Binney & Spergel (19828) Binney, J., Spergel, D. 1982. The Astrophysical Journal, 252, 308.
  * Binney & Spergel (1984) Binney, J., Spergel, D. 1984. MNRAS, 206(1), 159.
  * Bonnor & Sackfield (1968) Bonnor, W. B., Sackfield, A. 1968. Commun. Math. Phys. 8, 338
  * Carpintero & Aguilar (1998) Carpintero, D.D., Aguilar, L.A. 1998. MNRAS, 298, 1
  * Caranicolas (1996) Caranicolas, N. D. 1996. Astrophys. Space Sci. 246, 15
  * Caranicolas & Papadopoulos (2003) Caranicolas, N. D., Papadopoulos, N. J. 2003. Astron. Astrophys. 399, 957
  * Contopoulos (1979) Contopoulos, G.: in Stochastic behavior in classical and quantum Hamiltonian systems, G. Casati and J. Ford Eds., p. 1-17 (1979)
  * Dubeibe et al. (2017) Dubeibe, F. L., Lora-Clavijo, F. D., González, G. A. 2017. Physics Letters A, 381, 563.
  * González & Letelier (2004) González, G. A., Letelier, P. S. 2004. Phys. Rev. D. 69, 044013.
  * Greiner (1987) Greiner, J. 1987. Cel. Mech. 40, 171.
  * Häring & Rix (2004) Häring, N., Rix, H. W. 2004. The Astrophysical Journal Letters, 604, L89.
  * Hernquist (1990) Hernquist, L. 1990. The Astrophysical Journal, 356, 359.
  * Jaffe (1983) Jaffe, W. 1983. MNRAS, 202, 995.
  * King (1966) King, I.R. 1966. Astron J, 71, 64
  * Laskar (1993) Laskar, J. 1993. Physica D: Nonlinear Phenomena, 67, 257.
  * Lees & Schwarzschild (1992) Lees, J. F., Schwarzschild, M. 1992. The Astrophysical Journal, 384, 491.
  * Lemos & Letelier (1993) Lemos, J.P.S., Letelier, P.S. 1993. Class. Quantum Grav, 10, 6
  * Liu et al. (2019) Liu, J., Zhang, H., Howard, A. W., Bai, Z., Lu, Y., et al. 2019. Nature, 575, 7784
  * Manos & Athanassoula (2011) Manos, T., Athanassoula, E. 2011. MNRAS, 415, 1
  * Manabe (1979) Manabe, S. 1979. Publications of the Astronomical Society of Japan, 31, 369.
  * Martinet & Mayer (1975) Martinet, L., Mayer. F. 1975. Astron. Astrophys., 44, 45
  * Martinez-Medina et al. (2016) Martinez-Medina, L. A., Pichardo, B., Moreno, E., Peimbert, A., Velazquez, H. 2016. The Astrophysical Journal Letters, 817(1), L3.
  * Morgan & Morgan (1969) Morgan, T., Morgan, L. 1969. Phys Rev, 183, 5
  * Miyamoto & Nagai (1975) Miyamoto, M., Nagai, R. 1975. Publications of the Astronomical Society of Japan, 27, pp. 533-543 (1975)
  * Nagler (2004) Nagler, J. 2004. Phys. Rev. E, 69, 066218
  * Nagler (2005) Nagler, J. 2005. Phys. Rev. E, 71, 026227
  * Paczyński & Wiita (1980) Paczyński, B., Wiita, P.J. 1980, A&A, 88, 23
  * Plummer (1911) Plummer, H.C. 1911. MNRAS, 71, 5
  * Press et al. (1992) Press, H.P., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. 1992. Numerical Recipes in FORTRAN 77, 2nd edn. Cambridge University Press, Cambridge
  * Rees (1984) Rees, M. J. 1984. Black hole models for active galactic nuclei. Annual review of astronomy and astrophysics, 22, 471.
  * Semerak (2002) Semerak, O. 2002. Following the Prague Inspiration, to Celebrate the 60th Birthday of Jiri Bicák. World Scientific, Singapore.
  * Semerak (2004) Semerak, O. 2004. Class. Quantum Grav. 21, 2203.
  * Semerak & Karas (1999) Semerak, O., Karas, V. 1999. Astron Astroph, 343, 2
  * Skokos (2001) Skokos, C. 2001. Journal of Physics A, 34, 10029
  * Skokos et al. (2004) Skokos, C., Antonopoulos, C., Bountis, T. C., Vrahatis, M. N. 2004. Journal of Physics A: Mathematical and General, 37, 6269.
  * Sukova & Semerak (2013) Sukova, P., Semerak, O. 2013. MNRAS, 436(2), 978.
  * Steklain & Letelier (2006) Steklain, A.F., Letelier, P.S. 2006. Phys Lett A, 352, 4
  * Toomre (1963) Toomre, A. 1963. Astrophys J, 138, 385.
  * Ujevic & Letelier (2004) Ujevic, M., Letelier, P.S. 2004. Phys Rev D, 70, 8
  * Vogt & Letelier (2005) Vogt, D., Letelier, P.S. 2005. MNRAS, 363, 1
  * Wolfram (2003) Wolfram, S. 2003. The Mathematica Book, 5th edn. Wolfram Media, Champaign
  * de Zeeuw (1985) de Zeeuw, T. 1985. MNRAS, 216, 273.
  * Zhao (1996) Zhao, H. 1996. MNRAS, 278, 488.
  * Zotos (2012) Zotos, E.E. 2012. New Astronomy, 17, 6
  * Zotos (2014) Zotos, E.E. 2014. Mech. Res. Commun., 62, 102
  * Zotos (2016) Zotos, E.E. 2016. Baltic Astron., 25, 139
  * Zotos & Caranicolas (2013) Zotos, E.E., Caranicolas, N.D. 2013. A&A, 560, A110
  * Zotos & Carpintero (2013) Zotos, E.E., Carpintero, D.D. 2013. Celest. Mech. Dyn. Astron., 116, 417
  * Zotos et al. (2018) Zotos, E.E., Dubeibe, F.L., González, G.A. 2018. MNRAS, 477, 5388
  * Zotos et al. (2019) Zotos, E. E., Dubeibe, F. L., Nagler, J., Tejeda, E. 2019. MNRAS, 487, 2340.
  * Zotos & Steklain (2019) Zotos, E.E., Steklain, A.F. 2019. Astroph. Space Science, 364, 10