© ESO 2020

1. Introduction

The gyroscopic nature of radio pulsars (PSRs) allows them to be used as exceptionally stable and regular clocks. By measuring the times of arrival (ToAs) of the pulsar pulses, an observer can determine information on the rotational and astrometric properties of the PSR, in addition to information on any binary companion (e.g. Liu et al. 2011, 2018; Desvignes et al. 2016; Lazarus et al. 2016; Parthasarathy et al. 2019). In the Galactic centre, there is expected to be a large population (up to ∼10⁴ within the central parsec) of millisecond pulsars (MSPs; Wharton et al. 2012; Rajwade et al. 2017), whilst the number of pulsars in globular clusters is observed to be greater by a factor of 10³ per unit mass compared to the Galactic disk (Freire et al. 2013). In addition to MSPs, the Galactic centre is known to host a supermassive black hole (BH) with a mass of M ∼ 4 × 10⁶ M_⊙ (Gillessen et al. 2009; Boehle et al. 2016), whilst consilient strands of evidence from stellar kinematics (e.g. Gebhardt et al. 2002; van der Marel & Anderson 2010; Feldmeier et al. 2013), radiative accretion signatures (e.g. Ulvestad et al. 2007) and pulsar dynamics (Kızıltan et al. 2017) all suggest that globular clusters could host intermediate mass black holes (IMBHs), although definitive “smoking gun” evidence is still lacking, and some evidence is disputed (Miller-Jones et al. 2012; Baumgardt 2017). Although no MSPs have currently been detected in the Galactic centre or at the centre of globular clusters, the increased sensitivity of next generation of radio facilities – such as the Square Kilometer Array (SKA, Keane et al. 2015), the Five-hundred-meter Aperture Spherical Telescope (FAST, Nan et al. 2011), or the NASA Deep Space Network (DSN, Majid et al. 2019) – will make this population accessible.

Radio timing the orbit of an MSP around a much more massive BH has been identified as an excellent apparatus with which to probe general relativity (GR) and astrophysics in a highly strong-field, non-linear regime (Kramer et al. 2004; Liu et al. 2012; Psaltis et al. 2016). These systems are known as extreme (mass ratio q ∼ 10⁶) or intermediate (q ∼ 10³ − 10⁴) mass ratio binaries (E/IMRBs). Via radio timing of an MSP-E/IMRB, the essential parameters of the central BH – the mass, spin, and quadrupole moment (M, S, Q respectively) – can be determined to remarkable precision (e.g. Liu et al. 2012; Psaltis et al. 2016). Once the BH parameters have been established, it is then possible to explore fundamental questions of GR such as the cosmic censorship conjecture (Penrose 1969) and no hair theorem (Israel 1967). From the perspective of astrophysics, one could then establish whether the central dark nuclei of star clusters are indeed astrophysical BHs as described by the Kerr solution, or else something more exotic such as a boson star (Kleihaus et al. 2012) or a BH described by some deviation from the Kerr solution (“bumpy” black holes; Yagi & Stein 2016). Moreover, if an MSP could be detected orbiting an IMBH, this would also naturally settle the debate on the existence of intermediate mass, astrophysical BHs (Singh et al. 2014).

In order to determine the system parameters via radio timing, it is necessary to have a timing model which maps the time at which the pulse was received in the observers frame t to the time at which the pulse was emitted in the “clock frame” τ (i.e. the frame comoving with the pulsar). For pulsars in binary systems, there are four key delays which must be accounted for in the timing model:

$$\begin{array}{c}\hfill t=\tau +{\mathrm{\Delta }}_{\mathrm{E}}+{\mathrm{\Delta }}_{\mathrm{R}}+{\mathrm{\Delta }}_{\mathrm{S}}\end{array}$$(1)

where Δ_R is the Romer delay due to orbital motion of the pulsar, Δ_S is the Shaprio delay due to gravitational time dilation, and Δ_E is the Einstein delay, the relativistic time dilation (see Lorimer & Kramer 2004, for further details). The relationship between the timing delays (Δ_R, Δ_S, Δ_E) and the system parameters is based on the “Post-Keplerian” (PK) approximation (Blandford & Teukolsky 1976; Damour & Deruelle 1986) and is derived under the assumption that the gravitational field is weak. GR is subsequently treated perturbatively as an expansion about Newtonian gravity. Whilst this considerably simplifies the previously non-linear equations, such an approximation may not be appropriate for relativistic orbits close to massive BHs. The strength of the gravitational potential can be quantified via the parameter

$$\begin{array}{c}\hfill ϵ=\frac{\mathit{GM}}{r{c}^2}.\end{array}$$(2)

Tests of GR in the solar system which employ the PK framework are undertaken in an environment where typically ϵ ∼ 10⁻⁸. In binary pulsar systems ϵ ≲ 10⁻⁶. Conversely, for a pulsar at periapsis on a P = 0.1 year eccentric (e = 0.9) orbit around Sgr A^*, ϵ = 0.01. Due to external perturbations on the PSR orbit (Merritt et al. 2010), the cleanest information on the BH parameters will be obtained from the PSR at periapsis, in precisely the region of strongest gravity. Photons propagating from the pulsar to an observer may pass sufficiently close to the BH so as to traverse a spacetime where ϵ approaches unity (see e.g. discussion on the use of pulsars to probe quantum effects and the BH event horizon Estes et al. 2017). It therefore seems reasonable to question the validity of the PK weak-field approximation in a region where the gravitational potential is ≳10⁴ times stronger than the environments in which it is typically applied.

Pulsar timing delays given by the PK approximation in strong-field regimes was partly addressed in Hackmann & Dhani (2019) who investigate the propagation delays (Roemer and Shapiro) from circular orbits in a Schwarzschild spacetime. The influence of the central black hole spin on the measurement of pulsar orbital parameters was investigated in Zhang & Saha (2017), whilst Li et al. (2019) explored the impact of relativistic spin dynamics on the PSR emission. In this work, we further develop these previous studies by deploying a consistent, relativistic framework (Kimpson et al. 2019a,b) to accurately calculate the pulse ToA from a spinning pulsar orbiting a Kerr BH. We can then compare this “true” GR solution with the “approximate” PK solution and so determine the validity of the PK approach for BH parameter estimation at the Galactic centre and the centre of globular clusters.

This paper is organised as follows. In Sect. 2 we review the construction of a relativistic pulsar timing model via numerical methods. In Sect. 3 we overview the PK parameters that are used for tests of GR and compare the PK solution with the GR solution. Discussion and concluding remarks are made in Section 4. We adopt the natural units, with c = G = ℏ = 1, and a ( − , + , + , + ) metric signature. Unless otherwise stated, a c.g.s. Gaussian unit system is used in the expressions for the electromagnetic properties of matter. The gravitational radius of the black hole is r_g = M and the corresponding Schwarzschild radius is r_s = 2M, where M is the black-hole mass. A comma denotes a partial derivative (e.g. ​x_,r), and a semicolon denotes the covariant derivative (e.g. ​x_;r).

2. Relativistic MSP timing model

We now briefly review the methods of Kimpson et al. (2019b) for simulating the MSP timing signal. This has two primary components: the spin-orbital dynamics of a spinning MSP around a spinning BH, and the propagation of the MSP radio beam to the observer. We refer the reader to Kimpson et al. (2019b) for details regarding the algorithmic combination of these two effects to generate a consistent ToAs. For this work we ignore chromatic plasma effects and work solely in vacuum.

2.1. Kerr spacetime

The spacetime around a rotating BH is described by the Kerr solution. In Boyer-Lindquist coordinates, the metric interval is given by,

$$\begin{array}{cc}\hfill \mathrm{d}{s}^2& =-(1-\frac{2Mr}{\mathrm{\Sigma }})\mathrm{d}{t}^2-\frac{4aMr{sin}^2\theta }{\mathrm{\Sigma }}\mathrm{d}t\mathrm{d}\varphi +\frac{\mathrm{\Sigma }}{\mathrm{\Delta }}\mathrm{d}{r}^2+\mathrm{\Sigma }\mathrm{d}{\theta }^2\hfill \\ \hfill & +\frac{{sin}^2\theta }{\mathrm{\Sigma }}[{(r^2+a^2)}^2-\mathrm{\Delta }{a}^2{sin}^2\theta ]\mathrm{d}{\varphi }^2,\hfill \end{array}$$(3)

where Σ = r² + a²cos²θ, Δ = r² − 2Mr + a² and a is the black-hole spin parameter. Generally we set the BH mass M = 1 such that the lengthscale is normalised to the gravitational lengthscale r_g. The Kerr spacetime is both stationary -which leads to the conservation of energy E -and axisymmetric – associated with the conservation of the projection of the angular momentum, L_z. In additional to these two conserved quantities which are associated with the Killing vectors of the spacetime, the Carter constant Q (Carter 1968) is also conserved, which is associated with a Killing tensor. The Carter constant does not have an exact physical interpretation, but is related to the inclination of an orbiting body, and this is further discussed in De Felice & Preti (1999), Rosquist et al. (2009).

2.2. Spin-orbital dynamics

Astrophysical pulsars are not in reality point particles but extended spinning objects. This spin couples to the background spacetime and as a consequence a pulsar will not follow a geodesic of the Kerr metric. Instead, in limit of extreme mass ratio the MSP spin orbital dynamics can be described by the background gravitational field and the dynamical spin interaction with this field. The equations of motion of the pulsar are then determined via the MPD formulation (Mathisson 1937; Papapetrou 1951; Dixon 1974). Via this approach one constructs the “gravitational skeleton” of the body by performing a multipole expansion of the stress energy tensor, T^μν. Since the general equation of motion is defined via,

$$\begin{array}{c}\hfill {T^{\mu \nu }}_{;\nu }=0,\end{array}$$(4)

then to dipole order the equations of motion are,

$$\begin{array}{c}\hfill \frac{D{p}^{\mu }}{\mathrm{d}\tau }=-\frac{1}{2}{R^{\mu }}_{\nu \alpha \beta }{u}^{\nu }{s}^{\alpha \beta },\end{array}$$(5)

$$\begin{array}{c}\hfill \frac{D{s}^{\mu \nu }}{\mathrm{d}\tau }=p^{\mu }{u}^{\nu }-p^{\nu }{u}^{\mu }\end{array}$$(6)

(Mathisson 1937; Papapetrou 1951; Dixon 1974), where p^μ is the 4-momentum (mass monopole), s^μν the spin dipole, u^ν is the MSP 4-velocity and ${R^{\mu }}_{\nu \alpha \beta }$ the Riemann curvature tensor. We use the proper time τ to parameterise the MSP worldline and D/dτ denotes a covariant derivative. To make the system of equations determinate one must specify the reference worldline with which to define the multipole expansion, equivalent to choosing an observer with which with respect to which the centre of mass is defined, since in GR the centroid is observer-dependent. We take the centroid to be that measured in zero 3-momentum frame:

$$\begin{array}{c}\hfill s^{\mu \nu }{p}_{\nu }=0\end{array}$$(7)

(Tulczyjew 1959; Dixon 1964), the so-called TD condition. For a discussion on alternative centroid choices see Filipe Costa & Natário (2015), Babak et al. (2014).

Given that the Moller radius of the MSP if much less than the physical radius, and the MSP mass much less than the BH mass, the pole-dipole interaction is much stronger than the dipole-dipole interaction and the equations of motion reduce to

$$\begin{array}{c}\hfill \frac{D{u}^{\mu }}{\mathrm{d}\tau }=-\frac{1}{2m}{R^{\mu }}_{\nu \alpha \beta }{u}^{\nu }{s}^{\alpha \beta },\end{array}$$(8)

$$\begin{array}{c}\hfill \frac{D{s}^{\mu \nu }}{\mathrm{d}\tau }\approx 0\end{array}$$(9)

(Chicone et al. 2005; Mashhoon & Singh 2006). The complete set of ODEs is then

$$\begin{array}{c}\hfill \frac{\mathrm{d}{p}^{\alpha }}{\mathrm{d}\tau }=-{{\mathrm{\Gamma }}^{\alpha }}_{\mu \nu }{p}^{\mu }{u}^{\nu }+\left(\frac{1}{2m}{R^{\alpha }}_{\beta \rho \sigma }{ϵ}_{\mu \nu }^{\rho \sigma }{s}^{\mu }{p}^{\nu }{u}^{\beta }\right),\end{array}$$(10)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{s}^{\alpha }}{\mathrm{d}\tau }=-{{\mathrm{\Gamma }}^{\alpha }}_{\mu \nu }{s}^{\mu }{u}^{\nu }+\left(\frac{1}{2m^3}{R}_{\gamma \beta \rho \sigma }{ϵ}_{\mu \nu }^{\rho \sigma }{s}^{\mu }{p}^{\nu }{s}^{\gamma }{u}^{\beta }\right)p^{\alpha },\end{array}$$(11)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{x}^{\alpha }}{\mathrm{d}\tau }=-\frac{p^{\delta }{u}_{\delta }}{m^2}[p^{\alpha }+\frac{1}{2}\frac{(s^{\alpha \beta }{R}_{\beta \gamma \mu \nu }{p}^{\gamma }{s}^{\mu \nu })}{m^2+(R_{\mu \nu \rho \sigma }{s}^{\mu \nu }{s}^{\beta \sigma }/4)}]\end{array}$$(12)

(Singh 2005; Mashhoon & Singh 2006), where s^μ is the spin 4-vector contraction of the spin tensor,

$$\begin{array}{c}\hfill s_{\mu }=-\frac{1}{2m}{ϵ}_{\alpha \beta \mu \nu }{p}^{\nu }{s}^{\alpha \beta }.\end{array}$$(13)

Solving Eqs. (10)–(12) numerically via a 5th order Runge-Kutta-Fehlberg algorithm (Press 1977) then completely specifies the spin-orbital dynamics of the pulsar, accounting for the background spacetime and the dynamical interaction of the MSP spin.

As discussed, in Kerr geometry E, L_z and Q of an MSP are conserved along its orbit. For astronomical purposes it is useful to map these constants to orbital parameters. In this way we can specify the sort of Keplerian orbit that we want to describe. In particular, the E, L_z, Q can be mapped to the orbital elements 𝒫, e, ι, the semi-latus rectum, eccentricity and inclination respectively. In the gravitational strong field the semi-latus rectum and the eccentricity are effectively a re-parameterisation in turn of the apsidal approaches (periapsis distance r_p and apoapsis r_a),

$$\begin{array}{c}\hfill r_{\mathrm{p}}=\frac{\mathcal{P}}{1+e};r_{\mathrm{a}}=\frac{\mathcal{P}}{1-e},\end{array}$$(14)

whilst in the weak field limit 𝒫, e, correspond exactly to the usual Keplerian notion of orbital parameters. The inclination angle ι is defined (Glampedakis et al. 2002)

$$\begin{array}{c}\hfill cos\iota =\frac{L_z}{\sqrt{Q+L_z}}.\end{array}$$(15)

We also define i – as distinct from ι- as the inclination of the orbit with respect to the observer,

$$\begin{array}{c}\hfill cosi={\widehat{\mathit{L}}}_z·\widehat{\mathit{O}}\end{array}$$(16)

and $\widehat{\mathit{O}}$ is the line of sight unit vector between the observer and the BH. The observer is at some distant radius with polar and azimuthal coordinates Θ, Φ respectively. The relevant geometry is described in Fig. 1. The mapping between conserved quantities and orbital parameters is defined in (Schmidt 2002; Barausse et al. 2007) as

$$\begin{array}{c}\hfill E=\sqrt{\frac{\kappa \rho +2ϵ\sigma -2D\sqrt{\sigma (\sigma {ϵ}^2+\rho ϵ\kappa -\eta {\kappa }^2)}}{{\rho }^2+4\eta \sigma }},\end{array}$$(17)

$$\begin{array}{c}\hfill L_z=-\frac{g_{1}E}{h_1}+\frac{D}{h_1}\sqrt{g_1^2{E}^2+(f_1{E}^2-{\text{d}}_1)h_1},\end{array}$$(18)

$$\begin{array}{c}\hfill Q=z_{-}[a^2(1-E^2)+\frac{L_z^2}{1-z_{-}}],\end{array}$$(19)

Fig. 1. Schematics of the orbital quantities used in this work. The direction of the BH spin axis is S and the observer lies in direction $\widehat{O}$. The orbital plane of the pulsar is inclined at an angle ι with respect to the BH axis and the angle between the observer and the orbital angular momentum $\widehat{L}$ is given by i. The longitude of ascending mode is given by Ω.

where

$$\begin{array}{c}\hfill z_{-}={sin}^2\iota \end{array}$$(20)

and D = ±1 denotes prograde and retrograde orbits. In turn, the functions are defined as,

$$\begin{array}{c}\hfill f(r)\equiv r^4+a^2[r(r+2)+z_{-}\mathrm{\Delta }]\prime ,\end{array}$$(21)

$$\begin{array}{c}\hfill g(r)\equiv 2ar,\end{array}$$(22)

$$\begin{array}{c}\hfill h(r)\equiv r(r-2)+\frac{\mathrm{\Delta }{z}_{-}}{1-z_{-}},\end{array}$$(23)

$$\begin{array}{c}\hfill \mathrm{d}(r)\equiv \mathrm{\Delta }(r^2+a^2{z}_{-}),\end{array}$$(24)

and for the eccentric orbits considered in this work,

$$\begin{array}{cc}& (f_1,g_1,h_1,d_1)=f(r_p),g(r_p),h(r_p),\text{d}(r_p),\hfill \end{array}$$(25)

$$\begin{array}{cc}& (f_2,g_2,h_2,d_2)=f(r_a),g(r_a),h(r_a),\mathrm{d}(r_a),\hfill \end{array}$$(26)

and

$$\begin{array}{cc}\hfill \kappa & \equiv {\mathrm{d}}_1{h}_2-{\mathrm{d}}_2{h}_1,\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill \epsilon & \equiv {\mathrm{d}}_1{g}_2-{\mathrm{d}}_2{g}_1,\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill \rho & \equiv f_1{h}_2-f_2{h}_1,\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill \eta & \equiv f_1{g}_2-f_2{g}_1,\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill \sigma & \equiv g_1{h}_2-g_2{h}_1.\hfill \end{array}$$(31)

An MSP with particular orbital parameters 𝒫, e will have a corresponding semi-major axis a_* = 𝒫/(1 − e²). As a shorthand we sometimes refer not to the system semi-major axis, but instead to the orbital period P where the two are related via Kepler’s 3rd law, $P^2\propto a_{*}^3$. We use this only as a convenient re-parameterisation with the understanding that this relation only hold exactly in the weak-field, zero-spin limit and that generally the orbital period will be different due to contributions from BH spin.

2.3. Radio beam propagation

We describe the propagation of the pulsar radio signal to the observer in terms of geometrical optics. In this case the general vacuum Hamiltonian is,

$$\begin{array}{c}\hfill H(x^{\mu },k_{\mu })=\frac{1}{2}{g}^{\mu \nu }{k}_{\mu }{k}_{\nu }=0,\end{array}$$(32)

and the corresponding equations of motion, Hamilton’s equations are

$$\begin{array}{c}\hfill \dot{x^{\mu }}=\frac{\partial H}{\partial k_{\mu }},{\dot{k}}_{\mu }=-\frac{\partial H}{\partial x^{\mu }},\end{array}$$(33)

where g_μν is the spacetime metric, k^μ the contravariant 4-momenta, x^μ the spacetime coordinates, and an overdot denotes the derivative with respect to some affine parameter. For a Kerr spacetime and the associated conserved quantities the equations of motion reduce to a problem of quadratures whereby we have four ordinary differential equations $(\dot{t},\dot{r},\dot{\theta },\dot{\varphi })$ and four associated constants of motion E, L_z, Q_p, and H. The system of equations is integrable. It follows that the complete set of equations of motion is given, via Hamilton’s equations, as,

$$\begin{array}{cc}& \dot{t}=E+\frac{2r(r^2+a^2)E-2ar{L}_z}{\mathrm{\Sigma }\mathrm{\Delta }};\hfill \end{array}$$(34)

$$\begin{array}{cc}& \dot{r}=\frac{p_r\mathrm{\Delta }}{\mathrm{\Sigma }};\hfill \end{array}$$(35)

$$\begin{array}{cc}& \dot{\theta }=\frac{p_{\theta }}{\mathrm{\Sigma }};\hfill \end{array}$$(36)

$$\begin{array}{cc}& \dot{\varphi }=\frac{2arE+(\mathrm{\Sigma }-2r)L_z{csc}^2\theta }{\mathrm{\Sigma }\mathrm{\Delta }};\hfill \end{array}$$(37)

$$\begin{array}{cc}& {\dot{k}}_{\theta }=\frac{1}{2\mathrm{\Sigma }}[-2a^2{E}^{2}sin\theta cos\theta +2L_z^{2}cot\theta {csc}^2\theta ];\hfill \end{array}$$(38)

$$\begin{array}{cc}& {\dot{k}}_r=\frac{1}{\mathrm{\Sigma }\mathrm{\Delta }}[-\kappa (r-1)+2r(r^2+a^2)E^2-2aEL]-\frac{2p_r^2(r-1)}{\mathrm{\Sigma }}\hfill \end{array}$$(39)

where $\kappa =p_{\theta }^2+E^2{a}^2{\sin }^2\theta +L_z^2{\csc }^2\theta$. The two integration constants, E (energy at infinity) and L_z (the azimuthal component of the angular momentum at infinity), can be determined by the initial conditions using the following relations:

$$\begin{array}{cc}\hfill E^2& =(\mathrm{\Sigma }-2r)(\frac{{\dot{r}}^2}{\mathrm{\Delta }}+{\dot{\theta }}^2+)+\mathrm{\Delta }{\dot{\varphi }}^2{sin}^2\theta ;\hfill \end{array}$$(40)

$$\begin{array}{cc}\hfill L_z& =\frac{(\mathrm{\Sigma }\mathrm{\Delta }\dot{\varphi }-2arE){sin}^2\theta }{\mathrm{\Sigma }-2r}.\hfill \end{array}$$(41)

We also integrate this system of equations numerically so as to determine the evolution of the photon ray subject to all general relativistic effects.

3. Comparison with Post-Keplerian approximation

Given that we can now describe the relativistic orbital dynamics of the pulsar and the propagation of light between pulsar and observer, we are now in a position to investigate the pulsar timing delays for a MSP-E/IMRB. We can compare the fully general relativistic solution with the timing delays given by the PK approximation. For a given quantity X we can determine four “levels” of solution

1. X_PK, i – solution given by the i-th order PK approximation

2. X_{GEO, a = 0} – GR solution for zero BH spin (a = 0) and zero MSP spin (i.e. Schwarzchild geodesic)

3. X_{GEO, a ≠ 0} – GR solution for non-zero BH spin and zero MSP spin (i.e. Kerr geodesic)

4. X_MPD – GR solution for spinning BH and spinning MSP (i.e. full MPD solution)

We investigate the difference between each of these levels of solution, defining the respective differences as,

$$\begin{array}{c}\hfill {\delta }_{\alpha ,i}(X)=X_{\mathrm{GEO},a=0}-X_{\mathrm{PK},i}\end{array}$$(42)

$$\begin{array}{c}\hfill {\delta }_{\beta }(X)=X_{\mathrm{GEO},a\ne 0}-X_{\mathrm{GEO},a=0}\end{array}$$(43)

$$\begin{array}{c}\hfill {\delta }_{\gamma }(X)=X_{\mathrm{MPD}}-X_{\mathrm{GEO},a\ne 0}.\end{array}$$(44)

It is worth taking a moment to consider how we compare different quantities in different spacetimes. For δ_β(X) and δ_γ(X), the comparison is straightforward; in all cases we are working in Boyer-Lindquist coordinates, with the Schwarzschild coordinates of the Schwarzschild metric simple being a special (a = 0) case. However, solutions derived within the post-Keplerian framework are derived from the PK metric which typically uses a different coordinate system. In particular, by writing the Schwarzschild metric in harmonic coordinates (such that $g^{\mu \nu }{{\mathrm{\Gamma }}^{\alpha }}_{\mu \nu }=0$) where the radial harmonic coordinates r_h is related to the Boyer Lindquist coordinate as,

$$\begin{array}{c}\hfill r=r_h+1\end{array}$$(45)

and by taking the weak field (large r_h) limit, the metric components become,

$$\begin{array}{c}\hfill g_{00}=-1+2/r_h-\frac{1}{2}{(2/r_h)}^2+\cdots \end{array}$$(46)

$$\begin{array}{c}\hfill g_{\mathit{jk}}=(1+2/r_h){\delta }_{\mathit{jk}}+\frac{1}{4}{(2/r_h)}^2({\delta }_{\mathit{jk}}+n_j{n}_k)+\cdots \end{array}$$(47)

where n^j = x^j/r_h and n_j = δ_jkn^k (Poisson & Will 2014). At first order the metric interval is then

$$\begin{array}{c}\hfill \mathrm{d}{s}^2=-(1-\frac{2}{r_h})\mathrm{d}{t}^2+(1+\frac{2}{r_h}){(\mathrm{d}\overline{\rho })}^2\end{array}$$(48)

where ${(\mathrm{d}\overline{\rho })}^2=\mathrm{d}{x}^2+\mathrm{d}{\mathit{y}}^2+\mathrm{d}{z}^2$. This is the first order PK metric and is the one typically used (e.g. Blandford & Teukolsky 1976) for calculating the pulsar timing delays. Explicitly, the second order metric solution is:

$$\begin{array}{cc}\hfill \mathrm{d}{s}^2& =-(1-\frac{2}{r_h}+\frac{2}{r_h^2})\mathrm{d}{t}^2+(1+\frac{2}{r_h}+\frac{1}{r_h^2}){(\mathrm{d}\overline{\rho })}^2\hfill \end{array}$$(49)

$$\begin{array}{cc}& +\frac{1}{r_h^2}{n}_j{n}_k\mathrm{d}{x}^j\mathrm{d}{x}^k.\hfill \end{array}$$(50)

The primary search area for MSP-E/IMRBs is the Galactic centre. Although no MSPs have yet been detected in this region previous surveys can be shown to be insensitive to the population (Rajwade et al. 2017) and so the development of advanced radio facilities with increased sensitivities such as the SKA or the dedicated search methods of NASA DSN, for example – along with the development of appropriate relativistic search algorithms -should render this population detectable. Subsequently, to use an MSP in the Galactic centre as a precision probe of strong-field GR (e.g. Liu et al. 2012; Psaltis et al. 2016) it is important to time a pulsar that is solely affected by the gravitational field of the central BH. However there do exist astrophysical perturbations to the pulsar orbit e.g. gravitational perturbations due to the presence of other masses (Merritt et al. 2010). For orbits with periods ≲0.3 year external perturbations are expected to be negligible (Liu et al. 2012). Observations close to periapsis of some longer, sufficiently eccentric orbits may also be viable. Consequently in this work we consider eccentric systems with appropriately short orbital periods ≲0.3 year. We take the mass of the Galactic centre BH to be M ∼ 4.3 × 10⁶ M_⊙ (Gillessen et al. 2009). The PK timing delays that we investigate are,

1. Einstein delay, Δ_E – due to gravitational and relativistic time dilation of the MSP.

2. Roemer delay, Δ_R – due to the orbital motion of the MSP

3. Shapiro delay, Δ_S – due to the gravitational time dilation of light

We now explore each of these pulsar timing delays in turn

3.1. Einstein delay, Δ_E

The Einstein delay describes the Doppler shift and gravitational time delay (redshift), or the difference between the MSP proper time and the time measured by a distant observer,

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{E}}=t-\tau .\end{array}$$(51)

Determining the Einstein delay from the numerical relativistic solutions is straightforward, since at each integration time step we have the MSP proper time and the observer coordinate time by solving Eq. (12). Within the PK approximation, one can derive an analytical expression for the Einstein delay (Blandford & Teukolsky 1976; Zhang & Saha 2017) as,

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{E}}=\frac{\mathit{Pe}}{\pi a_{\ast }}(sinE\prime -sinE{\prime }_0)+\frac{3}{2}\frac{t}{a_{\ast }}\end{array}$$(52)

where a_* is the semi major axis (as distinct from spin parameter a) and E′ the eccentric anomaly, an angular parameter given by

$$\begin{array}{c}\hfill cosE\prime (t)=\frac{1}{e}(1-\frac{r_h(t)}{a_{\ast }}).\end{array}$$(53)

Evidentially the Einstein delay is then composed of an oscillatory part, of magnitude γ = Pe/πa_* and a linear term. Rather than computing the analytical evolution of Δ_E by e.g. solving the Kepler equation for ℰ, we employ a different approach. From Eq. (48) it follows that,

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}t}{\mathrm{d}\tau }\right)}^2=\frac{1+(1+\frac{2}{r_h})v^2}{1-\frac{2}{r_h}}\end{array}$$(54)

where ${\mathit{v}}^2=|{\mathrm{d}}^2\overline{\rho }/\mathrm{d}{\tau }^2|$. At second order,

$$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}t}{\mathrm{d}\tau }\right)}^2=\frac{1+\frac{1}{r_h^2}{n}_j{n}_k\frac{\mathrm{d}{x}^j}{\mathrm{d}\tau }\frac{\mathrm{d}{x}^k}{\mathrm{d}\tau }+(1+\frac{2}{r_h}+\frac{1}{r_h^2})v^2}{1-\frac{2}{r_h}+\frac{2}{r_h^2}}.\end{array}$$(55)

We can express the squared velocity magnitude in terms of BL coordinates by the straightforward coordinate transformation of a general contravariant vector,

$$\begin{array}{c}\hfill Y^a=\frac{\partial x^a}{\partial x^b}{Y}^b.\end{array}$$(56)

Explicitly, in conjunction with Eq. (45), the fact that dr_h/dr = 1 and the transformation between (spherical polar) harmonic and Cartesian coordinates then,

$$\begin{array}{cc}\hfill v^2& =v_x^2+v_y^2+v_z^2\hfill \\ \hfill & ={(sin\theta cos\varphi u^r+r_{h}cos\theta cos\varphi u^{\theta }-r_{h}sin\theta sin\varphi u^{\varphi })}^2\hfill \\ \hfill & +{(sin\theta sin\varphi u^r+r_{h}cos\theta sin\varphi u^{\theta }-r_{h}sin\theta cos\varphi u^{\varphi })}^2\hfill \end{array}$$(57)

$$\begin{array}{cc}\hfill & +{(cos\theta u^r-r_{h}sin\theta u^{\theta })}^2\hfill \end{array}$$(58)

where the u^α terms are the components of the MSP 4-velocity as given by Eq. (12). This gives us a differential equation for t as determined by the PK metric. It is equivalent to the 0-th term of the relativistic solution described by Eq. (12). As we integrate the complete relativistic solution numerically, we can then at each integration timestep also compute v and r_h and solve Eq. (54). This then gives us t (and hence Δ_E) as given by the PK solution for an MSP which has the same spatial orbital trajectory as described by the full relativistic solution. The Einstein delay given by the respective solutions for a Galactic centre MSP with orbital period P = 0.1 years, in the equatorial plane is presented in Fig. 2 at eccentricities e = 0, 7, 0.8, 0.9. The evolution of the Einstein delay for all the orbits considered has general oscillatory contribution, illustrated by the first term of Eq. (52) and an additional linear term. Rapid variations are observed as the MSP passes through periapsis and is moving at its fastest. After a single orbital period the Einstein delay has accumulated to ∼2 hr for all the solutions and all the eccentricities. Variations in the Einstein delay induced by relativistic effects are seen to occur at smaller scales: at the first level, δ_α, 1(Δ_E) is of the order 10¹ seconds, whilst δ_α, 2(Δ_E)∼10⁻² − 10⁻¹ seconds, with greater magnitude errors for more eccentric orbits which move faster at periapsis and probe a stronger gravitational field. The variations induced by the BH spin for a = +0.6, δ_β(Δ_E) introduce additional variations in the timing delay on the order of 0.1–1 s, again with more eccentric orbits exhibiting larger magnitude variations. The influence of “switching on” the MSP spin and the associated couplings is particularly interesting. Although subdominant to the δ_α, β variations, the MSP spin introduces periodic contributions to the Einstein delay. As before, these corrections peak in magnitude at periapsis and for e = 0.9, δ_γ(Δ_E)∼1 ms, whilst for e = 0.7 the maximal variations are of magnitude ∼30 μs. Radio facilities such as FAST or SKA are able to determine pulse ToA to precisions < 100 ns (Hobbs et al. 2014; Stappers et al. 2018), especially for stable MSPs which are free of glitches. Consequently, small variations in the timing residuals 𝒪 (μs – ms) will be both detectable and important for a consistent, long-term timing solution.

Fig. 2. Top panel: Einstein delay given by the 4 solutions (1 PK + 3 GR) from an MSP orbiting the Galactic centre BH in the equatorial plane at a range of difference eccentricities. The BH mass is set as 4.31 × 106 M⊙, with spin parameter a + 0.6. The delay accumulates up to 2 hrs after one orbital period, with rapid variations as the pulsar goes through periapsis. Bottom panel: differences between the respective solutions δα, 1(ΔE) (blue),δα, 2(ΔE) (red), δβ(ΔE) (green), δγ(ΔE) (orange) at eccentricities e = 0.9, 0.8, 0.7 (solid, dashed, dotted lines respectively).

3.2. Propagation delays: Shapiro and Roemer, Δ_S, Δ_R

In PK timing models for pulsar binaries, the photon propagation time between observer and source is the sum of the Roemer delay (Δ_R, i.e. time delay due to the finite speed of light propagating in a flat spacetime) and the Shapiro delay (Δ_S, time delay due to first order curvature corrections). This can be seen by decomposing the spacetime geodesic in terms of a Minkowskian (0th order) background plus a higher-order curvature perturbation (see e.g. Carroll 2003) i.e.

$$\begin{array}{c}\hfill x^{\mu }=x^{(0)\mu }+x^{(1)\mu }+{\mathit{x}}^{\mathbf{(}\mathbf{2}\mathbf{)}}\mathit{\mu }.\end{array}$$(59)

The background and perturbation tangent vectors are,

$$\begin{array}{c}\hfill k^{\mu }\equiv \frac{\mathrm{d}{x}^{(0)\mu }}{\mathrm{d}\lambda };l^{\mu }\equiv \frac{\mathrm{d}{x}^{(1)\mu }}{\mathrm{d}\lambda };q^{\mu }\equiv \frac{\mathrm{d}{x}^{(2)\mu }}{\mathrm{d}\lambda }.\end{array}$$(60)

By considering the perturbed geodesic equation, it can be shown that,

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\lambda }{l}^{\mu }={-}_{(1)}{{\mathrm{\Gamma }}^{\mu }}_{\rho \sigma }{k}^{\rho }{k}^{\sigma }\end{array}$$(61)

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}\lambda }{q}^{\mu }={-}_{(2)}{{\mathrm{\Gamma }}^{\mu }}_{\rho \sigma }{k}^{\rho }{k}^{\sigma }{-}_{(1)}{{\mathrm{\Gamma }}^{\mu }}_{\rho \sigma }(k^{\rho }{l}^{\sigma }+k^{\sigma }{l}^{\rho })\end{array}$$(62)

where ${}_{(i)}{{\mathrm{\Gamma }}^{\mu }}_{\rho \sigma }$ denotes the −-th order terms of the Christoffel connection. It then follows that the photon propagation time is,

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{prop}}& ={\displaystyle \int }\frac{\mathrm{d}{x}^0}{\mathrm{d}\lambda }\mathrm{d}\lambda \hfill \\ \hfill & ={\displaystyle \int }{k}^0\mathrm{d}\lambda +{\displaystyle \int }{l}^0\mathrm{d}\lambda +{\displaystyle \int }{q}^0\mathrm{d}\lambda \hfill \\ \hfill & ={\mathrm{\Delta }}_{\mathrm{R}}+{\mathrm{\Delta }}_{\mathrm{S}}+{\mathrm{\Delta }}_{\mathrm{S}}^{(2)}\hfill \end{array}$$(63)

where ${\mathrm{\Delta }}_{\mathrm{S}}^{(2)}$ labels the higher (second) order curvature time delay. Within the PK framework, one can derive analytical expressions for the evolution of Δ_R, Δ_S in terms of the system orbital parameters as (Shapiro 1964; Blandford & Teukolsky 1976; Damour & Deruelle 1986; Zhang & Saha 2017),

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{R}}=\stackrel{\sim }{\alpha }(cos{E}^{\prime }-e)+\stackrel{\sim }{\beta }sin{E}^{\prime }\end{array}$$(64)

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{S}}=2ln\left[\frac{1+ecosf}{1-sinisin(\omega +f)}\right]\end{array}$$(65)

where i is the inclination with respect to the observer, defined via Eq. (16), ω is the angle of periapsis and f the true anomaly, with

$$\begin{array}{c}\hfill \stackrel{\sim }{\alpha }=a_{\star }sinisin\omega \end{array}$$(66)

$$\begin{array}{c}\hfill \stackrel{\sim }{\beta }={(1-e^2)}^{1/2}{a}_{\star }sinicos\omega .\end{array}$$(67)

Similar to deriving the Einstein delay in the previous section, we do not use these expressions explicitly but instead proceed via an equivalent route to express the time delays in terms of coordinates rather than orbital parameters. The propagation time between two points r₁, r₂ can be calculated by first determining the temporal components of the tangent vectors by integrating Eqs. (61), (62)

$$\begin{array}{c}\hfill l^0=-2k\mathrm{\Phi }\end{array}$$(68)

$$\begin{array}{c}\hfill q^0=2k{\mathrm{\Phi }}^2\end{array}$$(69)

and so the propagation time is given from Eq. (63) as,

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{prop}}& ={\displaystyle \int }\mathrm{d}s-2{\displaystyle \int }\mathrm{\Phi }\mathrm{d}s+2{\displaystyle \int }{\mathrm{\Phi }}^2\mathrm{d}s\hfill \\ \hfill & =|{\mathit{r}}_2-{\mathit{r}}_1|+2{\displaystyle \int }\frac{1}{|\mathit{r}|}\mathrm{d}s+2{\displaystyle \int }\frac{1}{{|\mathit{r}|}^2}\mathrm{d}s\hfill \end{array}$$(70)

where ds denotes the unperturbed background path. Considering the second integral term (i.e. the Shapiro delay, Δ_S), we can parameterise in terms of t such that,

$$\begin{array}{c}\hfill \mathit{r}={\mathit{r}}_1+\widehat{\mathit{u}}t\end{array}$$(71)

where $\widehat{\mathit{u}}$ is the unit vector,

$$\begin{array}{c}\hfill \widehat{\mathit{u}}=\frac{{\mathit{r}}_2-{\mathit{r}}_1}{|{\mathit{r}}_2-{\mathit{r}}_1|}\end{array}$$(72)

and t ∈ [0, |r₂ − r₁|]. With this parameterisation, $\mathrm{d}s=|\widehat{\mathit{u}}|\mathrm{d}t$ and so,

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{S}}& ={\displaystyle \int }\frac{1}{|{\mathit{r}}_1+\widehat{\mathit{u}}t|}\mathrm{d}t\hfill \\ \hfill & ={\displaystyle \int }\frac{1}{\sqrt{|{\mathit{r}}_1{|}^2+t^2+2t{\mathit{r}}_1·\widehat{\mathit{u}}}}\mathrm{d}t.\hfill \end{array}$$(73)

Performing the integral gives,

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{S}}={[ln(t+{\mathit{r}}_1·\widehat{\mathit{u}}+|{\mathit{r}}_1+\widehat{\mathit{u}}t|)]}_{t=0}^{t=|{\mathit{r}}_2-{\mathit{r}}_1|},\end{array}$$(74)

and so,

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{S}}=ln\left(\frac{|{\mathit{r}}_2-{\mathit{r}}_1|+{\mathit{r}}_1·\widehat{\mathit{u}}+|{\mathit{r}}_2|}{{\mathit{r}}_1·\widehat{\mathit{u}}+|{\mathit{r}}_1|}\right).\end{array}$$(75)

Using the same parametrization for the quadratic term, it follows that

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{S}}^{(2)}& ={[\frac{1}{\zeta }arctan\left(\frac{t+{\mathit{r}}_1·\widehat{\mathit{u}}}{\zeta }\right)]}_{t=0}^{t=|{\mathit{r}}_2-{\mathit{r}}_1|}\hfill \end{array}$$(76)

$$\begin{array}{cc}& =\frac{1}{\zeta }arctan\left(\frac{|{\mathit{r}}_2-{\mathit{r}}_1|+{\mathit{r}}_1·\widehat{\mathit{u}}}{\zeta }\right)-\frac{1}{\zeta }arctan\left(\frac{{\mathit{r}}_1·\widehat{\mathit{u}}}{\zeta }\right)\hfill \end{array}$$(77)

with $\zeta =\sqrt{|{\mathit{r}}_1{|}^2-{({\mathit{r}}_1·\widehat{\mathit{u}})}^2}$ The propagation time along a straight line between two points in the PK framework at second order is then,

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{\mathrm{prop}}& =|{\mathit{r}}_2-{\mathit{r}}_1|+2ln\left(\frac{|{\mathit{r}}_2-{\mathit{r}}_1|+{\mathit{r}}_1·\widehat{u}+|{\mathit{r}}_2|}{{\mathit{r}}_1·\widehat{u}+|{\mathit{r}}_1|}\right)\hfill \end{array}$$(78)

$$\begin{array}{cc}& =\frac{1}{\zeta }arctan\left(\frac{|{\mathit{r}}_2-{\mathit{r}}_1|+{\mathit{r}}_1·\widehat{\mathit{u}}}{\zeta }\right)-\frac{1}{\zeta }arctan\left(\frac{{\mathit{r}}_1·\widehat{\mathit{u}}}{\zeta }\right),\hfill \end{array}$$(79)

where the first term corresponds to the Roemer delay and the second and third terms to the Shapiro delay at first and second order respectively. A linear decomposition of the propagation time in this way is not generally possible for the relativistic solutions, instead being found by solving Eq. (34). Therefore going forward, we consider the Roemer and Shapiro terms together in terms of the general propagation delay Δ_prop.

To explore the propagation delay according to the PK and relativistic solutions, we take a discrete sample of position points of an eccentric orbit. The orbit is calculated by solving the Kerr orbital geodesic equations numerically. The propagation delay between each sample point (i.e. the pulsar emitter) and the image plane of a distant observer (with Boyer Lindquist coordinates R = 10⁴, Θ = π/2, Φ = 0) can be calculated via the PK solution (Eq. (79)) and the relativistic solutions (Eq. (34).) Obviously, since the MSP spin does not influence the photon propagation, δ_γ(Δ_prop) = 0 (however, the MSP spin does generally influence the orbital evolution of the pulsar which in turn will influence the evolution of Δ_R. We return to this point later.) We consider the orbit of a Galactic centre MSP with semi-major axis a_* = 825r_g, ι=, ω = Ω = π/4 and ι = 20° at eccentricities e = 0.7, 0.8, 0.9. The spin parameter is a = +0.6. Since the observers image plane is in the positive x-direction, this corresponds to an inclination with respect to the observer of i = 76°. The propagation delay, along with the associated orbital trajectory, for this system is shown in Fig. 3. Since generally pulsar timing can only detect variations in the light travel time, rather than the absolute travel time itself, we normalise Δ_prop with respect to its initial value. The propagation delay evolves on the scale of hours over the course of the orbit, with more eccentric orbits having a slightly greater magnitude of variation. The error in the PK solutions at first and second order compared to the Schwarzschild relativistic solutions, δ_α, i(Δ_prop) evolve over the orbit and reaches an extrema when the pulsar is on the far side of the BH. In this case the photon tray traverses a more strongly curved spacetime and is subject to a greater degree of strong-field time dilation and gravitational lensing, both of which drive the PK solutions away from the relativistic solution. For this orientation, more eccentric orbits have greater extrema than less eccentric orbits since their far side points are at shorter orbital radii. The additional propagation delays induced by the BH spin, δ_β(Δ_prop), follow a similar general evolution to δ_α(Δ_prop), with greatest variations found for eccentric orbits on the far side of the BH. Variations here are of the order ∼0.1 s. Whilst for both δ_α(Δ_prop) and δ_β(Δ_prop) are greatest for far side sample points, the error remains significant over the entire orbit, especially given the high timing precision of pulsar timing (10 − −100 ns).

Fig. 3. Left panel: orbital trajectories of a Galactic centre MSP with a* = 825rg, i = 76°, ω = Ω = π/4 and eccentricities e = 0.7, 0.8, 0.9 (blue, green, orange lines respectively). The scatter points are included as a reference for the right panels and label the initial, final, middle and superior conjunction sample points (turquoise, purple, brown, red respectively). The black cross labels the BH singularity. Right panels: propagation timing delay for each of the eccentric orbits, along with the errors between the PK and relativistic solutions (δα, i(Δprop), δβ(Δprop)). For the middle two panels the dashed lines show the residual when the PK lensing correction is included (see also Fig. 5). Significant timing residuals are present between both δα, i(Δprop), δβ(Δprop) relative to typical radio pulsar timing precision.

Clearly the exact propagation time is going to depend upon the ray path. To first order in the PK solution rays travel in straight lines. Such an approximation may be appropriate for longer orbits and emission on the near side of the BH, but for more compact orbits where the pulsar is on the far side of the BH it is also important to account for the gravitational lensing which is naturally included by the relativistic solution. Lensing introduces two additional timing delays with the PK framework. The first is a geometrical delay Δ_geo due to the different spatial path traversed given by Rafikov & Lai (2006),

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{geo}}=2\left(\frac{\mathrm{\Delta }{R}_{\pm }}{R_E}\right)\end{array}$$(80)

where R_E is the Einstein radius and ΔR_± = |R_± − R_s| for image position R_± and “flat space position” R_s (see Rafikov & Lai 2006, for a further explanation of these terms). The second delay is due to a correction to the Shapiro delay since the lensed ray traverses a different gravitational potential than the unlensed ray. Indeed, for edge on orbits Eq. (65) diverges when f = π/2 since this ray path involves going directly through the BH singularity, which is evidently unphysical. At first order within the PK framework, the Shaprio delay is modified due to lensing as (Rafikov & Lai 2006),

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{S},\mathrm{L}}=2ln\left[\frac{a_{\ast }(1-e^2)}{\sqrt{r_{||}^2+R_{\pm }^2}-r_{||}}\right].\end{array}$$(81)

The net propagation delay the becomes

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{prop}}={\mathrm{\Delta }}_{\mathrm{R}}+{\mathrm{\Delta }}_{\mathrm{S},\mathrm{L}}+{\mathrm{\Delta }}_{\mathrm{geo}}.\end{array}$$(82)

The linear weak field description of lensing involves decomposing the ray trajectory into two parts. The first is the radial ray which travels from the image plane to the lensing plane. Upon reaching the lensing lane the ray is deflected and travels again in a straight line towards the emitter (see Fig. 4). The total propagation time is then the sum of the propagation time along these two linear trajectories. We can therefore the describe the PK time propagation delay taking in account lensing by evaluating Eq. (79) along the two relevant trajectories and taking the sum. The effect of the inclusion of PK lensing corrections on δ_α, i(Δ_prop)is shown in Figs. 3 and 5. The inclusion of lensing corrections significantly improves the PK solution, with the peak magnitude of δ_α(Δ_prop) reducing to 𝒪(0.1) for both the first and second order solutions, with the second order solution generally exhibiting lower magnitudes than the first. However, this error is still far in excess of the timing precisions that can be achieved via pulsar timing. It is also worth noting that whilst lensing is most prominent for points on the far side of the BH with respect to the observer, due to the BH spin the rays are deflected from their Minkowski straight line geodesics even on the near side, if sufficiently close to the BH.

Fig. 4. Schematic of weak-field, linear gravitational lensing.

Fig. 5. Enhanced view of the two middle right panels of Fig. 3. Error between the PK and Schwarzschild propagation delays δα, i(Δprop), for each of the eccentric systems when lensing corrections are included (dashed lines) and without. Lensing significantly improves the solution, but marked residuals remain. The typical magnitude of δα, 2(Δprop) are less than δα, 1(Δprop), but even in the case where we consider the second order solution with lensing, significant timing discrepancies are present.

As mentioned previously, the MSP spin will generally couple to the background spacetime and so influence the orbital evolution, driving it away from geodesic motion. These positional variations will in turn influence the evolution of the Roemer delay. The Roemer delay is generally given as,

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{R}}(\tau )=\widehat{\mathit{O}}·\mathit{x}(\tau )\end{array}$$(83)

where $\widehat{\mathit{O}}$ is the position unit vector of the observer and x(τ) the location of the orbiting pulsar. In this way positional variations will imprint on the pulsar timing residuals. The difference in the Roemer delay between a non-spinning and a spinning pulsar, i.e. δ_γ(Δ_R), is presented in Fig. 6. Generally the MSP spin-curvature coupling manifests in two ways: an additional component perpendicular to the orbital plane (e.g. Singh et al. 2014) and also a contribution to the precession of periastron (e.g. Li et al. 2019) (which is itself a PK parameter). We consider the same eccentric systems and observer orientation as used previously (a_* = 825r_g, i = 76° etc.). For these orbital parameters, δ_γ(Δ_R) varies periodically, with rapid, large magnitude variations as the pulsar goes through periapsis. For the most eccentric systems (e = 0.9) the extrema of δ_γ(Δ_R) are of magnitude a few μs, whilst for less eccentric systems the magnitude of these variations drops markedly (e.g. |δ_γ(Δ_R)| ∼ 0.1 μs for e = 0.8.) This is on account of more eccentric orbits probing a stronger gravitational field as they pass through periapsis and so the curvature coupling with the MSP spin becomes more significant. Whilst the magnitude of these spin-induced variations are much smaller than the lower order (δ_α, β) differences there are a few points to consider. Firstly, for the most eccentric orbits these variations are within the timing precision of the next generation of radio facilities. Consequently, in order to use these systems for precision parameter estimation of the central BH such small variation may need to be included in the residuals timing solution. Moreover, whilst we have focused on systems with P ∼ 0.1 years, for more compact orbits or different observer orientations this spin induced variation will be more pronounced. For stronger gravitational fields – precisely the region in which we wish to test GR – such effects will be more significant. Since the error between the two solutions grows with time upon repeated periapsis passages, systems which are observed over longer timescales may also need to correct for this effect.

Fig. 6. Error in the Roemer delay between the MPD and Kerr solutions. We consider a Galactic centre MSP with a* = 825rg, i = 76°, ω = Ω = π/4 and eccentricities e = 0.7, 0.8, 0.9 (blue, green, orange lines respectively). Rapid, large magnitude variations are observed as the pulsar goes through periapsis, with the most significant discrepancies – 𝒪(μs) – for the most eccentric systems which probe the strongest gravitational fields.

4. Discussion

Estimates based on PK frameworks suggest that via pulsar timing the parameters such as the mass and spin of the BH at the Galactic centre can be determined to a precision ∼10⁻⁵ − 10⁻³ (Liu et al. 2012; Psaltis et al. 2016). Some authors consider consistent timing campaigns (Liu et al. 2012) whilst others (Psaltis et al. 2016) consider instead dense observations close to periapsis. Both approaches have merit, however observations close to periapsis seems highly desirable for the purposes of minimising external perturbations (Merritt et al. 2010) and maximising the value of the gravitational potential ϵ which can be investigated.

However, our results highlight that these regions that are most scientifically attractive (cf. tests of astrophysics and fundamental physics) are precisely the regions where the accuracy and appropriateness of the PK framework is lacking. Any pulsar timing model is fundamentally a mapping between the proper time of the pulsar and the time recorded by the observer, τ → t. In these regions of high spacetime curvature, the convolution of spin, relativistic and gravitational effects become more prominent and the naive PK framework fails to account for these higher-order effects and the τ → t mapping does not hold. In particular:

• For a P = 0.1 year orbit at the Galactic centre, the error in the Einstein delay according to the PK framework is of the order ∼10 s of seconds compared to the relativistic (Schwarzschild) solution. This error grows with time; after 5 orbits with eccentricity e = 0.9, δ_α(Δ_E) = 55 s. Higher order corrections in the time delay due to the BH spin and the pulsar spin are of the order 10⁻¹ and 10⁻⁴ seconds respectively (see Fig. 2). The most severe errors are seen at periapsis when the pulsar is moving fastest and the gravitational potential is strongest.

• The propagation delay for a similar eccentric Galactic centre system, inclined with respect to the observer at i = 76° varies on the order of hours over the course of a single orbit. The error in the PK solution is δ_α(Δ_prop)∼𝒪(s) whilst effects due to the BH spin introduce errors δ_α(Δ_prop)∼𝒪(0.1s) (Fig. 3). For this particular orbital configuration, these errors are most drastic for the more eccentric orbits and the largest magnitude variations are seen when the pulsar is on the far side of the BH; in this case the photon ray must traverse the most strongly curved spacetime.

• The introduction of first order linear lensing corrections which account for the geometrical time delay Δ_geo and the modified Shapiro delay Δ_S, L dramatically improves the PK propagation solution, reducing δ_α(Δ_prop) to ∼𝒪(0.1s) (see Fig. 4). However, this error is still much greater than typical timing precisions enjoyed by pulsar radio timing. Consequently higher order corrections will be needed for strong field systems if viewed from a sufficiently “edge on” angle. Even for emission on the near side of the BH, the BH spin will cause an additional lensing effect.

• The coupling of the pulsar spin to the background spacetime will cause the orbital dynamics ot divert from that of pure geodesic motion. Such a variation on the coordinate variables will manifest in a change in the evolution of the Roemer delay. This effect is small, with δ_γ(Δ_R)∼1μs for a system with P = 0.1 years and e = 0.9, observer at i = 76°, but may still be significant over longer timescales, for more compact systems, or dense observations close to periapsis.

Evidently in the gravitational strong-field there is significant discrepancy between the PK and GR solutions across all of the parameters considered in this work. The timing delay according to the PK solution disagrees with the GR solution by at least a few μs in all parameters, up to discrepancies as larges as a few s and naturally this error is more severe as relativistic effects become more prominent; i.e. for highly spinning BHs, measurements near periapsis or rays which traverse a spacetime which is highly curved. Consequently, whilst the quoted precisions (e.g. BH mass to precision 10⁻⁵) may be achievable, this work suggests that in EMRB systems which probe the strong-field, higher-order relativistic effects should be included in a complete timing model to accurately determine the prospects of using PSR to test GR. Indeed, the EMRBs discussed in this work are the progenitors of the key LISA target Extreme Mass Ratio Inspirals (EMRIs), where weak-field methods are recognised as being entirely inadequate (Barack & Pound 2019).

The disagreement between the PK and GR solutions in strong-field regimes motivates the necessity of a consistent, relativistic, frequency-dependent timing model (e.g. Kimpson et al. 2019b). In developing this work, it would be highly desirable to use such a strong-field PSR-timing model to generate an accurate, relativistic, set of mock data. The ToAs generated by this relativistic model can then be used for a consistent covariance analysis (e.g. Liu et al. 2012), fed into a pulsar timing software package (e.g. TEMPO2, Hobbs et al. 2006), allowing the current pulsar data analysis algorithms to be assessed in the strong-field (cf. lack of detections in the Galactic Centre) and allow for quantitative estimates on the precision to which the BH parameters can be recovered. We defer this study for future work. In addition, whilst these results are indicative of the appropriateness or otherwise of a PK timing model for MSP-E/IMRBs, for a more definitive work it is necessary to consider a consistent phase-connected solution. We have also not thoroughly explored the whole parameter space instead just considering some typical example system; determining for what sort of set of orbital parameters relativistic effects start to influence the timing solution would also be a worthwhile enterprise.

Beyond the failure of the PK framework in the strong-field, there also exist additional corrections that will influence the photon ToA and pulse profile, such as spatial and temporal dispersion (Kimpson et al. 2019a) precession of the spin axis and relativistic aberration. Being able to accurately model the photon t − ν signal in relativistic regimes is essential to maximise the scientific return of the detection of pulsar systems (including conincident multimessenger detections, e.g. Kimpson et al. 2020a), realistically estimate the scientific prospects (e.g. to what precision can we hope to determine the BH mass?) and compare pulsar timing as a GR probe with complementary methods (e.g. Event Horizon telescope, stellar orbits, gravitational wave astronomy). It is worth noting however that the strengths of the PK method lie beyond the timing accuracy; the advantage of PK is their computational simplicity and the straightforward way in which deviations from GR can be incorporated. The numerical methods discussed in this work are naturally much more computationally demanding than the analytical PK solutions, however it is conceptually straightforward to develop the numerical approach beyond the GR solution and include beyond GR effects via the use of some alternative (non-Kerr) metric e.g. the quasi-Kerr metric of Glampedakis & Babak (2006) and used in e.g. Kimpson et al. (2020b).

In addition to correctly modelling strong-field PSR-EMRBs, another important consideration is the best environment to attempt to use PSRs to undertake precision, strong-field GR tests. Whilst much of the cited literature focuses on the Sgr A* BH, owing to its high mass-to-distance ratio, observations of the Galactic centre are made all the more difficult since it is not a clean environment. Material along the line of sight can cause scattering which results in temporal smearing of the pulse profile, whilst ionized gas can decrease the received flux density. Indeed, since the GC pulsar population is expected to be dominated by MSPs, such environmental effects mean that previous GC searches have been insensitive to the pulsar population (Macquart & Kanekar 2015; Rajwade et al. 2017). Perhaps unsurprisingly then, only one pulsar – the magnetar PSR J1745-29 (Kennea et al. 2013) – has been detected within 10′ of Sgr A*.

Whilst the GC remains an important target, potentially more fruitful hunting grounds are globular clusters and dwarf elliptical galaxies. These systems have exceptionally high stellar densities – up to 10⁶ stars per cubic parsec in the central regions (Freire et al. 2013). Due to the effects of mass segregation and dynamical friction, heavy objects like neutron stars are expected to sink to the centre of these stellar clusters. It is observed that per unit mass there are up to 10³ more pulsars in globular clusters than in the Galactic disk (Freire et al. 2013). Millisecond pulsars are thought to evolve from low mass X-ray binaries and globular clusters are known hosts of abundant LMXB populations. Consequently MSPs comprise a significant fraction of the globular cluster pulsar population (Camilo et al. 2005; Ransom et al. 2008). Some notable clusters are Terzan 5 with 37 pulsars (Cadelano et al. 2018) and 47 Tucanae which is known to host 25 pulsars, all of which have spin periods less than 8 ms (Freire et al. 2017). Observations of X-ray emission from globular clusters led to the first suggestion that globular clusters could host central BHs of intermediate mass (10³ − 10⁵ M_⊙, Silk & Arons 1975; Colbert & Miller 2006; Lin et al. 2018). Moreover, extrapolation to the low-mass end of the “M-σ relation” (Ferrarese & Merritt 2000; Gebhardt et al. 2000) also suggests that globular clusters should host IMBHs (Sadoun & Colin 2012; Graham & Soria 2019). Suggestive evidence for IMBHs is put forth via observational radiative accretion signatures (e.g. Ulvestad et al. 2007) and stellar kinematics (e.g. Gebhardt et al. 2002; van der Marel & Anderson 2010; Feldmeier et al. 2013), but other studies dispute this as evidence for IMBH (e.g. Miller-Jones et al. 2012; Baumgardt 2017). For a complete review and discussion on the current observational evidence for IMBH, see Mezcua (2017).

If globular clusters and dwarf galaxies do host IMBH, then the large stellar densities and expected pulsar population make these ideal, clean environments to search for PSR-BH systems from the perspectives of testing GR. As discussed, there is expected to be a large fraction of MSP which are highly desirable from the perspective of precision astronomy, given their enhanced spin period (i.e. the clock is more precise) and their long term stability (i.e. the clock is more accurate) owing to the lack of glitches – indeed the canonical pulsars are dominated by timing noise and rotational instability, whilst MSP timing models are white noise dominated (Verbiest et al. 2009) and so stand most to gain from the increased sensitivity of the next generation of radio telescopes. Moreover, whilst pulsars in globular clusters are more distant and so more faint, they are also localised in central regions, making it easier for deep radio observations, whilst the increased sensitivity of the next generation of radio telescopes will make it easier to detect their flux.

5. Conclusion

In summary, we have calculated the Einstein delay and the propagation delays for a typical MSP-E/IMRB system using a fully general relativistic method which accounts for both strong-field curvature effects on the photon spacetime trajectory and the impact of extended spinning body dynamics on the pulsar orbital motion. We have then compared this with the linear, weak-field PK solution which is typically used in pulsar studies. We have shown that there exists a significant discrepancy between the post-Keplerian timing delays and the fully general relativistic solutions in strong-field regimes of gravity. Accounting for this discrepancy is essential for both detecting MSP-E/IMRB systems and using them as probes of strong-field GR.

Acknowledgments

We thank the anonymous referee for constructive comments which served to much improve the original manuscript. TK acknowledges the support of a PhD studentship grant ST/N504488/1 awarded by the UK Science & Technology Facilities Council (UKRI-STFC).

References

All Figures

	Fig. 1. Schematics of the orbital quantities used in this work. The direction of the BH spin axis is S and the observer lies in direction $\widehat{O}$. The orbital plane of the pulsar is inclined at an angle ι with respect to the BH axis and the angle between the observer and the orbital angular momentum $\widehat{L}$ is given by i. The longitude of ascending mode is given by Ω.
In the text

Fig. 2. Top panel: Einstein delay given by the 4 solutions (1 PK + 3 GR) from an MSP orbiting the Galactic centre BH in the equatorial plane at a range of difference eccentricities. The BH mass is set as 4.31 × 106 M⊙, with spin parameter a + 0.6. The delay accumulates up to 2 hrs after one orbital period, with rapid variations as the pulsar goes through periapsis. Bottom panel: differences between the respective solutions δα, 1(ΔE) (blue),δα, 2(ΔE) (red), δβ(ΔE) (green), δγ(ΔE) (orange) at eccentricities e = 0.9, 0.8, 0.7 (solid, dashed, dotted lines respectively).

In the text

Fig. 3. Left panel: orbital trajectories of a Galactic centre MSP with a* = 825rg, i = 76°, ω = Ω = π/4 and eccentricities e = 0.7, 0.8, 0.9 (blue, green, orange lines respectively). The scatter points are included as a reference for the right panels and label the initial, final, middle and superior conjunction sample points (turquoise, purple, brown, red respectively). The black cross labels the BH singularity. Right panels: propagation timing delay for each of the eccentric orbits, along with the errors between the PK and relativistic solutions (δα, i(Δprop), δβ(Δprop)). For the middle two panels the dashed lines show the residual when the PK lensing correction is included (see also Fig. 5). Significant timing residuals are present between both δα, i(Δprop), δβ(Δprop) relative to typical radio pulsar timing precision.

In the text

	Fig. 4. Schematic of weak-field, linear gravitational lensing.
In the text

Fig. 5. Enhanced view of the two middle right panels of Fig. 3. Error between the PK and Schwarzschild propagation delays δα, i(Δprop), for each of the eccentric systems when lensing corrections are included (dashed lines) and without. Lensing significantly improves the solution, but marked residuals remain. The typical magnitude of δα, 2(Δprop) are less than δα, 1(Δprop), but even in the case where we consider the second order solution with lensing, significant timing discrepancies are present.

In the text

Fig. 6. Error in the Roemer delay between the MPD and Kerr solutions. We consider a Galactic centre MSP with a* = 825rg, i = 76°, ω = Ω = π/4 and eccentricities e = 0.7, 0.8, 0.9 (blue, green, orange lines respectively). Rapid, large magnitude variations are observed as the pulsar goes through periapsis, with the most significant discrepancies – 𝒪(μs) – for the most eccentric systems which probe the strongest gravitational fields.

In the text