Pulsar timing in extreme mass ratio binaries: a general relativistic approach

TL;DR

This paper develops a comprehensive relativistic framework to model pulsar signals in extreme mass ratio binaries, enabling precise tests of general relativity in strong gravitational fields.

Contribution

It introduces a fully relativistic, no-weak-field approximation model for pulsar timing in EMRBs, incorporating all higher-order gravitational effects and complex spin dynamics.

Findings

Framework accounts for spacetime curvature, frame dragging, and light bending.

Includes plasma effects and relativistic aberration in signal modeling.

Provides a basis for future observational tests of strong-field GR.

Abstract

The detection of a pulsar (PSR) in a tight, relativistic orbit around a supermassive or intermediate mass black hole - such as those in the Galactic centre or in the centre of Globular clusters - would allow for precision tests of general relativity (GR) in the strong-field, non-linear regime. We present a framework for calculating the theoretical time-frequency signal from a PSR in such an Extreme Mass Ratio Binary (EMRB). This framework is entirely relativistic with no weak-field approximations and so able to account for all higher-order strong-field gravitational effects, relativistic spin dynamics, the convolution with astrophysical effects and the combined impact on the PSR timing signal. Specifically we calculate both the spacetime path of the pulsar radio signal and the complex orbital and spin dynamics of a spinning pulsar around a Kerr black hole, accounting for spacetime curvature and frame dragging, relativistic and gravitational time delay, gravitational light bending, temporal and spatial dispersion induced by the presence of plasma along the line of sight and relativistic aberration. This then allows for a consistent time-frequency solution to be generated. Such a framework is key for assessing the use of PSR as probes of strong field GR, helping to inform the detection of an EMRB system hosting a PSR and, most essentially, for providing an accurate theoretical basis to then compare with observations to test fundamental physics.

Equations

$$H(x^{\mu },k^{\mu })=\frac{1}{2}{g}_{\mu \nu }{k}^{\mu }{k}^{\nu },$$

$$\dot{x^{\mu }}=\frac{\partial H}{\partial k_{\mu }},{\dot{k}}_{\mu }=-\frac{\partial H}{\partial x^{\mu }},$$

$$H(x^{\mu },k^{\mu })=\frac{1}{2}\left[g_{\mu \nu }{k}^{\mu }{k}^{\nu }+{\omega }_{\mathrm{p}}(x^{\nu }{)}^2\right],$$

$$\mathrm{d}{s}^2=-\left(1-\frac{2Mr}{\Sigma }\right)\mathrm{d}{t}^2-\frac{4aMr{\sin }^2\theta }{\Sigma }\mathrm{d}t\mathrm{d}\varphi +\frac{\Sigma }{\Delta }\mathrm{d}{r}^2+\Sigma \mathrm{d}{\theta }^2$$

$$+\frac{{\sin }^2\theta }{\Sigma }\left[(r^2+a^2{)}^2-\Delta a^2{\sin }^2\theta \right]\mathrm{d}{\varphi }^2,$$

$${\xi }^{(\mu ;\nu )}=0,$$

$${\omega }_{\mathrm{p}}^2=C\frac{f(r)+g(\theta )}{\Sigma (r,\theta )},$$

$$\dot{t}$$

$$\dot{r}$$

$$\dot{\theta }$$

$$\dot{\varphi }$$

$${\dot{k}}_{\theta }$$

$${\dot{k}}_r$$

$$\displaystyle-\frac{Cf(r)_{,r}\Delta}{2}-C(r-1)f(r)\biggl{]}-\frac{2p_{r}^{2}(r-1)}{\Sigma}$$

$$E^2=(\Sigma -2r)\left(\frac{{\dot{r}}^2}{\Delta }+{\dot{\theta }}^2+\frac{{\omega }_{\mathrm{p}}^2}{\Sigma }\right)+\Delta {\dot{\varphi }}^2{\sin }^2\theta ;$$

$$L_z=\frac{(\Sigma \Delta \dot{\varphi }-2arE){\sin }^2\theta }{\Sigma -2r}.$$

$$T_{;\nu }^{\mu \nu }=0,$$

$$\frac{D{p}^{\mu }}{d\tau }=-\frac{1}{2}{R}_{\nu \alpha \beta }^{\mu }{u}^{\nu }{s}^{\alpha \beta },$$

$$\frac{D{s}^{\mu \nu }}{d\tau }=p^{\mu }{u}^{\nu }-p^{\nu }{u}^{\mu },$$

$$s^{\mu \nu }{p}_{\nu }=0,$$

$$\frac{D{u}^{\mu }}{d\tau }=-\frac{1}{2m}{R}_{\nu \alpha \beta }^{\mu }{u}^{\nu }{s}^{\alpha \beta },$$

$$\frac{D{s}^{\mu \nu }}{d\tau }\approx 0,$$

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

$$\frac{d{s}^{\alpha }}{d\tau }=-{\Gamma }_{\mu \nu }^{\alpha }{s}^{\mu }{u}^{\nu }+\rho \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 },$$

$$\frac{d{x}^{\alpha }}{d\tau }=-\frac{p^{\delta }{u}_{\delta }}{m^2}\left[p^{\alpha }+\frac{1}{2}\frac{\rho (s^{\alpha \beta }{R}_{\beta \gamma \mu \nu }{p}^{\gamma }{s}^{\mu \nu })}{m^2+\rho (R_{\mu \nu \rho \sigma }{s}^{\mu \nu }{s}^{\beta \sigma }/4)}\right],$$

$$x_{\mathrm{rad}}^i=R_z(S_{\varphi })R_y(S_{\theta })R_{PSR}\left[\begin{array}{c}\sin (\psi )\cos (\chi )\\ \sin (\psi )\sin (\chi )\\ \cos (\psi )\end{array}\right]+x_{\mathrm{pulsar}}^i,$$

$${\mathbf{h}}_{j+1}={\mathbf{g}}_{j+1}+{\gamma }_j{\mathbf{h}}_j,$$

$${\gamma }_j=\frac{{\mathbf{g}}_{j+1}\cdot {\mathbf{g}}_{j+1}}{{\mathbf{g}}_j\cdot {\mathbf{g}}_j},$$

$${\mathbf{x}}_{j+1}={\mathbf{x}}_j+{\delta }_j{\mathbf{h}}_j,$$

$$\hat{\omega }=\arccos \left(\frac{n^i{\hat{k}}_i}{|n||k|}\right),$$

$$A^{(\beta )}={\eta }^{(\alpha )(\beta )}{g}_{\mu \nu }{e}_{(\alpha )}^{\mu }{A}^{\nu },$$

$$e_{(t)}^{\mu }=u^{\mu },$$

$$e_{(r)}^{\mu }=(u_r{u}^t,-u_t{u}^t-u_{\varphi }{u}^{\varphi },0,u_r{u}^{\varphi })/N_r,$$

$$e_{(\theta )}^{\mu }=(u_{\theta }{u}^t,u_{\theta }{u}^r,1+u_{\theta }{u}^{\theta },u_{\theta }{u}^{\varphi })/N_{\theta },$$

Full text

††pubyear: 2018††pagerange: Pulsar timing in extreme mass ratio binaries: a general relativistic approach–B.2

1 Introduction

The Galactic centre is thought to host a central black hole (BH) of mass $4.3\times 10^{6}M_{\odot}$ (Gillessen et al., 2009), and be surrounded by up to $10^{3}-10^{4}$ pulsars (PSRs) within the central parsec (Wharton et al., 2012; Rajwade et al., 2017). Moreover - if the “$M$-$\sigma$ relation” (Ferrarese & Merritt, 2000) holds universally - the centre of globular clusters should each contain a nuclear black hole with an ‘intermediate’ mass (IMBH, $M\sim 10^{3}-10^{5}\;\!{\rm M}_{\odot}$, see e.g. Feng & Soria 2011 for discussion of IMBH as ultraluminous X-ray sources and Perera et al. 2017 for the identification of a IMBH candidate in the globular cluster NGC 6624, based on PSR timing observations). The pulsar rate in Globular clusters is thought to be enhanced by a factor of $10^{2}-10^{3}$ per unit mass compared to the Galactic disk (Freire, 2013). Such systems where a pulsar attains a compact orbit around a massive ($\gtrsim 10^{3}M_{\odot}$) BH are known as Extreme Mass Ratio Binaries (EMRBs).

The detection and timing of PSR-EMRB systems presents several challenges. Currently no radio pulsars have been detected within 1 parsec of Sgr A* (Macquart et al., 2010; Rajwade et al., 2017). Neither have we detected pulsars in globular clusters at sufficiently compact orbits. A key challenge is the construction of a timing model which can be applied to the highly relativistic regimes which these systems inhabit. An accurate and precise model of the time-frequency behaviour of a PSR-EMRB is essential for both their detection - especially given the large data flow and real-time processing of the next generation radio telescopes like SKA (Norris, 2011) - and also in order to utilise these systems as a scientific apparatus, and to best exploit the vastity of physical information which may be encoded in observational data.

The main problem is that, in order to be realistic (and feasible for a concrete application to timing data), an accurate theoretical model of the radio signal of a PSR in an EMRB must consider a series of effects, which can be broadly classified into two categories:

Effects that influence the behaviour of the light through a modification of the geodesic followed by a radio pulsar light ray while traveling through a curved spacetime. In this category we can list, e.g., gravitational bending, gravitational and relativistic time dilation, relativistic Doppler frequency shift (Fuerst & Wu, 2004; Saxton et al., 2016), velocity induced intensity boost and also some astrophysical effects (e.g. the temporal and spatial dispersion induced by the interaction with line-of-sight material: interstellar medium, accretion disks, stellar winds Fuerst & Wu, 2007; Younsi et al., 2012; Kimpson et al., 2019). 2. 2.

Effects related to the relativistic orbital dynamics of a spinning pulsar around a spinning black hole, e.g. spin coupling, curvature coupling and time dilations (gravitomagnetic, relativistic, gravitational).

The two categories are naturally linked: the orbital dynamics of the pulsar provides varying initial conditions for the light ray trajectory to the observer, and both effects work in tandem to influence the key observables, i.e. the photon time of arrival (ToA) and photon frequency.

This work provides a first step in this direction. By using a ray tracing method, we develop a fully general relativistic approach to calculate the signal from a PSR in an EMRB, accounting for both spin and orbital motion. The method naturally includes phase-dependent relativistic effects (which are important for the relatively longer orbital periods of EMRB systems, see discussion in Zhang & Saha, 2017) and spin precession. We work in an approximated framework inasmuch we restrict our study to the extreme mass ratio of EMRB systems and so do not consider PSRs in stellar-mass black hole binaries with finite mass ratios (e.g. Blanchet, 2014; Liu et al., 2014). Also, in this work we do not address the task of how to perform mock data analysis nor consider external Newtonian perturbations due to the gravitational foreground from other bodies(see § 6 for a discussion of both of these extensions).

Beside being important “per se”, the discovery of a PSR in an EMRB may provide us with a powerful laboratory to investigate fundamental aspects of general relativity (GR) in regimes of large gravitational potential, $\epsilon={GM}/{rc^{2}}$, and spacetime curvature, $\xi_{c}=GM/r^{3}c^{2}$ (see e.g. Kramer et al., 2004; Wang et al., 2009a, b; Liu et al., 2012; Remmen & Wu, 2013; Nampalliwar et al., 2013; Singh et al., 2014; Kramer, 2016; Saxton et al., 2016; Li et al., 2018). This may be achievable thanks to both, the extreme gyroscopic stability of pulsars (e.g. Verbiest et al., 2009) and the high-precision measurements that are attainable with radio pulsar timing observations (e.g. Liu et al., 2011; Desvignes et al., 2016; Lazarus et al., 2016; Liu et al., 2018). Several aspects of GR have been tested experimentally, with electromagnetic observations of the solar system (Will, 2014), binary pulsar systems (Lorimer, 2008) and also with the recent breakthrough observation via gravitational radiation of binary BH-BH and BH-Neutron Star coalescing systems (Abbott et al., 2017a, b). However there remain open questions in GR, particularly concerning the non-uniqueness of the Einstein-Hilbert action (Psaltis, 2008), the existence of singular matter densities and curvatures (Pachner, 1970), in addition to problems regarding spinning objects in gravity, (see e.g. Plyatsko & Fenyk, 2016) and the spacetime dynamics associated with multiple spinning objects revolving around each other (Obukhov et al., 2009; d’Ambrosi et al., 2015). All existing observations based on electromagnetic radiation have probed the weak-field regime, where $\epsilon\lesssim 10^{-6},\xi_{c}\lesssim 10^{-28}$. The exploration of a more extreme parameter space is achievable through gravitational wave observations of binary black hole and binary BH-neutron stars mergers, or pulsar timing of extreme (and still undetected) systems, as those considered in this work. In particular, the expected precision attainable with pulsar timing potentially allows for system parameters to be determined with a remarkable accuracy: for a pulsar in a 0.3 year period orbit around Sgr A∗, determination of the post-Keplerian parameters allows the mass and spin to be determined to a precision of $10^{-5}$ and $10^{-3}$, respectively (see Liu et al., 2014). At shorter orbital radii, relativistic effects would become more pronounced, potentially allowing for greater precision tests and probing deeper into the strong-field regime.

If PSRs in an EMRB will be discovered, their precision timing with the next generation of radio telescopes such as SKA (Shao et al., 2015) or FAST (Smits et al., 2009) would then allow for an unrivaled probing of physics in an extreme gravity environment and precision measurements of the black hole parameters. This in turn may open a vast array of opportunities: with the determination of the mass and spin of the black hole, the spacetime is wholly determined and key questions of GR can then be investigated; the Cosmic Censorship Conjecture and the No Hair Theorem can then be possibly tested to $\leq 1\%$ precision (see e.g. Kramer et al., 2004; Liu et al., 2012, 2014; Wex & Kopeikin, 1999; Eatough et al., 2015). By extension, the nature of the central massive compact object can also be investigated, i.e. whether the Kerr solution is a true astrophysical solution or if parametrized deviations from the Kerr solution exist, so called ‘bumpy’ black holes (Yagi & Stein, 2016), or even whether the massive central compact object is some more exotic object like a boson star (Kleihaus et al., 2012). Pulsar timing of EMRBs can also be used to test alternative theories of gravity such as scalar-tensor theories (Liu et al., 2014), search for quantum gravitational effects (Yagi & Stein, 2016; Estes et al., 2017) and constrain the cosmological constant (Iorio, 2018). Moreover, EMRB systems are also of great interest in their own right, as they are progenitors of the extreme-mass-ratio-inspiral (EMRI) gravitational-wave sources, a major class of targets for the future LISA observations (Gair et al., 2010; Gair et al., 2017; Babak et al., 2017). Furthermore, such a system can be used for astrophysical purposes, e.g. precisely determining the mass of black holes at the centre of Local Group galaxies and globular clusters in order to constrain the low end of the $M-\sigma$ relation (Ferrarese & Merritt, 2000) and establishing the existence of otherwise of intermediate mass black holes (Singh et al., 2014).

In this work we compute the signal emitted by a spinning PSR in orbit around a Kerr black hole, accounting for the interaction between the radiation and the material close to the black hole, and by using a ray tracing method. The orbital dynamics of the spinning pulsar is computed considering the coupling between the pulsar spin and the spacetime curvature. We derive the theoretical time-frequency signal from a PSR in an EMRB, and discuss the key factors which influence the signal.

This paper is structured as follows. In Section 2 we present the equations for ray tracing a photon that propagates through a non-gravitating plasma in the Kerr spacetime. In Section 3 we review how to determine the orbital dynamics of an extended spinning body (i.e. a PSR) in the extreme mass ratio limit. In Section 4 we then combine the ray tracing and orbital dynamics calculations, and present a method for constructing the complete time-frequency signal from the pulsar. In Section 5 we apply this method to investigate the various relativistic factors which influence the photon time of arrival. Discussion and conclusions follow in Section 6.

We adopt the natural units, with $c=G=\hbar=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_{\rm g}=M$ and the corresponding Schwarzschild radius is $r_{\rm 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 Ray Tracing of light under gravity and influence of plasma

2.1 Formulation

The Hamiltonian for the photon propagation in vacuo is

[TABLE]

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

[TABLE]

where $g_{\mu\nu}$ is the spacetime metric, $k^{\mu}$ the contravariant 4-momenta, $x^{\mu}$ the spacetime coordinates, and an overdot denotes the derivative with respect to some affine parameter.

For photon propagation through a cold, non-magnetized electron-proton plasma, the Hamiltonian has an additional terms, expressed in terms of the electron plasma frequency $\omega_{\rm{p}}(x^{\mu})$:

[TABLE]

in the geometrical optics approximation (Synge, 1960). The plasma frequency of a cold non-magnetised plasma is $\omega^{2}_{\rm{p}}=C\;\!n$, where $C=4\pi e^{2}/m_{e}$ with $m_{e}$ the electron mass, $e$ the electron charge, and $n$ is the electron number density.

Due to the conservation of angular momentum, astrophysical black holes are thought to be rotating. The vacuo spacetime is then described by the Kerr metric, with spacetime interval,

[TABLE]

where $\Sigma=r^{2}+a^{2}\cos^{2}\theta$, $\Delta=r^{2}-2Mr+a^{2}$, $a$ is the black-hole spin parameter and we work in Boyer-Lindquist coordinates. Going forward, we normalize the black hole mass $M=1$ and so the spacetime interval is then a lengthscale in terms of the gravitational radius $r_{\rm g}(=1)$.

A vector $\xi^{\mu}$ which satisfies the Killing equation,

[TABLE]

is described as a Killing vector. Killing vectors are associated with the existence of spacetime symmetries since the inner product of a Killing vector with a tangent vector is conserved along a geodesic, i.e. if $\mathcal{K}_{\mu}=\xi^{\mu}k_{\mu}$ then $\dot{\mathcal{K}_{\mu}}=0$. The Kerr spacetime possesses two Killing vectors $\xi^{t},\xi^{\phi}$. The temporal Killing vector related with the stationarity of the spacetime is associated with the conservation of energy, $k_{t}=-E$. The azimuthal Killing vector is related with the asymmetry of the spacetime and associated with the angular momentum, $k_{\phi}=L_{z}$, where $L_{z}$ is the projection of the particle angular momentum along the black hole spin axis. In addition to the conservation of $E,L_{z}$, the particle rest mass rest mass is also conserved ($H=0$ for photons). Moreover, the Kerr spacetime possesses an additional rank-2 Killing tensor which in the vacuum case can be associated with Carter’s Constant, $Q$, (Carter, 1968), related to the separability of the Hamiltonian in $r$ and $\theta$ terms. The exact physical meaning of $Q$ is discussed in De Felice & Preti (1999); Rosquist et al. (2009).

For the general plasma density distributions, where the plasma frequency, $\omega_{\rm{p}}(r,\theta)$, has a spatial dependence, the Hamiltonian is not separable in terms of the usual co-ordinate variables and the vacuum Carter constant is no longer a constant along the geodesic, (see e.g. Perlick & Tsupko, 2017; Kimpson et al., 2019). In order for the equations of motion to be integrable, it is necessary for $\omega_{p}$ to take the form,

[TABLE]

where $f(r)$ and $g(\theta)$ are two general functions of the coordinate variables $r,\theta$. This represents a plasma density distributions with independent radial and polar terms. From the standpoint of astrophysical modelling, the form of $\omega^{2}_{\rm{p}}$ has specific advantageous properties, allowing the description - with an appropriate choices of $f(r)$ and $g(\theta)$ - of the characteristic features of an axisymmetric density profile. We recognize that this form of $\omega^{2}_{\rm{p}}$ will fundamentally be an approximation to the true astrophysical plasma frequency profile. However, we adopt this functional form in our calculations to serve as a decent first-order approximation to the type of plasma density profiles we expect. This then allows for the development of useful insights into the impact of non-vacuum effects on the time-frequency model signal from PSR in EMRB systems and enables the construction of more advanced model calculations that allow us to extract information from observations to test fundamental aspect gravitational physics in the future. The separable form of Eq. 6 allows us to define an analogous Carter constant, $Q_{p}$, for the ‘plasmic’ Hamiltonian (that is, Eq. 3). The equations of motion therefore reduce to a problem of quadratures whereby we have four ordinary differential equations $(\dot{t},\dot{r},\dot{\theta},\dot{\phi})$ 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 (Kimpson et al., 2019),

[TABLE]

where $\kappa=p_{\theta}^{2}+E^{2}a^{2}\sin^{2}\theta+L_{z}^{2}\csc^{2}\theta+a^{2}\omega_{\rm p}^{2}\cos^{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:

[TABLE]

Whilst it is possible to express the vacuum equations in terms of elliptic integrals (see e.g. Dexter & Agol, 2009), for our case we directly integrate the set of equations numerically (e.g. Fig 1) along the ray using a 5th-order Runge-Kutta integrator (Press et al., 1996). The complete specification of the initial conditions is described in the Appendix.

3 Spinning object dynamics in GR

The most basic approximation in GR is that of the test particle which has no spin, no internal structure and is not subject to self-force effects. Such a particle then directly follows a geodesic of the spacetime metric. However, real astrophysical objects are not in actuality test particles and to obtain an accurate description of their dynamics, higher-order effects must be considered. For a system in which one spinning object (e.g. a PSR) orbits a much more massive one (e.g. Galactic centre Kerr black hole), it is possible to construct two Hamiltonians in the extreme mass ratio and non-relativistic limit, one for the spin-orbit interaction and one for the spin-spin coupling (e.g. Iorio, 2012). This approach is sufficient to reproduce the conventional Lense-Thirring precession effects, however it neglects the spin-curvature coupling between the two masses. This coupling occurs because the spin of the small mass itself will modify the Kerr metric of the spacetime. Consequently, in the presence of spin-curvature coupling, the PSR will not follow a geodesic through the Kerr spacetime. Instead we must model the PSR motion through an alternative means, known as the Mathisson-Papatrou-Dixon (MPD) formalism (Mathisson, 1937; Papapetrou, 1951; Dixon, 1974).

3.1 Spin-Interaction in MPD formalism

The most general equation of motion of the PSR is given by,

[TABLE]

where $T^{\mu\nu}$ is the energy-momentum tensor. The ‘gravitational skeleton’ can be constructed via a multipole expansion of the energy-momentum tensor, defined with respect to some reference world line $z^{\alpha}(\lambda)$. Since the black hole mass ($M$) is much greater than pulsar mass ($m$), the pulsar can be treated as a test mass and its motion entirely determined by the background black hole spacetime and the dynamical spin interaction with this field. In this extreme mass ratio limit and since the typical pulsar radius $R_{\rm PSR}\ll r_{\rm g}$, moments greater than the quadrupole can be neglected. The first two moments are the 0th mass moment, encoded in the 4-momentum $p^{\mu}$ and the 1st dipole moment, given by the spin tensor $s^{\mu\nu}$. The corresponding equations of motion are (Mathisson, 1937; Papapetrou, 1951; Dixon, 1974),

[TABLE]

[TABLE]

where $D/d\tau$ denotes a covariant derivative with respect to the proper time along the PSR worldline, $u^{\nu}$ is the PSR 4-velocity and $R^{\mu}_{\nu\alpha\beta}$ the Riemann curvature tensor. This system of equations is not determinate, since there exist more unknowns than equations. This is related to the uncertainty in choosing a reference world line for the multipole expansion. A natural choice of a representative point of the bulk motion of the body is the centre of mass. However in GR the centre of mass of a spinning body is not invariant. It is therefore necessary to specify a spin supplementary condition (SSC) which renders the system of equations determinate. This is equivalent to choosing an observer with respect to which the centre of mass is defined. We adopt the Tulczyjew-Dixon (TD) condition,

[TABLE]

(Tulczyjew, 1959; Dixon, 1964). This is equivalent to choosing the centre of mass as measured in the zero 3 momentum frame. This SSC is advantageous since it specifies a unique worldline, whilst other choices of SSC are infinitely degenerate (for discussion see Filipe Costa & Natário, 2014). Now, since $m\ll M$ and the pulsar Moller radius $R_{\rm Moller}$ (the radius of the disk of all possible centroids) is much less than the radius of the pulsar, the pole-dipole terms are much stronger than the dipole-dipole terms. Therefore, to first order the 4-velocity and 4 momentum are parallel, i.e. $p^{\mu}\approx mu^{\mu}$. The equations of motion then become,

[TABLE]

[TABLE]

(Chicone et al., 2005; Mashhoon & Singh, 2006). The ordinary differential equations to then be integrated are (Singh, 2005; Mashhoon & Singh, 2006):

[TABLE]

[TABLE]

[TABLE]

where $s^{\mu}$ is the spin 4-vector and the dimensionless parameter $\rho$ is used to label the terms which contribute to MPD spin-curvature coupling ($\rho=1$ includes spin-curvature coupling, for $\rho=0$ the coupling is omitted). In the $\rho\rightarrow 0$ limit the conventional spin-spin and spin-orbit couplings are recovered. Spin-curvature coupling causes additional complexities in the pulsar’s orbital motion; a particle initialized in the orbital plane will remain in the orbital plane without spin-curvature coupling, but lifts in the vertical $z$-direction in the presence of coupling e.g. Fig 2 (see also Singh et al., 2014). We integrate this set of ODE’s numerically via a standard 4th-order Runge-Kutta algorithm (Press et al., 1996). This completely determines the motion of a pulsar around a spinning black hole. At any integration step we can determine not only the pulsar’s spacetime coordinates, but also the tangent momentum 4-vector and the spin 4-vector. The convolution of these vectors leads to complex relativistic dynamics, spin-axis precession (e.g. Figs 3, 4), relativistic aberration, relativistic and time dilation and gravitational and relativistic energy shift. The influence on the pulse arrival time, profile and energetics will be discussed in Section 5.

4 Constructing the PSR timing signal

We now combine the ray tracing and MPD orbital dynamics tools to create a theoretical timing model. Whilst it is possible to use a ‘forwards-in-time’ method and integrate a ray from the pulsar outwards, this is computationally wasteful since it involves integrating many geodesics that do not hit the observer’s image plane. Instead we use a ‘backwards-in-time’ approach and integrate rays from the observer image plane towards the PSR-BH system, the method outlined in Section 2. For each integration step on the PSR orbit, we want to find the ray which intersects with this location and arrives at the distant observer. Were pulsars isotropic emitters, an intersection (to within some tolerance, e.g. $|x^{i}_{\rm pulsar}-x^{i}_{\rm photon}|<R_{\rm PSR}$, such that the ray strikes the PSR surface) would be a sufficient condition to declare that the pulsar is ‘seen’ by some distant observer. However, pulsars clearly do not emit isotropically but instead have characteristic beamed radiation. We therefore instead search for an intersection not with the pulsar’s centre of mass, but instead with a ‘radiation point’, $x^{i}_{\rm rad}$, taken as being the point which lies at a radius $R_{\rm PSR}$ from the PSR centre of mass, in the direction of the radiation beam $\mathbf{n}(\tau)$ (see Fig. 6a) . If the radiation axis is at angles $(\psi,\chi)$ from the pulsar spin axis $\mathbf{S}$, and the spin axis is at angles $S_{\theta},S_{\phi}$ from the coordinate $z$-axis (i.e. parallel to the black hole spin axis), then we can transform from the pulsar centre of mass $x^{i}_{\rm pulsar}$ to the radiation point, accounting for the misalignment of the spin and global vertical $z$-axis as,

[TABLE]

where $R_{z}$ and $R_{y}$ are 3-dimensional rotation matrices about the coordinate $z$ and $y$ axes respectively. Typically, $\psi$ is set to be some constant angle corresponding to the case where the latitude of the radiation beam does not evolve with respect to the spin axis. The phase $\chi=2\pi\tau/P$, for pulsar period $P$. The spin axis angles are time dependent, i.e. $S_{\theta}=S_{\theta}(\tau)$, $S_{\phi}=S_{\phi}(\tau)$ and are related to the precession and nutation of the spin axis with their evolution described by the MPD formalism. The evolution of $x^{i}_{\rm rad}$ and $x^{i}_{\rm pulsar}$ is illustrated in Fig. 5.

4.1 Finding an intersection

The problem of finding an intersection is effectively a two point boundary value problem; we know the desired final boundary value (i.e. radiation point $x^{i}_{\rm rad}$) and we want to find the initial boundary value, i.e. the values of the image plane coordinates $\alpha,\beta$. The ray then traces between these two limits. We solve this problem numerically using a shooting method (Press et al., 1996). The ray tracing is effectively a black box function $f(\alpha,\beta)$ which takes some $\alpha,\beta$ and returns the minimum distance along the ray from the target value $x^{i}_{\rm rad}$ for each time-step of the PSR orbital integration. We want to minimize $ds=f(\alpha,\beta)$, to within some tolerance $\epsilon$ by varying initial conditions $\alpha,\beta$. The adjustment of $\alpha,\beta$ proceeds via a non-linear conjugate gradient descent algorithm (Fletcher & Reeves, 1964; Press et al., 1996), where the new computed gradient is conjugate to previous gradients. The conjugate direction vector $\mathbf{h}_{j}$, is updated at each iteration step $j$ as,

[TABLE]

where $\mathbf{g}_{j}=-\nabla f(\alpha_{j},\beta_{j})$ and

[TABLE]

(Fletcher & Reeves, 1964; Press et al., 1996). The variables $\alpha,\beta$ are then updated at each iteration as,

[TABLE]

for vector $\mathbf{x}_{j}=(\alpha_{j},\beta_{j})$ and where $\delta_{j}$ is the variable stepsize, determined via an inexact line search. This method is advantageous over a vanilla implementation of gradient descent since for certain target points $x^{i}_{\rm rad}$, the function $f(\alpha,\beta)$ takes the form of an ill-conditioned narrow valley. In such an environment, gradient steepest descent becomes inordinately slow, since the direction of steepest descent is not, in general, in the direction of the minimum and the algorithm instead follows a ‘criss-cross’ pattern, oscillating between the sides of the valley. This issue is avoided by moving in a direction which is conjugate to previous directions. Since the analytical form of $f(\alpha,\beta)$ is unknown, the gradients necessary for the optimization are evaluated numerically via the difference quotient. The algorithm exhibits dependable and fast convergence and for this work we declare the algorithm to have converged - i.e. an intersection to be found - when $ds^{2}<10^{-19}$ natural units,working in quadruple precision. Taking the black hole mass to be $4\times 10^{6}M_{\odot}$, then this corresponds to a spatial distance of $\sim 1.9m$ or a light travel time of $\sim 6$ ns. The intersection of multiple rays with a section of pulsar orbit is illustrated in Fig 6b. More complex ray paths and intersections are presented in Section 5.

4.2 Pitch angle in the comoving local frame

Once an intersection with the radiation point is found, in order to determine if the beam is seen by a distant observer we require the pitch angle $\tilde{\omega}$, i.e. the angle the photon ray makes with the normal to the stellar surface. All the MPD and ray tracing calculations take place in the coordinate frame, and so any calculation will result in a pitch angle also measured in the coordinate frame. However the angle that determines whether or not the beam will be visible to an observer is the angle as measured in the orthonormal tetrad frame comoving with the star, $\hat{\omega}$. Due to relativistic aberration these two angles are not the same, $\tilde{\omega}\neq\hat{\omega}$. The pitch angle in the comoving frame is given by,

[TABLE]

where $n^{i}$ is the vector normal to the stellar surface and $\hat{k}_{i}$ is a vector tangent to the ray path, as measured in the comoving frame. Naturally the the covariant 4-momentum of the ray $k_{\mu}$, is tangent to the ray. We can transform a general vector $A^{\nu}$ from the coordinate basis to the tetrad basis as,

[TABLE]

where $\eta^{(\alpha)(\beta)}$ is the Minkowski metric, indices with braces (e.g. $A^{(\beta)}$) indicate the comoving basis and the basis 4-vectors are,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

(Krolik et al., 2005; Kulkarni et al., 2011; Shcherbakov & Huang, 2011; Dexter, 2016) where,

[TABLE]

[TABLE]

[TABLE]

for transformation to a frame moving with 4-velocity $u^{\mu}$ at position $(r,\theta,\phi)$. We can then transform the ray tangent vector to the comoving frame, $k_{\mu}\rightarrow k_{(\mu)}$. The comoving frame is locally flat and so indices are raised and lowered with the Minkowski metric $\eta^{(\mu)(\nu)}$. The vector normal is $n^{i}=x^{i}_{\rm rad}-x^{i}_{\rm pulsar}$ and considering only the spatial components of the comoving ray tangent vector, the pitch angle can be calculated via Eq. 29. Once $\hat{\omega}$ has been determined, we declare an observation if $\hat{\omega}<\omega_{c}$ where $\omega_{c}$ is some critical angle, i.e. the pulsar jet opening angle.

The variation in the pitch angle over 2 rotations of a millisecond pulsar is presented in Fig. 7a, for pitch angle calculated in both the global ($\tilde{\omega}$) and comoving ($\hat{\omega}$) frames. Whilst the periodicity of $\hat{\omega}$ and $\omega$ are the same, $\hat{\omega}$ is shifted in both time and amplitude with respect to $\omega$. If we change the latitude angle of the spin axis from $S_{\theta}=0$ to $S_{\theta}=\pi/4$ (Fig. 7b) the amplitude of the pitch angle reduces in both the global and comoving frames since the radiation beam is shifted further from the observer’s line of sight. Furthermore, the amplitude of the pitch angle in the comoving frame is now markedly greater than in the global frame.

5 Application

We have established the framework for determining the pulsar ray path and orbital dynamics, along with the algorithm to find the intersection of the ray and calculate the pitch angle. We can now apply this framework to investigate the impact of the convolution of relativistic and astrophysical effects on the signal from a PSR in an EMRB, in particular the photon time of arrival and time-frequency profile.

Clearly, any pulsar timing model is going to depend on both the BH mass and spin. The mass of the BH at the centre of the Milky Way is thought to be $\sim 4\times 10^{6}M_{\odot}$ (Gillessen et al., 2009). If intermediate mass black holes exist at the centre of globular clusters then the advanced capabilities of next-generation radio telescopes may allow for the detection of an extragalactic EMRB of the Local Group (Keane et al., 2015). In this case the mass would be more intermediary, $\sim 10^{3}-10^{5}M_{\odot}$ (Wrobel et al., 2018). In contrast to the mass determination, there remains considerable uncertainty in the spin parameter of astrophysical BHs. Astrophysical measurements of the nearest BH source, Sgr A*, range from $0.44-0.996$ (Aschenbach, 2010; Kato et al., 2010; Dokuchaev, 2014). The determination of the spin of BHs other than Sgr A* have been possible via continuum fitting and iron lines measurements, and again cover a wide range from $a<-0.2$ to $a>0.98$ (see Table 1 of Bambi et al., 2016). Going forward we set $M=4\times 10^{6}M_{\odot}$ and have the BH as rapidly spinning such that $a=0.998$. These parameters are adopted in our demonstrative calculations so as to illustrate the sorts of effects a strong field pulsar timing model must account for. Exploration of the effects of alternative spin parameters and masses and their consequent impact on the photon time of arrival will take place in a later paper.

5.1 Effects of relativistic spin orbital dynamics

Treating the pulsar as a rigid spinning object and accounting for the associated spin couplings causes substantial variation in the pulsar orbital dynamics. In particular, spin-curvature coupling causes the pulsar to deviate from geodesic motion in a Kerr spacetime and exhibit complex orbital dynamics. In particular there is noticeable motion in the vertical direction, out of the orbital plane, that would not be present in the $\rho=0$ case. The severity of the effects of spin-curvature coupling become more pronounced for lower mass ratios (e.g. IMBH in globular clusters), for faster spinning objects (e.g. millisecond pulsars) and retrograde orbits (see e.g Singh et al., 2014). These dynamical spin effects will influence the pulse time of arrival. Figure 8 illustrates the change in the light travel time, in vacuum, between a photon emitted at the same orbital phase in the $\rho=0$ and $\rho=1$ cases over a single orbit. We neglect any complications induced by secondary pulses from strong gravitational bending (see Sec 5.3). The light travel time varies on the order of $\sim\pm 10\mu$s (similar values are found in Singh et al., 2014) a variation which should be detectable with the next generation of radio telescopes; 10 minute integrations of millisecond pulsars at 1.4 GHz should achieve a ToA precision of the order $100$ ns (Liu et al., 2011). Furthermore, any variation in the light travel time due to the orbital dynamics will be further compounded by the consideration of additional factors (e.g. time dispersion, time dilation, spin axis precession etc. see subsequent sections).

An accurate description of the PSR orbital dynamics is essential, since the apparent or observed pulse frequency is related to the intrinsic frequency

[TABLE]

In order to determine the Doppler shift accurately, one must therefore have an appropriate description of the PSR velocity. Furthermore, for systems with short orbital periods, the line of sight velocity may not be constant over an integration period. The complex orbital dynamics of a spinning pulsar around a spinning black hole will also naturally influence both the acceleration and the change in acceleration (i.e. jerk) of the pulsar, depending on the orbital and observer configuration. Highly accelerated systems - the most scientifically interesting systems from the perspectives of testing GR - are difficult to search for when the integration time is comparable to the orbital period. Assuming the acceleration to be constant, the changing velocity of the pulsar causes the signal to drift into multiple frequency bins, with the number of bins drifted,

[TABLE]

(Eatough, 2009), for constant line-of-sight acceleration $a_{0}$ and integration time $T$. Search algorithms can correct for this effect (‘acceleration searches’), but the correction is typically computationally intensive. Spin-curvature coupling can cause additional acceleration as measured by the observer, depending on the line of sight, and will influence the number of frequency bins drifted. Furthermore, in the presence of spin-curvature dynamics the assumption of constant acceleration may need modification, since in highly relativistic regimes the jerk may be non-negligible. In this case the number of frequency bins drifted becomes time dependent,

[TABLE]

at time $t$ with jerk $j_{0}$, which may introduce additional complications in searching for highly relativistic systems.

5.2 Effects of spin axis precession

Since the motion of the pulsar is inextricably tied to its spin, the geodetic precession of the spin axis works in conjunction with the relativistic aberration to influence the observability of the pulsar signal, the pulse profile (e.g. Rafikov & Lai, 2006) and the pulse arrival time. Over the course of its orbital motion the pulsar spin axis oscillates (e.g. Fig. 3), with the timescale of oscillation shortening for more compact orbital radii where the gravitational curvature is greater. For some orbital configurations, this oscillation will leave the pulsar beam invisible at certain times (e.g. Istomin, 1991) as $\hat{\omega}>\omega_{c}$. The pulsar signal would consequently be intermittent, leading to additional complications in search observations and any consistent timing model.

To inspect the effects of spin axis precession, we consider variations in $S_{\theta}$ whilst holding $S_{\phi}$ constant at $S_{\phi}=0$. (Fig 9). We arbitrarily define the pulse width $W_{40}$ as the phase subtended when $\hat{\omega}$ varies from $40^{\circ}$, to the minimum, and back to $40^{\circ}$. Since the orientation of the spin axis affects the effective latitude on the radiation cone with which the line of sight intersects, the pulse width exhibits a clear variation with the orientation of the spin axis. The pulse width can be considered as a proxy for the observed pulse duration, and so precession of the spin axis over the course of the pulsar’s orbit will directly influence the time over which a pulse is viewable. Whilst the variations in the pulse width is due to how close the line of sight is to the edge of the emission cone, the asymmetry evident in the broadening is a result of the relativistic aberration. Consequently, due to spin effects the pulse profile will vary over the course of the orbit, and also vary between different kinds of orbit due to the difference in the 4-velocity and 4-spin vectors.

The asymmetry in the $\hat{\omega}-t$ profiles can be quantified by considering the phase angle $\chi$ at which the minimum of $\hat{\omega}$ occurs, which corresponds to the centroid of the arc subtended by the radiation beam across our line of sight (denoted by vertical dashed lines in Fig 9). Variation in $S_{\theta}$ in conjunction with relativistic aberration causes a shift in the centroid location $\Delta\chi$. This centroid shift in turn causes the pulse arrival times from a pulsar at a given location to vary with $S_{\theta}$ (Fig. 10), due to the extra time elapsed traversing the angle $\Delta\chi$. For a pulsar with a 1 ms spin period, the centroid of the profile from a pulsar with $S_{\theta}=\pi/12$ is delayed on the order of $\sim 10\mu\rm s$ compared to a pulsar with $S_{\theta}=\pi/4$ and typical delays are of the order $\sim(1-10)\mu$s, i.e. $0.1-1$ % of a MSP spin period. Since the centroid shift will occur smoothly over orbital timescales, during observations this effect should be traceable and measurable and used to inform the evolution of the pulsar spin axis. Furthermore, at steeper $S_{\theta}$ angles, the amplitude of $\hat{\omega}$ is also reduced since the radiation beam is shifted further from the observer’s line of sight. Approximating the pulsar beam as a Lambertian surface, precession of the spin axis will cause variations in the pulse intensity $dI$ as,

[TABLE]

Consequently, the precession of the spin axis can directly affect the observed intensity. Clearly this is an oversimple approximation and to correctly determine the pulse intensity also accounting for scattering, absorption and Doppler boosting would require covariant general relativistic radiative transfer (e.g. Younsi et al., 2012) along the geodesics determined via the framework presented in this paper. However, the Lambertian approximation serves as a adequate first order estimate to inspect the influence of spin precession on the pulse intensity. The peak intensity of a beam that originates from a pulsar with $S_{\theta}=\pi/4$ is $\lesssim 85\%$ that of a pulsar with $S_{\theta}=\pi/12$ for the observer at the same location (Fig. 10).

In addition to the orbital dynamics, spin precession can then also affect the pulse width (and hence duration), the pulse ToA due to the relativistic centroid shift and the observed pulse intensity. The precision that can be achieved with pulsar timing is a function of both the pulse intensity and the sharpness of the pulse profile, and so an accurate description of these quantities is key for realistic assessment of the prospects of the detection of a PSR-EMRB. Variations in the pulse profile width may provide a method to determine the PSR spin behaviour (e.g. Rafikov & Lai, 2006), and hence the orbital parameters of the system. For very compact orbits or periastron passages, the timescale of spin orientation variation can severely shorten, causing rapid changes in the pulse profile (width, intensity, ToA). Corrections for this change in the overall pulse profile shape will be necessary for both long term timing and - if the timescale of spin axis oscillation is comparable to the observation integration time - initial detection. This analysis has involved just varying $S_{\theta}$. Changes in $S_{\phi}$ over the course of the orbit adds an additional degree of freedom to influence the pulse profile, which are naturally included within our framework.

5.3 Gravitational Light Bending

The notion of a geodesic, i.e. a straight line in curved spacetime, means that the spatial path of the ray emitted from the pulsar can be bent due to the spacetime curvature induced by the central massive black hole. That is, the signal from a pulsar may not follow a purely radial path. A strongly deflected pulsar beam may propagate directly though the strong-field regime, considered in the analytical case for a Schwarzchild black hole in Wang et al. (2009a, b). Gravitational lensing by Sgr A* has been considered for emission for S-stars close to the Galactic centre (Bozza & Mancini, 2009; Bin-Nun, 2010; Bozza & Mancini, 2012) as observed by the next generation of Very Large Telescope Interferometer instrumentation (e.g. GRAVITY, Gravity Collaboration et al., 2017). The magnitude of the astrometric shift that results from gravitational bending is naturally dependent on the configuration of the system (i.e. relative alignment of source, lens, observer. See e.g. Bozza & Mancini, 2012, for a description of the relevant angles and geometry). For the S-stars, at alignment angles of $\sim 35^{\circ}$, the expected astrometric shift is of order $30\mu$as, whilst the contribution of post-Newtonian and spin corrections to gravitational bending are expected to appear at $\sim 5\mu$as, below current interferometric instrumental sensitivity (Bozza & Mancini, 2012).

Gravitational bending in a Kerr spacetime is naturally included in the ray tracing set up described in Sec. 2. In the weak deflection case, gravitational bending directly influences the photon ToA since the spatial trajectory of the ray is modified from flat Minkowski spacetime. Furthermore, in addition to the primary pulses received by the observer that suffer weak or no gravitational bending, there also exist secondary pulses that are associated with the primary pulses, but which are strongly bent (see e.g. Fig 11 and also Wang et al., 2009a). Whilst higher order pulses are possible, for the purposes of this paper we only consider the primary and secondary pulses, since the intensity of higher order pulses is naturally low. Our framework and algorithm as described in Sec. 4 is able to find both the primary and secondary pulses. Naturally these strong bending effects are most prominent when the pulsar is on the far side of the black hole. For certain orbital configurations, the pulsar ‘primary’ pulses may be invisible, and instead the pulsar is only visible via the strongly bent ‘secondary’ pulses. Whether strong bending of pulsar rays will occur astrophysically is highly dependent on both the system orbital configuration and the observer viewing angle, with the maximal bending occurring when the pulsar lies close to the central black hole and on the far side with respect to the observer. Analysis from Stovall et al. (2012) suggests that with current radio observation facilities, the probability of detecting strong bent pulsar beams is small but non-negligible. Regarding future radio facilities, optimistic estimates suggest that SKA should be able to detected strongly bent beams from multiple pulsars. The analysis of Stovall et al. (2012) is focused towards the Galactic centre where there are complications due to scattering which decreases the detection probability. Observations in globular clusters of the Local Group might provide more favourable hunting grounds. Furthermore, in estimating a probability Stovall et al. (2012) take as a prototypical model a pulsar at $r=10^{4}r_{\rm g}$. Taking this as the semi-major axis, this gives a Keplerian period of $\sim 4$ years. At such radii, as noted in Stovall et al. (2012), the degree of strong deflection is rather small and if an observer receives the primary beam it will also receive the secondary beams. However, for pulsars with smaller orbital periods, or at small orbital radii (e.g. at periastron) the probability of observing a strongly deflected beam would be increased. Extending the considered pulsar population to include not just pulsars of the Milky Way, but also nearby globular clusters would also increase this rate. Moreover, the large baselines offered by future radio interferometers like SKA will enable $\sim\mu$as astrometry (Fomalont & Reid, 2004; Smits et al., 2011). This provides a complementary pathway to timing for identifying secondary rays.

Both strong and weak gravitational bending have important implications for not only the photon ToA due to differing spatial paths, but also due to gravitational and relativistic time dilation, relativistic energy shift (see Sec. 5.4) and time-frequency behaviour due to temporal and spatial dispersion (Sec. 5.5).

5.4 Gravitational and relativistic time dilation

For a pulsar in a strong-field environment, the interplay of both gravitational time dilation (clocks run slower in potential wells) and the relativistic Doppler shift induced by the motion of the pulsar will have an impact on the observed radio signal; an apparent modulation in the observed pulsar period where the intrinsic pulsar rotation frequency is different to that recorded by some distant observer. This effect can be quantified as,

[TABLE]

where $\nu$ can be either the apparent pulsar rotation frequency or the photon frequency. This quantity is frame-invariant and we choose to evaluate it in the global coordinate frame. Since $\gamma$ is a function of both the gravitational time dilation and the relativistic motion of the pulsar, the net shift varies over the pulsar orbit (Figs 12, 13). Due to the gravitational bending of light rays, for some sections of the orbit the observer will also receive a secondary beam. When the pulsar is directly behind the black hole the primary ray cannot reach the observer. Gravitational light bending allows for secondary rays to reach the observer from sections of the orbit that would otherwise exist in the black hole shadow. Since this bent ray follows a different geodesic to the primary ray, it is subject to a different general relativistic time dilation and so will also exhibit a different value of $\gamma$.

The modulation of the observed pulsar rotation frequency must be considered on both long and short timescales. For the former, the variation in the observed rotation frequency will cause a pulse smearing, akin to that observed in relativistic binaries (e.g. Jouteux et al., 2002), but with the added contribution due to strong gravitational time dilation. This would be especially pronounced in the detection and timing of pulsars close to periastron and may be further complicated by the existence of secondary rays which in addition to having different light travel times, will have different values for $\gamma$. Over longer timescales the modulation in the observed pulsar frequency will need to be considered for any accurate timing model, especially one which is then used to test fundamental aspects of GR.

Taking into account these relativistic effects on the apparent pulsar rotation period is essential to accurately model the signal from a pulsar in an extreme gravity environment. The dependence of the net time dilation on the motion of the emitter also emphasizes the necessity to have an accurate description of the pulsar orbital dynamics - accounting for all general relativistic effects e.g. spin-curvature coupling - in modeling the pulsar signal.

5.5 Time-frequency signal

Typically, pulsar signals are dispersed in time due to the interstellar medium (ISM) and any signal needs correcting to account for this dispersion (e.g. Cordes et al., 2016). The time delay induced by the presence of plasma along the line of sight,

[TABLE]

for frequency $\nu$ and dispersion measure,

[TABLE]

where $n_{e}$ is the free electron number density. During pulsar searches, attempts are made to de-disperse the signal at a number of trial DMs, since the true DM is unknown, before searching for periodicities in Fourier space. Once an accurate DM is determined, long-term precision timing measurements can then take place. Dispersion can also be induced due to relativistic fluids close to the black hole (e.g. accretion flow, Psaltis, 2012). Dispersive effects due to the interaction with plasma along the ray path are fully accounted for within our framework, as all electromagnetic effects are included within the $\omega_{\rm p}^{2}$ term of the Hamiltonian. As discussed, due to strong gravitational bending we receive two radio signals from a PSR; the primary minimally bent ray and the ray which undergoes strong bending due to the curvature of the space time (illustrated in Fig 11). The secondary ray is typically retarded with respect to the first ray, since it follows a longer spatial path. The alternative trajectory ensures that the secondary ray is subject to a different gravitational time dilation and also the DM along the two rays will be different. Higher order time variations also occur since the pulsar must rotate by a certain angle between the primary and secondary emission, during which time the pulsar is also traveling, but these variations are a fraction of the pulsar period and dominated by the difference in the ray path. All these factors conspire to affect the pulse arrival time and the degree of dispersion across the frequency bandwidth within that pulse (Fig 14).

Gravitational light bending means that the apparent position of the pulsar - the position which lies on the asymptote of the tangent line to the ray that converges at the observer - is distinct from the true pulsar position. This introduces complexities in timing observations since two pulses which arrive at a similar time could originate from different times of the orbit and so be are subject to different time dilations, energy shifts, spin and orbital effects and have different dispersion measures and pulse profiles. Since the DM and $t-\nu$ profile will be different for primary and secondary rays this may provide a mechanism to distinguish between the two in any timing solution.

In addition to temporal dispersion, the combination of strong-field curvature and electron plasma induces a spatial dispersion such that the trajectory followed by the ray after being deflected by the black hole is frequency dependent (Kimpson et al., 2019). This has several implications for rays which are gravitationally bent. Firstly, some rays may no longer be visible in specific frequency bins since if the spatial dispersion is sufficiently severe the ray path is bent such that is does not hit the observer’s image plane. Moreover, each ray which does reach the observer has followed a different spatial trajectory and so suffers from a time delay due to the alternative ray path, in addition to different DM’s along that path and encounters a different spacetime. Any corrections for gravitational light bending are frequency dependent; the apparent position of the pulsar can be related to multiple true positions, depending on the ray frequency, each necessitating a different DM correction to be applied in each frequency bin. The total received signal is then not some function that varies smoothly with the pulsar orbital phase, but instead the convolution of different energy rays emitted at different orbital phases and consequently subject to differing relativistic and line-of-sight effects. This may, depending on the orbital configuration, result in additional difficulties in detecting signals from PSR-EMRB systems.

6 Discussion and Conclusions

We have presented the principles and framework for calculating the radio signal from a PSR in an EMRB. We restrict our study to the extreme mass ratio of EMRB systems and so do not consider PSRs in stellar-mass black hole binaries with finite mass ratios (e.g. Blanchet, 2014; Liu et al., 2014). We account for both relativistic and astrophysical effects and the convolution between the two. This includes gravitational and relativistic time dilation and energy shift, gravitational light bending, complex orbital dynamics induced by spin couplings, temporal variation and distortion of the pulse profile due to spin axis precession and relativistic aberration, 2nd order pulses due to gravitational bending, and dispersions (temporal and spatial) induced by the material along the line of sight. We have demonstrated that within our framework we are able to determine the time-frequency behaviour accounting for all these effects. The framework also applies for any orbital configuration, e.g. we are not restricted to orbital motion in the equatorial plane or beaming confined to the orbital plane. The methods used are entirely covariant and general relativistic, rather than working under any post-Newtonian approximation and so are inherently more accurate. Indeed, the post-Newtonian method is an explicitly weak-field method, and the validity of its application to strong-field dynamical regimes is unclear (Will, 2011). Whilst working explicitly in the Kerr metric means that we are unable to independently probe either alternative gravitational theories or extensions to Kerr (e.g. Kerr spacetime with an arbitrary mass quadrupole, see Bini et al., 2009), our framework provides the basis for a theoretical timing model which can then be compared with observations for tests of strong field GR. We approximate the PSR body as a perfect sphere. However due to the spin of the PSR the true shape is more oblate. This will ultimately influence the pitch angle of the ray with the neutron star surface. This effect is considered to be minor, but the method could easily be extended to account for this oblateness (see Nättilä & Pihajoki, 2018). We neglect the effects of hydrodynamic drag due to the plasma that surrounds that black hole, since at compact radii ( $\lesssim 10^{4}r_{g}$) the gravitational and relativistic effects dominate (Psaltis, 2012). We also do not take account of any potential Newtonian perturbations on the motion of the pulsar (e.g. Merritt et al., 2011) due to the presence of other masses (e.g. stars, other compact objects etc.) since these factors are likely negligible for the orbital periods considered in this work ($\lesssim 0.3$ years, Liu et al., 2012). Indeed, the potential for external perturbations to hamper tests of strong-field GR necessitates that an ideal PSR-EMRB systems should have orbital periods on the order of $0.1$ years (or better), or else observations should be taken close to periapsis (see discussion in Psaltis et al., 2016). These are precisely the regions where the spacetime curvature and orbital acceleration is greatest, further stressing the importance of a strong-field timing model. We also neglect any influence of gravitational radiation on the orbit or the ray trajectory. The neglect of gravitational radiation is justified since in the extreme mass ratio limit, the timescale for orbital decay due to gravitational wave emission is (Misner et al., 1973),

[TABLE]

where $M$ is the mass of the black hole, $m$ the pulsar mass and $r$ the orbital separation. The eccentricity function is,

[TABLE]

If we take the PSR orbital period $P$ to be Keplerian, then for a pulsar with mass $1.4M_{\odot}$ on an eccentric ($e=0.8$), $P=0.1$ year orbit around a BH with mass $4.3\times 10^{6}M_{\odot}$ ,

[TABLE]

and so the effects of gravitational radiation can be neglected. Even for smaller radii and more eccentric orbits the spacetime is well approximated as stationary (e.g. $\tau_{GW}/P\ \sim 10^{5}$ for $e=0.9$, $r=100M$ ). Whilst the effects of gravitational radiation are then not important for a single orbit, for observations over longer periods of time the effect of gravitational emission on the orbit and hence the timing solution will need to be considered. The PSR may also emit a gravitational wave burst during passage through periastron (Berry & Gair, 2013a, b). The influence of this gravitational radiation on both the PSR trajectory and the photon ToA is highly non-trivial and not considered here.

In this work we do not address the task of how to use our calculation to perform mock data analysis and extract orbital parameters from simulated PSR-EMRB timing data. This is achievable, for instance, by using a software package such as TEMPO (Hobbs et al., 2006) to determine the parameters of the timing model from the simulated data (e.g. Liu et al., 2014), or/and by performing a Markov Chain Monte Carlo fitting to investigate constraints on the orbital parameters (e.g. Zhang & Saha, 2017). In principle, the method outlined in this work can be used to generate time-frequency data which could then be analysed by the aforementioned methods, but this is beyond the scope of this paper. We also caveat that, due to the high stellar density of the Galactic centre, timing data can be influenced by external Newtonian perturbations (from e.g. stars, stellar mass compact objects Merritt et al., 2011). Any consistent timing solution should therefore provide a method to correct for these gravitational foreground disturbances (e.g. Angélil & Saha, 2014; Zhang & Saha, 2017). The handling of such perturbations is not explicitly considered in this work, but we restrict our analysis to pulsars on orbits of $P\lesssim 0.1$ years where such perturbations are more likely to be negligible (Liu et al., 2012). The analysis of timing data taken from orbital systems with longer periods would need a method to remove these effects (another possibility is to analyse only data taken close to periapsis, where the magnitude of these perturbations is expected to be less severe).

Whilst we have started the theoretical basis for timing observations of a strong-field PSR, there are a number of further potential developments of this work. With the ray tracing solution, we can then perform general relativistic radiative transfer along the rays (e.g. Fuerst & Wu, 2004, 2007; Younsi et al., 2012) so as to determine the effects of line-of-sight material on the beam intensity. Understanding the received intensity and consequent S/N ratio is essential for accurately exploring the prospect of using PSR as probes of strong-field GR. As mentioned, it would also be of interest to investigate how well the weak-field post-Keplerian parametrization can describe strong-field effects. This is important both for determining the types of PSR we require to test GR and for creating an accurate model to then compare with observations. A coherent $t-\nu$ model which accounts for all relativistic and line of sight effects could then also be potentially used to inform detections; e.g. are acceleration searches necessary for the detection of MSP close to the Galactic centre? As noted by Faucher-Giguère & Loeb (2011), the high stellar density in the Galactic Centre may allow for the creation of some rare binaries (e.g. triple systems). The subsequent dynamics (e.g. Remmen & Wu, 2013) and impact on the PSR signal would be another interesting pursuit.

To summarize, accurately modelling the time-frequency behaviour from a radio PSR in the strong-field regime leads to a number of higher-order effects which will influence the photon ToA. These include:

• •

Consideration of spin curvature coupling can lead to variations of order $\pm 10\mu$s in the photon ToA compared to when spin-curvature coupling is neglected. Lower order estimates to the orbital motion (e.g. pure Keplerian dynamics) will further exacerbate the discrepancy, whilst accurate determination of the orbital motion (i.e. including spin-spin, spin-orbit, spin-curvature couplings) is essential for accurately modeling the frequency modulation.

• •

Precession of the spin axis in conjunction with relativistic aberration influences both the pulse profile, pulse duration, the pulse arrival time and the pulse intensity. Aperiodicity in the spin precession may introduce additional complications in the detection of pulsars, whilst severe precession could leave the pulsar signal intermittent.

• •

Gravitational bending causes deviation from a simple Minkowski geodesic which naturally influences the photon ToA. Strong gravitational bending can cause multiple (primary/secondary) pulses to be received by the observer, emitted from a pulsar at approximately the same location. Each ray follows a distinct spacetime path and so suffers different transfer effects, both gravitational (e.g. time dilation) and those due to interaction with material along the line of sight (i.e. temporal dispersion with the DM different for the primary and secondary rays). The interplay of gravitational bending with spatial dispersion further compounds these effects.

• •

Gravitational and relativistic time dilation causes substantial shift in the observed pulse frequency over the orbit. The magnitude of the dilation varies over the orbit and is more pronounced for more eccentric orbits. Secondary, highly bent rays traverse a different curvature of spacetime and so the frequency modulation is different from that of the primary rays.

• •

The presence of material along the ray path causes a temporal dispersion in the photon ToA, whilst a coupling with the strong spacetime curvature causes a spatial dispersion of the rays. Since each ray of a given frequency follows a different path - and so is subject to different gravitational and line of sight effects - the ToA will vary.

To conclude, in order to undertake precision tests of GR it is necessary to have a coherent, accurate theoretical model with which to compare observations. In this work we present a framework for calculating the time-frequency behaviour from a pulsar signal in an entirely general relativistic context, including the effects of spin axis precession and nutation, relativistic aberration, relativistic and gravitational time shift, relativistic energy shifts, spatial and temporal dispersion and gravitational light bending. The convolution of these effects will all ultimately influence the photon arrival time. Such a framework is the first step toward the creation of an accurate timing model of a PSR signal in the strong-field regime.

Acknowledgments

We thank the reviewer Prasenjit Saha for useful critique and suggestion which helped to improve the paper, and in particular his recommendations towards making explicit what the timing model does and does not include. We also thank Ellis Owen and Ziri Younsi for careful reading of the manuscript and useful helpful suggestions towards its improvement. TK acknowledges support from a UK Science and Technology Facilities Council studentship. This research has made use of NASA’s Astrophysics Data Systems.

Appendix A Ray Tracing Initial Conditions

The differential equations are integrated ‘backwards-in-time’ from the observer image plane to the black hole, using a fifth-order Runge-Kutta-Fehlberg algorithm with adaptive step-size (see Press et al., 1996). The centre of the observer’s image plane is defined at some location $r_{\rm obs},\theta_{\rm obs}$, where $r_{\rm obs}$ is the distance from the black hole center and $\theta_{\rm obs}$ the angle from the positive black hole $z$-axis. Since the Kerr metric is axisymmetric we can set $\phi_{\rm obs}=0$. $r_{\rm obs}$ is chosen so as to be sufficiency large such that the observer’s grid can be considered as a Euclidean grid with zero spacetime curvature, and all rays are perpendicularly incident on this grid. The coordinates $\alpha,\beta$ denote the location of a point on this grid location, where for $\theta_{\rm obs}=\pi/2$, $\beta$ is parallel to the positive black hole $z$-axis. The black hole coordinate system is right-handed with axes $\mathbf{x}=(x,y,z)^{\rm T}$. The observer’s coordinate system is left-handed with axes $\mathbf{x^{\prime}}=(x^{\prime},y^{\prime},z^{\prime})^{\rm T}$ and the $z^{\prime}$ axis oriented towards the centre of the black hole. Since we are interested in the plane coordinates $\alpha,\beta$ we can say $x^{\prime}=\alpha,y^{\prime}=\beta,z=0$. We determine the initial conditions of rays starting on the observer’s grid. To do this we transform from $\mathbf{x^{\prime}}\rightarrow\mathbf{x}$ via the methods outlined in Younsi (2013); Pu et al. (2016):

Rotate clockwise by $(\pi-\theta_{\rm obs})$ about the $x^{\prime}$-axis ($R_{x^{\prime}}$) 2. 2.

Rotate clockwise by $(2\pi-\phi_{\rm obs})$ about the $z^{\prime}$-axis ($R_{z^{\prime}}$). 3. 3.

Reflect in the plane $y^{\prime}=x^{\prime}$ ($A_{y^{\prime}=x^{\prime}}$). 4. 4.

Translate $\mathbf{\bar{x^{\prime}}}$ so that the origins of both coordinate systems coincide ($T_{x^{\prime}\rightarrow x}$)

The net transformation is then

[TABLE]

where $\mathcal{D}=(\sqrt{r_{\rm obs}^{2}+a^{2}}-z^{\prime})\sin\theta_{\rm obs}-y^{\prime}\,\cos\theta_{\rm obs}$. We then transform from Cartesian to Boyer-Lindquist coordinates,

[TABLE]

where $w=x^{2}+y^{2}+z^{2}-a^{2}$. This defines the initial $(r,\theta,\phi)$ for a photon on the observer grid.

The initial velocities of the ray can then be determined. Since each ray arrives perpendicular to the image plane, $(\dot{x^{\prime}},\dot{y^{\prime}},\dot{z^{\prime}})=(0,0,1)$. Consequently, the velocity components in the black hole frame are given by

[TABLE]

Converting to Boyer-Lindquist coordinates gives expressions for $(\dot{r},\dot{\theta},\dot{\phi})$ in the black hole frame:

[TABLE]

where $\mathcal{R}=\sqrt{r^{2}+a^{2}}$ and $\Phi=\phi-\phi_{\rm obs}$. This completely defines our initial conditions.

Appendix B MPD Spin-Interaction

B.1 Momentum-Velocity relation in the Tulczyjew-Dixon condition

With the TD condition the pulsar mass is given by,

[TABLE]

The spin vector is given by,

[TABLE]

We can express the spin tensor in terms of the spin-vector,

[TABLE]

and so it follows that,

[TABLE]

which is also a constant of the motion. Contracting Eq. 18 with $u_{\nu}$ gives an expression for the momentum as,

[TABLE]

B.2 Initialisation of MPD orbits

We can initialize the $p^{\alpha},s^{\alpha},x^{\alpha}$ 4-vectors of our particle, which are related to the conserved quantities of the Kerr spacetime, the mass $m$ of the pulsar, its energy $E$, angular momentum $L_{z}$ and Carter constant $Q$. More usefully, the conserved parameters $E,L_{z},Q$ can be mapped to geometrical orbital parameters $p,e,z_{-}$ with the transformation described in Barausse et al. (2007) where $p$ is the semi-latus rectum, $e$ the eccentricity, $z_{-}=\cos^{2}\theta_{\rm min}$ and $\theta_{\rm min}$ is the minimum angle reached by the pulsar. This framework is fundamentally an approximation since it does not include spin effects and so both the eccentricity and semi-latus rectum are not constant, as for a Keplerian orbit, but evolve with times (see e.g. Singh et al., 2014). Nevertheless these variations are typically small they provide a decent first order approximation to the sorts of orbits that we want to model.