# Classical and general relativistic post-Keplerian effects in binary pulsars hosting fast rotating main sequence stars

Lorenzo Iorio, Michel Rieutord, Jean-Pierre Rozelot, Armando Domiciano de Souza

<figure><img src="image_1.png" /></figure>

## To cite this version:

Lorenzo Iorio, Michel Rieutord, Jean-Pierre Rozelot, Armando Domiciano de Souza. Classical and general relativistic post-Keplerian effects in binary pulsars hosting fast rotating main sequence stars. European Physical Journal C: Particles and Fields, 2019, 79 (8), pp.690. 10.1140/epjc/s10052-019-7194-5. hal-02147941

HAL Id: hal-02147941

https://hal.science/hal-02147941v1

Submitted on 28 May 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Classical and general relativistic post-Keplerian effects in binary pulsars hosting fast rotating main sequence stars

Lorenzo Iorio<fnref n="1" />,ª, Michel Rieutord<fnref n="2" />,ᵇ, Jean-Pierre Rozelot<fnref n="3" />,ᶜ, Armando Domiciano de Souza<fnref n="4" />,⁵

Received: 8 July 2019 / Accepted: 1 August 2019 / Published online: 16 August 2019
© The Author(s) 2019

**Abstract** We consider a binary system composed of a pulsar and a massive, fast rotating, highly distorted main sequence star of mass *M*, spin angular momentum *S*, dimensionless mass quadrupole moment *J*₂, equatorial and polar radii *R*ₑ, *R*ₚ, flattening *v* = (Rₑ − Rₚ)/Rₑ, and ellipticity *ε* = √(1 − Rₚ<fnref n="2" />/Rₑ<fnref n="2" />) as a potential scenario to dynamically put to the test certain post-Keplerian effects of both Newtonian and post-Newtonian nature. We numerically produce time series of the perturbations Δ(*δτ*) of the Rømer-like, orbital component of the pulsar's time delay *δτ* induced over 10 years by the pN gravitoelectric mass monopole (Schwarzschild, *GMc*⁻<fnref n="2" />), quadrupole (*GMR*₂*J*₂*c*⁻<fnref n="2" />), gravitomagnetic spin dipole (Lense–Thirring, *GSc*⁻<fnref n="2" />) and octupole (*GSR*₂*ε*<fnref n="2" />*c*⁻<fnref n="2" />) accelerations along with the Newtonian quadrupolar (*GMR*₂*J*₂) one. We do not deal with the various propagation time delays due to the travelling electromagnetic waves. It turns out that, for a Be-type star with *M* = 15 M☉, *R*ₑ = 5.96 R☉, *v* = 0.203, *S* = 3.41 × 10⁴⁵ J s, *J*₂ = 1.92 × 10⁻<fnref n="3" /> orbited by a pulsar with an orbital period *P*ₚ ≈ 40–70 days, the classical oblateness-driven effects are at the ≲ 4 − 150 s level, while the pN shifts are of the order of ≲ 1.5 − 20 s (*GMc*⁻<fnref n="2" />), ≲ 10 − 40 ms (*GMR*₂*J*₂*c*⁻<fnref n="2" />), ≲ 0.5 − 6 ms (*GSc*⁻<fnref n="2" />), ≲ 5 − 20 µs (*GSR*₂*ε*<fnref n="2" />*c*⁻<fnref n="2" />), depending on their orbital configuration. The root-mean-square (rms) timing residuals στ of almost all the existing non-recycled, non-millissecond pulsars orbiting massive, fast rotating main sequence stars are ≲ ms. Thus, such kind of binaries have the potential to become interesting laboratories to measure, or, at least, constrain, some Newtonian and post-Newtonian (*GMc*⁻<fnref n="2" />, *GMJ*₂*c*⁻<fnref n="2" />,

and, perhaps, *GSc*⁻<fnref n="2" /> as well) key features of the distorted gravitational fields of the fast rotating stars hosted by them.

## 1 Introduction

In its weak-field and slow-motion approximation, general relativity predicts that, in addition to the time-honored post-Newtonian (pN) gravitoelectric and gravitomagnetic precessions induced by the mass *M* (Schwarzschild) and the spin angular momentum *S* (Lense–Thirring) of the central body acting as source of the gravitational field, other pN gravitoelectric and gravitomagnetic orbital effects related to its oblateness arise as well [9,18,19,24,46,48,57,59–61,67]. So far, they have never been put to the test in any astronomical and astrophysical scenarios, despite some recent preliminary investigations pertaining the planet Jupiter in our solar system [23,25]; for some embryonic thoughts about an Earth-spacecraft scenario, see Iorio [23,24]. To the pN level, the oblateness of astronomical bodies modifies also the propagation of the electromagnetic waves in their deformed spacetime. About the perspectives of measuring the resulting deflection due to Jupiter with astrometric techniques, see, e.g., Crosta and Mignard [11], Kopeikin and Makarov [32], Le Poncin-Lafitte and Teyssandier [37], Abbas et al. [1], and references therein.

An analysis of the analytical expressions of the pN gravitoelectric and gravitomagnetic orbital precessions due to the asphericity of the primary [24,25] shows that the key ingredients needed to enhance their magnitude are a strongly distorted central body, and a highly eccentric and close orbit of the moving particle.

Binaries composed by a pulsar and a main-sequence star [66] may offer, in principle, interesting natural laboratories to try to investigate such little known pN effects. Indeed,

<footnote n="1">¹ Ministero dell'Istruzione, dell'Università e della Ricerca (M.I.U.R.)-Istruzione, Viale Unità di Italia 68, 70125 Bari, BA, Italy</footnote>

<footnote n="2">² IRAP, Université de Toulouse, CNRS, UPS, CNES, 14, avenue Edouard Belin, 31400 Toulouse, France</footnote>

<footnote n="3">³ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Nice & 77, Chemin des Basses Moulières, 06130 Grasse, France</footnote>

<footnote n="4">⁴ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, UMR7293 Lagrange, 28 Av. Valrose, 06108 Nice Cedex 2, France</footnote>

<footnote>ª e-mail: lorenzo.iorio@libero.it</footnote>

<footnote>ᵇ e-mail: Michel.Rieutord@irap.omp.eu</footnote>

<footnote>ᶜ e-mail: Jean-Pierre.ROZELOT@univ-cotedazur.fr</footnote>

<footnote>ᵈ e-mail: Armando.Domiciano@oca.eu</footnote>
<table><caption>**Table 1** Binary pulsars with a massive, fast rotating main sequence star companion discovered so far. For each of them, we list the companion, the distance (when available), the mass *M*, the equatorial and polar radii *R*<sub>e</sub>, *R*<sub>p</sub> (when available), the orbital period *P*<sub>b</sub>, the orbital eccentricity *e*, the spin period *P*, the rms timing residuals *σ*<sub>τ</sub>, the data analysis time</caption><thead><tr><th>Pulsar</th><th>Companion</th><th>Distance (kpc)</th><th><i>M</i> (M<sub>☉</sub>)</th><th><i>R</i><sub>e</sub> (R<sub>☉</sub>)</th><th><i>R</i><sub>p</sub> (R<sub>☉</sub>)</th><th><i>P</i><sub>b</sub> (d)</th><th><i>e</i></th><th><i>P</i> (ms)</th><th><i>σ</i><sub>τ</sub> (ms)</th><th>Δ<i>T</i> (yr)</th></tr></thead><tbody><tr><td>PSR J0045-7319</td><td>B1V star</td><td>In SMC</td><td>8.8 ± 1.8</td><td>6.4 ± 0.7</td><td>—</td><td>51.17</td><td>0.808</td><td>930</td><td>7.4</td><td>2</td></tr><tr><td>PSR J1740-3052</td><td>Main sequence star</td><td>—</td><td>> 11</td><td>—</td><td>—</td><td>231</td><td>0.578</td><td>570</td><td>0.8</td><td>2.29</td></tr><tr><td>PSR B1259-63</td><td>B2e star LS 2883</td><td>2.75</td><td>≈ 30</td><td>≈ 9.7</td><td>≈ 8.1</td><td>1237</td><td>0.870</td><td>48</td><td>0.46</td><td>13</td></tr><tr><td>PSR J1638-4725</td><td>Be star</td><td>—</td><td>> 4</td><td>—</td><td>—</td><td>1940</td><td>0.95</td><td>764</td><td>5.3</td><td>4.35</td></tr><tr><td>PSR J2032 + 4127</td><td>Be star MT91 213</td><td>1.7</td><td>≈ 15</td><td>—</td><td>—</td><td>8578</td><td>0.93</td><td>143</td><td>0.5–1</td><td>6</td></tr></tbody></table>

they are systems composed of a neutron star regularly emitting electromagnetic radio pulses orbiting an usually more massive main sequence star, which, in most cases, is highly distorted due to its fast rotation, along a generally elliptical orbit. A very accurate observable quantity in binary pulsars is represented by the measurement of the times of arrival (TOAs) τ of the pulses emitted by the neutron star which, in case it has a gravitationally bound companion, exhibit a regular variation δτ due to, among other things, the Keplerian motion about the common centre of mass: it is the Rømer-like time delay. The full variation of the pulses' times of arrival is due to several other effects connected, e.g., with the propagation of the electromagnetic waves through the deformed spacetime of the system [66].

The first binary pulsar hosting a main sequence star to be discovered was PSR B1259-63 [27,58]; it is characterized by a highly eccentric orbit (e = 0.870) with an orbital period of $P_b = 1237$ days = 3.38 yr. The pulsar's non-degenerate companion is the fast spinning Be star LS 2883, whose equatorial velocity $V_e$ is about 280 km s⁻<fnref n="1" /> corresponding to ~70% of its break-up velocity [52]. Its mass *M* and equatorial radius *R*<sub>e</sub> amounts to about 30 Solar masses (M<sub>☉</sub>) and 9.7 Solar radii (R<sub>☉</sub>) [47]. Later, the eccentric (e = 0.808) binary PSR J0045-7319 was discovered [29]. To date, it is the fastest orbiting system since it is $P_b = 51.17$ days. Its primary is a main sequence B-star spinning close to its break-up velocity [28,36]. PSR J1638-4725, having an orbital period of $P_b = 1940$ days = 5.3 yr and e = 0.95, was found by Lorimer et al. [39]. Its stellar companion should be a rapidly rotating Be star. PSR J1740-3052 [63], with $P_b = 231$ days and e = 0.578, hosts most likely a B-type main-sequence star [43,64]. The most recently discovered main-sequence-star binary pulsar is the highly eccentric (e = 0.93) PSR J2032+4127 [42] characterized by $P_b = 8578$ days = 23.5 yr. The companion of the neutron star is the massive Be star MT91 213 with $M \approx 15 M_{\odot}$.

For the sake of completeness, we mention also a few other binary pulsars hosting a non-degenerate star, although they are not relevant for our purposes in view of the nature of

span Δ*T*. For the values of the listed parameters, see the references cited in the text and the online pulsar catalog at http://www.atnf.csiro.au/research/pulsar/psrcat/. For PSR B1256-63 we give the equatorial *R*<sub>e</sub> and polar *R*<sub>p</sub> radii as estimated by Negueruela et al. [47]

their non massive and fast-rotating partners. They are PSR J1903+0327 [4], whose companion is a F5V-GOV ~ 1 M<sub>☉</sub> star moving in $P_b = 95$ days along a rather eccentric orbit with e = 0.44, the transitional millisecond pulsar PSR J1023 + 0038 [3] orbiting a low-mass (0.2 M<sub>☉</sub>) companion star in a circular path with $P_b = 4.75$ h. The rms timing residuals of the aforementioned binary pulsars are all of the order of ~0.5 ms, apart from PSR J1903 + 0327 which is at the ~1 µs level; more specifically, they are ~0.46 ms over 13 year for PSR B1259-63 [65], 7.4 ms over 2 year for PSR J0045-7319 [29], ~5.3 ms over 4.35 year for PSR J1638-4725 [39], ~0.8 ms over 2.29 year for PSR J1740-3052 [63], ~0.5-1 ms over about 6 year for PSR J2032 + 4127 [42], ~1 µs over about 3 yr for PSR J1903 + 0327 [4], 0.1 ms over about 4 yr for PSR J1023 + 0038 [2]. Table 1 summarizes the key data for the binary pulsars hosting massive, fast rotating main sequence stars.

About the achievable accuracy level in timing residuals, in the case of the pulsars orbiting a main sequence star, their timing seems doomed to stay at the ~ms level. It is so because, for evolutionary reasons, they are not fully recycled [62]. Thus, their spinning periods *P* are not at the millisecond level, and their TOAs are not measured with a precision of the order of ~µs, and their timing is often contaminated by timing noise [21]. The timing of the non-recycled pulsars is almost always less accurate than for the fully recycled pulsars. Indeed, PSR J1903+0327 rotates with a period *P* = 2.5 ms, and its rms timing residuals are as little as ~1 µs. The spin period of PSR J1023+0038 is *P* = 1.7 ms, and its rms timing residuals are 0.1 ms. Incidentally, we mention the fact that, according to Table 2 of Arzoumanian et al. [4], the rms timing accuracy of some fully recycled pulsars with *P* = 1.5–10 ms, isolated or with a white dwarf as companion, is of the order of 0.1–0.2 s. It is expected that future instrumental improvements may push the rms timing accuracy of some of them to the<fnref n="1" /> ~10 ns level over time spans some yr long. For binary pulsars hosting another neutron star, the rms timing accuracy is of the order of 1–100 µs.

<footnote n="1">¹ A. Possenti, personal communication to L.I., April 2019.</footnote>
Here, we will preliminarily investigate the size and the temporal patterns of the perturbations $\Delta (\delta\tau)$ induced on the Rømer-like orbital part of the pulsar's time delay $\delta\tau$ by both the standard (Schwarzschild and Lense-Thirring) and the oblateness-driven pN accelerations felt by a fictitious neutron star orbiting a highly distorted, fast rotating B-type main sequence star in view of a possible detection in new binaries that may eventually be discovered in the future. However, caution is in order before inferring too optimistic conclusions from a straightforward comparison of our simulated time series with the rms timing residuals listed in Table 1. Even if the size of some pK signatures were to be larger than the ≃ ms level, it does not necessarily mean that such effects will be measurable in the actual processing of the real observations. Indeed, careful, dedicated simulated data reductions and covariance analyses should be performed by explicitly modeling the signals of interest, estimating its characteristic parameters and inspecting the resulting correlations with the other parameters usually estimated. It should be kept in mind that, in principle, an unmodeled effect may be removed from the post-fit residuals, at least to a certain extent, being “absorbed” in the estimated values of the other parameters determined in the data reduction. Thus, our investigation should be deemed just as a sensitivity analysis able to preliminarily explore the potential of the scenarios considered. We will assume the validity of general relativity throughout the paper, which is organized as follows.

In Sect. 2, we discuss the magnitude of the angular momentum $S$ and the dimensionless quadrupole mass moment $J_2$ of typical fast rotating massive B-type stars. Section 3 is devoted to the numerical calculation of the perturbations $\Delta (\delta\tau)$ induced on the pulsar's Rømer time delay by some post-Keplerian (pK) classical and pN accelerations. We do not calculate the propagation delays accounting, e.g., for the effects on the pulsar's travelling electromagnetic waves through the deformed spacetime of the B-type star. Section 4 summarizes our findings and contains our conclusions. Once again, we stress the preliminary nature of our sensitivity investigation; we do not perform a full covariance analysis implying, e.g., the simulation of the pulsar's TOAs and their reduction along with parameter estimation.

## 2 Quadrupole mass moment, flattening, and angular momentum of fast rotating main sequence B-type massive stars

The binary pulsars on which we shall focus presumably own a main sequence B-type massive star. Such stars are mostly fast rotators Levato and Grosso [38] and are hence distorted by the centrifugal acceleration. Such a distortion takes the mass distribution away from spherical symmetry and endows these stars with a quadrupolar and higher grav-

itational moments. In the case of Be stars (namely B stars with emission lines), rotation is believed to be almost critical, namely the rotational velocity is taken close (> 70%) to the Keplerian velocity at equator (e.g. [55]). When critical rotation is reached, the surface distortion is maximum and the flattening $v \doteq (R_e - R_p)/R_e$, expressed in terms of the equatorial and polar radii $R_e, R_p$, respectively, is close to one third. For such stars the computation of the gravitational quadrupolar moment cannot be done perturbatively as it is the case for the Sun, which is a slow rotator (e.g. [56]). Modelling these stars requires two-dimensional models. Fortunately, self-consistent 2D-models have recently been achieved with the ESTER code Espinosa Lara and Rieutord [16], Rieutord et al. [54]. Compared to previous 2D-models, ESTER models include self-consistently the baroclinicity of the stellar envelopes and can thus predict the associated differential rotation. They therefore provide unambiguously the total angular momentum of the star given, for instance, its equatorial velocity.

With the ESTER code, we computed the parameters of three stellar models of 10, 15, and 30 $M_\odot$ as they can represent the companions of PSR J0045-7319, PSR J2032+4127 and PSR B1256-63 respectively. Since the evolutionary status of the stars is unknown but presumably on or close to the main sequence, we computed their steady state at ZAMS (Zero-Age Main Sequence) and at half-main sequence to get an idea of the effects of evolution. We use standard galactic metallicity $Z = 0.02$ with the solar mixture (which may be approximate for PSR J0045-7319 which is in the SMC, known to be less metallic than the Galaxy). Results of ZAMS and evolved models are displayed in Tables 2 and 3 respectively. There we give the total spin angular momentum $S$ and the dimensionless quadrupole mass moment $J_2$ along with other bulk parameters of the models. We recall that multipole gravitational moments $J_\ell$ are defined by the multipole expansion of the gravitational potential of a mass $M$, namely

$$U(\mathbf{r}) = -\frac{\mu}{r} \left[ 1 - \sum_{\ell} J_{\ell} \left( \frac{R_e}{r} \right)^{\ell} \mathcal{P}_{\ell}(\xi) \right], \quad (1)$$

<table><caption>**Table 2** ESTER models for Zero-Age Main Sequence (ZAMS) stars with an equatorial angular velocity at 70% of the critical angular velocity. $R_e$ and $V_e$ are the equatorial radius and velocity, respectively, $S$ is the total spin angular momentum, $v$ is the flattening, and $J_2$ is the dimensionless mass quadrupole moment. The metallicity is $Z=0.02$ and the hydrogen mass fraction is $X = 0.7$. In our simulations, we use the parameters of the star with 15 $M_\odot$ listed here</caption><thead><tr><th>M (M<sub>☉</sub>)</th><th>R<sub>e</sub> (R<sub>☉</sub>)</th><th>V<sub>e</sub> (km s<sup>-1</sup>)</th><th>S (×10<sup>44</sup> J s)</th><th>v</th><th>J<sub>2</sub> (×10<sup>-3</sup>)</th></tr></thead><tbody><tr><td>10</td><td>4.74</td><td>444</td><td>15.3</td><td>0.201</td><td>1.63</td></tr><tr><td>15</td><td>5.96</td><td>485</td><td>34.1</td><td>0.203</td><td>1.92</td></tr><tr><td>30</td><td>8.89</td><td>562</td><td>125.</td><td>0.210</td><td>2.17</td></tr></tbody></table>
<table><caption>**Table 3** Same as in Table 2 but for ESTER models of stars at mid-main-sequence, namely when the hydrogen mass fraction in the convective core is half of the initial one</caption><thead><tr><th><i>M</i> (M<sub>☉</sub>)</th><th><i>R</i><sub>e</sub> (R<sub>☉</sub>)</th><th><i>V</i><sub>e</sub> (km s<sup>-1</sup>)</th><th><i>S</i> (×10<sup>44</sup> J s)</th><th><i>ν</i></th><th><i>J</i><sub>2</sub> (×10<sup>-3</sup>)</th></tr></thead><tbody><tr><td>10</td><td>6.93</td><td>367</td><td>13.3</td><td>0.203</td><td>0.816</td></tr><tr><td>15</td><td>9.07</td><td>393</td><td>28.0</td><td>0.207</td><td>0.788</td></tr><tr><td>30</td><td>14.8</td><td>435</td><td>870.</td><td>0.227</td><td>0.495</td></tr></tbody></table>

where $\mu \doteq GM$ is the star's gravitational parameter, $G$ is the Newtonian constant of gravitation, $\xi \doteq \hat{\mathbf{S}} \cdot \hat{\mathbf{r}}$ is the cosine of the angle $\theta$ between the directions of the body's spin axis and of an external point at $\mathbf{r}$, $\mathcal{P}_{\ell}(\xi)$ is the Legendre polynomial of degree $\ell$. We consider the mass distribution of the stellar models to be symmetric with respect to equator thus making the $J_{\ell}$ of odd order all vanish. The remaining $J_{2p}$ can be computed with the integral expression

$$J_{2p} = -\frac{1}{MR_e^{2p}} \int_V r^{2p} \mathcal{P}_{2p}(\xi) \rho(\mathbf{r}) d^3\mathbf{r} \quad (2)$$

where the integration is over the volume $V$ of the star. The same expression is given in, for instance, Helled et al. [20]. Here we are especially interested in $J_2$ and $S$, namely in

$$J_2 = -\frac{1}{MR_e^2} \int_V r^2 \mathcal{P}_2(\xi) \rho(\mathbf{r}) d^3\mathbf{r} \quad \text{and} \\ S = \int_V r^2(1-\xi^2) \Xi(\mathbf{r}) \rho(\mathbf{r}) d^3\mathbf{r} \quad (3)$$

which are directly computed from the ESTER models; $\Xi(\mathbf{r})$ is the local angular speed.

We choose an angular rotation rate of 70% of the actual critical (Keplerian) angular velocity of the star. Such a rotation rate is typical of the nearby fast rotating stars that have been measured by interferometry. Their flattening is typically ~0.2 (e.g. [14]), as our models.

From Tables 2 and 3, we clearly see that as evolution proceeds, namely as the hydrogen content of the core decreases, $J_2$ decreases as expected from the resulting contraction of the convective core. From the work of James [26] we can compute $J_2$ for a polytrope of index $n = 3$ with a similar flattening as the ESTER models. We find that $J_2 \simeq 2.1 \times 10^{-3}$ which is quite similar to the ZAMS models. For evolved models one should use polytopes with a higher polytropic index, as they are more centrally condensed, and typically $n = 3.43$ matches the ESTER evolved models as far as $J_2$ is concerned. This result may be useful for simulating the orbital evolution of binary pulsars since polytropic models are much easier to compute.

### 3 The perturbations of the Römer-type pulsar's time delay due to some pK Newtonian and pN accelerations

Here, we will assume a coordinate system centered in the binary's barycenter whose reference z-axis is directed along the line of sight from the binary to the observer, while the reference $\{x, y\}$ plane spans the plane of the sky.

#### 3.1 The pulsar as a structureless, pointlike particle

We will, first, consider the following pK accelerations experienced by a test particle moving with velocity $\mathbf{v}$ in the external field of an oblate body of mass $M$, equatorial and polar radii $R_e$, $R_p$, ellipticity $\epsilon \doteq \sqrt{1 - R_p^2/R_e^2}$, angular momentum $\mathbf{S}$ and dimensionless quadrupole mass moment $J_2$. In Sect. 3.2, we will discuss the limits of validity of the point-particle approximation.

To the Newtonian level, the external potential of the distorted star at the position $\mathbf{r}$ is, from Eq. (1),

$$U(\mathbf{r}) = U_0 + \Delta U_2 = -\frac{\mu}{r} \left[ 1 - \left(\frac{R_e}{r}\right)^2 J_2 \mathcal{P}_2(\xi) \right], \quad (4)$$

where $\mathcal{P}_2(\xi) = (3\xi^2 - 1)/2$ is the Legendre polynomial of degree 2. The classical acceleration due to $J_2$ is

$$\mathbf{A}^{\text{NJ}_2} = -\nabla \Delta U_{J_2} = \frac{3\mu J_2 R_e^2}{2r^4} \left[ (5\xi^2 - 1) \hat{\mathbf{r}} - 2\xi \hat{\mathbf{S}} \right]. \quad (5)$$

The 1pN gravitoelectric Schwarzschild-like acceleration affecting the motion of a test particle in the static field of a nonrotating, spherically symmetric body is [51]

$$\mathbf{A}^{\text{1pNM}} = \frac{\mu}{c^2 r^2} \left[ \left( \frac{4\mu}{r} - v^2 \right) \hat{\mathbf{r}} + 4v_r \mathbf{v} \right], \quad (6)$$

where $c$ is the speed of light in vacuum, and $v_r \doteq \mathbf{v} \cdot \hat{\mathbf{r}}$ is the radial velocity of the test particle. Equation (6) is responsible for the formerly anomalous perihelion precession of Mercury whose explanation was the first empirical confirmation of general relativity [15].

The 1pN gravitomagnetic Lense–Thirring acceleration in the stationary field due to the rotating primary is [51]

$$\mathbf{A}^{\text{1pNS}} = \frac{2GS}{c^2r^3} \left[ 3\xi \hat{\mathbf{r}} \times \mathbf{v} + \mathbf{v} \times \hat{\mathbf{S}} \right]. \quad (7)$$

The gravitomagnetic field of the Earth was unambiguously measured for the first time by the Gravity Probe B (GP-B) mission [17]. Tests of the Lense–Thirring orbital precessions with some terrestrial geodetic satellites are ongoing; see, e.g. Renzetti [53], and Lucchesi et al. [40] and references therein for comprehensive reviews.
The 1pN gravitoelectric acceleration felt by a test particle in the field of an oblate body is [9, 24, 60, 61, 67]

$$A^{1pNM} J_2 = \frac{\mu J_2 R_e^2}{c^2 r^4} \left\{ \frac{3}{2} \left[ (5\xi^2 - 1) \hat{\mathbf{r}} - 2\xi \hat{\mathbf{S}} \right] \left( v^2 - \frac{4\mu}{r} \right) - 6 \left[ (5\xi^2 - 1) v_r - 2\xi v_S \right] \mathbf{v} - \frac{2\mu}{r} (3\xi^2 - 1) \hat{\mathbf{r}} \right\}, \quad (8)$$

where $v_S \doteq \mathbf{v} \cdot \hat{\mathbf{S}}$ is the component of the particle's velocity along the direction of the primary's spin. Note that the parameter $J_2$ in Eq. (8) is the same entering Eq. (5), as per Equations 1 to 2 of Soffel et al. [60].

The 1pN gravitomagnetic acceleration imparted to a test particle by the spin octupole moment<fnref n="2" /> of a uniformly rotating homogenous oblate spheroid [25, 46, 48] can be cast into the compact form

$$A^{1pNS} J_2 = \frac{3GSR_e^2\epsilon^2}{7c^2r^5} \mathbf{v} \times \left\{ 5\xi(7\xi^2 - 3)\hat{\mathbf{r}} + 3(1 - 5\xi^2)\hat{\mathbf{S}} \right\}. \quad (9)$$

The pK accelerations of Eqs. (5)–(9) perturb the otherwise Keplerian motion of the binary causing a change $\Delta(\delta\tau)$ of the regular variation $\delta\tau$ of the TOAs due to the relative orbital motion of the pulsar and the massive companion. It can be modeled as the ratio of the projection of the barycentric orbit of the pulsar onto the line of sight to $c$ [12, 31]. Thus, $\Delta(\delta\tau)$ can be calculated by looking at the perturbations $\Lambda z$ induced by Eqs. (5)–(9) on the $z$-component of the pulsar's barycentric orbital motion. We do that by numerically integrating the equations of motion of a fictitious pulsar with mass $M_p = 1.4 M_\odot$ having as a companion a Be-type star with $M = 15 M_\odot$, $R_e = 5.96 R_\odot$, $v = 0.203$, $S = 3.41 \times 10^{45}$ J s, $J_2 = 1.92 \times 10^{-3}$, as per Table 2, for different values of its orbital configuration, determined by the initial values of the semimajor axis $a$, the eccentricity $e$, the orbital inclination $I$ to the plane of the sky, the longitude of the ascending node $\Omega$, the argument of periastron $\omega$, the true anomaly at epoch $f_0$, and of the stellar spin axis characterized by its inclination $i$ to the line of sight, and the longitude $\phi$ of the projection of the stellar spin onto the plane of the sky. For a chosen pK acceleration $A^{pK}$, in order to compute its perturbation $\Delta(\delta\tau) = \Delta z(t)/c$ over, say, 10 yrs, we perform two runs sharing the same initial conditions with and without $A^{pK}$, calculate the resulting time series of $z(t)/c$ and take their difference. Figures 1, 2, 3, 4 and 5 depict our results for a given orbital configuration and the aforementioned Be-type main sequence star; in the panels of each Figure, we vary the Keplerian orbital elements and the orientation of $\hat{\mathbf{S}}$ in order to investigate the sensitivity to the parameter space

of the adopted binary system. For the sake of a comparison, the Keplerian delay for the adopted reference orbital configuration lies in the range $-250$ s $\lesssim \delta\tau \lesssim 50$ s over one orbital revolution.

Figure 1 displays the Newtonian signatures due to the star's $J_2$. It can be noted that the resulting signal is strongly dependent on the initial value of the true anomaly, spanning the range from 50 s to $-150$ s. Instead, for the given value $f_0 = 228$ deg, the sensitivity to the other orbital parameters is rather modest, amounting to about 10 s.

The 1pN gravitoelectric Schwarzschild-like signatures due to the stellar mass monopole are reproduced in Fig. 2. Also in this case, the initial value of the true anomaly induces a marked variability in the decadal time series which ranges from to $-5$ to $25$ s. The other orbital parameters have less impact since the resulting variation of the signals is of the order of about 2 s.

The 1pN gravitoelectric time series induced by the quadrupole mass moment $J_2$ of the star are the subject of Fig. 3. In this case, the ranges of variation due to all the orbital parameters are rather similar, amounting to about 20–50 ms.

Figure 4 shows the 1pN gravitomagnetic Lense-Thirring signatures due to the spin dipole moment of the star. They are mainly sensitive to the orbital inclination $I$ and to the spin's inclination $i$ to the line of sight which induces a variability as large as 10–12 ms. The ranges of variation induced by the other parameters are, instead, of the order of 1–5 ms.

The 1pN gravitomagnetic time series caused by the spin octupole moment of the star are depicted in Fig. 5. The widest range of variability, 30–50 µs, is due to the inclination $I$, the node $\Omega$, and the longitude $\phi$ of the spin's projection onto the plane of the sky. We also checked that, for a pulsar orbiting in $\simeq 20–30$ days a star as massive as those in the last lines of Tables 2 and 3 with periastron distances of $r_{\min} \simeq 1.2–1.1 R_{\text{eq}}$, the magnitude of the timing signatures would reach the $\simeq 1–10$ ms level.

## 3.2 The effects of the mass and spin multipoles of the pulsar

Until now, we have considered the neutron star as a structureless, point particle moving around its more massive companion. In fact, a pulsar is an extended body with its own mass and spin multipole moments which, at least in principle, may have an impact on its orbital motion in a full two-body framework.

The modification of Eq. (6) for two finite bodies of masses $M_A, M_B$ is [61]

$$A^{1pNM} = \frac{\mu_{\text{tot}}}{c^2 r^2} \left\{ \left[ (4 + 2\zeta) \frac{\mu_{\text{tot}}}{r} - (1 + 3\zeta) v^2 + \frac{3}{2} \zeta v_r^2 \right] \hat{\mathbf{r}} + (4 - 2\zeta) v_r \mathbf{v} \right\}, \quad (10)$$

<footnote n="2">² It will be shown that its effects are small enough to justify order-of-magnitude calculations, without need of detailed stellar models.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Numerically integrated time series of the timing shift $\Delta (\delta\tau)$ due to Eq. (5), in s, for variations of the parameter space of a fictitious binary pulsar, characterized by $P_b = 50$ days, $e = 0.8$, $I = 50$ deg, $\Omega = 140$ deg, $\omega = 149$ deg, $f_0 = 228$ deg, orbiting a Be-type star with $M = 15 M_{\odot}$, $R_e = 5.96 R_{\odot}$, $J_2 = 1.92 \times 10^{-3}$, $i = 60$ deg, $\phi = 217$ deg</figcaption></figure>

where $\mu_{\text{tot}} = GM_{\text{tot}} = G(M_{\text{A}} + M_{\text{B}})$, and $\zeta \doteq M_{\text{A}}M_{\text{B}}/M_{\text{tot}}^2$. By taking the standard value $M_{\text{p}} = 1.4 M_{\odot}$ for the mass of the pulsar, it is $\zeta = 0.078$ for the Be-star assumed in Sect. 3.1. It turns out that the introduction of $\zeta$ in our numerical code changes the size of the time series of Fig. 2 by $\simeq 0.1-0.2$ s, while their temporal patterns remain essentially unchanged. Given the current level of accuracy in the timing

residuals, such a discrepancy might be significative, and Eq. (10) should be used instead of Eq. (6).

In regard to the angular momentum, the spin $\mathbf{S}$ of the Be-star in Eq. (7) should be replaced by the sum [5]

$$ \mathbf{S} \rightarrow \left( 1 + \frac{3}{4} \frac{M_{\text{p}}}{M} \right) \mathbf{S} + \left( 1 + \frac{3}{4} \frac{M}{M_{\text{p}}} \right) \mathbf{S}_{\text{p}}, \quad (11) $$
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Numerically integrated time series of the timing shift $\Delta(\delta\tau)$ due to Eq. (6), in s, for variations of the parameter space of a fictitious binary pulsar, characterized by $P_b = 50$ days, $e = 0.8$, $I = 50$ deg, $\Omega = 140$ deg, $\omega = 149$ deg, $f_0 = 228$ deg, orbiting a Be-type star with $M = 15 M_\odot$</figcaption></figure>

where $S_p$ is the angular momentum of the pulsar. By assuming, e.g., $S_p = 3 \times 10^{40}$ J s as for PSR J0737-3039A [10,30,33,34,41], it turns out that the magnitude of the second term in Eq. (11) amounts to about $\simeq 7 \times 10^{-5}$ of that of the first one. Moreover, it is $(3/4) (M_p/M) = 0.07$. Thus, as far as the 1pN gravitomagnetic Lense-Thirring effect is concerned, our scenario can be well approximated by a restricted two-body system with a spinning primary, and Eq. (7) is substantially adequate. Indeed, it turns out that rescaling the star's angular momentum in our numerical simulations as dictated by Eq. (11) slightly modifies the size of the time series of Fig. 4 by just $\simeq 1 - 1.5$ ms.

The effect of the pulsar's quadrupole mass moment $J_2^p$ can be accounted for by the replacement $M \rightarrow M_{\text{tot}}$ in Eq. (5) and by writing in its right-hand-side another term for $J_2^p$ analogous to that for the stellar oblateness $J_2$. As a result, by introducing the dimensional quadrupole mass moment $Q_2 \doteq -J_2 M R_e^2$, the two terms in the right-hand-side of the resulting modified version of Eq. (5) are weighed by [5]

$$Q_2 = \left(1 + \frac{M_p}{M}\right) Q_2, \qquad (12)$$

$$Q_2^p = \left(1 + \frac{M}{M_p}\right) Q_2^p. \qquad (13)$$

For a neutron star, it is [7,8,35,50]

$$Q_2^p = -q \frac{M_p^3 G^2}{c^4}, \qquad (14)$$

with $0.07 \lesssim q \lesssim 3.507$ for a variety of Equations of State (EOSs). Thus, we have

$$Q_2 = -9.8 \times 10^{47} \text{ kg m}^2, \qquad (15)$$

$$Q_2^p = -q \ 1.2 \times 10^{37} \text{ kg m}^2, \qquad (16)$$

so that $Q_2^p \simeq 10^{-10} Q_2$. On the other hand, $M_p/M = 0.09$ in Eq. (12). Thus, for the Newtonian signature of the quadrupole mass moment, the restricted two-body scenario with an oblate primary is adequate in the present case, and the use
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Numerically integrated time series of the timing shift $\Delta (\delta\tau)$ due to Eq. (8), in ms, for variations of the parameter space of a fictitious binary pulsar, characterized by $P_b = 50$ days, $e = 0.8$, $I = 50$ deg, $\Omega = 140$ deg, $\omega = 149$ deg, $f_0 = 228$ deg, orbiting a Be-type star with $M = 15 M_{\odot}$, $R_{\odot} = 5.96 R_{\odot}$, $J_2 = 1.92 \times 10^{-3}$, $i = 60$ deg, $\phi = 217$ deg</figcaption></figure>

of Eq. (5) is justified provided that the stellar quadrupole moment $Q_2$ is replaced by Eq. (12). Indeed, it turns out that the introduction of $Q_2$ in the numerical integration changes the size of the time series in Fig. 1 by $\simeq 1$ s which may not be neglected, given the current level in $\sigma_{\tau}$.

Despite Eqs. (8) and (9) were derived so far only for the motion of a test particle around a spinning, oblate mass, there

are no doubts that they are adequate to the scenario considered here.

Actually, general relativity predicts that, in general, there is also a self-force due to the spin angular momentum of an extended rotating body in motion in an external gravitational field which acts on it modifying its orbit through a spin-orbit coupling [6, 13, 22, 44, 45, 49]. In order to quickly evaluate the
<figure><img src="image_6.png" /><figcaption>**Fig. 4** Numerically integrated time series of the timing shift $\Delta(\delta\tau)$ due to Eq. (7), in ms, for variations of the parameter space of a fictitious binary pulsar, characterized by $P_b = 50$ days, $e = 0.8$, $I = 50$ deg, $\Omega = 140$ deg, $\omega = 149$ deg, $f_0 = 228$ deg, orbiting a Be-type star with $M = 15 M_{\odot}$, $S = 3.41 \times 10^{45}$ J s, $i = 60$ deg, $\phi = 217$ deg</figcaption></figure>

possible impact of such an effect in our case, let us proceed as follows. To the 1pN level, the precession $\dot{\Omega}_{S_p}$ of, say, the node $\Omega$, averaged over one orbital revolution of the spinning pulsar in its motion around the massive Be-type star, assumed nonrotating, is [22]

$$ \dot{\Omega}_{S_p} = \frac{3\mu\sigma_p \csc I (\hat{\sigma}_p \cdot \hat{\mathbf{m}})}{2c^2 a^3 (1 - e^2)^{3/2}}. \quad (17) $$

In this expression, $\sigma_p = S_p/M_p$ is the pulsar's spin angular momentum per unit mass, while $\hat{\mathbf{m}} = \{-\cos I \sin \Omega, \cos I \cos \Omega, \sin I\}$ is a unit vector in the orbital plane perpendicular to the line of the nodes. Let us compare Eq. (17) to the analogous precession induced by some of the effects previously considered in which the pulsar was treated as a point particle. The Lense-Thirring node rate, induced by Eq. (7) which is responsible for the $\simeq$ ms time series of Fig. 4, is [22]
<figure><img src="image_7.png" /><figcaption>**Fig. 5** Numerically integrated time series of the timing shift $\Delta(\delta\tau)$ due to Eq. (9), in µs, for variations of the parameter space of a fictitious binary pulsar, characterized by $P_b = 50$ days, $e = 0.8$, $I = 50$ deg, $\Omega = 140$ deg, $\omega = 149$ deg, $f_0 = 228$ deg, orbiting a Be-type star with $M = 15$ M$_\odot$, $R_e = 5.96$ R$_\odot$, $\nu = 0.203$, $S = 3.41 \times 10^{45}$ J s, $i = 60$ deg, $\phi = 217$ deg</figcaption></figure>

$$ \dot{\Omega}_{\text{LT}} = \frac{2GS \csc I (\hat{\mathbf{S}} \cdot \hat{\mathbf{m}})}{c^2 a^3 (1 - e^2)^{3/2}}. \quad (18) $$

It turns out that, in the case of a pulsar like, e.g., PSR J0737-3039A orbiting the Be-type star of the Figs. 1, 2, 3, 4 and 5, it is $|\dot{\Omega}_{S_p}/\dot{\Omega}_{LT}| \simeq 10^{-5}$. Thus, we conclude that the spin-orbit self-force experienced by the pulsar is completely negligible

in our scenario. In fact, there is also a further self-force acting on an extended rotating body moving in an external gravitational field due to the spin-spin coupling between the angular momenta of the source and of the orbiter itself [6,22]. By following the same strategy for the spin-orbit coupling, the results in Iorio [22] allow to conclude that, in our case, such an
effect is even smaller than the previous one, being of the order of $\simeq 10^{-7}$ of, say, the $\simeq$ ms-level Lense-Thirring signal.

## 4 Summary and conclusions

We preliminarily explored the possibility of putting to the test several pK features of motion of Newtonian and pN origin in binaries hosting a pulsar and a massive, fast rotating, highly distorted main sequence star characterized by mass *M*, angular momentum *S*, equatorial and polar radii $R_e$, $R_p$, flattening $v$, ellipticity $\epsilon$ and dimensionless quadrupole moment $J_2$. Indeed, in addition to the usual 1pN Schwarzschild and Lense-Thirring effects due to the mass monopole ($\propto GMc^{-2}$) and spin dipole ($\propto GSc^{-2}$) moments, respectively, of the distorted stellar field, there are also other 1pN orbital effects, induced by the mass quadrupole ($\propto GM R_c^2 J_2 c^{-2}$) and spin octupole ($\propto G S R_c^2 \epsilon^2 c^{-2}$) moments, whose magnitudes may, perhaps, lie above the sensitivity threshold of the pulsar timing residuals in yet-to-be-discovered close binaries. However, the Newtonian perturbations due to $J_2$ are larger than the pN ones.

In order to perform a preliminary sensitivity analysis, we numerically investigated the orbital shifts $\Delta(\delta\tau)$ induced over 10 years by all of such pK perturbations on the otherwise Keplerian Rømer-type delay $\delta\tau$ in the pulsar's TOAs for a Be-type main sequence star characterized by $M = 15 M_\odot$, $R_e = 5.96 R_\odot$, $v = 0.203$, $S = 3.41 \times 10^{45}$ J s, $J_2 = 1.92 \times 10^{-3}$ orbited by a pulsar with an orbital geometry compatible ($P_b \simeq 40-70$ days) with some of the tightest binaries of this kind out of those discovered so far. We also investigated the sensitivity of the pK timing shifts to the whole system's parameter space by varying both the orientation of the stellar spin axis and the orbital elements of the pulsar's orbit. It turns out that the magnitude of the Newtonian signature due to $J_2$ can be as large as $\lesssim 4-150$ s, while the pN gravitoelectric quadrupolar one is $\lesssim 10-40$ ms. The pN gravitoelectric (Schwarzschild-like) signal due to the stellar mass monopole can be as large as $\lesssim 1.5-20$ s level. The pN gravitomagnetic shifts due to the spin dipole (Lense-Thirring) and octupole moments, for the evaluation of whose size the knowledge of the stellar spin angular momentum *S* and ellipticity *ε* is crucial, are of the order of $\lesssim 0.5-6$ ms, and $\lesssim 5-20$ µs, respectively. The rms of the timing residuals of all the non-recycled, non-millisecond pulsars like those having a fast rotating main sequence star companion discovered so far are $\lesssim$ ms over 2-13 years. It seems difficult that they can be substantially improved in the future. This implies that, in principle, all the pK effects considered fall within the potential measurability domain, except the pN spin octupole which is $\simeq 1-2$ orders of magnitude weaker, at least for the orbital configurations and the star considered in this study. The pN Lense-Thirring signatures are just at the $\simeq$ ms level. The different tempo-

ral patterns characteristic of the signals investigated may be helpful in separating them. It turns out that the tiniest pN effect may reach the $\simeq 1-10$ ms level only for a very tight, eccentric binary ($P_b \simeq 20-30$ days, $r_{min} \simeq 1.1-1.2 R_e$) hosting a Be-star with $30 M_\odot$, $S = 125-870 \times 10^{44}$ J s.

Finally, we stress the preliminary nature of our sensitivity analysis. To this aim, we remark that we did not compute the other kinds of time delay connected, e.g., with the propagation of the electromagnetic waves in the deformed spacetime. Moreover, we did not perform a full covariance analysis implying a simulation of the pulsar's TOAs, their reduction, and parameter estimation.

**Acknowledgements** We are grateful to A. Possenti for useful information. MR and ADS also acknowledge the strong support of the French Agence Nationale de la Recherche (ANR), under Grant ESRR (ANR-16-CE31-0007-01).

**Data Availability Statement** This manuscript has no associated data or the data will not be deposited. [Authors' comment: This is a theoretical study and no experimental data has been listed.]

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Funded by SCOAP³.

## References

1. U. Abbas, B. Bucciarelli, M.G. Lattanzi, MNRAS 485, 1147 (2019)

2. A.M. Archibald, V.M. Kaspi, J.W.T. Hessels, B. Stappers, G. Janssen, A. Lyne (2013), arXiv:1311.5161

3. A.M. Archibald et al., Science 324, 1411 (2009)

4. Z. Arzoumanian et al., Astrophys. J. Suppl. 235, 37 (2018)

5. B.M. Barker, R.F. O'Connell, Phys. Rev. D 12, 329 (1975)

6. B.M. Barker, R.F. O'Connell, Gen. Relativ. Gravit. 11, 149 (1979)

7. M. Bauböck, E. Berti, D. Psaltis, F. Özel, Astrophys. J. 777, 68 (2013)

8. E. Berti, N. Stergioulas, MNRAS 350, 1416 (2004)

9. V.A. Brumberg, *Essential Relativistic Celestial Mechanics* (Adam Hilger, Bristol, 1991)

10. M. Burgay et al., Nature 426, 531 (2003)

11. M.T. Crosta, F. Mignard, Class. Quantum Gravity 23, 4853 (2006)

12. T. Damour, G. Schaefer, Phys. Rev. Lett. 66, 2549 (1991)

13. W.G. Dixon, Philos. Trans. R. Soc. Lond. A 277, 59 (1974)

14. S.A. de Domiciano et al., A&A 569, A10 (2014)

15. A. Einstein, (Sitzber. Preuss. Akad., 1915), p. 831

16. L.F. Espinosa, M. Rieutord, A&A 552, A35 (2013)

17. C.W.F. Everitt et al., Phys. Rev. Lett. 106, 221101 (2011)

18. F. Frutos-Alfaro, M. Soffel, R. Soc. Open Sci. 5, 180640 (2018)

19. J. Heimberger, M. Soffel, H. Ruder, Celest. Mech. Dyn. Astron. 47, 205 (1990)

20. R. Helled, J.D. Anderson, M. Podolak, G. Schubert, Astrophys. J. 726, 15 (2011)

21. G. Hobbs, A.G. Lyne, M. Kramer, MNRAS 402, 1027 (2010)

22. L. Iorio, Gen. Relativ. Gravit. 44, 719 (2012)

23. L. Iorio, Class. Quantum Gravity 30, 195011 (2013)

24. L. Iorio, Int. J. Mod. Phys. D 24, 1550067 (2015)
25. L. Iorio, MNRAS 484, 4811 (2019)

26. R.A. James, Astrophys. J. 140, 552 (1964)

27. S. Johnston, R.N. Manchester, A.G. Lyne, M. Bailes, V.M. Kaspi, G. Qiao, N. D'Amico, Astrophys. J. 387, L37 (1992)

28. V.M. Kaspi, M. Bailes, R.N. Manchester, B.W. Stappers, J.F. Bell, Nature 381, 584 (1996)

29. V.M. Kaspi, S. Johnston, J.F. Bell, R.N. Manchester, M. Bailes, M. Bessell, A.G. Lyne, N. D'Amico, Astrophys. J. 423, L43 (1994)

30. M. S. Kehl, N. Wex, M. Kramer, K. Liu, in *Proceedings of the MG14 Meeting on General Relativity The Fourteenth Marcel Grossmann Meeting*, eds. M. Bianchi, R. Jantzen, R. Ruffini (World Scientific, Singapore, 2017), pp. 1860–1865

31. M. Konacki, A.J. Maciejewski, A. Wolszczan, Astrophys. J. 544, 921 (2000)

32. S.M. Kopeikin, V.V. Makarov, Phys. Rev. D 75, 062002 (2007)

33. M. Kramer, in *The Twelfth Marcel Grossmann Meeting. Proceedings of the MG12 Meeting on General Relativity*, eds. T. Damour, R. Jantzen, R. Ruffini (World Scientific, Singapore, 2012), pp. 241–260

34. M. Kramer, in *IAU Symposium, Vol. 337, Pulsar Astrophysics the Next Fifty Years*, eds. P. Weltevrede, B.B.P. Perera, L.L. Preston, S. Sanidas, (Cambridge University Press, Cambridge 2018), pp. 128–133

35. W.G. Laarakkers, E. Poisson, Astrophys. J. 512, 282 (1999)

36. D. Lai, L. Bildsten, V.M. Kaspi, Astrophys. J. 452, 819 (1995)

37. C. Le Poncin-Lafitte, P. Teyssandier, Phys. Rev. D 77, 044029 (2008)

38. H. Levato, M. Grosso, PASP 125, 1191 (2013)

39. D.R. Lorimer et al., MNRAS 372, 777 (2006)

40. D.M. Lucchesi, L. Anselmo, M. Bassan, C. Pardini, R. Peron, G. Pucacco, M. Visco, Class. Quantum Gravity 32, 155012 (2015)

41. A.G. Lyne et al., Science 303, 1153 (2004)

42. A.G. Lyne, B.W. Stappers, M.J. Keith, P.S. Ray, M. Kerr, F. Camilo, T.J. Johnson, MNRAS 451, 581 (2015)

43. E.C. Madsen et al., MNRAS 425, 2378 (2012)

44. B. Mashhoon, D. Singh, Phys. Rev. D 74, 124006 (2006)

45. M. Mathisson, Gen. Relativ. Gravit. 42, 1011 (2010)

46. J. Meichsner, M.H. Soffel, Celest. Mech. Dyn. Astron. 123, 1 (2015)

47. I. Negueruela, M. Ribó, A. Herrero, J. Lorenzo, D. Khangulyan, F.A. Aharonian, Astrophys. J. Lett. 732, L11 (2011)

48. M. Panhans, M.H. Soffel, Class. Quantum Gravity 31, 245012 (2014)

49. A. Papapetrou, Philos. Trans. R. Soc. Lond. A 209, 248 (1951)

50. G. Pappas, T.A. Apostolatos, Phys. Rev. Lett. 108, 231104 (2012)

51. G. Petit, B. Luzum et al., IERS Tech. Note 36, 1 (2010)

52. J.M. Porter, MNRAS 280, L31 (1996)

53. G. Renzetti, Open Phys. 11, 531 (2013)

54. M. Rieutord, L.F. Espinosa, B. Putigny, J. Comput. Phys. 318, 277 (2016)

55. T. Rivinius, A.C. Carciofi, C. Martayan, A&A Rev. 21, 69 (2013)

56. J.P. Rozelot, S. Godier, S. Lefebvre, Solar Phys. 198, 223 (2001)

57. M. Schanner, M. Soffel, Celest. Mech. Dyn. Astron. 130, 40 (2018)

58. R.M. Shannon, S. Johnston, R.N. Manchester, MNRAS 437, 3255 (2014)

59. M. Soffel, F. Frutos, J. Geodesy 90, 1345 (2016)

60. M. Soffel, R. Wirrer, J. Schastok, H. Ruder, M. Schneider, Celest. Mech. Dyn. Astron. 42, 81 (1987)

61. M.H. Soffel, *Relativity in Astrometry, Celestial Mechanics and Geodesy* (Springer, Heidelberg, 1989)

62. G. Srinivasan, New Astron. Rev. 54, 93 (2010)

63. I.H. Stairs et al., MNRAS 325, 979 (2001)

64. C.R. Tam, I.H. Stairs, S. Wagner, M. Kramer, R.N. Manchester, A.G. Lyne, F. Camilo, N. D'Amico, MNRAS 406, 1848 (2010)

65. N. Wang, S. Johnston, R.N. Manchester, MNRAS 351, 599 (2004)

66. N. Wex, MNRAS 298, 67 (1998)

67. C.M. Will, Phys. Rev. D 89, 044043 (2014)