# The effect of environment in the timing of a pulsar orbiting SgrA∗

Amodio Carleo<EMAIL_ADDRESS>Dipartimento di Fisica “E.R Caianiello”,
Università degli Studi di Salerno,
Via Giovanni Paolo II, 132 - 84084 Fisciano (SA), Italy  Istituto Nazionale di
Fisica Nucleare - Gruppo Collegato di Salerno,
Via Giovanni Paolo II, 132 - 84084 Fisciano (SA), Italy.  INAF, Osservatorio
Astronomico di Cagliari, Via della Scienza 5, 09047 Selargius (CA), Italy
Bilel Ben-Salem<EMAIL_ADDRESS>ZARM, Universität Bremen,
Am Fallturm 2, 28359 Bremen, Germany  Fakulät für Physik, Universität
Bielefeld, Postfach 100131, 33501 Bielefeld, Germany

###### Abstract

Pulsars are rapidly rotating neutron stars emitting intense electromagnetic
radiation that is detected on Earth as regular and precisely timed pulses. By
exploiting their extreme regularity and comparing the real arrival times with
a theoretical model (pulsar timing), it is possible to deduce many physical
information, not only concerning the neutron star and its possible companion,
but also the properties of the interstellar medium as well as tests of General
Relativity. In light of recent works according to which dark energy may have
an astrophysical origin, in this paper we investigate the effect of ’matter’
on the propagation time delay of photons emitted by a pulsar orbiting a
spinning black hole using the rotational Kiselev metric. We deduce an
analytical formula for the time delay from geodesic equations, showing how it
changes as the type of matter around the black hole (radiation, dust or dark
energy) varies with respect to previous results, where matter has not been
taken into account. It turns out that while the spin $a$ only induces a shift
in the phase of the maximum delay without increasing or decreasing the delay,
the effect of matter surrounding the black hole results in a noticeable
alteration of it. Our results show that dark energy would give the strongest
effect and that, interestingly, when the pulsar is positioned between the
observer and the black hole a slightly lower pulse delay than in the no-matter
case appears. We estimated these effects for SGR J1745-2900, the closest
magnetar orbiting SgrA*.

##  I Introduction

Pulsars are the result of the explosion of massive stars showing a repeated
emission of radio waves which we detect as an extremely regular series of
pulses. Since they are remarkably precise clocks, pulsars can be used to
investigate many different aspects of physics, like testing theories of
gravity, studying the magnetic field of the Galaxy and the interior of neutron
stars, investigating the effect of the interstellar medium, and, last but not
least, the detection of the gravitational waves background (GWB). Most
applications of pulsars involve a technique called “pulsar timing”, i.e. the
measurement of the time of arrival (ToA) of photons emitted by the pulsar,
which are then compared with a theoretical model. With a collection of ToAs in
hand, it becomes possible to fit a model of the pulsar’s timing behaviour,
accounting for every rotation of the neutron star. Depending on whether the
pulsar is isolated or binary, the multi-parameter fit gives several important
parameters (so-called ephemeris), like period, period derivative, orbital
period (if binary), position in the sky, eccentricity, ecc. The first hint of
the power of this method was the case of the binary system PSR B1913+16, whose
orbital decay agree with their predicted values to better than 0.5% Taylor and
Weisberg [1]: the observed accumulated shift of periastron is in excellent
agreement with the General theory of Relativity, leaving little room for
alternative theories of gravity. Over the years there have been numerous
studies on possible violations of General Relativity (GR) using the timing of
the pulsars [e.g. 2, 3], revealing that pulsars are a great research tool in
this field, offering also the possibility to test the no-hair theorem as well
as the cosmic censorship conjecture [4, 5].
More recently, pulsar timing is also used to detect and characterize the low-
frequency gravitational wave universe through timing an array of approximately
100 millisecond pulsars (MSP) using the largest radio telescopes in the world.
Indeed, gravitational waves will cause changes in the travel times of pulses
between pulsars and the Earth, detectable as perturbations in pulsar time-of-
arrival measurements. Most importantly, this signature will show a
characteristic sky correlation, predicted by Einstein’s theory of General
Relativity, detectable by correlating the data from all of the pulsars in the
array. Any unmodelled effects will appear in the timing residuals, and the
timing model is revised and/or extended accordingly: trends in the residuals
are indicative of non-optimized parameters, while white noise residuals
suggest a good timing model. Generally, the timing models just requires the
precise modelling of the pulsar’s rotation, orbital motion and the signal’s
propagation in space, and not the details of the radiation’s physics or
emission mechanism. In the case of GWB, the unmodelled residulas would be a
red noise at nHz frequencies [6]. Even if to date there has not been a reveal
of such background, the road looks promising [7, 8].
Usually, the timing model to predict ToAs is based on Damour and Deruelle’s
approach using a post-Newtonian expansion to treat the relativistic two body
problem [9], and its corresponding relativistic effects are described by a set
of post-Keplerian parameters, see e.g. [10]. However, the validity of the
post-Newtonian approximation (that assumes a weak field) it is no longer
guaranteed for a pulsar orbiting closely a supermassive black hole (SMBH), in
particular if pulsar, black hole and observer are (nearly) aligned. When the
weak field approximation holds, then three distinct delay terms arise: the
Roemer delay, the Shapiro delay and the Einstein delay. The Roemer delay is
the difference of arrival time between the case in which the Earth is “on the
same side” of the pulsar and the one in which is on the other side during the
orbit around the sun. In order to avoid the modulation induced by Earth’s
orbit, ToAs are referred to the Solar System barycenter (SSB), where
coordinate time is defined as
$t_{SSB}=t_{em}+(1/c)|\mathbf{r}_{p}-\mathbf{r}_{b}|$, where $t_{em}$ is the
time of photon emission, and $\mathbf{r}_{p,b}$ is the position of the pulsar
or SSB (usually calculated using distant quasars). The Shapiro delay is the
(always positive) additional delay to take into account the deviation of light
caused by the gravitational field of the Solar System and it is easily
obtained by solving the geodetic equations in the weak field approximation, as
the case of the Solar System. Finally, the Einstein delay is due to the
difference between coordinate time $t$ of the pulses and the proper time
$\tau$ of observation, which is different because of a relative motion w.r.t.
the pulsar. In the case of a binary pulsar, there will be a Roemer, Shapiro
and Einstein delays also for the binary system 111In the case of a binary
pulsar, the time of emission is localized at the binary barycenter., with the
difference that, in the case of a binary pulsar, general relativistic effects
are much more important than in the Solar System, being a relativistic system.
This implies that the expressions for the above three different delays become
more complex, since we have to treat a full general relativistic two-body
problem. The two-body problem does not have an exact solution in GR and hence
solving it requires a perturbative post-Newtonian expansion, whose
coefficients (post-Keplerian parameters) can be inferred by pulsar timing.
They are linked not only to the shape of the orbit (like eccentricity and
semi-major axis), but also to other non-Newtonian effects, such as the above
mentioned Shapiro and Einstein delays or the decay of the orbit due to
gravitational radiation labeled with $\dot{P}$. First post-newtonian order
terms are usually sufficient, but if data set is accurate, higher order
contributions can in principle be added in order to get even more information
on the binary system (for a review see e.g. [11]).
Among the various non-Keplerian parameters, recently the effect of dragging
was studied: a compact companion induce a delay in the ToAs due to the frame
dragging on the pulsar’s radiation, if the latter rotates close enough. In
[12], a comparison between an exact analytical result for the frame dragging
time delay and two post-Newtonian derivations was analysed. The exact formula
was derived as the difference between exact geodesic solution for a Kerr black
hole (the compact companion) and the equivalent solution for a Schawarzchild
[13]. The conclusion is that post-Newtonian based treatments overestimate the
frame dragging effect on the lightlike signals, in particular around and after
superior conjunction, hence the analytical solution provide a more reliable
and accurate approach, especially in extreme mass ratio binary configurations,
as the case of a pulsar orbiting Sagittarius A∗ (SgrA∗). In this type of
setup, where the pulsar is considered as a test particle, the effect of the
spin is not negligible and analytical formulas for the propagation delay are
necessary, since fully non-linear numerical approaches suggest that easier
post-Newtonian treatments may not be very accurate [14, 15].
In the wake of these results and since recently evidence has emerged that dark
energy could have an astrophysical origin [16, 17], in this paper we study the
further effect of the presence of matter (radiation, dust or dark energy) in
the surrounding of the central massive black hole. The aim of this work is to
investigate how the propagation time delay is affected by the type of matter,
whose presence is certainly not negligible in the case of SgrA∗. This could
serve, in principle, to understand the environment in which the binary system
is located by exploiting the timing of the pulsars, provided that this effect
is separable from the others. The outline of the paper is as follow. In
Section 2 we derive the equations of motion for lightlike geodesic using a
rotational Kiselev metric, and characterize the corresponding parameter space
$(\lambda,q)$ where $\lambda$ is the adimensional angular momentum and $q$ is
the Carter constant. In Section 3 we analytically solve them using the Mino
time and the elliptic functionals, while in Section 4 we apply the results to
three different type of matter (radiation, dust and dark energy; see later)
and compare them to simpler Kerr case. Finally, we close the manuscript with a
summary and an outlook in Section 5.

##  II Geodesic Equations

###  II.1 The Kiselev metric

Figure 1: Curves, in the parameter space ($\lambda$,$q$), of double roots
outside of the horizons for the radial potential $R(r)$ in the case of
radiation (orange), dust (gray) and dark energy (blue) orbiting a central
black hole, when $a=0.9$ and $c=0.1$ ($M=1$). Kerr black hole (black) with no
matter is also shown. When $q=0$, the corresponding $\lambda$ values,
$\lambda_{\pm}$, on the double roots curve, coincide with the impact parameter
of the unstable circular photon orbits. Above each curve, photons enter from
infinity, reach a turning point outside of the horizons and returns to
infinity (flyby); below each curve, photons move directly between infinity and
the horizon, with no turning points (direct trajectory). The latter behavior
also occurs when $q<0$ (not shown in figure).

The solution of Einstein’s field equation for a Schwarzchild black hole
surrounded by quintessence (a type of dark energy) has been obtained in Ref.
[18, 19]. Even if in these works only quintessence is considered, the Kiselev
solution contemplates any type of energy-matter, once a state parameter has
been established. Indeed, a rotational Kiselev black hole looks like [20].

$\begin{array}[]{ccc}ds^{2}=-\Big{(}1-\frac{2Mr+cr^{1-3\omega}}{\Sigma^{2}}\Big{)}dt^{2}+\frac{\Sigma^{2}}{\Delta}dr^{2}\\\\[6.0pt]
-\frac{2a\sin^{2}\theta\left(2Mr+cr^{1-3\omega}\right)}{\Sigma^{2}}d\varphi
dt+\Sigma^{2}d\theta^{2}\\\\[6.0pt]
+\sin^{2}\theta\left(r^{2}+a^{2}+a^{2}\sin^{2}\theta\frac{2Mr+cr^{1-3\omega}}{\Sigma^{2}}\right)d\varphi^{2}\par\end{array}$
(1)

where we defined

$\Delta=r^{2}-2Mr+a^{2}-cr^{1-3w},\;\;\;\Sigma^{2}=r^{2}+a^{2}\cos^{2}{\theta}.$

$M$ is the mass of the black hole and $a$ is the spin parameter. Moreover, $c$
is the strength parameter and $w$ defines the EoS, $p=w\rho$. Eq. (1) is the
rotational symmetry solution for a black hole wrapped in any kind of energy-
matter definable by the EoS. In general, for dark energy, we would expect
$w<0$. In the following, we will investigate dust ($w=0$) and radiation
($w=1/3$), as well as a dark energy like component with $w=-1/3$ [21].
The number of horizons depends on the value of $w$. For $-1\leq w<-1/3$,
$\Delta=0$ has three positive solutions, corresponding to a Cauchy horizon, an
event horizon and a cosmological horizon. For $w=\pm 1/3$ and $w=0$, the
cosmological horizon disappears and only two horizons, $r_{\pm}$, exist. As an
upper limit on the strength parameter $c$ for the cases we want to study, we
choose $c<1$, compatibly with [22]. We decided not to make a perturbation for
small values $c\ll 1$, since the presence of matter around a supermassive
black hole might not be negligible, as also shown by the image of the SMBH in
the center of M87∗ [23]. Notice, finally, that the drag effect is higher at
higher $|w|$ values, and, when $w$ is fixed, is stronger at higher $c$ values;
when $c\ll 1$, the dependence on $c$ is approximately linear. Therefore, these
effect could combine with the non-zero spin parameter $a$ and create a
degeneracy between different parameters.

###  II.2 Geodesic Equations

Using the metric components $g_{\mu\nu}$ from the line element $ds$, geodesic
equations are given by

$g_{\mu\nu}\dfrac{dx^{\mu}}{d\lambda}\dfrac{dx^{\nu}}{d\lambda}=\epsilon$ (2)

where $\lambda$ is an affine parameter and $\epsilon=0,-1$. While in
Schawarzchild spacetime orbits starting in the plane (for example
$\theta=\pi/2$) remain planar, in Kerr and hence Kiselev metric this is not
true and an additional motion constant is needed. Here, to find equations of
motion (EOM), we adopt the well-known Hamilton-Jacobi equation and Carter
constant separable method 222A similiar treatment has been done for a specific
state parameter value ($w=-2/3$) in [24].  [25]. Therefore, action $S$ must
satisfy:

$\dfrac{\partial S}{\partial\lambda}=\dfrac{1}{2}k\lambda-
Et+L\varphi+S_{r}(r)+S_{\theta}(\theta)$ (3)

where $p_{\alpha}\doteq\partial S/\partial x^{\alpha}$, $L\doteq p_{\varphi}$,
$E\doteq-p_{t}$, $k=0$ for photons and $k=-m^{2}$ for massive particles. With
the inverse metric components

$g^{\mu\nu}=\left(\begin{array}[]{cccc}g^{tt}&0&0&-\dfrac{a(2Mr+cr^{1-3w})}{\Sigma\Delta}\\\
0&\dfrac{\Delta}{\Sigma}&0&0\\\ 0&0&\dfrac{1}{\Sigma}&0\\\
-\dfrac{a(2Mr+cr^{1-3w})}{\Sigma\Delta}&0&0&g^{\varphi\varphi}\end{array}\right)$
(4)

where

$\begin{array}[]{ll}g^{tt}=\dfrac{a^{2}\Delta\sin^{2}{\theta}-(a^{2}+r^{2})^{2}}{\Delta\Sigma}\\\
g^{\varphi\varphi}=\dfrac{\Delta-a^{2}}{\sin^{2}{\theta}}{\Sigma\Delta\sin^{2}{\theta}}.\end{array}$
(5)

After rearranging, Eq. (3) comes down to two separated equations, namely

$\begin{array}[]{lll}E^{2}a^{2}{\cos^{2}{\theta}}+ka^{2}{\cos^{2}{\theta}}-\left(\dfrac{\partial
S}{\partial\theta}\right)^{2}-L^{2}{\cot^{2}{\theta}}=\mathcal{C}\\\
kr^{2}-(aE-L)^{2}-\Delta\left(\dfrac{\partial S}{\partial
r}\right)^{2}+\dfrac{1}{\Delta}\left[E(r^{2}+a^{2})-aL\right]^{2}=\mathcal{C}\end{array}$
(6)

where $\mathcal{C}$ is a separation constant. From definitions of energy $E$
and angular momentum $L$, the following relations hold

$\begin{array}[]{ll}\dot{t}=-\frac{E}{g_{tt}}\left[1+g_{t\varphi}\left(\lambda+\frac{g_{t\varphi}}{g_{tt}}\right)\left(g_{\varphi\varphi}-\frac{g_{t\varphi}^{2}}{g_{tt}}\right)^{-1}\right]\\\
\dot{\varphi}=E\left(\lambda+\frac{g_{t\varphi}}{g_{tt}}\right)\left(g_{\varphi\varphi}-\frac{g_{t\varphi}^{2}}{g_{tt}}\right)^{-1}\,\
,\end{array}$ (7)

where $\lambda\doteq L/E$. Hence, the geodesic equations are

$\dot{t}=\dfrac{(r^{2}+a^{2})(r^{2}+a^{2}-a\lambda)}{\Delta}-a(a-\lambda)+a^{2}\cos^{2}{\theta}$
(8)
$\dot{\varphi}=\dfrac{a(r^{2}+a^{2}-a\lambda)}{\Delta}-a+\dfrac{\lambda}{\sin^{2}{\theta}}$
(9)
$\dot{\theta}^{2}=q+\cos^{2}{\theta}\left[\left(1+\dfrac{k}{E^{2}}\right)a^{2}-\dfrac{\lambda^{2}}{\sin^{2}{\theta}}\right]=\Theta(\theta)$
(10)
$\dot{r}^{2}=-\Delta\left[q-\dfrac{k}{E^{2}}r^{2}+(\lambda-a)^{2}\right]+\left(r^{2}+a^{2}-\lambda
a\right)^{2}=R(r).$ (11)

Here, a dot means derivative w.r.t. the so called Mino time $\gamma$, which
satisfies the condition $dx^{\mu}=(\Sigma/E)p^{\mu}d\gamma$, while
$q\doteq\mathcal{C}/E^{2}$. Quantities $m$, $E$, $L$, $\mathcal{C}$ are
constants of motion for the EOMs (8)-(11). Notice that when $c=0$ Eqs.(8)-(11)
are equal to the Kerr case [12], although they formally remain similar even
when $c$ is non-zero. Furthermore, as for a Kerr metric, $\Theta(\pi/2)=0$
$\longleftrightarrow$ $q=0$, i.e. a geodetic lies entirely in the equatorial
plane if and only if $q=0$. One difference is the passage for a point with
$r=0$ and $\theta\not=\pi/2$ when $w=1/3$. Indeed, the condition $R(0)\geq 0$
for $w=1/3$ implies

$q\leq\dfrac{c(a-\lambda)^{2}}{(a^{2}-c)}\geq 0$ (12)

where we assumed $c<<1$. In this regime, $r=0$ is approachable for both
negative and positive values of $q$, while for $w\in\\{-1,-2/3,-1/3,0\\}$ this
happens only for negative values ($q<0$), as in Kerr metric. On the other
hand, from the positivity condition $\Theta(\theta)\geq 0$, motion is allowed
only when [26]

$q\geq\left\\{\begin{array}[]{ll}0&|\lambda|\geq a\\\
-(|\lambda|-a)^{2}&|\lambda|\leq a\end{array}\right.$ (13)

which border an allowed region for the parameter space $(q,\lambda)$ and is
independent of the new parameter $c$.

###  II.3 Roots of radial potential

The calculation and classification of the roots of the radial potential $R(r)$
proceed in a similar way to the Kerr case [26]. We consider here the most
interesting cases for our scopes, i.e. $w\in\\{\pm 1/3,0\\}$.
When $w=-1/3$, the solution of equation $R(r)=0$ is given by

$\displaystyle r_{1}$
$\displaystyle=-z-\sqrt{-\frac{A}{2}-z^{2}+\frac{B}{4}z^{2}}$ (14)
$\displaystyle r_{2}$ $\displaystyle=-z+\sqrt{-\frac{A}{2}-z^{2}+\frac{B}{4}}$
$\displaystyle r_{3}$ $\displaystyle=z-\sqrt{-\frac{A}{2}-z^{2}-\frac{B}{4z}}$
$\displaystyle r_{4}$ $\displaystyle=z+\sqrt{-\frac{A}{2}-z^{2}-\frac{B}{4z}}$

where

$\displaystyle A$ $\displaystyle=a^{2}-q-\lambda^{2}+c\chi,\;\;\;\;\;B=2M\chi$
$\displaystyle C$
$\displaystyle=-a^{2}q,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\chi=q+(\lambda-a)^{2}\geq
0.$

A very similarly computation is required for the remaining two values of $w$
that we consider, hence we omit it for brevity. As we expect, the discrepancy
between the Kerr and Kiselev radial potential roots increases as the strength
parameter $c$ increases. The effect is more evident for negative values of
$w$. When $w=0$, the additional term appearing always competes with the mass
term $Mr$, being therefore almost negligible.
In order to find quadruple roots of the potential, we impose the form
$R(r)=(r-r_{0})^{4}$ and after comparing with its explicit expression, we
obtain

$q=0,\;\;\;\;\lambda=a$ (15)

besides that $r_{0}=0$. Here, we assumed $a\not=0$. Triple solutions satisfy
the relations $R(r)=R^{\prime}(r)=R^{\prime\prime}(r)=0$, whose solution is

$q=\dfrac{8r^{3}}{2M}-(\lambda-a)^{2},\;\;\;\;\lambda=\dfrac{M(a^{2}+3r^{2})-2r^{3}(1-c)}{aM}$
(16)

and hence

$r=\dfrac{M-\left[M(M^{2}-(1-c)a^{2})\right]^{1/3}}{(1-c)}.$ (17)

Finally, double solutions ($R^{\prime}(r)=R(r)=0$) occur when

$\displaystyle q$
$\displaystyle=\dfrac{r^{3}}{a^{2}}\left[\dfrac{4M\Delta}{\left(M-r(1-c)\right)^{2}}-r\right],$
(18) $\displaystyle\lambda$
$\displaystyle=a+\dfrac{r}{a}\left[r-\dfrac{2\Delta}{r(1-c)-M}\right].$

We then evaluate these relations on the border of the admissible region (13),
i.e. $q=0$ and $q=-(\lambda\pm a)^{2}$. In the first case, three real
solutions for $r$ are possible:

$r=e\cos\left[\dfrac{1}{3}\arccos\left(\dfrac{2a^{2}(1-c)-M^{2}}{M^{2}}\right)+\dfrac{2\pi
k}{3}\right]+e$ (19)

with $k=0,1,2$ and $e=2M/(1-c)$. When $q=-(\lambda+a)^{2}$, four real and
distinct solutions for $r$ are allowed; their structure is similar to (14) and
we omit them for brevity. Finally, in the last case, $q=-(\lambda-a)^{2}$, in
addition to the Kerr case analogous solution, i.e. $r_{\pm}$, other options
appear, namely $r_{N}=\pm a/\sqrt{c}$. Since double solution curves separates
the parametric $(q,\lambda)$ regions with two or four real solutions, once the
values of mass $M$, spin $a$ and strength of matter $c$ have been fixed,
different behaviours of the roots are delimited by such curves, which, indeed,
border the range of radial coordinate $r$ for which (13) holds, i.e. motion is
allowed. In particular, when $q=0$, the corresponding $\lambda$ values on the
double roots curve, $\lambda_{\pm}$, coincide with the impact parameters of
the so-called unstable circular photon orbits (counter-rotating and co-
rotating). Radii of such (equatorial) orbits are simply given by the largest
roots (19), and in the case $w=-1/3$ are

$r_{ph}^{\pm}=e\cos\left[\dfrac{1}{3}\arccos\left(\dfrac{2a^{2}(1-c)-M^{2}}{M^{2}}\right)+k_{\pm}\dfrac{4\pi}{3}\right]+e$

where $k_{+}=0$, $k_{-}=2$. When $c=0$, they reduce to the well-known Kerr
results. Notice that with this notation $\lambda_{+}\equiv\lambda(r_{ph}^{+})$
and $\lambda_{-}\equiv\lambda(r_{ph}^{-})$, it will be
$\lambda_{+}<\lambda_{-}$. A possible trend for the three different types of
matter is showed in Fig. (1), where each point of the curves corresponds to a
choice of parameters ($\lambda$,$q$) for which two roots for the potential
$R(r)$ coincide ($r_{3}=r_{4}>r_{+}$). Above the curve, four real roots of the
potential are allowed, two of which outside of the horizons. Conversely, below
the curve ($\lambda_{+}<\lambda<\lambda_{-}$ when $q=0$) no real roots outside
of the horizons exist. Since $R(r)>0$ at infinity, radial potential is
positive in the ranges $r<r_{1}$, $r_{2}<r<r_{3}$ and $r>r_{4}$. Therefore,
only in the first case photons enter from infinity, reach a turning point
outside of the horizons (at $r=r_{4}$), and returns to infinity 333In addition
to this motion (flyby), for parameter space points above the curve a bound
orbit is also possible, but that is not interesting for our purposes.; in the
second case, photons move directly between infinity and the horizon (there may
be turning points inside the horizon). The presence of matter also reduces the
parametric region corresponding to four real roots all inside the horizons
(not shown in Fig. (1)). As we expect, the effect is more evident for larger
value of $c$, but it is more pronounced in the presence of radiation rather
than dark energy.
In what follows, we will deal with a beam of photons from a pulsar orbiting a
black hole in the equatorial plane. In this case, the trajectory followed can
be of only two types: either the photons move directly from the pulsar to the
observer, or they first encounter a turning point outside of the horizons
($r_{4}$) and then move from it to the observer. The parametric region
corresponding to the first case (direct trajectory) is made up of all points
below the curve (even for $q<0$), while the second case (flyby) describes
trajectories for photons with motion parameters above the curve. Overall, the
presence of matter attached to a rotating black hole increases the region of
parameter space in which a direct trajectory towards a distant observer is
possible, at the expense of a flyby. More precisely, it is easy to find out
the behaviour of the emitted photons from considerations on the potential
shape and its positivity ranges [26].
Computation of quadruple, triple, and double roots for the remaining two cases
$w=0$, $w=1/3$ are reported in Appendix A.

##  III Time Delay

In this section we analytically solve Eqs. (8)-(11) using elliptic integrals
in the Legendre form, following the strategy used in [12]. As in the previous
section, for brevity we show here only the case $w=-1/3$. Final results for
different cases are however summarized in Appendix B.

###  III.1 The exact time delay

Combining Eqs. (8), (10) and (11), one get an equation for $t$ in integral
form, i.e.

$\mathrm{c}(t_{a}-t_{e})=\int_{\gamma_{r}}\dfrac{G(r)}{\Delta\sqrt{R(r)}}dr+\int_{\gamma_{\theta}}\dfrac{a^{2}\cos^{2}{\theta}}{\sqrt{\Theta(\theta)}}d\theta$
(20)

where $\mathrm{c}$ is the speed of light (not to be confused with $c$) and we
have defined

$G(r)=r^{2}(r^{2}+a^{2}+ac(a-\lambda))+2Mra(a-\lambda).$ (21)

We notice that the term proportional to $c$ is the novelty with respect to the
Kerr case (as well as the different definition of $\Delta$). On the other
hand, the angular integral remains unchanged. The integral path $\gamma_{r}$
starts at the radial point of emission $r_{e}$ and either runs directly to
infinity (direct trajectory) or first decreases in radius towards a turning
point outside of the horizons ($r_{4}$) and then return to infinity (flyby),
according to the motion parameters of photons (see Fig. (1)). Therefore, we
split radial integral as

$\left(\int_{r_{4}}^{\infty}\pm\int_{r_{4}}^{r_{e}}\right)\dfrac{G(r)}{\Delta\sqrt{R(r)}}dr$
(22)

where we choose the minus sign for a direct trajectory 444Notice that, in this
case, $r_{4}$ is just a generic point between the emission and the observer
position, which we fixed at infinity, i.e. $r_{a}=\infty$. Therefore
$r_{e}<r_{4}<\infty$. and the plus sign for a flyby motion. Similarly, the
angular integral in Eq. (20) can be written as [12]

$\left(\int_{0}^{u_{e}}\pm\int_{0}^{u_{a}}\right)\dfrac{a^{2}u}{2\sqrt{U(u)}}du$
(23)

where we defined

$U(u)=u(q+u(a^{2}-\lambda^{2}-q)-a^{2}u^{2})$

and we changed the integration variable to $u=\cos^{2}\theta$. Hence, $u_{e}$
and $u_{a}$ represent the emission and the observer latitudinal positions,
respectively. Here, we choose the plus sign if the equatorial plane is crossed
and the negative sign else 555We adopt the usual convention that $\theta=0$
corresponds to the north pole. In this way, $\theta=\pi/2$ correspond to the
equatorial plane.. The above equation strictly only holds in absence of
latitudinal turning points $u_{\pm}$ (the non-zero roots of $U(u)$). If
latitudinal turning points are encountered, we have to add complete integrals
in the form $\int_{0}^{u_{+}}$ if $q>0$ or $\int_{u_{-}}^{u_{+}}$ if $q<0$.
However, below we will assume that we are not in these cases.
The integral (22) can then be solved exactly in terms of elliptic integrals.
We notice that similar expressions for the time evolution of lightlike
geodesics in Kerr metric have been derived before, in slightly different ways
than the one used here and in [12]. For example, [27] give an expression using
partly Carlson’s elliptic integrals. Expressions in terms of Weierstrass
functions have been derived in [13], while in [26] Jacobi elliptic integrals
in ’Jacobi form’ (instead of Legendre form employed here) were used. Here, we
are not concerned with the optimal choice; what is certain is that having an
analytical solution allows to avoid divergences which would be impossible to
avert in a purely numerical calculation. The result is (details in Appendix B)

$\begin{array}[]{ccc}\mathrm{c}(t_{a}-t_{e})=&T_{r}(\infty,\lambda_{e},q_{e})\pm
T_{r}(r_{e},\lambda_{e},q_{e})\\\ &+|T_{u}(u_{e},\lambda_{e},q_{e})\pm
T_{u}(u_{a},\lambda_{e},q_{e})|\end{array}$ (24)

with the definitions

$\begin{array}[]{l}T_{u}(u,\lambda,q)=\int_{0}^{u}\frac{a^{2}}{2}\frac{udu}{\sqrt{U_{\lambda,q}}}\\\
=\frac{a}{\sqrt{u_{+}-u_{-}}}\left[u_{-}F(v,w)+\left(u_{+}-u_{-}\right)E(v,w)-\frac{u_{+}v\sqrt{\mathrm{1}-v^{2}}}{\sqrt{1-w^{2}v^{2}}}\right]\end{array}$
(25)

and

$\begin{array}[]{l}T_{r}(r,\lambda,q;c)=\delta\cdot\Big{[}F(x,k)\cdot\Big{(}4M^{2}\gamma^{-3}-a^{2}c^{2}\gamma^{-2}-ac\lambda\gamma^{-1}\\\
+2Mr_{3}\gamma^{-2}+\frac{1}{2}\gamma^{-1}\cdot\big{[}r_{1}\left(r_{3}-r_{4}\right)+r_{3}\left(r_{3}+r_{4}\right)\big{]}+\frac{B_{+}l}{l_{+}}+\frac{B_{-}l}{l_{-}}\Big{)}\\\
+E\big{(}x,k\big{)}\cdot\Big{(}-\frac{1}{2}\gamma^{-1}\cdot\big{(}r_{4}-r_{2}\big{)}\big{(}r_{3}-r_{1}\big{)}\Big{)}\\\
+\Pi(x,l,k)\cdot\Big{(}2Mr_{4}\gamma^{-2}-2Mr_{3}\gamma^{-2}\Big{)}\\\
+\Pi(x,l_{+},k)\cdot\Big{(}B_{+}-\frac{lB_{+}}{l_{+}}\Big{)}+\Pi\big{(}x,l_{-},k\big{)}\cdot\Big{(}B_{-}-\frac{lB_{-}}{l_{-}}\Big{)}\Big{]}\\\
+\frac{\sqrt{R(r)}}{r-r_{3}}\gamma^{-1}\\\ \end{array}$ (26)

where $r_{1,..,4}$ are the radial potential roots (13), $u_{\pm}$ the non-zero
roots of $U,$ $r_{\pm}$ are the horizons, and for brevity, we defined (for
$v$, $w$, $x$, $k$ and $B_{\pm}$ see Appendix B)

$\delta=\frac{2}{\sqrt{\left(r_{4}-r_{2}\right)\left(r_{3}-r_{1}\right)}},\;\;\;\;\;\;\;\gamma=1-c,$
$l=\dfrac{r_{1}-r_{4}}{r_{1}-r_{3}},\;\;\;\;\;\;\;l_{\pm}=\dfrac{l(r_{3}-r_{\pm})}{r_{4}-r_{\pm}}.$

The functions $F$, $E$ and $\Pi$ appearing in Eq. (26) are the well-known
elliptic functions of first, second and third kind, respectively (see again
Appendix B). To obtain Eq. (26) we also used the vanishing rule
$\sum_{i=1}^{4}r_{i}=0$. The above equations reduce to Kerr ones [12] when
$c\rightarrow 0$. As in Kerr, some divergences appear in Eq.(26), specifically
in $\Pi(x,l,k)$ and in the last term, which is not clear how to handle in
fully numerical calculations. In Eq. (24), the constants of motion $\lambda$
and $q$ only depend on the emission point, since we have fixed the arrival
point at infinity. Furthermore, in the case of equatorial orbits
($\theta=\pi/2$) the angular integral in Eq. (20), and hence $T_{u}$ in Eq.
(25), are vanishing.
To avoid divergences, it is usual to subtract the time delay w.r.t. a fixed
reference point from the actual time delay (24), namely

$\begin{array}[]{l}\Delta
t_{ex}\left(t_{e},\varphi_{e},u_{e};c\right)=\left(t_{a}-t_{e}\right)-\left(t_{a}-t_{\text{ref
}}\right)\\\
=\frac{1}{\mathrm{c}}\Big{[}T_{r}\left(\infty,\lambda_{e},q_{e}\right)\pm
T_{r}\left(r_{e},\lambda_{e},q_{e}\right)\\\
+|T_{u}\left(u_{e},\lambda_{e},q_{e}\right)\pm
T_{u}\left(u_{a},\lambda_{e},q_{e}\right)|\Big{]}\\\
-\frac{1}{\mathrm{c}}\Big{[}T_{r}\left(\infty,\lambda_{\text{ref
}},q_{\text{ref }}\right)\pm T_{r}(r_{\text{ref }},\lambda_{\text{ref
}},q_{\text{ref }})\\\ +|T_{u}\left(u_{\text{ref }},\lambda_{\text{ref
}},q_{\text{ref }}\right)\pm T_{u}\left(u_{a},\lambda_{\text{ref
}},q_{\text{ref }}\right)|\Big{]},\end{array}$ (27)

where $\lambda_{ref},q_{ref}$ are the angular momentum and Carter constants at
the reference point and $c$ in the r.h.s. means that the expression is
evaluated in presence of matter 666In the following, we simply denote Eq. (27)
with $\Delta t_{ex}(c)$. This implies that $\Delta t_{ex}(c=0)$ reproduce same
results of [12], i.e. in absence of surrounding matter.. The ascending node
w.r.t. the plane of the sky ($\phi=-w$) is used as the reference point, since
the time delay is zero for photons leaving the pulsar in such a position.
Actually, the Roemer, first-order Shapiro (when $e=0$) and geometric delays
all vanish at the ascending node ($\phi=-w$), but other delays, like Einstein
delay, second-order Shapiro delay, as well as the exact formula (27) do not
vanish at the ascending node, and the addition of individual offset to each
last types of delay is usually needed to have a vanishing delay point.

###  III.2 Orbital parameters

In order to fully calculate Eq. (24), we need the coordinates of the emission
point ($r_{e},\theta_{e},\varphi_{e},t_{e}$) on the pulsar orbit as well as
the observer latitude $\theta_{a}$ ($r_{a}=\infty$ and we assume, for
simplicity, $\varphi_{a}=0$). The position of the emission point follows the
pulsar’s orbit around the black hole (see Fig. (2)); for each point of the
orbit, we need the motion parameters ($\lambda_{e},q_{e}$) of geodesics
connecting to points to infinity. As pointed out by [12], there is no general
analytical solution to such problem (emitter-observer problem). To obtain
$\lambda$ and $q$, in general, one needs to numerically solve both the
equations

$\int_{\gamma_{r}}\dfrac{dr}{\sqrt{R(r)}}=\int_{\gamma_{\theta}}\dfrac{d\theta}{\sqrt{\Theta(\theta)}},$
(28)
$\varphi_{a}-\varphi_{e}=\int_{r}\dfrac{2Mra-a^{2}\lambda+acr^{2}}{\Delta\sqrt{R(r)}}dr+\int_{\gamma_{\theta}}\dfrac{\lambda}{\sin^{2}(\theta)\sqrt{\Theta(\theta)}}d\theta$
(29)

where the first integral in Eq. (28) is still convergent, since given any
$r_{4}<s<\infty$ (no problems for $r\rightarrow\infty$), then

$\int_{r_{4}}^{s}\dfrac{2Mra-a^{2}\lambda+acr^{2}}{\Delta\sqrt{R(r)}}dr\leq\delta
MF\big{(}x_{s},k\big{)}<\infty$

where $M$ is the maximum of the function $(2Mra-a^{2}\lambda+acr^{2})/\Delta$
on the range of the integral and $x_{s}$ is the variable $x$ (see Appendix B)
evaluated in $r=s$. Proceeding similarly to what was done for Eq. (26) and
assuming $\varphi_{a}=0$ and $q=0$ for both pulsar and observer (pulsar, black
hole and observer lay on the equatorial plane, from Eq. (29) we get

$\begin{array}[]{l}-\varphi_{e}=\dfrac{ac}{\gamma}\delta\cdot\left[F(\infty,k)\pm
F(r_{e},k)\right]\\\ +\dfrac{\delta
a(a^{2}c-2Mr_{+}+a\lambda\gamma)}{\gamma^{2}(r_{-}-r_{+})(r_{4}-r_{+})}\Big{[}\Big{(}1-\dfrac{l}{l_{+}}\Big{)}\Pi(\infty,l_{+},k)\\\
+\Big{(}\pm
1\mp\dfrac{l}{l_{+}}\Big{)}\Pi(r_{e},l_{+},k)+\dfrac{l}{l_{+}}F(\infty,k)\pm\dfrac{l}{l_{+}}F(r_{e},k)\Big{]}\\\
+\dfrac{\delta
a\Big{(}2Mr_{-}-a^{2}c-a\lambda\gamma\Big{)}}{\gamma^{2}(r_{-}-r_{+})(r_{4}-r_{-})}\Big{[}\Big{(}1-\dfrac{l}{l_{-}}\Big{)}\Pi(\infty,l_{-},k)\\\
+\Big{(}\pm
1\mp\dfrac{l}{l_{-}}\Big{)}\Pi(r_{e},l_{-},k)+\dfrac{l}{l_{-}}F(\infty,k)\pm\dfrac{l}{l_{-}}F(r_{e},k)\Big{]}\\\
+\dfrac{2\lambda[F(\infty,k)\pm
F(r_{e},k)]}{\sqrt{(r_{4}-r_{2})(r_{3}-r_{1})}}\end{array}$ (30)

where we choose upper or lower sign for flyby or direct trajectory,
respectively. We also note that, in general, the angular integral in Eq. (29)
can be computed using

$\int_{0}^{u_{e,a}}\dfrac{\lambda}{\sin^{2}(\theta)\sqrt{\Theta(\theta)}}d\theta=\pm\dfrac{\lambda}{\sqrt{u_{-}}}\Pi\left(\sqrt{\dfrac{u_{e,a}}{u_{+}}},u_{+},\sqrt{\dfrac{u_{+}}{u_{-}}}\right),$

but, in the assumption of equatorial orbit (i.e. $\theta=\pi/2$) and using Eq.
(28), it simply reduces to a radial integral

$\int_{\gamma_{\theta}}\dfrac{\lambda}{\sin^{2}(\theta)\sqrt{\Theta(\theta)}}d\theta=\int_{\gamma_{r}}\dfrac{\lambda}{\sqrt{R(r)}_{q=0}}dr$

and this gives the last term in Eq. (30).
Generally, the particular case in which everything is restricted to the
equatorial plane is simpler and, at the same time, more interesting, as this
case corresponds to the strongest relativistic effects. Before schematizing
the procedure we have followed, let us first explain how to relate metric
coordinates to the coordinates of a pulsar orbiting a black hole. To this end,
we will adopt the geometrical set-up already introduced in [13]. Due to the
big difference in mass between pulsar and (supermassive) black hole, we may
consider the pulsar as a test particle, i.e. the center of mass coincide with
the BH center. Of course, the pulsar will not remain, in general, in a fixed
plane due to the frame dragging, and that the same orbit would rotate because
of the relativistic precession of the periapsis. Nevertheless, we will assume
a keplerian orbit in first approximation, with (almost all) relativistic
effects encoded in post-Newtonian orbital parameters, as usual in pulsar
timing models. However, here we are not interested in individual post-
Keplerian effects, since Eq. (27) is a full relativistic formula for time
delay; what we need to do is just to compare formula (27) with the analogue
one without the surrounding matter ($c=0$, i.e. Kerr), in order to highlight
the effect of the latter. This is a difference with [12], where the
relativistic effect of frame dragging has been studied and compared to weak
field post Newtonian approximations.

Figure 2: Orientation of the black hole-pulsar system in the sky with respect
to an observer sitting at infinity. The blue shaded region around the black
hole reflects the possible environment of dark energy, radiation or dust
investigated in this work.

In this case we can express the coordinate $(x,y,z)$ as follows:
$x=r_{e}\cos(\omega+\phi)$, $y=r_{e}\sin(\omega+\phi)$, $z=0$, where $\omega$
is the argument of the periastron. A rotation around the $x$-axis by the
inclination angle $i$ suffices to transform to the $(X,Y,Z)$ system. For the
case of an edge-on equatorial pulsar orbit, that we will discuss later in the
paper, the desired angle between pulsar and observer is then given by the
angle $\vartheta$ in spherical coordinates $X=r\cos\psi\sin\vartheta$,
$Y=r\cos\psi\sin\vartheta$, $Z=r\cos\vartheta$. In the common plane of pulsar
and observer, the angle $\phi_{e}$ is then determined by
$\varphi_{e}=\vartheta$ with $\cos\vartheta=-\sin\textit{i}\sin(\omega+\phi)$
and therefore

$\cos\varphi_{e}=-\sin\textit{i}\sin(\omega+\phi)\,.$ (31)

##  IV Results

First of all, we assume an extreme binary system of a pulsar orbiting a
supermassive black hole with a mass of $M=4\times 10^{6}M_{\odot}$ (solar
masses), i.e. $GM_{\odot}/c^{2}=1476M$. Therefore, the propagation time delay
will be expressed in seconds and the corresponding adimensional value can be
recovered by dividing by a factor $GM/c^{3}\approx 19.7s$. Different black
hole masses (say $M_{2}$) lead to time delays multiplied by a factor
$M_{2}/M$.

First, we notice that by choosing $c\approx 0$, the effects of matter are
negligible and we are able to reproduce the results of [12] from our Eq. (27).
Just to give an example, we show a single plot for this case in Fig. (3),
where we also show the cases of a rotating and a not-rotating Kiselev black
holes in presence of dark energy. It turns out that while the spin $a$ induces
a shift in the phase of the maximum delay without increasing or decreasing the
delay, the effect of matter surrounding the black hole results in a noticeable
alteration of the delay, which increases as the strength parameter $c$
increases. Interestingly, when the pulsar is positioned between the observer
and the black hole ($\theta\approx 0$) we predict a slightly lower pulse delay
than in the no-matter case. Difference with Schwarzchild and Kerr cases are
less evident (but still appreciable) in presence of dust (see Fig. (4)) and
practically imperceptible in the case of pure radiation (see Fig. (5)).

We also may isolate the only effect of matter by making the subtraction

$\Delta t_{matter}=\Delta t_{ex}(c)-\Delta t_{ex}(c=0)$ (32)

where $\Delta t_{ex}(c=0)$ is obtained by putting $c=0$ in Eq. (27). Here a
comment is necessary. If one identifies $\Delta t_{ex}(c=0)$ in Eq. (32) with
the exact time delay in Kerr metric, then, in order to compare results derived
in different space-times (Kiselev vs Kerr), we first need to identify a
physical invariant. The common idea is to fix the circumference of a circle,
which is an invariant characteristic. In particular, in Kerr space-time, such
a circumference is given by
$\mathcal{C}_{kerr}=2\pi\sqrt{r^{2}_{kerr}+a^{2}+2Ma^{2}/r_{kerr}}$, while our
metric (1) leads to the (equatorial) circumference

$\mathcal{C}_{kis}=2\pi\sqrt{r^{2}_{kis}+a^{2}+a^{2}\left(\dfrac{2M}{r_{kis}}+cr_{kis}^{-1-3w}\right)}.$
(33)

Then the equality is achieved with

$r_{kerr}=2\sqrt{-Q}\cos\left(\dfrac{\theta}{3}\right)$ (34)

where $Q$ depends on the type of matter

$Q_{de}=-\dfrac{r^{3}_{kis}+a^{2}(2M+cr_{kis})}{3r_{kis}},\;\;\;\;Q_{dust}=-\dfrac{r^{3}_{kis}+a^{2}(2M+c)}{3r_{kis}}$
$Q_{rad}=-\dfrac{r^{4}_{kis}+a^{2}(2Mr_{kis}+c)}{3r^{2}_{kis}},\;\;\theta=\arccos\Big{(}-\dfrac{Ma^{2}}{\sqrt{-Q^{3}}}\Big{)}.$

However, to make things easier, one could decide to use only Eq. (27) to
calculate the matter delay (32), for both the contributions $\Delta
t_{ex}(c)$, $\Delta t_{ex}(c=0)$. This is also justified by the fact that Eq.
(34) gives very convergent $r_{kerr}$, $r_{kis}$ values as long as $c$ is far
from 1, as we are assuming.

###  IV.1 Effect of the black hole environment on the propagation time delay

Figure 3: The exact time delay $\Delta t_{\text{ex}}-\Delta t_{\text{R}}$ for
a Schwazschild black hole (black line), a Kerr black hole with $a=0.9$ (dashed
blue line), a non-rotating Kiselev black hole surrounded by dark energy with
$c=0.01$ (dashed green line) and a rotating one with $a=0.9$ and $c=0.01$ (red
line) for a circular edge-on orbit with a Schwarzschild radius
$r_{\text{S}}=100\,M$. Bottom figures are a zoom of the top. At this distance,
the additional time delay when the pulsar is behind the BH (superior
conjunction) is more than 30 seconds, while it is reduced by a few seconds
when the pulsar is in front of the BH (inferior conjunction).

Given the Roemer delay $\Delta t_{\text{R}}$ which is given by [28] :

$\Delta
t_{\text{R}}=\frac{A(1-\textit{e}^{2})\sin\textit{i}\sin(\omega+\phi)}{c(1+\textit{e}\cos\phi)}$
(35)

where i is the inclination of the orbital plane with respect to the plane of
sky and $\omega$ is the argument of periapsis, $\phi$ is the argument of the
pulsar’s position and e is the eccentricity of the orbit, we show in figure 3,
the time delay $\Delta t_{\text{ex}}-\Delta t_{\text{R}}$ (see equation 27)
for a Schwarzschild and a Kerr black hole both in the case where they are in a
vacuum as well as surrounded by dark energy with $c=0.01$. We choose a simple
pulsar trajectory of an edge-on ($\mathit{i}=\pi/2$) circular orbit. The
ascending node with respect to the plane of sky is used as the referenece
point i.e $\varphi_{\text{ref}}=\pi/2$, which with $\omega=-\pi/2$ simplifies
to $\varphi_{\text{ref}}=\phi_{\text{ref}}=\pi/2$. Note that all exact
propagation time delays $\Delta t_{\text{ex}}$ includes considerable offset.
Therefore we correct them by adding global constants to the individual delays
such that they exactly vanish at $\phi_{\text{ref}}=\pi/2$.
As expected for a circular edge-on orbit in spherical symmetry, the curve
representing the Schwarzschild black hole ($a=0,c=0$) is symmetric with
respect to the superior conjunction at $\phi=\pi$. Once the black hole rotates
this symmetry is broken with a slight shift in the top of the curve which
corresponds to the switch from contra-rotating to co-rotating lightlike
geodesics emitted from the pulsar to the observer which was already discussed
in details in [12]. However, if the dark energy is present ($w=-1/3$), a
considerably additional amount of time delay (more than 30 seconds if the
pulsar is at a distance $r_{s}=100$ M) is induced if the pulsar comes behind
the black hole along its trajectory. On the other side, dark energy
contributes negatively to the time delay (with an advance of a few seconds)
once the pulsar is at the front of the observer. This behaviour is seen both
for a Schwarzschild and a Kerr black hole as it is illustrated in the bottom
figure 3.

Surprisingly, in the case the dust is present at the vicinity of the black
hole ($w=0$), the same above feature is observed as in the case of dark energy
but with a very small amplitude (see figure 4). This finding appears to
contradict the expected behaviour, in which the dust would induce a constant
time delay for a circular pulsar orbit, regardless of its position with
respect to the black hole.

Figure 4: The exact time delay $\Delta t_{\text{ex}}-\Delta t_{\text{R}}$ for
a Schwarzschild black hole (black line), a Kerr black hole with $a=0.9$
(dashed blue line), a non-rotating Kiselev black hole surrounded by dust with
$c=0.01$ (dashed green line) and a rotating one with $a=0.9$ and $c=0.01$ (red
line) for a circular edge-on orbit with a Schwarzschild radius
$r_{\text{S}}=100\,M$. Bottom figures are a zoom of the top.

Figure 5: The exact time delay $\Delta t_{\text{ex}}-\Delta t_{\text{R}}$ for
a Schwazschild black hole (black line), a Kerr black hole with $a=0.9$ (dashed
blue line), a non-rotating Kiselev black hole surrounded by radiation with
$c=0.01$ (dashed green line) and a rotating one with $a=0.9$ and $c=0.01$ (red
line) for a circular edge-on orbit with a Schwarzschild radius
$r_{\text{S}}=100\,M$. Bottom figures are a zoom of the top.

In figure 5, the propagation time delay of light geodesics emitted from the
pulsar in the presence of radiation ($w=1/3$) is investigated. As it is
expected, the radiation contributed positively to the propagation time delay
along the pulsar trajectory around the black hole. However this contribution
is maximal once the pulsar is at a superior conjecture and minimal at the
inferior conjecture.

##  V Conclusions

Most applications of pulsars involve a technique called “pulsar timing”, i.e.
the measurement of the time of arrival (ToA) of photons emitted by the pulsar,
which are then compared with a theoretical model. With a collection of ToAs in
hand, it becomes possible to fit a model of the pulsar’s timing behaviour,
accounting for every rotation of the neutron star. Depending on whether the
pulsar is isolated or binary, the multi-parameter fit gives several important
parameters (so-called ephemeris), like period, period derivative, orbital
period (if binary), position in the sky, eccentricity, ecc. Over the years
there have been numerous studies on possible violations of General Relativity
(GR) using the timing of the pulsars [e.g. 2, 3], revealing that pulsars are a
great research tool in this field, offering also the possibility to test the
no-hair theorem as well as the cosmic censorship conjecture [4, 5]. More
recently, pulsar timing is also used to detect and characterize the low-
frequency gravitational wave universe through timing an array of approximately
100 millisecond pulsars (MSP) using the largest radio telescopes in the world.
Usually, the timing model to predict ToAs is based on Damour and Deruelle’s
approach using a post-Newtonian expansion to treat the relativistic two body
problem [9], and its corresponding relativistic effects are described by a set
of post-Keplerian parameters. However, the validity of the post-Newtonian
approximation (that assumes a weak field) it is no longer guaranteed for a
pulsar orbiting closely a SMBH, in particular if pulsar, black hole and
observer are (nearly) aligned. This means that full general relativistic
computations are sometimes necessary.
In this paper we study, from an analytical point of view the further effect of
the presence of matter (radiation, dust, dark energy) in the surrounding of
the central massive black hole. Hence, the aim of this work was to investigate
if and how the propagation time delay is affected by different types of
matter, whose presence is certainly not negligible in the case of SgrA∗. We
have modeled such a presence with the (rotational) Kiselev metric, which in
addition to the spin parameter $a$, depends on the strength parameter $c$, as
well as on the state parameter $w$ depending on the type of matter ($w=\pm
1/3,0$). Recently, evidence has emerged that black holes could contribute to
the dark energy content of the universe, resulting cosmologically coupled [16,
17]; therefore considering the possibility of dark energy around SMBHs could
give clues in this direction.
We first computed geodesic equations in a general, rotational Kiselev metric
and characterized the corresponding parametr space $(\lambda,q$), where
$\lambda$ is the (adimensional) angular momentum and $q$ is the (adimensional)
Carter constant, finding that generally the presence of matter enlarge the
region of direct trajectories of photons moving from the pulsar and a far
observer (see Fig. 1), at the expense of a flyby. This effect is more evident
for "dark energy" and less evident for radiation, and the divergence with Kerr
case ($c=0$) increases as the value of $c$ increases, as we expected. Then, we
analytically solve motion equations using the Mino time and elliptic
integrals, following the strategy used in [12]. Our main formula, Eq. (26),
reproduce Kerr case when $c=0$ (regardless of the value of $w$). To estimate
the magnitude of the effect, we assumed a Keplerian orbit for the pulsar
orbiting a supermassive black hole in the equatorial plane ($q=0$), where
relativistic effect are strongest. Due to the big difference mass between
pulsar and BH, we considered the pulsar as a test particle. We did not deal
with individuals post-Keplerian effects, but just compared our model to the
Kerr case in order to isolate the effect of matter. We found that a deviation
is real and very pronounced in presence of dark energy ($w=-1/3$): while the
spin $a$ only induces a shift in the phase of the maximum delay without
increasing or decreasing the delay, the effect of matter surrounding the black
hole results in a noticeable alteration of the delay, which increases as the
strength parameter $c$ increases. For example, if the pulsar is at a distance
$r_{s}=100$ M, the additional time delay when the pulsar is behind the BH
(superior conjunction) would be more than 30 seconds. Interestingly, when the
pulsar is positioned between the observer and the black hole ($\phi\approx 0$,
i.e. inferior conjunction), we predict a slightly lower pulse delay (with an
advance of a few seconds) than in the no-matter case (see Fig. (3)).
Differences with Schwarzchild and Kerr cases are less evident (but still
appreciable) in presence of dust (see Fig. (4)) and practically imperceptible
in the case of pure radiation (see Fig. (5)).

Even if at the moment no pulsars close enough to SgrA* are known, their search
is attracting more and more efforts from the scientific community [29], due to
the enormous implications that such a discovery would have [30]. Since our
delay formulas do not depend on the energy of the emitted photons, they are
also valid for magnetars, which often show variability at higher energies than
radio waves. This allows us to estimate the delay of SGR J1745-2900, the
closest magnetar orbiting SgrA*. At a distance $r_{s}=0.1$ pc [31], i.e.
$r_{s}\simeq 5\times 10^{5}M$, and assuming an edge-on orbit 777More
observations are needed to determine the orientation of SGR J1745-2900’s
orbit., it turns out that the difference in delays between the no-matter case
(Kerr metric with $a=0.9$) and the dark-energy case (Kiselev metric with
$a=0.9$, $w=-1/3$ and $c=10^{-4}$) is approximately larger than $300$ seconds
at superior conjunction and shorter than $40$ seconds at inferior conjunction.
When $c=10^{-5}$ differences flatten significantly and are almost
imperceptible (see Fig. 6). For completeness, when $c=10^{-3}$ the delay in
presence of dark energy would be $\sim 4000$ seconds greater, a value so large
as to make us consider $c=10^{-3}$ unlikely at that distance from SgrA*.
Finally, we notice that, at the considered distance, the spinning effects are
almost completely absent almost completely absent and this makes the search
for matter effects even cleaner. By comparing the observed magnetar delays at
superior and inferior conjunctions with theoretical predictions, a signature
of the presence of matter (dark energy in particular) could appear. However,
at the moment the timing of SGR J1745-2900 has not yet reached sufficient
levels of precision for this purpose, mainly due to the intrinsic variability
of the source.

Figure 6: (Top) The exact time delay ($\Delta t_{\text{ex}}-\Delta
t_{\text{R}}$) for SGR J1745-2900, the closest magnetar orbiting SgrA*, for a
Schwarzschild black hole (black line), a Kerr black hole with $a=0.9$ (dashed
blue line), a non-rotating Kiselev black hole surrounded by dark energy with
$c=10^{-4}$ (dashed green line) and a rotating one with $a=0.9$ and
$c=10^{-4}$ (red line) for a circular edge-on orbit with a Schwarzschild
radius $r_{\text{S}}=10^{5}\,M$, where $M=4\times 10^{6}$ $M_{\odot}$.
(Bottom) As before, but with $c=10^{-5}$. By comparing the observed magnetar
delays at superior and inferior conjunctions with theoretical predictions, a
signature of the presence of dark energy could appear.

Beyond the information on any matter present and its effects on timing, the
advantage of this work, as compared to numrical ray-tracing methods, is the
possibility to isolate diverging terms providing at the same time non-
approximate results. Besides, our formulas could be integrated into a new
relativistic timing model that is suitable for extreme binary systems where
the presence of matter is not-negligible. However, this could be a not-easy
task, since it requires to know the environment around the black hole.
Conversely, one could exploit pulsar timing to constrain $w$ in the vicinity
of the black hole, if the pulsar is close enough. The difficulty in this case
would be isolating the matter effect from all other effects, which is
impossible in strong field situations, where the (non-linear) full
relativistic equations must be used.
Possible extensions of this work include the study of non-equatorial orbits,
the addition of a pulsar’s spin, the treatment of the pulsar as a timelike
geodesic. Parallel works, on the other hand, may concern the use of other
metrics, such as the Kerr-Newmann and Kerr-Sen ones. In principle, the
equations of motion for lighlike geodesics are solvable in terms of elliptic
or hyperelliptic integrals. Such investigations would probe the possibility to
test the no-hair theorem by predicting the time delays induced by the
additional parameters (axion, charge, ecc…). A degeneracy between different
parameters, however, could arise, invalidating an accurate measurement of
single parameters. In particular, in light of the results found here, the
presence of matter surrounding the black hole could affect accurate
measurements of the spin parameter $a$ by pulsar timing methods.
In conclusion, our results could serve, in principle, to better understand the
environment in which the binary system is located by exploiting the timing of
the pulsar, provided that this effect is separable from the others.

###### Acknowledgements.

A.C. would acknowledge the support by the Istituto Nazionale di Fisica
Nucleare (INFN) Iniziativa Specifica QGSKY. Furthermore, A.C. would express
his gratitude to Prof. Alberto Sesana for ospitality at University of Milano
’Bicocca’ and Prof. Andrea Possenti for hospitality at INAF-Cagliari, as well
as for their valuable advice.

##  Appendix A Multiple Roots

In this section we report quadruple, triple and double roots for the radial
potential $R(r)$ for a Kiselev metric, in presence of dust and radiation. We
will omit the computation of simple roots, since it is entirely equivalent to
Eq. (14). Furthermore, quadruple roots remain unchanged and, hence, equal to
the Kerr metric.
For radiation ($w=1/3$), triple roots are

$q=\dfrac{8r^{3}}{2M}-(a+\lambda)^{2},\;\;\;\;\;\lambda=a+\dfrac{r^{2}}{a}\left(3-\dfrac{2r}{M}\right)$

with the only non-zero real radial coordinate

$r=M-\left[M\left(M^{2}-a^{2}+c\right)\right]^{1/3}.$

Double solutions are obtained for

$q=\dfrac{r^{2}}{a^{2}(M-r)^{2}}\left[4a^{2}\left(Mr+c\right)-\left(3Mr-r^{2}+2c\right)^{2}\right],$
(36) $\lambda=a+\dfrac{r^{2}}{a}\left(1-\dfrac{2\Delta}{r(r-M)}\right).$

On the boundary of the allowable region (13), $q=0$ or $q=-(\lambda\pm
a)^{2}$. The first condition is realized for four real and distinct radial
values, two of which are always positive. The second condition
($q=-(\lambda+a)^{2}$) give

$r=\dfrac{3M^{2}+c}{3M}\cos\left[\dfrac{1}{3}\theta_{1}+\dfrac{2\pi
k}{3}\right]-\dfrac{A}{3}$

where $\theta_{1}$ and $A$ satisfy

$1+\cos\theta_{1}=\dfrac{54M^{4}(M^{2}-a^{2}+c)}{3M^{2}+c}^{3},\;\;\;\;\;A=\dfrac{c-3M^{2}}{2M}.$

As in the Kerr case, the last condition ($q=-(\lambda-a)^{2}$) implies that
$r=r_{\pm}$, where $r_{\pm}$ is the radial position of inner and outer
horizons.
In presence of dust ($w=0$), for triple roots we get

$q=\dfrac{r^{3}}{a^{2}}\left[\dfrac{8a^{2}(c+2M)-r(6M+3c-4r)^{2}}{(2M+c)^{2}}\right],$
$\lambda=a+\dfrac{r^{2}}{a}\left(\dfrac{6M+3c-4r}{2M+c}\right)$ (37)

as well as

$r=\dfrac{1}{2}\left[2M+c+n^{1/3}\right].$

with $n=\left(4a^{2}c-c^{3}+8a^{2}M-6c^{2}M-12cM^{2}-8M^{3}\right)$. Finally,
double solutions occur when

$q=\dfrac{a^{2}\left[16a^{2}\Delta-\left(4\Delta+r(2M-2r+c)\right)^{2}\right]}{a^{2}\left(2M-2r+c\right)^{2}},$
$\lambda=a+\dfrac{r^{2}}{a}-\dfrac{4r\Delta}{a(2r-2M-c)}$

which leads to two cubic equations when $q=0$ and $q=-(\lambda+a)^{2}$, as
well as $r=r_{\pm}$ when $q=-(\lambda-a)^{2}$. Their solutions are,
respectively

$r=2\sqrt{M(c+M)}\cos\left(\dfrac{\theta_{2}}{3}+\dfrac{2\pi
k}{3}\right)+2M+c,$
$r=\dfrac{c+2M}{2}\cos\left(\dfrac{\theta_{3}}{3}+\dfrac{2\pi
k}{3}\right)+\dfrac{2M+c}{4}$ (38)

where

$\displaystyle\cos\theta_{2}$
$\displaystyle=\dfrac{(c+2M)(2a^{2}+2c^{2}-cM-M^{2})}{\left(M^{2}+cM+c^{2}\right)^{3/2}}\;,$
$\displaystyle\cos\theta_{3}$ $\displaystyle=1-\dfrac{8a^{2}}{(c+2M)^{2}}.$

##  Appendix B Time delay integrals

In the following, we summarize the computation of the time delay integral (26)
in presence of a dark energy component. We also report the final formulas for
dust and radiation, the steps being completely analogous. For the angular
integral (25) we refer to [12].

###  B.1 Dark Energy

To solve analytically the integral (22), we start with the cahnge of variable:

$x^{2}=\dfrac{(r-r_{4})(r_{3}-r_{1})}{(r-r_{3})(r_{4}-r_{1})}$ (39)

where $r_{i}$ are the roots of the radial potential $R(r)$. This leads the
integral into the form

$\begin{array}[]{l}T(r,\lambda,q;c)=\displaystyle\int_{r_{4}}^{r}\dfrac{G(r)}{\Delta\sqrt{R(r)}}dr=\\\
\delta\cdot\displaystyle\int_{0}^{x(r)}\dfrac{G(x)}{\Delta(x)\sqrt{(1-x^{2})(1-k^{2}x^{2})}}dx\end{array}$
(40)

where $G(r)$ is defined in Eq. (21) and

$k^{2}=\dfrac{(r_{3}-r_{2})(r_{4}-r_{1})}{(r_{3}-r_{1})(r_{4}-r_{2})}.$

We then expand the quantity $\dfrac{G(x)}{\Delta}$ in partial fractions:

$\dfrac{G(x)}{\Delta}=\dfrac{1}{\gamma^{3}}\Big{[}4M^{2}+(2Mr-a^{2}c^{2})\gamma+(r^{2}-ac\lambda)\gamma^{2}\Big{]}+A(r)+B(r)$

where we have defined

$A(r)=\frac{8M^{3}r_{+}-4a^{2}M\left(M+cr_{+}\right)+a^{3}c\lambda-2\tilde{a}r_{+}M\lambda}{2\sqrt{M^{2}-a^{2}\gamma}\left(r-r_{+}\right)\gamma^{3}},$

$B(r)=\frac{-8M^{3}r_{-}+4a^{2}M\left(M+cr_{-}\right)-a^{3}c\lambda+2\tilde{a}M\lambda
r_{-}}{2\sqrt{M^{2}-a^{2}\gamma}\left(r-r_{-}\right)\gamma^{3}}.$

where $\tilde{a}=a(1-2c)$. We then use the following well-known closed
integrals

$F(x,k)=\displaystyle\int_{0}^{x}\dfrac{dx}{\sqrt{(1-x^{2})(1-k^{2}x^{2})}},$
$E(x,k)=\displaystyle\int_{0}^{x}\dfrac{\sqrt{1-k^{2}x^{2}}}{\sqrt{1-x^{2}}}dx,$
$\Pi(x,l,k)=\displaystyle\int_{0}^{x}\dfrac{dx}{(1-lx^{2})\sqrt{(1-x^{2})(1-k^{2}x^{2})}},$

having defined the constant $l=\dfrac{r_{1}-r_{4}}{r_{1}-r_{3}}$. Noting that

$\displaystyle\int_{0}^{x}\dfrac{-x^{2}}{(1-lx^{2})\sqrt{(1-x^{2})(1-k^{2}x^{2})}}=\dfrac{1}{l}\Big{[}F(x,k)-\Pi(x,l,k)\Big{]},$

and since [32]

$\displaystyle\dfrac{1}{(1-lx^{2})^{2}y(x)}=$ $\displaystyle
C_{1}\dfrac{d}{dx}\Big{[}\dfrac{xy(x)}{1-lx^{2}}\Big{]}+\dfrac{C_{2}}{y(x)}+\dfrac{C_{3}(1-k^{2}x^{2})}{y(x)}+\dfrac{C_{4}}{y(x)(1-lx^{2})}$

(to be integrated), one obtains Eq. (26) (for $w=-1/3$), where we also called

$B_{+}=\frac{8M^{3}r_{+}-4a^{2}M\left(M+cr_{+}\right)+a^{3}c\lambda-2\tilde{a}r_{+}M\lambda}{2\sqrt{M^{2}-a^{2}\gamma}\left(r_{4}-r_{+}\right)\gamma^{3}},$
(41)
$B_{-}=\frac{-8M^{3}r_{-}+4a^{2}M\left(M+cr_{-}\right)-a^{3}c\lambda+2\tilde{a}M\lambda
r_{-}}{2\sqrt{M^{2}-a^{2}\gamma}\left(r_{4}-r_{-}\right)\gamma^{3}}$ (42)

where

$r_{\pm}=\dfrac{M\pm\sqrt{M^{2}-a^{2}\gamma}}{\gamma}$

are the positions of the horizons in the case $w=-1/3$. The presence of
diverging parts for $r\rightarrow\infty$ in $\Pi(x,l,k)$ and
$\dfrac{\sqrt{R(r)}}{r-r_{3}}$ can be isolated by exploiting the relation

$\Pi(x,l,k)=F(x,k)-\Pi\Big{(}x,\frac{k^{2}}{l},k\Big{)}+\dfrac{\ln(Z)}{2P}$
(43)

with

$\displaystyle Z$
$\displaystyle=\frac{\sqrt{\left(1-x^{2}\right)\left(1-k^{2}x^{2}\right)}+Px}{\sqrt{\left(1-x^{2}\right)\left(1-k^{2}x^{2}\right)}-Px},$
(44) $\displaystyle P^{2}$
$\displaystyle=\frac{(l-1)(l-k^{2})}{l}=\frac{\left(r_{3}-r_{4}\right)^{2}}{\left(r_{4}-r_{2}\right)\left(r_{3}-r_{1}\right)}.$

Therefore, Eq. (26) becomes

$\begin{array}[]{l}T_{r}\left(r,\lambda,q;c\right)=\\\
\delta\cdot\Big{[}F(x,k)\cdot\Big{(}4M^{2}\gamma^{-3}-a^{2}c^{2}\gamma^{-2}-ac\lambda\gamma^{-1}+2Mr_{3}\gamma^{-2}\\\
+\frac{1}{2}\gamma^{-1}\big{[}r_{1}\big{(}r_{3}-r_{4}\big{)}+r_{3}\big{(}r_{3}+r_{4}\big{)}\big{]}+\frac{B_{+}l}{l_{+}}+\frac{B_{-}l}{l_{-}}+2M(r_{4}-r_{3})\gamma^{-2}\Big{)}\\\
+E\big{(}x,k\big{)}\Big{(}-\frac{1}{2}\gamma^{-1}\big{(}r_{4}-r_{2}\big{)}\big{(}r_{3}-r_{1}\big{)}\Big{)}-\\\
-\Pi\big{(}x,\frac{k^{2}}{l},k\big{)}\Big{(}2Mr_{4}\gamma^{-2}-2Mr_{3}\gamma^{-2}\Big{)}+\\\
+\Pi\big{(}x,l_{+},k\big{)}\Big{(}B_{+}-\frac{lB_{+}}{l_{+}}\Big{)}+\Pi\big{(}x,l_{-},k\big{)}\Big{(}B_{-}-\frac{lB_{-}}{l_{-}}\Big{)}\Big{]}+\\\
+2M\delta\gamma^{-2}\big{(}r_{4}-r_{3}\big{)}\dfrac{\ln(Z)}{2P}+\dfrac{\sqrt{R}}{r-r_{3}}\gamma^{-1}.\\\
\end{array}$ (45)

where the last two terms are divergent when $r\rightarrow\infty$:

$\dfrac{\sqrt{R(r)}}{r-r_{3}}=r+r_{3}+\mathcal{O}\left(\dfrac{1}{r}\right)$
(46)
$\ln(Z)=\ln\left(\dfrac{2}{r_{3}+r_{4}}\right)+\ln(r)+\mathcal{O}\left(\dfrac{1}{r}\right)$
(47)

For the angular integral (25) as well as for the constants $v$ and $w$ we
directly refer to Eqs. (A17-A19) in [12].

###  B.2 Radiation

In a very similar way we obtain for the radiation case ($w=1/3$)

$\begin{array}[]{l}T_{r}\left(r,\lambda,q;c\right)=\delta\cdot\Big{[}F(x,k)\Big{(}4M^{2}+c+2Mr_{3}\\\
+\dfrac{1}{2}\big{[}r_{1}(r_{3}-r_{4})+r_{3}(r_{3}+r_{4})\big{]}+2M(r_{4}-r_{3})\\\
+\dfrac{B_{+}l}{l_{+}}+\dfrac{B_{-}l}{l_{-}}\Big{)}+E(x,k)\Big{(}-\dfrac{1}{2}(r_{4}-r_{2})(r_{3}-r_{1})\Big{)}\\\
-\Pi\Big{(}x,\dfrac{k^{2}}{l},k\Big{)}\Big{(}2M(r_{4}-r_{3})\Big{)}\\\
+\Pi(x,l_{+},k)\Big{(}B_{+}-\dfrac{B_{+}l}{l_{+}}\Big{)}+\Pi(x,l_{-},k)\Big{(}B_{-}-\dfrac{B_{-}l}{l_{-}}\Big{)}\Big{]}\\\
+2M\delta(r_{4}-r_{3})\dfrac{\ln{Z}}{2P}+\dfrac{\sqrt{R(r)}}{r-r_{3}}\end{array}$
(48)

where

$\displaystyle R(r)=\Big{(}r^{2}-2M$ $\displaystyle
r+a^{2}-c\Big{)}\Big{(}-q-\lambda^{2}-a^{2}+2a\lambda\Big{)}$ (49)
$\displaystyle+\Big{(}r^{2}+a^{2}-\lambda a\Big{)}^{2}\;,$
$r_{\pm}=M\pm\sqrt{M^{2}-a^{2}+c}\;,$
$B(r)=\dfrac{c^{2}+4M^{2}(c-a^{2})+4cMr+8M^{3}r-a\lambda(c+2Mr)}{2(r_{4}-r)\sqrt{M^{2}-a^{2}+c}},$
$B_{+}=B(r_{+}),\;\;\;\;\;B_{-}=-B(r_{-}).$ (50)

Finally, the $\varphi$ coordinate is given in the integral form

$\varphi_{a}-\varphi_{e}=\int_{\gamma_{r}}\dfrac{2Mra-a^{2}\lambda+ac}{\Delta\sqrt{R(r)}}dr+\dfrac{2\lambda
F(x,k)}{\sqrt{(r_{4}-r_{2})(r_{3}-r_{1})}}.$ (51)

###  B.3 Dust

As done in the previous section, the time delay in presence of dust ($w=0$)
reads as

$\begin{array}[]{l}T_{r}\left(r,\lambda,q;c\right)=\delta\cdot\Big{[}F(x,k)\Big{(}4M^{2}+4cM+c^{2}+r_{3}(2M+c)\\\
+(2M+c)(r_{4}-r_{3})+\dfrac{1}{2}\big{[}r_{1}(r_{3}-r_{4})+r_{3}(r_{3}+r_{4})\big{]}+\dfrac{B_{+}l}{l_{+}}+\dfrac{B_{-}l}{l_{-}}\Big{)}\\\
+E(x,k)\Big{(}-\dfrac{1}{2}(r_{4}-r_{2})(r_{3}-r_{1})\Big{)}\\\
-\Pi\big{(}x,\dfrac{k^{2}}{l},k\big{)}(2M+c)(r_{4}-r_{3})\\\
+\Pi(x,l_{+},k)\Big{(}B_{+}-\dfrac{lB_{+}}{l_{+}}\Big{)}+\Pi(x,l_{-},k)\Big{(}B_{-}-\dfrac{lB_{-}}{l_{-}}\Big{)}\Big{]}\\\
+\delta(2M+c)(r_{4}-r_{3})\dfrac{\ln{Z}}{2P}+\dfrac{\sqrt{R(r)}}{r-r_{3}}\end{array}$
(52)

where now

$R(r)=-\Big{(}r^{2}-2Mr+a^{2}-cr\Big{)}\Big{(}q+\lambda^{2}+a^{2}-2a\lambda\Big{)}+\Big{(}r^{2}+a^{2}-\lambda
a\Big{)}^{2},$ (53)
$r_{\pm}=\dfrac{1}{2}\Big{[}c+2M\pm\sqrt{4M^{2}-4a^{2}+4cM+c^{2}}\Big{]},$
$B(r)=\dfrac{(c+2M)\Big{[}-a^{2}(c+2M)+(c+2M)^{2}r-a\lambda
r\Big{]}}{(r_{4}-r)\sqrt{(2M+c)^{2}-4a^{2}}},$
$B_{+}=B(r_{+}),\;\;\;\;\;B_{-}=-B(r_{-}).$ (54)

Finally, the $\varphi$ coordinate is given in the integral form

$\varphi_{a}-\varphi_{e}=\int_{\gamma_{r}}\dfrac{(2Ma+ac)r-a^{2}\lambda}{\Delta\sqrt{R(r)}}dr+\dfrac{2\lambda
F(x,k)}{\sqrt{(r_{4}-r_{2})(r_{3}-r_{1})}}.$ (55)

## References

  * Taylor and Weisberg [1989] J. H. Taylor and J. M. Weisberg, Further experimental tests of relativistic gravity using the binary pulsar psr 1913+ 16, The Astrophysical Journal 345, 434 (1989).
  * Stairs [2003] I. H. Stairs, Testing general relativity with pulsar timing, Living Reviews in Relativity 6, 10.12942/lrr-2003-5 (2003).
  * Freire [2022] P. C. C. Freire, Tests of gravity theories with pulsar timing (2022).
  * Liu  _et al._ [2012] K. Liu, N. Wex, M. Kramer, J. M. Cordes, and T. J. W. Lazio, PROSPECTS FOR PROBING THE SPACETIME OF WITH PULSARS, The Astrophysical Journal 747, 1 (2012).
  * Izmailov  _et al._ [2019] R. N. Izmailov, E. R. Zhdanov, A. Bhadra, and K. K. Nandi, Relative time delay in a spinning black hole as a diagnostic for no-hair theorem, Eur. Phys. J. C 79, 105 (2019).
  * Detweiler [1979] S. Detweiler, Pulsar timing measurements and the search for gravitational waves, Astrophys. J. 234, 1100 (1979).
  * Goncharov  _et al._ [2021] B. Goncharov, R. M. Shannon, D. J. Reardon, G. Hobbs, A. Zic, M. Bailes, M. Curyło, S. Dai, M. Kerr, M. E. Lower, R. N. Manchester, R. Mandow, H. Middleton, M. T. Miles, A. Parthasarathy, E. Thrane, N. Thyagarajan, X. Xue, X.-J. Zhu, A. D. Cameron, Y. Feng, R. Luo, C. J. Russell, J. Sarkissian, R. Spiewak, S. Wang, J. B. Wang, L. Zhang, and S. Zhang, On the evidence for a common-spectrum process in the search for the nanohertz gravitational-wave background with the parkes pulsar timing array, The Astrophysical Journal Letters 917, L19 (2021).
  * Chen  _et al._ [2021] S. Chen, R. N. Caballero, Y. J. Guo, A. Chalumeau, K. Liu, G. Shaifullah, K. J. Lee, S. Babak, G. Desvignes, A. Parthasarathy, H. Hu, E. van der Wateren, J. Antoniadis, A.-S. B. Nielsen, C. G. Bassa, A. Berthereau, M. Burgay, D. J. Champion, I. Cognard, M. Falxa, R. D. Ferdman, P. C. C. Freire, J. R. Gair, E. Graikou, L. Guillemot, J. Jang, G. H. Janssen, R. Karuppusamy, M. J. Keith, M. Kramer, X. J. Liu, A. G. Lyne, R. A. Main, J. W. McKee, M. B. Mickaliger, B. B. P. Perera, D. Perrodin, A. Petiteau, N. K. Porayko, A. Possenti, A. Samajdar, S. A. Sanidas, A. Sesana, L. Speri, B. W. Stappers, G. Theureau, C. Tiburzi, A. Vecchio, J. P. W. Verbiest, J. Wang, L. Wang, and H. Xu, Common-red-signal analysis with 24-yr high-precision timing of the european pulsar timing array: inferences in the stochastic gravitational-wave background search, Monthly Notices of the Royal Astronomical Society 508, 4970 (2021).
  * Damour and Deruelle [1986] T. Damour and N. Deruelle, General relativistic celestial mechanics of binary systems. II. The post-Newtonian timing formula., Ann. Inst. Henri Poincaré Phys. Théor 44, 263 (1986).
  * Edwards  _et al._ [2006] R. T. Edwards, G. B. Hobbs, and R. N. Manchester, TEMPO2, a new pulsar timing package - II. The timing model and precision estimates, mnras 372, 1549 (2006), arXiv:astro-ph/0607664 [astro-ph] .
  * Manchester [2015] R. N. Manchester, Pulsars and gravity, International Journal of Modern Physics D 24, 1530018 (2015).
  * Ben-Salem and Hackmann [2022] B. Ben-Salem and E. Hackmann, Propagation time delay and frame dragging effects of lightlike geodesics in the timing of a pulsar orbiting SgrA*, Monthly Notices of the Royal Astronomical Society 516, 1768 (2022).
  * Hackmann and Dhani [2019] E. Hackmann and A. Dhani, The propagation delay in the timing of a pulsar orbiting a supermassive black hole, General Relativity and Gravitation 51, 10.1007/s10714-019-2517-2 (2019).
  * Zhang and Saha [2017] F. Zhang and P. Saha, Probing the spinning of the massive black hole in the galactic center via pulsar timing: A full relativistic treatment, The Astrophysical Journal 849, 33 (2017).
  * Kimpson  _et al._ [2019] T. Kimpson, K. Wu, and S. Zane, Pulsar timing in extreme mass ratio binaries: a general relativistic approach, Monthly Notices of the Royal Astronomical Society 486, 360 (2019).
  * Farrah  _et al._ [2023a] D. Farrah, S. Petty, K. S. Croker, G. Tarlé , M. Zevin, E. Hatziminaoglou, F. Shankar, L. Wang, D. L. Clements, A. Efstathiou, M. Lacy, K. A. Nishimura, J. Afonso, C. Pearson, and L. K. Pitchford, A preferential growth channel for supermassive black holes in elliptical galaxies at z $\lesssim$ 2, The Astrophysical Journal 943, 133 (2023a).
  * Farrah  _et al._ [2023b] D. Farrah, K. S. Croker, M. Zevin, G. Tarlé , V. Faraoni, S. Petty, J. Afonso, N. Fernandez, K. A. Nishimura, C. Pearson, L. Wang, D. L. Clements, A. Efstathiou, E. Hatziminaoglou, M. Lacy, C. McPartland, L. K. Pitchford, N. Sakai, and J. Weiner, Observational evidence for cosmological coupling of black holes and its implications for an astrophysical source of dark energy, The Astrophysical Journal Letters 944, L31 (2023b).
  * Kiselev [2003] V. V. Kiselev, Quintessential solution of dark matter rotation curves and its simulation by extra dimensions, (2003), arXiv:gr-qc/0303031 .
  * Ghosh [2016] S. G. Ghosh, Rotating black hole and quintessence, Eur. Phys. J. C 76, 222 (2016), arXiv:1512.05476 [gr-qc] .
  * Toshmatov  _et al._ [2017] B. Toshmatov, Z. Stuchlík, and B. Ahmedov, Rotating black hole solutions with quintessential energy, Eur. Phys. J. Plus 132, 98 (2017), arXiv:1512.01498 [gr-qc] .
  * Melia [2015] F. Melia, The Cosmic Equation of State, Astrophys. Space Sci. 356, 393 (2015), arXiv:1411.5771 [astro-ph.CO] .
  * Carleo  _et al._ [2022] A. Carleo, G. Lambiase, and L. Mastrototaro, Energy extraction via magnetic reconnection in Lorentz breaking Kerr–Sen and Kiselev black holes, Eur. Phys. J. C 82, 776 (2022), arXiv:2206.12988 [gr-qc] .
  * Akiyama  _et al._ [2019] K. Akiyama _et al._ (Event Horizon Telescope), First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole, Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Pratap Singh [2017] B. Pratap Singh, Rotating charged black holes shadow in quintessence, arXiv e-prints , arXiv:1711.02898 (2017), arXiv:1711.02898 [gr-qc] .
  * Carter [1968] B. Carter, Global structure of the kerr family of gravitational fields, Phys. Rev. 174, 1559 (1968).
  * Gralla and Lupsasca [2020] S. E. Gralla and A. Lupsasca, Null geodesics of the kerr exterior, Physical Review D 101, 10.1103/physrevd.101.044032 (2020).
  * Dexter and Agol [2009] J. Dexter and E. Agol, A FAST NEW PUBLIC CODE FOR COMPUTING PHOTON ORBITS IN a KERR SPACETIME, The Astrophysical Journal 696, 1616 (2009).
  * Blandford and Teukolsky [1976] R. Blandford and S. A. Teukolsky, Arrival-time analysis for a pulsar in a binary system., Astrophys. J. 205, 580 (1976).
  * Torne  _et al._ [2021] P. Torne, G. Desvignes, R. P. Eatough, M. Kramer, R. Karuppusamy, K. Liu, A. Noutsos, R. Wharton, C. Kramer, S. Navarro, G. Paubert, S. Sanchez, M. Sanchez-Portal, K. F. Schuster, H. Falcke, and L. Rezzolla, Searching for pulsars in the galactic centre at 3 and 2 mm, Astronomy & Astrophysics 650, A95 (2021).
  * Pfahl and Loeb [2004] E. Pfahl and A. Loeb, Probing the Spacetime around Sagittarius A* with Radio Pulsars, Astrophys. J. 615, 253 (2004), arXiv:astro-ph/0309744 [astro-ph] .
  * Yusef-Zadeh  _et al._ [2015] F. Yusef-Zadeh, R. Diesing, M. Wardle, L. O. Sjouwerman, M. Royster, W. D. Cotton, D. Roberts, and C. Heinke, RADIO CONTINUUM EMISSION FROM THE MAGNETAR SGR j1745-2900: INTERACTION WITH GAS ORBITING sgr a*, The Astrophysical Journal 811, L35 (2015).
  * Boas Jr [1954] R. Boas Jr, Higher transcendental functions, vols. i and ii., Science 120, 302 (1954).