Gravity Tests with Radio Pulsars

Abstract

The discovery of the first binary pulsar in 1974 has opened up a completely new field of experimental gravity. In numerous important ways, pulsars have taken precision gravity tests quantitatively and qualitatively beyond the weak-field slow-motion regime of the Solar System. Apart from the first verification of the existence of gravitational waves, binary pulsars for the first time gave us the possibility to study the dynamics of strongly self-gravitating bodies with high precision. To date there are several radio pulsars known which can be utilized for precision tests of gravity. Depending on their orbital properties and the nature of their companion, these pulsars probe various different predictions of general relativity and its alternatives in the mildly relativistic strong-field regime. In many aspects, pulsar tests are complementary to other present and upcoming gravity experiments, like gravitational-wave observatories or the Event Horizon Telescope. This review gives an introduction to gravity tests with radio pulsars and its theoretical foundations, highlights some of the most important results, and gives a brief outlook into the future of this important field of experimental gravity.

1. Introduction

Already one week before Albert Einstein presented his final field equations of general relativity (GR) to the Prussian Academy of Science [1], he demonstrated that his new theory of gravity naturally explains the anomalous perihelion precession of Mercury [2]. After this—at least in hindsight—first experimental verification of GR, many more Solar System tests followed, confirming different predictions by GR with ever increasing precision [3]. However, there are many aspects of gravity that cannot be tested in the weak-field slow-motion environment of the Solar System, since they require higher velocities and/or strong gravitational fields. A well known example is the emission of gravitational waves (GW). The experimental verification of such radiative aspects of gravity requires large, highly accelerated masses (m), which are absent in the Solar System. Furthermore, aspects related to strong gravitational fields cannot be investigated in the Solar System, where gravitational potentials ($\mathsf{\Phi }$) and binding energies ($E_{\mathrm{B}}$) are small ($\left|\mathsf{\Phi }\right|/c^2\sim \left|E_{\mathrm{B}}\right|/m{c}^2\lesssim {10}^{-6}$, where c denotes the speed of light in vacuum). To test such gravity regimes, one has to go beyond the Solar System.

An important step in that direction was the discovery of radio pulsars, i.e., rotating neutron stars (NS) that emit coherent radio waves along their magnetic poles, by Jocelyn Bell Burnell and Antony Hewish in 1967 [4]. This discovery itself was already a mile stone for relativistic astrophysics. Once identified as NSs, it was clear that pulsars with their strong gravitational fields and enormous matter densities could only be fully understood on the basis of GR. As a comparison to the Solar System, for NSs one typically finds $\left|\mathsf{\Phi }\right|/c^2\sim \left|E_{\mathrm{B}}\right|/m{c}^2\sim 0.15$. The discovery of the first pulsar in a binary system by Russell Hulse and Joseph Taylor [5], then marked the beginning of precision gravity tests with radio pulsars. Since 1974, many more pulsars in binary systems have been discovered [6], several of them suitable for gravity tests. Binary pulsars provide unique gravity experiments that allow the exploration of various predictions by GR and its alternatives, depending on aspects like the nature of the companion and the size of the binary orbit. One of the most constraining tests, when it comes to alternative gravity theories, do come from a pulsar in a stellar triple system. All this will be discussed in details below.

In the next section we give a short introduction to pulsars and observational details relevant for this review. In Section 3 we review some of the theoretical foundations of current gravity experiments with pulsars, with a focus on GR and some general aspects of alternative theories of gravity. In Section 4 we present different tests of the gravitational interaction of strongly self-gravitating bodies based pulsars as members of double NS systems. This includes experiments that test the propagation of photons in the gravitational field of a NS. Section 5 covers the most important binary pulsar tests for the emission of dipolar GWs, a prediction of many alternatives to GR. In Section 6 we present pulsar tests of the universality of free fall, a key part of the strong equivalence principle. In Section 7 we summarize pulsar tests of various gravitational symmetries, some related to strong-field generalizations of different parameters of the parametrized post-Newtonian (PPN) formalism, others to a temporal and spatial constancy of the gravitational constant. In Section 8 we give an outlook into the future, with a particular focus on the prospects coming from new and future radio telescopes. A summary is given in Section 9.

On a final note, the field of gravity experiments with pulsars has grown considerably since its beginning, nearly half a century ago. Within the scope of this review, it is sheer impossible to give an overview that is even roughly complete. For that reason, it was necessary to be selective in the topics covered, a selection that is inevitably biased by personal preferences.

2. Pulsar Population and Pulsar Observations

To date, close to 3000 radio pulsars are known, ranging in their rotational periods from about 1.4 ms up to 23.5 s [6]. More than 10% of these “cosmic light-houses” are members of a binary system. As binary companion one finds various types of stars, like NSs, white dwarfs (WD), and main sequence stars, but also rather unexpected companions like planets [7]. There is also a pulsar known to be member of a stellar triple system1, which we will discuss in more details in Section 6.1. Most of the known pulsars populate the disc of our Galaxy or reside in Galactic globular clusters. A few pulsars are known in the large and small Magellanic clouds.

The population of radio pulsars is best presented in a diagram that gives the rotational period (P) on one axis and its derivative with respect to time ($\dot{P}$) on the other (see Figure 1). The slow-down in the rotational period is due to the loss of rotational energy due to magnetic dipole radiation, and for a given period P scales with the square of the surface magnetic field strength. Typical values for the surface magnetic fields for radio pulsars range from about ${10}^7$ to about ${10}^{14}$ Gauss. Fast rotating pulsars with small $\dot{P}$, the so-called millisecond pulsars (MSP; bottom left corner of Figure 1), appear to be particularly stable in their rotation. In terms of stability on long time-scales, some of them can compare with the best atomic clocks on Earth [8]. MSPs are the result of a previous, stable mass transfer from the companion, which leads to a “recycling” of an old pulsar, often spinning it up to millisecond periods [9]. For that reason, they are mostly found as members of a binary system (see Figure 1). As a result of the recycling process, fully-recycled binary pulsars ($P\lesssim 10$ ms) are usually found in very circular orbits.

The rotational stability of recycled pulsars makes them ideal tools for precision measurements, including experiments in gravitational physics. Most gravity tests with pulsars are based on the precise measurement of the time of arrival (TOA) of the radio signals at the telescope, while keeping track of the rotational phase of the pulsar. This technique is called pulsar timing.

2.1. Pulsar Timing

A TOA is generally not determined for an individual radio pulse of a pulsar. There are two reasons for that. First, in most cases individual pulses are too weak to be detected, even with the largest radio telescopes, and therefore it requires the coherent addition of a large number of pulses to obtain a good signal-to-noise ratio. Furthermore, individual pulses are rather erratic in their appearance and strongly fluctuate in shape and intensity. Only once a sufficient number of pulses (many hundreds or even thousands) have been added, one gets a stable pulse profile. This stable, high signal-to-noise profile is then cross-correlated with a standard profile (different for every pulsar) to obtain a TOA for a particular rotational phase of the pulsar (see e.g., [10] for details). For some of the MSPs, such an integrated pulse profile consists of several 100 thousand pulses, achieving a timing precision of well below 1 $\mathsf{\mu }$s [11].

To connect the TOA of a pulsar signal at the telescope with the rotational phase of the pulsar at the signal’s emission, one uses a so-called timing model that accounts for all the relevant contributions in the pulsar system, in the Solar System, and in the interstellar space which have an effect on the proper times of pulsar and observer and on the signal propagation. The most important contribution from interstellar space is the presence of free electrons which affect the propagation speed of the radio waves. This is a frequency dependent effect and can be determined and corrected for by having observations at (at least) two different radio frequencies. For further details we refer the reader to the standard literature of pulsar astronomy, like [10]. Once the influence of the Solar System on the TOA, i.e., the motion of the telescope in the Solar System, the time dilation at the observatory, and the signal propagation delays in the curved spacetime of the Solar System (Shapiro delay), have been corrected for, one has the Solar System barycenter (SSB) TOA ($t_{\mathrm{SSB}}$). For a binary pulsar, $t_{\mathrm{SSB}}$ is linked to the proper time $T_{\mathrm{p}}$ of the pulsar by an expression of the following form

$$t_{\mathrm{SSB}}-t_0=D^{-1}\times \left[T_{\mathrm{p}}+{\Delta }_{\mathrm{R}}\left(T_{\mathrm{p}}\right)+{\Delta }_{\mathrm{E}}\left(T_{\mathrm{p}}\right)+{\Delta }_{\mathrm{S}}\left(T_{\mathrm{p}}\right)+{\Delta }_{\mathrm{A}}\left(T_{\mathrm{p}}\right)\right].$$

(1)

The time $t_0$ is a (chosen) reference epoch. D is the Doppler factor related to the (generally unknown) radial velocity between the SSB and the pulsar system. Only temporal changes of D are of interest for gravity tests, and therefore $D\equiv 1$ can be chosen at a given epoch (see discussion in [12]). The expressions ${\Delta }_i$ correspond to different effects that influence the propagation time of the pulsar signal. Quite generally, ${\Delta }_{\mathrm{R}}$ models propagation delays related to changes in the pulsar position due to its orbital motion, ${\Delta }_{\mathrm{E}}$ describes the time dilation, ${\Delta }_{\mathrm{S}}$ encompasses effects that are related to the propagation of the pulsar signal through the curved spacetime of the binary system, and finally ${\Delta }_{\mathrm{A}}$ accounts for aberration effects, related to the fact that the pulsar acts as a moving “light-house” having to send signals in a specific direction, i.e., towards Earth. All these contributions depend on a number of so called Keplerian and post-Keplerian (PK) parameters, which are determined in a fit of the model to the TOAs. Keplerian parameters are parameters which are known from Newtonian celestial mechanics, like the orbital period $P_{\mathrm{b}}$ and eccentricity e, while PK parameters represent phenomenological descriptions of different relativistic contributions present in the timing signal. More specific details are given in the next section and can be found in [12,13].

The proper time of a pulsar is directly related to its rotational phase $\varphi$ via the rotational frequency of the pulsar, ${\nu }_{\mathrm{p}}$, and its derivatives [14]:

$$\frac{\varphi \left(T_{\mathrm{p}}\right)}{2\pi }=N_0+{\nu }_{\mathrm{p}}(T_{\mathrm{p}}-t_0)+\frac{1}{2}{\dot{\nu }}_{\mathrm{p}}{(T_{\mathrm{p}}-t_0)}^2+\cdots .$$

(2)

The quantity $N_0$ is an arbitrary constant. If model Equation (1) represents a complete description of all relevant contributions to SSB TOAs of a pulsar, then a fit of the Keplerian and PK parameters will result in a phase coherent solution that accounts for every single rotation of a pulsar via Equation (2). For some of the best pulsars this covers many years or even decades, and in particular for MSPs that corresponds to ∼${10}^{11}$ rotations. Covering such a long baseline with high timing precision, which in some cases is better than 100 ns for a single TOA, leads to an unprecedented precision in many of the Keplerian and PK parameters (see e.g., [15]). The phase coherent solution is therefore key to most of the precision gravity tests with pulsars.

To conclude this section, we would like to give some general statements about the actual measurability of the particularly relevant individual parameters in pulsar timing. Specific examples are discussed in some of the following sections. The spin parameters of Equation (2) are generally known to many significant digits, up to around 15 digits for ${\nu }_{\mathrm{p}}$ in some cases. The pulsar’s location and proper motion at the sky, i.e., its astrometry, are usually also well known, in particular for pulsars with high timing precision. Concerning the orbital parameters of pulsars in binary systems, it all depends on the orbital configuration which of the Keplerian and PK parameters can be measured, and what precision one obtains for them. For any orbit, the orbital period $P_{\mathrm{b}}$ and the projected semi-major axis of the pulsar orbit x (exact definition given in Section 3.2) are usually measured with very high precision. For (sufficiently) eccentric orbits, in addition the orbital eccentricity e, the longitude of the periastron $\omega$ (angle between ascending node and periastron), and a reference time of periastron passage $T_0$ are known with high precision. In terms of relativistic effects, in eccentric orbits it is generally the relativistic precession of periastron $\dot{\omega }$ that is detected first. Once the periastron has precessed for a sufficient amount, the amplitude of the time dilation effect ${\Delta }_{\mathrm{E}}$ becomes measurable. For systems with $P_{\mathrm{b}}\lesssim 1\mathrm{d}$, at some stage also the change in the orbital period due to the emission of gravitational waves ${\dot{P}}_{\mathrm{b}}$ is observable, although that parameter is often “contaminated” by external effects, which we will discuss in details further below. If a binary system is sufficiently close to the edge-on orientation, then the Shapiro delay ${\Delta }_{\mathrm{S}}$ results in a prominent feature in the timing data that cannot be absorbed in other timing parameters. The Shapiro delay gives then direct access to the mass of the companion and the inclination angle i of the orbit. For a pulsar in a (nearly) circular orbit the location of periastron is not well measured and therefore neither the advance of periastron nor ${\Delta }_{\mathrm{E}}$ can be measured in such systems. Many examples for the measurement precision achieved in state of the art pulsar timing can be found in [11].

2.2. Pulse-Structure Analysis

The precision measurement of the arrival time of pulsar signals at the radio telescope is not the only kind of pulsar observation that has been used in pulsar gravity experiments. Tests related to relativistic spin precession require the determination of the orientation of a pulsar and how it changes over time. Basically there are two features that can be utilized for this, the shape of the (integrated) pulse profile and the polarization of the radio signal. Both of them are expected to change in a characteristic way if precession changes the angle between the spin of the pulsar and the line of sight (LOS) [16,17,18,19,20].

While on the one hand it is generally difficult to convert pulse profile changes (shape and intensity) into quantitative measurements, since the geometry of the pulsar emission in many cases does not follow a simple intrinsic geometry, the polarization, on the other hand, is often a more promising tool to obtain quantitative information about the pulsar orientation. In terms of polarization, many pulsars show a characteristic “swing” in their polarization as the LOS cuts through the magnetosphere of the pulsar. In the simplest cases the linear polarization angle in the sky $\psi$ can be sufficiently well described by the rotating vector model [21]:

$$tan\left(\psi \left(\varphi \right)-{\psi }_0\right)=\frac{sin\alpha sin(\varphi -{\varphi }_0)}{cos\alpha sin\zeta -sin\alpha cos\zeta cos(\varphi -{\varphi }_0)},$$

(3)

with $\alpha$ being the angle between the pulsar spin and the (observed) magnetic pole, and $\zeta$ denoting the angle between the spin and the direction towards the observer. A precession of the pulsar will generally lead to a long-term change in ${\psi }_0$ and $\zeta$, which encode the information on the change in the spin orientation.

Periodic variations of ${\psi }_0$ and $\zeta$ can also result from the aberration caused by the orbital motion of the pulsar. They depend on the details of the orbit and the orientation of the pulsar with respect to the observer [13]. This in turn leads to orbital variations of the pulse profile and the polarization pattern, the latter is evident from Equation (3). Such short-term periodic changes can be used to constrain the geometry of the pulsar, in particular when combined with the long-term secular changes.

In [13] the phenomenological pulse-structure analysis has been worked out in great details, by introducing a set of 11 PK parameters that in principle can be extracted from pulse-structure data. Long-term changes in ${\psi }_0$ and how this can be used in tests of spin precession has been discussed in details in [22].

3. Theoretical Foundations

Pulsar timing experiments can be viewed as ranging experiments which are based on the comparison of two very precise clocks. On the receiver side there is the atomic clock at the radio telescope, measuring the proper time of the observer ${\tau }_{\mathrm{obs}}$. On the emitter side there is the “pulsar clock” whose proper time ${\tau }_{\mathrm{psr}}$ can be determined form the rotational phase of the pulsars, via Equation (2), where $T_{\mathrm{p}}\propto {\tau }_{\mathrm{psr}}$. Any constant factor between $T_{\mathrm{p}}$ and ${\tau }_{\mathrm{psr}}$ gets absorbed when fitting for the frequency ${\nu }_{\mathrm{p}}$ and its derivatives. Figure 2 illustrates a spacetime view of a binary pulsar experiment.

In the previous section we already gave a general discussion of the timing model that connects ${\tau }_{\mathrm{psr}}$ and ${\tau }_{\mathrm{obs}}$. In this section we will have a closer look at the theoretical aspects that enter Equation (1), as they are at the heart of gravity tests with pulsars. In terms of details, we will restrict our discussion to effects that already have been quantitatively tested with pulsars, and therefore are needed for the interpretation of the observational results in the following sections. In the picture laid out in Figure 2, the theoretical treatment of the timing model Equation (1) can be split into two parts, the pulsar’s world line, including the proper time of the pulsar, and the world line of the photon propagating from the pulsar to the observer. In other words, one has a combination of the orbital motion and the signal propagation. We first discuss these two aspects within GR, and then address the most relevant changes one expects from alternative gravity theories.

3.1. Orbital Dynamics

The orbital dynamics of a binary pulsar is described on the basis of the post-Newtonian (PN) approximation scheme, where in GR the problem is greatly simplified due to the effacement of the internal structure of the strongly self-gravitating NS (see [23] and references therein). At the first PN approximation the orbital dynamics of the relative motion in harmonic coordinates is given by the Lagrangian

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\mathrm{O}}& =& m_{\mathrm{p}}\left(\frac{{\overrightarrow{v}}_{\mathrm{p}}^2}{2}+\frac{{\overrightarrow{v}}_{\mathrm{p}}^4}{8c^2}\right)+m_{\mathrm{c}}\left(\frac{{\overrightarrow{v}}_{\mathrm{c}}^2}{2}+\frac{{\overrightarrow{v}}_{\mathrm{c}}^4}{8c^2}\right)\hfill \\ & & +\frac{G{m}_{\mathrm{p}}{m}_{\mathrm{c}}}{r}\left[1-\frac{(\overrightarrow{r}·{\overrightarrow{v}}_{\mathrm{p}})(\overrightarrow{r}·{\overrightarrow{v}}_{\mathrm{c}})}{2c^2{r}^2}-\frac{{\overrightarrow{v}}_{\mathrm{p}}·{\overrightarrow{v}}_{\mathrm{c}}}{2c^2}+\frac{3{({\overrightarrow{v}}_{\mathrm{p}}-{\overrightarrow{v}}_{\mathrm{c}})}^2}{2c^2}\right]-\frac{G^{2}M{m}_{\mathrm{p}}{m}_{\mathrm{c}}}{2c^2{r}^2},\hfill \end{array}$$

(4)

where $m_{\mathrm{p}}$ and $m_{\mathrm{c}}$ are the (inertial) masses of pulsar and companion respectively, $M\equiv m_{\mathrm{p}}+m_{\mathrm{c}}$, $\overrightarrow{r}$ denotes their (coordinate) separation vector with $r\equiv |\overrightarrow{r}|$, and ${\overrightarrow{v}}_{\mathrm{p}}$ and ${\overrightarrow{v}}_{\mathrm{c}}$ the respective (coordinate) velocities. The coordinate system assumed in Equation (4) is harmonic at lowest order [24]. At the first PN level energy and angular momentum are conserved, a consequence of Equation (4). A particularly elegant solution for this dynamics has been found by [24], where the orbital motion can be described by a simple quasi-Keplerian parameterization. In terms of structure, the quasi-Keplerian paramerization is similar to the Keplerian solution in Newtonian celestial mechanics, with two differences at the 1PN level. First, instead of one eccentricity there are three different eccentricities in the three equations (Kepler’s equation and the equations for polar coordinates r and $\varphi$). More importantly, there is a factor in the equation for $\varphi$, that describes the relativistic precession of periastron. That factor is one of the PK parameters and in GR is given by

$$k=\frac{3{\beta }_{\mathrm{O}}^2}{1-e^2},$$

(5)

where ${\beta }_{\mathrm{O}}\equiv {\left(GM{n}_{\mathrm{b}}\right)}^{1/3}/c$, with $n_b\equiv 2\pi /P_{\mathrm{b}}$. Concerning the eccentricity e in Equation (5), at the 1PN level there is no need to distinguish between the three different eccentricities. The secular change of the longitude of periastron $\omega$ is given by $\dot{\omega }=k{n}_{\mathrm{b}}$, while it is important to note that in the Damour & Deruelle (DD) solution, $\omega$ does not advance linearly in time but linearly in the true anomaly (angle between periastron and pulsar position) [24].

The mutual difference between the three eccentricities gives rise to two PK parameters. However, while one is (practically) unobservable due to covariances with other timing parameters, the other one, ${\delta }_{\theta }$, is extremely difficult to measure, even for highly eccentric relativistic binary pulsars. So far, ${\delta }_{\theta }$ has not been measured with high significance in any binary pulsar. It is important to note, that the actual observable eccentricity in the DD timing model, as it is implemented in the pulsar standard analysis software (e.g., [25]), is a combination of the three eccentricities [12], which we will simply denote by e in this review.

Although in GR the GW damping contributions to the orbital dynamics appear at the 2.5PN level ($\mathcal{O}\left({\beta }_{\mathrm{O}}^5\right)$) in the equations of motion, by now there are quite a few binary pulsars where this contribution needs to be accounted for. GWs carry away energy and angular momentum from the binary orbit [26,27]. In the DD timing this is incorporated as adiabatic changes of the Keplerian parameters, first and foremost a secular change in the orbital period, which for GR reads [27]

$${\dot{P}}_{\mathrm{b}}=-\frac{192\pi }{5}\frac{1+73e^2/24+37e^4/96}{{(1-e^2)}^{7/2}}\frac{m_{\mathrm{p}}{m}_{\mathrm{c}}}{M^2}{\beta }_{\mathrm{O}}^5.$$

(6)

This is the result of the well-known quadrupole formula [28], which also holds for binary systems containing strongly self-gravitating bodies (see [23] and references therein). In spite of its high PN level ($\mathcal{O}\left({\beta }_{\mathrm{O}}^5\right)$), ${\dot{P}}_{\mathrm{b}}$ resulting from GW damping is observable in several relativistic binary pulsars, as it leads to a secular evolution in the orbital phase that grows quadratic in time [12,13,14].

The orbital motion gives rise to a further PK effect, that is closely linked to the proper time of the pulsar. The transformation between $T_{\mathrm{p}}$ and the coordinate time t of the binary system is the so-called Einstein delay, denoted by ${\Delta }_{\mathrm{E}}$ in Equation (1). To leading order it accounts for a changing second-order Doppler effect due to the variation in the pulsar speed $|{\overrightarrow{v}}_{\mathrm{p}}|$, and a changing gravitational redshift in the gravitational field of the companion due to the orbital variation in r. It is only relevant for eccentric orbits, as constant factors drop out in pulsar timing. The periodic difference between $T_{\mathrm{p}}$ and t is (with sufficient accuracy) described by a single amplitude, a PK parameter which in GR is given by [14]

$${\gamma }_{\mathrm{E}}=\frac{P_{\mathrm{b}}}{2\pi }e\left(1+\frac{m_{\mathrm{c}}}{M}\right)\frac{m_{\mathrm{c}}}{M}{\beta }_{\mathrm{O}}^2.$$

(7)

Although the effect becomes stronger for wider orbits ($\propto P_{\mathrm{b}}$), it is generally only measurable in relativistic binary pulsars with short orbital periods. The reason for this is, that the Einstein delay is apriori covariant with the over all Rømer delay (${\Delta }_{\mathrm{R}}$ in Equation (1); see Section 3.2 for details), and can only be separated once the longitude of periastron $\omega$ has sufficiently advanced, which requires a sufficiently large $\dot{\omega }$ [12,14].

3.2. Signal Propagation

The strongest influence on the signal propagation time, related to the binary system, comes simply from the orbital motion of the pulsar. This is called the Rømer delay and denoted by ${\Delta }_{\mathrm{R}}$ in Equation (1). The strength of that effect is given by the projected semi-major axis of the pulsar orbit, defined as $x\equiv a_{\mathrm{p}}sini/c$, where $a_{\mathrm{p}}$ is the semi-major axis of the pulsar orbit and i is the orbital inclination, i.e., the angle between the orbital angular momentum and the direction from the observer to the pulsar. For relativistic binary pulsars, x is typically in the order of one to a few seconds.

The calculation for the Rømer delay is solely based on the orbital dynamics described in Section 3.1, and assumes a flat spacetime for the signal propagation. For binary systems close to edge-on, i.e., i close to 90${}^{°}$, however, the gravitational field of the companion cannot be ignored. While the radio signal propagates through the curved spacetime of the companion it is subject to the so-called Shapiro delay [29], a well tested effect in the Solar System [30]. In binary pulsars, to leading order the Shapiro delay depends on two PK parameters, the range $r_{\mathrm{S}}$ and the shape $s_{\mathrm{S}}$, and is given by

$${\Delta }_{\mathrm{S}}=2r_{\mathrm{S}}ln\left(\frac{1+ecosf}{1-s_{\mathrm{S}}sin(\omega +f)}\right),$$

(8)

where f is the true anomaly of the orbital motion [14]. In GR one finds

$$r_{\mathrm{S}}=\frac{G{m}_{\mathrm{c}}}{c^3}\mathrm{and}{s}_{\mathrm{S}}=sini.$$

(9)

For the Shapiro range one finds $r_{\mathrm{S}}\simeq 4.9255(m_{\mathrm{c}}/M_{\odot })\mathsf{\mu }$s, which for high inclinations, where $sini$ approaches 1, gets strongly enhanced around superior conjunction ($\omega +f={90}^{°}$).

Finally, since the time of emission of a pulsar signal is linked to the pulsar’s rotational phase, at least in principle, one also has to account for aberration. Because of aberration, the rotational phase of the pulsar at which a signal is sent towards Earth periodically changes as the pulsar moves around the center of mass of the binary system [12,17]. While aberration can lead to periodic changes in the pulse profile, as we have discussed in Section 2.2, aberration contributions to the TOAs, however, become relevant only if there is a change in the spin orientation of the pulsar. For a fixed spin direction, the aberration delay ${\Delta }_{\mathrm{A}}$ cannot be separated from the Rømer delay ${\Delta }_{\mathrm{R}}$ [12,13]. Having said that, for pulsars in relativistic binary orbits, one might expect a change in the spin orientation, as we discuss next.

3.3. Spin Precession

If the pulsar spin is misaligned with respect to the orbital motion, then relativistic spin-orbit coupling leads to the so-called geodetic precession of the spin around the orbital angular momentum $\overrightarrow{L}$. Since $\overrightarrow{L}$ is many orders of magnitude larger than the pulsar spin, it very nearly coincides with the total angular momentum of the system. The rate of precession for a binary system has first been worked out within GR in [31]. Averaged over one orbit, one finds

$${\mathsf{\Omega }}_{\mathrm{p}}^{\mathrm{SO}}=\frac{n_{\mathrm{b}}}{1-e^2}\left(2+\frac{3m_{\mathrm{c}}}{2m_{\mathrm{p}}}\right)\frac{m_{\mathrm{p}}{m}_{\mathrm{c}}}{M^2}{\beta }_{\mathrm{O}}^2,$$

(10)

For the test particle limit ($m_{\mathrm{p}}\to 0$), Equation (10) approaches the well known expression for the de Sitter-Fokker precession [32,33].

A change in the orientation of the pulsar spin with respect to the LOS has two major observable effects. First, it leads to changes (intensity and shape) in the pulse profile as the LOS cuts the active emission region of the pulsar differently at different epochs [16,18,20]. Second, it changes the shape of the characteristic swing of the linear polarization (Equation (3)), as well as the absolute value of the position angle of the polarization in the sky [13,16,22]. A test of the relativistic spin precession in binary pulsars is therefore generally based on a, to some extent model dependent, pulse structure analysis (see Section 2.2). A somewhat unique exception will be discussed in Section 4.2.

3.4. Testing General Relativity with Post-Keplerain Parameters

In most cases the two masses of a binary pulsar, $m_{\mathrm{p}}$ and $m_{\mathrm{c}}$, are apriori unknown. There are exceptions, for instance, where the companion is a bright WD which allows for high-resolution spectroscopy, giving access to the masses with the help of optical observations (see Section 5), but generally one needs the measurement of at least three PK parameters for a test of GR (see e.g., [13,34]). With two PK parameters the system is fully determined for a given gravity theory, as both of the masses can be computed based on that theory. Once a third PK parameter has been measured the system is over-determined, and therefore provides a test of the theory assumed. Only if all three PK parameters agree on a certain set of masses, i.e., a common region in the mass-mass plane, the theory used has passed the test. This will become clearer in Section 4.

3.5. Alternative Gravity Theories and Pulsar Timing

GR is based on two basic postulates [1,35,36]. The first one is the postulate of universal coupling between the matter fields ${\mathsf{\Psi }}_{\mathrm{mat}}$ and gravity, meaning that the matter part of the total action has the form ${\mathcal{S}}_{\mathrm{mat}}[{\mathsf{\Psi }}_{\mathrm{mat}};g_{\mu \nu }]$, where $g_{\mu \nu }$ is the spacetime metric. Theories that fulfil this principle are called metric theories of gravity [3]. The second postulate is the field equations that describe the dynamics of $g_{\mu \nu }$, and can be derived from the total action [37]

$$\mathcal{S}={\mathcal{S}}_{\mathrm{mat}}[{\mathsf{\Psi }}_{\mathrm{mat}};g_{\mu \nu }]+\frac{c^4}{16\pi G}\int R\left[g_{\mu \nu }\right]\sqrt{-g}{d}^{4}x,$$

(11)

where $g\equiv det\left(g_{\mu \nu }\right)$ and R is the curvature/Ricci scalar. Under the assumptions stated in Lovelock’s theorem [38,39] (see also [40]), GR emerges as the unique theory of gravity. However, there are various ways to build a viable gravity theory that circumvents Lovelock’s theorem (see Figure 1 in [41]). One of the most popular approaches is the introduction of additional dynamical fields in the second term of the right hand side of Equation (11), while keeping the postulate of universal coupling. Amongst these theories, scalar-tensor theories are presumably the most well studied alternatives to GR [42].

Introducing additional fields in Equation (11), generally has a significant impact on the structure and motion of compact bodies (see a detailed discussion in [3]). The “gravitational charges” of bodies linked to the additional gravitational fields are generally structure dependent, and can be very different for a NS, as compared to Solar System bodies or WDs. At the same time, these additional gravitational fields can have a significant impact on the motion of a binary pulsar and the propagation of its radio signals in the spacetime of the binary system. Pulsar experiments are therefore sensitive to deviations from GR resulting from the highly relativistic interior of NSs. In this sense, pulsar experiments are strong-field gravity tests, even though their orbital velocities are small compared to the speed of light [3,43]. In other words, pulsars probe the mildly relativistic strong-field regime of gravity [44].

Structure dependent modifications of the gravitational interaction enter already at the Newtonian level, by introducing a body dependent effective gravitational constant ${\widehat{G}}_{ab}$ ($={\widehat{G}}_{ba}$) for the interaction between mass a and mass b. For a wide range of fully conservative theories of gravity, including scalar-tensor theories, the orbital dynamics up to first post-Newtonian approximation is given by the modified Einstein-Infeld-Hoffmann formalism (mEIH) [45]. The 1PN equations of motion for two strongly self-gravitating bodies can be derived from the following Lagrangian

$$\begin{array}{ccc}\hfill {\mathcal{L}}_{\mathrm{O}}^{\mathrm{mEIH}}& =& m_{\mathrm{p}}\left(\frac{{\overrightarrow{v}}_{\mathrm{p}}^2}{2}+\frac{{\overrightarrow{v}}_{\mathrm{p}}^4}{8c^2}\right)+m_{\mathrm{c}}\left(\frac{{\overrightarrow{v}}_{\mathrm{c}}^2}{2}+\frac{{\overrightarrow{v}}_{\mathrm{c}}^4}{8c^2}\right)\hfill \\ & & +\frac{{\widehat{G}}_{\mathrm{pc}}{m}_{\mathrm{p}}{m}_{\mathrm{c}}}{r}\left[1-\frac{(\overrightarrow{r}·{\overrightarrow{v}}_{\mathrm{p}})(\overrightarrow{r}·{\overrightarrow{v}}_{\mathrm{c}})}{2c^2{r}^2}-\frac{{\overrightarrow{v}}_{\mathrm{p}}·{\overrightarrow{v}}_{\mathrm{c}}}{2c^2}+(1+2{\widehat{\gamma }}_{\mathrm{pc}})\frac{{({\overrightarrow{v}}_{\mathrm{p}}-{\overrightarrow{v}}_{\mathrm{c}})}^2}{2c^2}\right]\hfill \\ & & +\left(1-2\frac{m_{\mathrm{c}}}{M}{\widehat{\beta }}_{\mathrm{cc}}^{\mathrm{p}}-2\frac{m_{\mathrm{p}}}{M}{\widehat{\beta }}_{\mathrm{pp}}^{\mathrm{c}}\right)\frac{{\widehat{G}}_{\mathrm{pc}}^{2}M{m}_{\mathrm{p}}{m}_{\mathrm{c}}}{2c^2{r}^2},\hfill \end{array}$$

(12)

with the three structure dependent parameters ${\widehat{\gamma }}_{\mathrm{pc}}$, ${\widehat{\beta }}_{\mathrm{cc}}^{\mathrm{p}}$, and ${\widehat{\beta }}_{\mathrm{pp}}^{\mathrm{c}}$ at the 1PN level [34,45,46].2 Sub- and superscripts “p” and “c” correspond to pulsar and companion respectively. In the weak-field limit, these parameters approach their corresponding body-independent parameters of the PPN formalism, i.e., ${\widehat{\gamma }}_{\mathrm{pc}}={\gamma }_{\mathrm{PPN}}$ and ${\beta }_{\mathrm{cc}}^{\mathrm{p}}={\beta }_{\mathrm{pp}}^{\mathrm{c}}={\beta }_{\mathrm{PPN}}$. In this sense, they can be seen as the strong-field analogues of the Eddington parameters [47]. Depending on the details of the theory under investigation, in the presence of NSs these 1PN strong-field parameters can be very different from their PPN counterparts. This is particularly prominent in the presence of phenomena like spontaneous scalarization [48], where the weak-field PPN parameters can be identical to GR but ${\widehat{\gamma }}_{\mathrm{pc}}$ and ${\widehat{\beta }}_{\mathrm{pc}}$ can have large deviations. In GR, where the internal structure of a body is effaced, one has ${\widehat{\gamma }}_{\mathrm{pc}}={\widehat{\beta }}_{\mathrm{cc}}^{\mathrm{p}}={\widehat{\beta }}_{\mathrm{pp}}^{\mathrm{c}}=1$.

The structural similarity between Equations (4) and (12) results in the fact, that for both the solution to the equations of motion can be cast into the same simple quasi Keplerian structure of [24]. It is then only the parameters as function of energy, angular momentum, and the masses which are different. For pulsar timing this means, that the DD timing formula can be applied as a phenomenological model, where the PK parameters as functions of the Keplerian parameters and the (inertial) masses are theory dependent [13]. For instance, for the precession of periastron, Equation (5) changes to

$$k=\left(2+2{\widehat{\gamma }}_{\mathrm{pc}}-\frac{m_{\mathrm{c}}}{M}{\widehat{\beta }}_{\mathrm{cc}}^{\mathrm{p}}-\frac{m_{\mathrm{p}}}{M}{\widehat{\beta }}_{\mathrm{pp}}^{\mathrm{c}}\right)\frac{{\widehat{\beta }}_{\mathrm{O}}^2}{1-e^2},$$

(13)

where ${\widehat{\beta }}_{\mathrm{O}}\equiv {\left({\widehat{G}}_{\mathrm{pc}}M{n}_{\mathrm{b}}\right)}^{1/3}/c$. We will demonstrate below that the statement of the (rather) general applicability of the DD model also holds for the Einstein and the Shapiro delay.

The damping of the orbital motion and therefore the change in the orbital period, i.e., Equation (6), also gets generally modified in alternatives to GR, as the additional (dynamical) gravitational fields contribute to the loss of orbital energy (and angular momentum). In asymmetric binaries, where pulsar and companion have different amount of “gravitational charges”, the orbital motion gives rise to a time-varying gravitational dipole which is the source of dipolar GWs [3,49]. The GW damping related to the dipolar GWs generally enters the equations of motion already at the 1.5PN level, i.e., order ${\widehat{\beta }}_{\mathrm{O}}^3$ (see e.g., [50]). This contribution is therefore (formally) many orders of magnitude stronger than the quadrupolar GW damping of GR. However, the actual amount of dipolar radiation depends quite sensitively on the difference in the gravitational charges. For that reason pulsar-WD systems are particularly interesting to test this aspect of GW generation, as we show in Section 5. An explicit equation for the orbital period change due to dipolar GW will be given further below, within a specific class of gravity theories.

Time dilation, i.e., Equation (6), gets modified as well, since both the orbital velocity of the pulsar and the gravitational redshift caused by the companion are usually different from GR. Quite generally, one finds the same orbital behaviour as in GR but a change in the amplitude of the Einstein delay according to

$${\gamma }_{\mathrm{E}}=\frac{P_{\mathrm{b}}}{2\pi }e\left(\frac{{\widehat{G}}_{0\mathrm{c}}}{{\widehat{G}}_{\mathrm{pc}}}+\frac{m_{\mathrm{c}}}{M}\right)\frac{m_{\mathrm{c}}}{M}{\widehat{\beta }}_{\mathrm{O}}^2,$$

(14)

where ${\widehat{G}}_{0\mathrm{c}}$ is the effective gravitational constant between the companion and a test particle, i.e., the effective gravitational constant that enters the $g_{00}$ component of the companion’s spacetime metric (cf. Equation (16)) [13,45].

In alternative gravity theories, the local gravitational constant $G_{\mathrm{loc}}$ as seen by the pulsar when moving in the gravitational field of the companion is generally expected to change as a function of the distance r [49]. As a consequence, the pulsar’s compactness and therefore its moment of inertia $I_{\mathrm{p}}$ change along the orbit, leading to a periodic variations in the rotational velocity of the pulsar while its angular momentum (i.e., its spin $S_{\mathrm{p}}=2\pi {\nu }_{\mathrm{p}}{I}_{\mathrm{p}}$) is conserved. To leading order one generally has $G_{\mathrm{loc}}=G_0(1-{\widehat{\eta }}_{\mathrm{c}}{m}_{\mathrm{c}}/c^{2}r)$, where ${\widehat{\eta }}_{\mathrm{c}}$ is a body-dependent parameter related to the gravitational field of the companion [45]. As a consequence, the variation in the time of emission of the pulsar signal due to the change of $G_{\mathrm{loc}}$ has the same orbital dependency as the time dilation effect. It is therefore part of the observed Einstein delay amplitude, i.e., ${\gamma }_{\mathrm{E}}^{\mathrm{new}}={\gamma }_{\mathrm{E}}+\delta {\gamma }_{\mathrm{E}}$ with

$$\delta {\gamma }_{\mathrm{E}}=\frac{P_{\mathrm{b}}}{2\pi }e\left(\frac{{\kappa }_{\mathrm{p}}{\widehat{\eta }}_{\mathrm{c}}}{{\widehat{G}}_{\mathrm{pc}}}\right)\frac{m_{\mathrm{c}}}{M}{\widehat{\beta }}_{\mathrm{O}}^2,$$

(15)

where ${\kappa }_{\mathrm{p}}\equiv -\partial (ln{I}_{\mathrm{p}})/\partial (lnG)$ measures the change of the moment of inertia due to a change in the gravitational constant [51]. In GR, where there is no variation of the local gravitational constant, one has $\delta {\gamma }_{\mathrm{E}}=0$.

The modification of the gravitational field of the companion also affects the signal propagation, i.e., the Shapiro delay Equation (8). If the spacetime metric of the companion in the first post-Minkowskian approximation has the form ($i,j=1,2,3$)

$$g_{00}=-1+2\frac{{\widehat{G}}_{0\mathrm{c}}{m}_{\mathrm{c}}}{c^{2}r},g_{0j}=0,g_{ij}=\left(1+2{\widehat{\gamma }}_{0\mathrm{c}}\frac{{\widehat{G}}_{0\mathrm{c}}{m}_{\mathrm{c}}}{c^{2}r}\right){\delta }_{ij},$$

(16)

as it is the case for many alternatives to GR, then the functional form of the Shapiro delay as given in Equation (8) remains unchanged [13,45]. Also the Shapiro shape parameter $s_{\mathrm{S}}$ can still be identified with $sini$. It is only the range of the Shapiro delay in Equations (9) that gets modified according to

$$r_{\mathrm{S}}=\frac{1+{\widehat{\gamma }}_{0\mathrm{c}}}{2}\frac{{\widehat{G}}_{0\mathrm{c}}{m}_{\mathrm{c}}}{c^3}.$$

(17)

Similarly to ${\widehat{\gamma }}_{\mathrm{pc}}$, ${\widehat{\gamma }}_{0\mathrm{c}}$ corresponds to ${\gamma }_{\mathrm{PPN}}$ in the weak-field limit. Note, while $s_{\mathrm{S}}$ is still the sine of the orbital inclination i, its connection to the masses via Kepler’s third law (with ${\widehat{G}}_{\mathrm{pc}}$) and the observed projected semi-major axis of the pulsar orbit, x, gets however modified with respect to GR: $s_{\mathrm{S}}=x{n}_{\mathrm{b}}M/m_{\mathrm{c}}{\widehat{\beta }}_{\mathrm{O}}$.

Finally, in [13] a generic expression for the geodetic precession of a pulsar is given, which reads

$${\mathsf{\Omega }}_{\mathrm{p}}^{\mathrm{SO}}=\frac{n_b}{1-e^2}\left[\frac{{\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}}{{\widehat{G}}_{\mathrm{pc}}}+\left(\frac{{\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}}{{\widehat{G}}_{\mathrm{pc}}}-\frac{1}{2}\right)\frac{m_{\mathrm{c}}}{m_{\mathrm{p}}}\right]\frac{m_{\mathrm{p}}{m}_{\mathrm{c}}}{M^2}{\widehat{\beta }}_{\mathrm{O}}^2,$$

(18)

where ${\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}$ is the body-dependent coupling function that enters the Lagrangian for the spin-orbit interaction. In GR one has ${\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}=2G$, and in the weak-field (PPN) limit one finds ${\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}=(1+{\gamma }_{\mathrm{PPN}})G$. Although the latter expression suggests ${\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}=(1+{\widehat{\gamma }}_{\mathrm{pc}}){\widehat{G}}_{\mathrm{pc}}$, in a fully generic approach ${\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}$ can be seen as an independent strong-field parameter, where even ${\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}\ne {\widehat{\mathsf{\Gamma }}}_{\mathrm{c}}^{\mathrm{p}}$ [13].

For a given (conservative) gravity theory, the strong-field parameters ${\widehat{G}}_{\mathrm{pc}}$, ${\widehat{G}}_{0\mathrm{c}}$, ${\widehat{\gamma }}_{\mathrm{pc}}$, ${\widehat{\gamma }}_{0\mathrm{c}}$, ${\widehat{\beta }}_{\mathrm{cc}}^{\mathrm{p}}$, ${\widehat{\beta }}_{\mathrm{pp}}^{\mathrm{c}}$, ${\kappa }_{\mathrm{p}}$, ${\widehat{\eta }}_{\mathrm{c}}$, ${\widehat{\mathsf{\Gamma }}}_{\mathrm{p}}^{\mathrm{c}}$, and ${\widehat{\mathsf{\Gamma }}}_{\mathrm{c}}^{\mathrm{p}}$ can be calculated from first principles as a function of the masses $m_{\mathrm{p}}$ and $m_{\mathrm{c}}$. Since they generally depend on the structure of the bodies, one also needs to specify the equation of state (EOS) for the body under consideration. In practice, there is some uncertainty for the EOS of NSs, which needs to be accounted for when testing a theory.

To illustrate the above, with the help of a specific gravity theory, we will use the mono-scalar-tensor theory of Damour & Esposito-Farèse (called DEF gravity in this review; see [48,52,53] for the details given below). This class of alternatives to GR is well studied and shows a number of effects that quite generally illustrate how gravity could deviate from GR, in particular in the presence of strongly self-gravitating bodies. Besides a spacetime metric $g_{\mu \nu }$, it contains a mass-less scalar field $\phi$, with asymptotic value ${\phi }_0$ at spatial infinity. The field equations of DEF gravity can be derived from the (Einstein frame) action

$$\mathcal{S}={\mathcal{S}}_{\mathrm{mat}}[{\mathsf{\Psi }}_{\mathrm{mat}};{\tilde{g}}_{\mu \nu }]+\frac{c^4}{16\pi G_{*}}\int (R\left[g_{\mu \nu }\right]-2g^{\mu \nu }{\partial }_{\mu }\phi {\partial }_{\nu }\phi )\sqrt{-g}{d}^{4}x,$$

(19)

where all matter fields couple universally to the physical (Jordan) metric ${\tilde{g}}_{\mu \nu }\equiv g_{\mu \nu }exp[2{\alpha }_0(\phi -{\phi }_0)+{\beta }_0{(\phi -{\phi }_0)}^2]$. The Newtonian gravitational constant, as measured in a Cavendish-type experiment, is related to the bare gravitational constant $G_{*}$ by $G=G_{*}(1+{\alpha }_0^2)$. The parameters ${\alpha }_0$ and ${\beta }_0$ define the two-dimensional space of DEF gravity. Jordan-Fierz-Brans-Dicke gravity [54,55,56] corresponds to ${\beta }_0=0$, and GR corresponds to ${\alpha }_0={\beta }_0=0$. The quantities (“gravitational form factors”) of a body with mass $m_a$ and moment of inertia $I_a$ that enter the PK parameters are

$${\alpha }_a\equiv \frac{\partial ln{m}_a}{\partial {\phi }_0},{\beta }_a\equiv \frac{\partial {\alpha }_a}{\partial {\phi }_0},{\mathcal{K}}_a\equiv -\frac{\partial ln{I}_a}{\partial {\phi }_0},$$

(20)

where the number of baryons is kept fixed when taking the partial derivatives. The quantity ${\alpha }_a$ is the effective scalar coupling of the body and gives its specific scalar charge. For weakly self-gravitating masses ${\alpha }_a$ approaches ${\alpha }_0$. In the parameter space where spontaneous scalarization does occur for NSs (${\beta }_0\lesssim -4.5$), ${\alpha }_a$ can be of order unity even if ${\alpha }_0=0$. The effective gravitational constant for the interaction of two bodies is given by ${\widehat{G}}_{ab}=G_{*}(1+{\alpha }_a{\alpha }_b)$, and for the strong-field PPN parameters one finds

$${\widehat{\gamma }}_{ab}=1-\frac{2{\alpha }_a{\alpha }_b}{1+{\alpha }_a{\alpha }_b}\mathrm{and}{\widehat{\beta }}_{bc}^a=1+\frac{{\beta }_a{\alpha }_b{\alpha }_{\mathrm{c}}}{2(1+{\alpha }_a{\alpha }_b)(1+{\alpha }_a{\alpha }_c)}.$$

(21)

Consequently, the PK parameters read

$$\begin{array}{ccc}\hfill k& =& \left(\frac{3-{\alpha }_{\mathrm{p}}{\alpha }_{\mathrm{c}}}{1+{\alpha }_{\mathrm{p}}{\alpha }_{\mathrm{c}}}-\frac{m_{\mathrm{p}}{\alpha }_{\mathrm{p}}^2{\beta }_{\mathrm{c}}+m_{\mathrm{c}}{\alpha }_{\mathrm{c}}^2{\beta }_{\mathrm{p}}}{2M{(1+{\alpha }_{\mathrm{p}}{\alpha }_{\mathrm{c}})}^2}\right)\frac{{\widehat{\beta }}_{\mathrm{O}}^2}{1-e^2},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\gamma }_{\mathrm{E}}& =& \frac{P_{\mathrm{b}}}{2\pi }e\left(\frac{1+{\mathcal{K}}_{\mathrm{p}}{\alpha }_{\mathrm{c}}}{1+{\alpha }_{\mathrm{p}}{\alpha }_{\mathrm{c}}}+\frac{m_{\mathrm{c}}}{M}\right)\frac{m_{\mathrm{c}}}{M}{\widehat{\beta }}_{\mathrm{O}}^2,\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill r_{\mathrm{S}}& =& \frac{G_{*}{m}_{\mathrm{c}}}{c^3},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill s_{\mathrm{S}}& =& \frac{x{n}_{\mathrm{b}}}{{\widehat{\beta }}_{\mathrm{O}}}\frac{M}{m_{\mathrm{c}}}.\hfill \end{array}$$

(25)

Above we did not provide an expression for ${\mathsf{\Omega }}_{\mathrm{p}}^{\mathrm{SO}}$, since so far this quantity has not played any role in constraining DEF gravity. For the orbital period decay due to dipolar GWs one has

$${\dot{P}}_{\mathrm{b}}^{\mathrm{dipole}}=-2\pi \frac{1+e^2/2}{{(1-e^2)}^{5/2}}\frac{m_{\mathrm{p}}{m}_{\mathrm{c}}}{M^2}\frac{{({\alpha }_{\mathrm{p}}-{\alpha }_{\mathrm{c}})}^2}{1+{\alpha }_{\mathrm{p}}{\alpha }_{\mathrm{c}}}{\widehat{\beta }}_{\mathrm{O}}^3+\mathcal{O}\left({\widehat{\beta }}_{\mathrm{O}}^5\right).$$

(26)

Apart from the dipole contribution, there are also monopole and quadrupole contributions related to the scalar field. Both of them start at the 2.5PN level ($\mathcal{O}\left({\widehat{\beta }}_{\mathrm{O}}^5\right)$) and are usually subdominant to ${\dot{P}}_{\mathrm{b}}^{\mathrm{dipole}}$. Detailed expressions can be found in [52].

There are various theories that deviate from above assumptions. For instance, semi-conservative theories that admit a global preferred frame of reference for the gravitational interaction, like Einstein-Æther and Bekenstein’s TeVeS. We will not go into details here, and refer the interested reader to [3,57], and references therein. Some of these alternatives have recently been tightly constrained or even been ruled out by the multi-messenger observations of the double NS merger GW170817 [58], which confirmed the speed of the tensorial GW modes to be practically equal to c within about one part in ${10}^{15}$ (see e.g., [59]). Note, however, this is not the case for DEF-like gravity theories, where the tensorial GW modes travel at the speed of light (if the scalar field is mass-less, this also holds for the scalar mode) [3].

4. Gravitational Interaction of Strongly Self-Gravitating Bodies

Binary pulsars where the companion is likewise a NS are in their orbital dynamics driven by the interaction of two strongly self-gravitating bodies. If the binary orbit is sufficiently compact ($P_{\mathrm{b}}\lesssim 1$ day), in order to show prominent relativistic effects in the orbital motion of the pulsar, then this should allow a number of unique tests of GR and its alternatives (see Section 3). Furthermore, if the system is seen sufficiently close to edge-on, meaning its orbital inclination i being close to 90${}^{°}$, then this allows to probe how the radio signals from the pulsar travel in a spacetime geometry sourced by the strongly self-gravitating companion.3 Currently, out of the ∼300 binary pulsars there are about 20 which have (most likely) a NS companion, half of which have an orbital period of less than one day [6]. Quite a few of the latter group has turned out to be excellent gravity laboratories. Notably, the very first binary pulsar discovered, PSR B1913+16, belongs to them.

4.1. Hulse-Taylor Pulsar and the Existence of Gravitational Waves

PSR B1913+16, which by now called the Hulse-Taylor pulsar, was discovered in summer 1974 by Russell Hulse and Joseph Taylor [5]. It is a mildly recycled 59-ms pulsar in a quite eccentric ($e=0.62$) 7.8-h orbit with an unseen companion, which almost certainly is a NS as well. Already a few months after its discovery, a significant change in the longitude of periastron was measured, amounting to $\dot{\omega }=4.2\mathrm{deg}{\mathrm{yr}}^{-1}$ [60]. This relativistic precession of periastron is 35,000 times larger than the GR contribution to the precession of the Mercury perihelion. At that stage it was therefore already clear that the discovery of the Hulse-Taylor pulsar would be a milestone in experimental gravity.

Since the two masses, $m_{\mathrm{p}}$ and $m_{\mathrm{c}}$ were apriori unknown, Taylor and his group had to wait a few years to measure two more PK parameters, in order to conduct a GR test (cf. Section 3.4). By the end of 1978, besides a greatly improved measurement of $\dot{\omega }$, the Einstein delay amplitude ${\gamma }_{\mathrm{E}}$ and the change in the orbital period ${\dot{P}}_{\mathrm{b}}$ had been measured [61]. At that time, the measurement error ${\dot{P}}_{\mathrm{b}}$ was still about 20%, which nevertheless marked the first confirmation of the existence of GWs as predicted by GR, seen through their back reaction onto the emitting system. Over the time the precision on these three PK parameters has greatly improved [62,63,64,65], leading to a precise measurement of the masses ($m_{\mathrm{p}}=1.438\pm 0.001M_{\odot }$ and $m_{\mathrm{c}}=1.390\pm 0.001M_{\odot }$). The observed change in the orbital period is by now measured with a precision of 0.04%, i.e., ${\dot{P}}_{\mathrm{b}}^{\mathrm{obs}}=(-2.423\pm 0.001)\times {10}^{-12}$ [65]. This, however, is not the intrinsic change of the orbital period purely due to GW emission. ${\dot{P}}_{\mathrm{b}}^{\mathrm{obs}}$ is “contaminated” by an apparent orbital period change due to the Shklovskii effect (an apparent acceleration due to the transverse motion of the pulsar system [66]) and the Galactic differential acceleration between the SSB and the pulsar system [67]. In order to correct for these contributions, one needs to know the distance to the Hulse-Taylor pulsar. Unfortunately, there is a large uncertainty in the estimation of that distance ($9\pm 3$ kpc [68]), which dominates the error in the determination of the intrinsic ${\dot{P}}_{\mathrm{b}}$. As a result, on gets for the ratio between the intrinsic and the predicted orbital period change ${\dot{P}}_{\mathrm{b}}^{\mathrm{intr}}/{\dot{P}}_{\mathrm{b}}^{\mathrm{GR}}=0.9983\pm 0.0016$ [65]. The GR prediction has been computed from Equation (6), based on the masses determined from Equations (5) and (7). In spite of the limitation due to the distance uncertainty, currently the Hulse-Taylor pulsar is the second best test for the quadrupole formula of GW emission.

Apart from the $\dot{\omega }$-${\gamma }_{\mathrm{E}}$-${\dot{P}}_{\mathrm{b}}$ test, the Hulse-Taylor pulsar in the meantime allowed for the detection of a Shapiro delay in that system, and there is a hint on the presence of the relativistic deformation of the orbit ${\delta }_{\theta }$ [65].

4.2. The Double Pulsar: A Wealth of Relativistic Effects

The year 2003 saw the discovery of a truly remarkable system, the so-called Double Pulsar (PSR J0737−3039A/B) [69,70]. The Double Pulsar consists of two active radio pulsars, a mildly recycled 23 ms pulsar (pulsar A) and a non-recycled slow pulsar with 2.8 s rotational period (pulsar B). They move around each other in a mildly eccentric ($e=0.089$) orbit in a bit less than 2.5 h. The system is therefore clearly more relativistic than the Hulse-Taylor pulsar (Section 4.1), showing an advance of periastron $\dot{\omega }$ of nearly 17 degrees per year. From the simultaneous timing observations of A and B one has the projected semi-major axes of both of the orbits, which directly gives the mass ratio of the system, since $q\equiv m_{\mathrm{A}}/m_{\mathrm{B}}=x_{\mathrm{B}}/x_{\mathrm{A}}=1.0714\pm 0.0011$ [70,71]. In fact, this ratio is correct up to 1PN order for any Lorentz-invariant gravity theory (see discussion in [34]).

In terms of relativistic effects, the Double Pulsar shows, to begin with, the ones already known from the Hulse-Taylor pulsar [71]: $\dot{\omega }$, ${\gamma }_{\mathrm{E}}$, and ${\dot{P}}_{\mathrm{b}}$. The GW damping test from ${\dot{P}}_{\mathrm{b}}$ by now has reached a level of well below 0.1% [15], which is currently the best test for the quadrupole formula of GW generation in GR. Compared to the Hulse-Taylor pulsar, the Double Pulsar has the big advantage of being comparably close by (distance ∼1 kpc), where Shklovskii and Galactic corrections to ${\dot{P}}_{\mathrm{b}}$ are small. Furthermore, there are direct ways of measuring its distance via very long baseline interferometry (VLBI) [72], and eventual via pulsar timing [22].

It is interesting to compare this with tests of GW emission by the LIGO/Virgo collaboration [73,74]. At the 0PN level in the radiation reaction force (corresponding to the 2.5PN level in the equations of motion), the Double Pulsar is about three orders of magnitude more constraining. This is also the case for the dipolar term (−1PN in the radiation reaction force, 1.5PN in the equations of motion). One has to keep in mind, however, that in particular the LIGO/Virgo tests based on binary black hole mergers probe the GW generation of a completely different class of objects, and it therefore depends on the details of a gravity theory how these limits can be compared. Needless to say, that the Double Pulsar becomes very quickly unconstraining as one moves to higher PN orders, due to its comparably small orbital velocity ($v\sim 0.002c$).

The Double Pulsar is also unique in its orientation. It is the most highly inclined binary pulsar known to date. The tilt of the orbital plane with respect to the LOS is so small, that at every superior conjunction of pulsar A, for about 30 s its signals get (periodically) blocked by the rotating plasma-filled magnetosphere of pulsar B [75,76]. From modelling these eclipses, [77] found an inclination of $i={89.3}^{°}\pm {0.1}^{°}$. Such a high inclination leads to a very prominent Shapiro delay of about 130 $\mathsf{\mu }$s, allowing the measurement of the Shapiro range, $r_{\mathrm{S}}$, and the Shapiro shape, $s_{\mathrm{S}}$, parameters. The shape parameter, which can be identified with $sini$ within GR and a broad class of other gravity theories (see Section 3), has been measured with a precision of $0.05\%$ already a few years after the discovery of the Double Pulsar [71]. The resulting (timing) inclination agrees well with the inclination derived from the modelling of the eclipses of A [77]. The range $r_{\mathrm{S}}$ is somewhat less precise (a few %), but agrees well with the mass of B derived from other PK parameters. This is the most precise test for the propagation of photons in the spacetime of a strongly self-gravitating object. When looking at the maximum spacetime curvature probed, this is many orders of magnitude stronger than in Solar System signal-propagation experiments or in Event Horizon Telescope (EHT) observations [78], since the signals of A come as close as ∼10,000 km to B (see Section 9 for more details).

Finally, due to its nearly edge-on orientation, the Double Pulsar allowed for a quite novel test of relativistic spin precession (Equations (10) and (18)), by using the eclipses of pulsar A at its superior conjunction caused by the magnetosphere of pulsar B [79]. More details will be given in Section 4.3. Interestingly, the spin-precession of B around the orbital angular momentum, which is about 5 degrees per year, has changed the orientation of B such that since 2008 B is no longer visible from Earth [80]. This, however, does not pose a significant limitation to our gravity tests with the Double Pulsar system, since most of these tests do come from the anyway much more precise timing of A. It is not expected that pulsar A will suffer the same fate as B, since in the case of A the rotational axis is closely aligned with the orbital angular momentum [81].

Figure 3 gives the mass-mass diagram for the Double Pulsar, which combines all the GR tests discussed here. As one can see, GR passes all these tests with flying colors.

With this wealth on PK parameters and the mass ratio at hand, one can even derive some generic limits on the strong-field PK parameter of the modified Einstein-Infeld-Hoffmann formalism in Section 3.5 (see [22,79] for details).

4.3. Spin Precession

Shortly after the discovery of the Hulse-Taylor pulsar, the idea was put forward that this system could allow for a test of the spin precession resulting from relativistic spin-orbit coupling (geodetic precession) [82]. While eventually changes in the pulse profile could be observed and modelled under the assumption of geodetic precession [18,20,83], till today the Hulse-Taylor pulsar has not allowed for a quantitative test of the precession rate. Quantitative tests of the relativistic spin precession of a pulsar came then from other systems, some of them discovered decades later. Although the precision of these tests is moderate (order few %) compared to the test of geodetic precession with Lunar laser ranging (LLR) [84] or Gravity Probe B (GPB) [85], there are two qualitatively different aspects in binary pulsar tests. Firstly, it tests spin-orbit coupling in systems where both masses are of similar size, therefore testing terms beyond the test-particle limit. Secondly, it is a test with a strongly self-gravitating “gyroscope”, and therefore also providing constraints on strong-field effects in spin-orbit coupling. In the following we will summarize the three binary pulsars that so far have allowed for a good quantitative test of geodetic precession.

4.3.1. PSR B1534+12

PSR B1534+12 was the second pulsar with a NS companion discovered in the Galactic disk [86]. It is a mildly recycled 38 ms pulsar in an eccentric ($e=0.27$) 10.1-h orbit. By now the PK parameters $\dot{\omega }$, ${\gamma }_{\mathrm{E}}$, and ${\dot{P}}_{\mathrm{b}}$ have been measured with high precision. Unfortunately there is a large uncertainty in the distance to PSR B1534+12 leading to a large uncertainty in the Shklovskii correction for ${\dot{P}}_{\mathrm{b}}$. As a consequence, the GW damping in this system is not well determined. The PSR B1534+12 system has an orbital inclination of 78${}^{°}$, leading to a prominent Shapiro delay and the measurement of the PK parameters $r_{\mathrm{S}}$ and $s_{\mathrm{S}}$. Details on the latest timing results for this pulsar can be found in [87].

The spin of PSR B1534+12 is tilted by about $\delta ={27}^{°}$ (or 153${}^{°}$) with respect to the orbital angular momentum, as a result of the second supernova explosion in that system [87]. Consequently, orbital modulations of the pulse profile and the polarization along the orbit are expected, as well as long-term secular changes in these as the spin precesses with a rate of 0.514 ${\mathrm{deg}\mathrm{yr}}^{-1}$ (GR prediction) around the orbital angular momentum. And indeed, as a result of the good signal-to-noise ratio obtained for PSR B1534+12 with the 305-m William E. Gordon Arecibo radio telescope, such changes were first reported by [88], providing the first quantitative test of geodetic precession in a binary pulsar. While the uncertainties in [88] were still comparably large, a refined analysis based on a longer observing time span by [87] lead to a measurement of ${\mathsf{\Omega }}_{\mathrm{p}}^{\mathrm{SO}}={0.59}_{-0.08}^{+0.12}{\mathrm{deg}\mathrm{yr}}^{-1}$, in agreement with the GR prediction. Modelling the precession of PSR B1534+12 also lead to the actual measurement of the tilt $\delta$ of its spin (within the $\delta \to {180}^{°}-\delta$ ambiguity) obtaining a precision of about three degrees, with the value for $\delta$ already given above. Considerations from binary evolution and arguments with respect to asymmetric supernova explosions and resulting kicks onto the newborn NS suggest that the smaller misalignment angle, i.e., $\delta <{90}^{°}$, is more likely [89].

4.3.2. PSR J1906+0746

The so far best test for the geodetic precession of a pulsar comes from the non-recycled 144 ms pulsar PSR J1906+0746, which is in a relativistic 4.0-hour orbit with an unseen companion, that most likely is also a NS. In terms of certain properties, like the eccentricity ($e=0.085$) and the masses ($m_{\mathrm{p}}=1.29\pm 0.01M_{\odot }$ and $m_{\mathrm{c}}=1.32\pm 0.01M_{\odot }$), it shares similarities with the Double Pulsar, but here we only seem to see the non-recycled component [90].

For PSR J1906+0746 the (present) orientation of the spin and the magnetic axis is such that the LOS passes close to both of the magnetic poles while the pulsar rotates. Consequently, observations revealed the presence of two polarized components in the pulse profile about half a period apart and well described by Equation (3) [91]. Having this information from both of the emission regions allowed a precise determination of the viewing geometry. Regular monitoring with the 305-m William E. Gordon Arecibo radio telescope since 2012, revealed a gradual change in the orientation due to geodetic precession. In combination with earlier data from the Arecibo and Nançay telescopes, it was possible to determination of the precession rate with good precision: ${\mathsf{\Omega }}_{\mathrm{p}}^{\mathrm{SO}}=2.17\pm 0.11\mathrm{deg}{\mathrm{yr}}^{-1}$, in excellent agreement with the GR value of 2.23 deg yr${}^{-1}$ [91].

4.3.3. PSR J0737–3039B

A quite unique approach to test the spin precession of a pulsar was taken in the Double Pulsar. As mentioned in Section 4.2, pulsar A gets eclipsed by the magnetosphere by pulsar B, every 2.45 h when it moves through its superior conjunction [75]. The eclipse lasts for about 30 s. Soon it had been noticed in high signal-to-noise observations from the 100-m Robert C. Byrd Green Bank Telescope that the eclipse is not continuous during these 30 s [92]. In fact, as the “doughnut-shaped” magnetosphere of B rotates with B’s period of 2.8 s, there are phases where the signals from A can get through to Earth. The spacing between these “bright” phases during the eclipse encodes the orientation of B’s magnetosphere as it rotates and therefore the orientation of B’s spin. A simple geometrical model was proposed in [79], able to quantitatively describe the eclipse pattern over a period of several years during which the spin precession had significantly altered the direction of B’s spin, before it eventually precessed away from the LOS [80]. The model allowed for a fit of the precession rate, leading to ${\mathsf{\Omega }}_{\mathrm{B}}^{\mathrm{SO}}={4.77}_{-0.65}^{+0.66}$$\mathrm{deg}{\mathrm{yr}}^{-1}$ which agrees with the GR prediction of 5.073 $\mathrm{deg}{\mathrm{yr}}^{-1}$ (see also Figure 3). The result is somewhat less precise than the one from PSR J1906+0746 (Section 4.3.2), but comes from a completely different method. Moreover, the fact that the Double Pulsar gives a (generally) theory independent access to the mass ratio (see Section 4.2) allowed for a generic constraint on the body-dependent spin-orbit coupling function in Equation (18): ${\widehat{\mathsf{\Gamma }}}_{\mathrm{B}}^{\mathrm{A}}/2{\widehat{G}}_{\mathrm{AB}}=0.95\pm 0.11$ [22,79]. Note, in GR ${\widehat{\mathsf{\Gamma }}}_{\mathrm{B}}^{\mathrm{A}}/2{\widehat{G}}_{\mathrm{AB}}=1$.

5. Dipolar Gravitational Radiation

In the previous section we have seen that pulsars with NS companions provide the so far best test for the emission of quadrupolar gravitational waves as predicted by GR. However, as discussed in Section 3, in most alternatives to GR that predict additional gravitational charges (e.g., scalar charges) for NSs, the leading order contribution to the GW damping of the orbit comes from dipolar GWs, which enter the equations of motion already at the 1.5PN level (see Section 3). Yet, a significant emission of dipolar GWs generally requires a sufficient asymmetry in the compactness (respectively in the fractional binding energy) of the two masses in a binary system.4 For that reason, the double NS systems of Section 4, where $m_{\mathrm{p}}\sim m_{\mathrm{c}}$, provide only weak limits on dipolar GWs, in spite of the high precision obtained in the measurement of their orbital period decay. Sufficient asymmetry in compactness can be found in binary pulsars with WD companions. While the compactness $\mathcal{C}\equiv Gm/R{c}^2$ for typical pulsar (mass $m\sim 1.4M_{\odot }$) is of the order of 0.2, for WDs one finds, due to their comparably large radii R, $\mathcal{C}\lesssim {10}^{-3}$. That makes such systems certainly interesting for dipolar radiation tests.

The large majority of the ∼300 binary pulsars have a WD as a companion, and quite a few of them are in compact orbits with $P_{\mathrm{b}}\lesssim 1$ day. Some of the most precise “pulsar clocks” are found in such systems. They are fully recycled MSPs, which then generally means that the binary orbits have very low eccentricities (see Section 2). As a consequence, often neither the advance of periastron ($\dot{\omega }$) nor the Einstein delay (${\gamma }_{\mathrm{E}}$) is observable in such systems. And if the system is not seen sufficiently edge on, meaning that no Shapiro delay is detected, then the orbital period change ${\dot{P}}_{\mathrm{b}}$ is often the only PK parameter measurable. Still, a small number of these systems provide some of the best tests of dipolar GWs. In fact, the current best limits on dipolar GW damping do come from PSR J1738+0333, a MSP where the only relativistic effect measured from timing is the orbital period decay ${\dot{P}}_{\mathrm{b}}$.

5.1. PSR J1738+0333: Constraining Dipolar Gravitational Radiation

PSR J1738+0333 is a 5.9 ms pulsar with a Helium WD companion. Pulsar and WD orbit each other in a nearly circular orbit ($e<4\times {10}^{-7}$) in just 8.5 h (see [94] for the latest timing results). The orbital inclination of the system is only 33${}^{°}$, and a Shapiro delay therefore undetectable, as it is completely covariant with the Rømer delay (cf. [95]). Fortunately, the Helium WD companion is seen optically, with prominent Balmer lines, which eventually allowed the determination of the WD mass ($m_{\mathrm{c}}={0.181}_{-0.005}^{+0.007}{M}_{\odot }$) and, in combination with the timing observations, the mass ratio ($q\equiv m_{\mathrm{p}}/m_{\mathrm{c}}=8.1\pm 0.2$) [96]. With this, and the well measured Keplerian parameters, the system is fully determined.

Extensive timing observations, in particular with the 305-m William E. Gordon Arecibo radio telescope, allowed the measurement of the distance to the system ($1.5$ kpc), and more importantly the observation of the orbital period decay due to GW emission. After correcting for the Shklovskii effect and the Galactic differential acceleration, [94] obtained for the (intrinsic) change in the orbital period

$${\dot{P}}_{\mathrm{b}}^{\mathrm{intr}}=-25.9\pm 3.2\mathrm{fs}{\mathrm{s}}^{-1},$$

(27)

which agrees well with the GW damping predicted by GR (see Figure 4).

Compared to the <$0.1$% GW test with the Double Pulsar (Section 4.2), the limit (27) appears to be quite weak and of no interest. However, in the dipolar radiation test the lack of precision is more than compensated by the high asymmetry of the system. To illustrate this, one can do a comparison of the compactness between pulsar and companion. For the Double Pulsar $\Delta \mathcal{C}={\mathcal{C}}_{\mathrm{A}}-{\mathcal{C}}_{\mathrm{B}}\approx {10}^{-2}$, while for the PSR J1738+0333 system one finds $\Delta \mathcal{C}={\mathcal{C}}_{\mathrm{p}}-{\mathcal{C}}_{\mathrm{c}}\simeq {\mathcal{C}}_{\mathrm{p}}\approx 0.2$ (these numbers depend to some extent on the chosen EOS). Keeping in mind that the figure of merit for dipolar radiation tests is proportional to $\Delta {\mathcal{C}}^2$, it becomes clear why PSR J1738+0333 can still provide better tests on certain aspects of GWs. Considerations based on $\Delta \mathcal{C}$ give only a rough estimate on the limits one can obtain from PSR J1738+0333. The actual limits on dipolar radiation depend on the specifics of the gravity theory under consideration. For DEF gravity Equation (26) applies, and one obtains the following constraint on the scalar dipole:

$$|{\alpha }_{\mathrm{p}}-{\alpha }_{\mathrm{c}}|<2\times {10}^{-3}(95\%\mathrm{C}.\mathrm{L}.).$$

(28)

Since ${\alpha }_{\mathrm{c}}\simeq {\alpha }_0$, this directly limits the (specific) scalar charge of the pulsar. As ${\alpha }_{\mathrm{p}}$ depends on the two parameters ${\alpha }_0$ and ${\beta }_0$ of the scalar coupling strength of DEF gravity, Equation (28) can be converted into constraints in the ${\alpha }_0$−${\beta }_0$ parameter space (details on this will be given in a comparison with other experiments in Section 6.1).

The limit Equation (28) is actually more general. To obtain the mass of the WD only well tested Newtonian physics was needed, and the mass ratio is also quite theory independent (see discussion in Section 4.2). Hence, the limit Equation (28) on the scalar dipole in DEF gravity can be generically seen as a limit on any gravitational dipole linked to dynamical gravitational fields (cf. [97]).

5.2. Other Systems for Dipolar Radiation Tests

Although, formally speaking, PSR J1738+0333 currently provides the best pulsar limit on the existence of a gravitational dipole, this does not necessarily mean that it generally leads to the best pulsar constraints on alternatives to GR in the radiative sector. This is even the case within the class of DEF gravity theories. In the highly non-linear area of DEF gravity, where one also finds the phenomenon of spontaneous scalarization (${\beta }_0\lesssim -4$) [48], the scalar charge of a NS depends crucially on the mass of the NS. On top of that, the scalarization behaviour is quite sensitive to the EOS for NS matter, which currently is not very well constrained. Figure 5 illustrated this for a specific choice of ${\alpha }_0$ and ${\beta }_0$ in DEF gravity. For that reason, other binary pulsar experiments become important, even if they provide a weaker limit for the difference in the scalar charges, i.e., the scalar dipole. In particular pulsars with a higher mass than PSR J1738+0333 are important to probe such potential strong-field deviations from GR. The timing observations of PSR J0348+0432, a two Solar mass pulsar in a 2.5-h orbit [98] was an important step in that direction, as one can see in Figure 5.

Depending on the parameters ${\alpha }_0$ and ${\beta }_0$ and the EOS for NS matter, PSR J1738+0333 and PSR J0348+0432 limits can still be circumvented (see [103] and EOS ENG in Figure 5). Including more systems with masses between 1.46 and 2.0 Solar masses is therefore necessary, and will eventually help to close this “mass gap” for spontaneous scalarization [104,105].

DEF gravity serves here only as an example on how nature can deviate in the presence of strongly self-gravitating masses. Quite generally, when it comes to gravity tests with pulsar it is important to obtain a good coverage in terms of NS masses. We will revisit this point in Section 8.

Finally, for gravity theories that invoke screening mechanisms to suppress deviations from GR in strong gravitational fields, the less compact WD might be more relevant than the pulsar, and therefore the criteria for gravity tests are more closely linked to the WD properties (see for instance [106]).

6. Universality of Free Fall

The universality of free fall (UFF) of (electrically neutral) test bodies stood right at the beginning of the genesis of GR, when in 1907 Einstein had the “most fortunate thought” of his life [107] that eventually led him to the insight that gravity is a manifestation of curved spacetime. To date the UFF for test masses, also known as the weak equivalence principle (WEP) [3], has been confirmed to a precision of ∼${10}^{-14}$ with the MICROSCOPE satellite experiment [108]. For that reason, many alternatives to GR are so called metric theories of gravity, which by design fulfill the WEP. In fact, metric theories of gravity fulfill, as a result of their postulate of universal coupling (see Section 3.5), the Einstein equivalence principle (EEP), which in addition to the WEP demands the fulfillment of the local position invariance (LPI) and the local Lorentz invariance (LLI) for the outcome of any non-gravitational experiment. In a sense, the WEP plays insofar a leading role here, as Schiff’s conjecture states that the WEP implies the full EEP for any consistent theory of gravity (see e.g., Section 2.4 in [3]).

While all metric theories of gravity agree on the EEP, they generally do not fulfill the strong equivalence principle (SEP), which extends the EEP to (local) gravitational experiments. In fact, GR might be the only valid gravity theory in four spacetime dimensions that fully embodies the SEP [3,109]5. This explains the importance of SEP tests for gravitational physics.

In the SEP, the WEP gets promoted to the gravitational WEP (GWEP), which states that in an external gravitational field all bodies fall the same, independent of their amount of gravitational self-energy. Solar System experiments have tightly constrained the violation of the GWEP, i.e., the presence of a so-called Nordtvedt effect. For instance, [84] find a confirmation of the GWEP at the ≲${10}^{-13}$ level in their analysis of LLR data. One has to keep in mind, however, that experiments in the Solar System can probe the GWEP only in the weak field at leading (linear) order in the fractional binding energy $ϵ$. To give an example, for the Earth $ϵ\approx -4.6\times {10}^{-10}$, which makes order ${ϵ}^2$ terms totally inaccessible to current and (foreseeable) future LLR experiments.

In [111] it has been suggested to use binary pulsars with WD companions for a test of the UFF of strongly self-gravitating bodies. In this experiment the pulsar ($ϵ\sim -0.15$) and its WD companion ($\left|ϵ\right|\lesssim {10}^{-3}$) fall in the external gravitational field of our Galaxy. A violation of the GEWP by the strong field of the NS would then lead to a gravitational Stark effect on the binary system, causing a characteristic change of its orbital eccentricity [111]. On the upside, the fractional binding energy of the pulsar is many orders of magnitude larger than that of any Solar System body. On the downside, the Galactic gravitational filed at the location of the pulsar is comparably weak, typically of the order of $2\times {10}^{-8}\mathrm{cm}{\mathrm{s}}^{-2}$. As a result, such UFF tests have not lead to any tight constraints on alternatives to GR (see [112] and references therein). For that reason, the discovery of a MSP in a stellar triple system by [113] was a real game changer (see also discussion in [114]).

6.1. The Pulsar in a Stellar Triple System

PSR J0337+1715 is a 2.7 ms pulsar in a hierarchical triple system with two WDs and orbital periods of 1.63 and 327 days [113]. The inner orbit consists of the pulsar and a 0.2 $M_{\odot }$ Helium WD (optically identified). The outer WD has a mass of 0.4 $M_{\odot }$, and is the source of the external gravitational field in the UFF test. The acceleration of the inner orbit due to the outer WD is about 0.17 cm s${}^{-2}$. Therefore the polarizing force, in case of a violation of GWEP, is ${10}^7$(!) times stronger than that for binary pulsars falling in the gravitational field of our Galaxy.

A violation of the GWEP manifests itself already at the Newtonian level in the three-body dynamics. It is best described by the presence of a body-dependent gravitational constant ${\widehat{G}}_{ab}$ (cf. Section 3.5). The acceleration of a mass $m_a$ in the gravitational field of the other masses is then determined according to

$${\ddot{\overrightarrow{x}}}_a=-\sum _{b\ne a}{\widehat{G}}_{ab}{m}_b\frac{{\overrightarrow{x}}_a-{\overrightarrow{x}}_b}{|{\overrightarrow{x}}_a-{\overrightarrow{x}}_b{|}^3}+\mathcal{O}\left(c^{-2}\right),$$

(29)

where ${\overrightarrow{x}}_i$ denotes the (coordinate) position of the ith body (see the modified Einstein-Infeld-Hoffmann equations for a N-body system in [13,45]). We use the definition ${\widehat{G}}_{ab}\equiv G(1+{\Delta }_{ab})$ with G denoting the Newtonian gravitational constant, as measured in a Cavendish-type experiment. A violation of GWEP in the hierarchical triple of PSR J0337+1715 leads to a non-vanishing

$${\delta }_{\mathrm{GWEP}}\equiv \frac{{\widehat{G}}_{\mathrm{AC}}-{\widehat{G}}_{\mathrm{BC}}}{({\widehat{G}}_{\mathrm{AC}}+{\widehat{G}}_{\mathrm{BC}})/2}\simeq {\Delta }_{\mathrm{AC}}-{\Delta }_{\mathrm{BC}},$$

(30)

where A denotes the pulsar, B the inner and C the outer WD.6 Since B and C are both weakly self-gravitating bodies, one quite generally finds that ${\Delta }_{\mathrm{BC}}\simeq 0$ can be assumed in Equation (30). This however depends on the details of the theory. For theories with screening mechanisms, for instance, the situation can be different (cf. discussion at the end of Section 5.2).

Already a few years after its discovery, a tight limit of order few times ${10}^{-6}$ on a violation of GWEP has been obtained from PSR J0337+1715 [115]. The so far best limit comes from regular timing observations with the Nançay radio telescope [116]:

$${\delta }_{\mathrm{GWEP}}=(0.5\pm 1.8)\times {10}^{-6}(95\%\mathrm{C}.\mathrm{L}.).$$

(31)

This limit is three orders of magnitude better than the best limit from a binary pulsar in the Galactic field [117].

A comparison of limit Equation (31) with limits from the Solar System is not at all straight forward, since these are two different regimes of gravity [112]. While for the Solar System calculations can be done within the generic PPN formalism, the presence of the strongly self-gravitating pulsar in the triple system requires the full non-linearity of the field equations of a theory to consistently calculate its motion [118]. In an expansion of ${\delta }_{\mathrm{GWEP}}\simeq {\Delta }_{\mathrm{AC}}$ with respect to the fractional binding energy of the pulsar $ϵ$, like in [111], one has

$${\delta }_{\mathrm{GWEP}}=\eta ϵ+{\eta }^{\prime }{ϵ}^2+\cdots$$

(32)

where $\eta$ is the so-called Nordtvedt parameter [119]. Solar System experiments already provide tight constraints for $\eta$ of the order of ${10}^{-4}$ [84,120]. The higher order terms, as mentioned already above, are out of reach for Solar System tests due to the small fractional binding energy of Solar System bodies. Then again, for PSR J0337+1715 with its $1.44M_{\odot }$ one finds (within GR) a typical of $ϵ\approx -0.13$, which makes that test sensitive to higher order contributions in Equation (32). If we ignore higher order contributions, then the limit Equation (31) converts into $\left|\eta \right|\lesssim 2\times {10}^{-5}$, which is significantly tighter than the Solar System limit. Ignoring higher order contributions, however, means rather restrictive assumptions about the strong-field behaviour of a gravity theory.

An alternative approach to compare pulsar experiments with tests in other regimes, like the Solar System or LIGO/Virgo observations, is the use of a theory-dependent framework [34]. For the triple-system test, such a comparison has been done on the basis of DEF gravity [115,116]. In DEF gravity (cf. Section 3.5)

$${\widehat{G}}_{ab}=G_{*}(1+{\alpha }_a{\alpha }_b)=G\frac{1+{\alpha }_a{\alpha }_b}{1+{\alpha }_0^2},$$

(33)

and consequently

$${\Delta }_{ab}=\frac{{\alpha }_a{\alpha }_b-{\alpha }_0^2}{1+{\alpha }_0^2}\simeq {\alpha }_a{\alpha }_b-{\alpha }_0^2,$$

(34)

where the experimental fact that ${\alpha }_0\ll 1$ has been used in the last step.7 Since ${\alpha }_{\mathrm{B}}\simeq {\alpha }_{\mathrm{C}}\simeq {\alpha }_0$ in the triple system, it is the quantity

$${\delta }_{\mathrm{GWEP}}\simeq {\alpha }_0({\alpha }_{\mathrm{A}}-{\alpha }_0)$$

(35)

that is being constrained by the limit Equation (31). Figure 6 shows the limit on the parameter space of DEF gravity, in comparison with other limits. Currently, PSR J0337+1715 provides the most stringent limits on DEF gravity for most of the parameter space. In particular for positive ${\beta }_0$, it will be difficult to surpass these limits in the near future with other experiments, including merger events observed with gravitational wave observatories [121]. This is not restricted to DEF gravity, but should be the case for many alternatives to GR that violate GWEP due to the presence of additional “gravitational charges” (see e.g., [122,123]).

6.2. UFF Towards Dark Matter

As mentioned above, limits on a violation of the GWEP coming from pulsar-WD binaries, falling in the gravitational field of our Galaxy, are generally in no way competitive to the ∼${10}^{-6}$ test conducted with the pulsar in the stellar triple (Section 6.1). There is an exception, however, which is the UFF towards dark matter (DM) [124]. If the attraction towards the DM distribution in our Galaxy is different between a NS and a WD, then a pulsar-WD system would experience a gravitational Stark effect, leading to a polarization of the binary orbit [111]. The result of such a polarization would be a temporal change of the orbital eccentricity. Wide orbits are particularly sensitive to such an effect. As it turns out, the MSP PSR J1713+0747 is particularly suitable for such a test, for several reasons. It is in a wide 68-day orbit with a WD companion. It is one of the most precisely timed pulsars with more than 20 years of high-precision timing data, leading to a firm limit on any change in the eccentricity. And from precise timing one could infer with high precision its location in Galaxy and the spatial orientation of the orbit. The latter is important to calculate how a potential polarizing force acts on the orbit. For details on that system and the latest timing parameters we refer the reader to [117]. Based on the results of [117] and a model of the DM distribution in the Galaxy, in [124] a limit on the difference in the accelerations $a_{\mathrm{DM}}$ toward DM was derived:

$$\left|\frac{a_{\mathrm{DM}}\left(\mathrm{PSR}\right)-a_{\mathrm{DM}}\left(\mathrm{WD}\right)}{[a_{\mathrm{DM}}\left(\mathrm{PSR}\right)+a_{\mathrm{DM}}\left(\mathrm{WD}\right)]/2}\right|<0.004(95\%\mathrm{C}.\mathrm{L}.).$$

(36)

There are two aspects which make that experiment particularly interesting in comparison with UFF experiments related to DM and conducted in the Solar System. On the one hand, it is the large material difference between the NS and the WD, where the composition of the NS is dominated by neutrons. On the other hand, it is the large (fractional) gravitational binding energy of the NS. A comparison with other experiments of the UFF towards DM can be found in [124].

7. Gravitational Symmetries

Apart from the universality of free fall (see Section 6, there are are two symmetries that enter the SEP [3,125], the local position invariance (LPI) and the local Lorentz invariance (LLI). For SEP, these symmetries get extended to the gravitational sector. LPI therefore states that the outcome of any local gravitational experiment is independent of when and where the experiment is performed. And LLI demands that a gravitational experiment does not depend on the velocity of the system. More specifically, the existence of a (global) preferred frame for the gravitational interaction would violate the LLI in the gravitational sector. Experiments with pulsars allow to test if these gravitational symmetries also hold in the presence of strongly self-gravitating bodies.

In the following we will summarize some of the pulsar tests for a spatial and temporal change in the local gravitational constant, and experiments related to symmetries that are linked to PPN parameters, more specifically their strong-field generalization.

Pulsars can also be used to test certain parameters of the standard model extension (SME; [126]). We will not cover this important field of experimental gravity here, but refer the reader to [104,127,128,129].

7.1. Variation of the Gravitational Constant

For alternative theories of gravity that violate the SEP, it is generally expected that the locally measured gravitational constant G changes with the expansion of the Universe. A classical example are scalar-tensor theories where the gravitational constant, as measured in a Cavendish-type experiment, depends on the cosmological value of the scalar field(s) [56,130]. Tight constraints on a (present day) variation of the gravitational constant do come from Solar System experiments, many of them limiting $|\dot{G}|/G$ to the order of ${10}^{-13}{\mathrm{yr}}^{-1}$, and even below [3]. The limit given in [120], based on the NASA MESSENGER mission to Mercury, is |$\dot{G}|/G<4\times {10}^{-14}{\mathrm{yr}}^{-1}$.

The orbital motion of binary pulsars gets affected by a temporal variation of G in two ways [131]. First, like for planetary orbits, it directly modifies the size of the orbit and consequently the orbital period of the pulsar due to a change in the universal gravitational coupling parameter, i.e., $G\left(t\right)$. Secondly, because of the compactness of the pulsar (and its companion if it is also a NS) a change in G has a significant impact on the mass of the pulsar, which in turn as well changes the gravitational attraction between pulsar and companion. More specifically, a change in G changes the gravitational binding energy of a self-gravitating body, and therefore it leads to a change of the pulsar and companion mass, i.e., $m_{\mathrm{p}}$ and $m_{\mathrm{c}}$. In total one has [104,131]

$$\frac{{\dot{P}}_{\mathrm{b}}}{P_{\mathrm{b}}}=-2\frac{\dot{G}}{G}\left[1-\left(1+\frac{m_{\mathrm{c}}}{2M}\right)s_{\mathrm{p}}-\left(1+\frac{m_{\mathrm{p}}}{2M}\right)s_{\mathrm{c}}\right],$$

(37)

where the “sensitivity”

$$s_a\equiv -{\left.\frac{\partial (ln{m}_a)}{\partial (lnG)}\right|}_N,$$

(38)

measures how the mass of body a changes with a change of the local gravitational constant G, for a fixed baryon number N (see e.g., Ref. [45] for details). For a given mass, $s_a$ depends on the specifics of the gravity theory and the nature of the star. WDs have very small sensitivities (comparable to their fractional binding energy ∼${10}^{-5}$⋯∼${10}^{-3}$), while for a NS $|s_a|$ can in principle be very large (≫1). For that reason, binary pulsar experiments for $\dot{G}$ either have to be interpreted in a theory specific context, or if generic have to make assumptions about the sensitivity. But it is also obvious that the test aspects of a change in G are related to the strong internal fields of NSs.

The best limit of $\dot{G}$ from binary pulsars currently comes from PSR J1713+0747, a pulsar in a wide orbit ($P_{\mathrm{b}}\simeq 67.8\mathrm{d}$) with a 0.3 $M_{\odot }$ WD companion [117], already introduced in Section 6.2. For a typical sensitivity of $s_p\sim 0.16$ for the 1.33 $M_{\odot }$ pulsar, in [117]

$$\dot{G}/G=(-1\pm 9)\times {10}^{-13}{\mathrm{yr}}^{-1}(95\%\mathrm{C}.\mathrm{L}.).$$

(39)

is found. This limit is still an order of magnitude weaker than the Solar System limits, however it comes from a very different technique with different experimental challenges, and it is an experiment with a strongly self-gravitating mass. The importance of the latter becomes clear in two ways. For instance in theories with $\dot{G}=0$, a NS could still change its mass with the cosmic expansion. For instance, in Barker’s constant-G theory [132] or in DEF gravity with ${\beta }_0=-(1+{\alpha }_0^2)$, a change of the background scalar field leads to a change of the NS mass and therefore a change in the orbital period of a binary pulsar, while the Solar System remains practically unaffected, since $\dot{G}=0$. Another regime, where binary pulsars can be superior to Solar System tests, is the strongly non-linear regime, for instance for ${\beta }_0\lesssim -4$ in DEF gravity. In such a regime the sensitivities in Equation (37) can become large. Moreover, in Equation (37) the gravitational constant needs to be extended by a body-dependent parameter $K_{pc}\left(t\right)$ [133], which has been neglected so far. Depending on the details of the theory, such an “effective gravitational constant” $G\left(t\right)K_{\mathrm{pc}}\left(t\right)$ ($={\widehat{G}}_{\mathrm{pc}}$ in Section 3.5) can respond much stronger to a change of, for instance, the background scalar field, than only $G\left(t\right)$ [112], making binary pulsars (with suitable masses) a more sensitive probe to a temporal variation of the gravitational constant.

Finally, a modification of the local gravitational constant does not only change the mass of a pulsar. In Section 3.5 we have discussed the change of the moment of inertia of a pulsar due to a local gravitational constant that depends on the distance to the companion, i.e., $G_{\mathrm{loc}}\left(r\right)$. For that reason, a pulsar on an eccentric orbit around its companion would periodically speed up and slow down in its rotation. As shown in Section 3.5, this can to first order be incorporated in the PK parameter ${\gamma }_{\mathrm{E}}$ of the Einstein delay. In this sense, every binary pulsar experiment that includes the PK parameter ${\gamma }_{\mathrm{E}}$ tests a spatial variation of the gravitational constant (cf. Section 4).

7.2. Strong-Field Generalizations of PPN Parameters

The parametrized post-Newtonian (PPN) formalism is one of the most successful frameworks in the interpretation of gravity experiments. In the standard PPN gauge, the framework contains 10 dimensionless PPN parameters in the metric components as coefficients of various gravitational potentials. These 10 parameters describe in a general way how metric theories of gravity can deviate from GR at the post-Newtonian level. The PPN parameters take different values in different gravity theories. We refer the reader to [3] and references therein for further details on the PPN formalism. Some of the PPN parameters are directly linked to a violation of LPI ($\xi$) or LLI (${\alpha }_1$, ${\alpha }_2$, ${\alpha }_3$). The PPN parameter ${\alpha }_3$ combines a violation of LLI with a violation of the conservation of the total momentum. Another PPN parameter of interest in the context of pulsar experiments is ${\zeta }_2$, which is also linked to a violation of the conservation of total momentum. For GR one has $\xi ={\alpha }_1={\alpha }_2={\alpha }_3={\zeta }_2=0$. A fairly up to date summary on experimental constraints for PPN parameters can be found in [3].

In the following, we quickly summarize the most important observable effects for pulsars related to the PPN parameters listed above. A non-vanishing Whitehead parameter $\xi$ leads to an anisotropy of gravitational interaction, induced by the gravitational field of the Galaxy. As a result of this, the spin of a solitary pulsar should precess around a direction that approximately coincides with the direction to the Galactic center [134,135]. The most important effect of a non-vanishing ${\alpha }_1$ is a polarization of the orbit (gravitational Stark effect) of a binary which moves with respect to a preferred frame of reference for the gravitational interaction (most naturally the rest frame of the cosmic microwave background) [136]. Similarly to $\xi$, a non-vanishing ${\alpha }_2$ leads to a spin precession, this time around the direction of motion of the pulsar with respect to the preferred frame [134]. In the presence of an ${\alpha }_3$, a pulsar experiences a self-acceleration, which is perpendicular to the pulsar spin and the motion of the pulsar with respect to the preferred frame. The most important consequence, for testing that parameter, is the polarization of a binary pulsar orbit [137]. Finally, ${\zeta }_2$ leads to a self-acceleration of an asymmetric ($m_{\mathrm{p}}\ne m_{\mathrm{c}}$) eccentric binary pulsar system in the direction of periastron [138].

For binary pulsars, the weak-field PPN parameters are expected to have modifications related to the strong gravitational self-field of a NS [136]. See for instance [139], who have explicitly calculated ${\widehat{\alpha }}_1$ and ${\widehat{\alpha }}_2$ (corresponding to ${\alpha }_1$ and ${\alpha }_2$ respectively in the weak field) for NSs in Einstein-Æther and khronometric gravity. Quite generally, one expects a dependence on the (fractional) gravitational binding energy $ϵ$ of the pulsar and its companion. In the absence of non-perturbative phenomena, one could think of an expansion in terms of ${ϵ}_{\mathrm{p}}$ and ${ϵ}_{\mathrm{c}}$. For a generalized strong-field PPN parameter $\widehat{\mathcal{X}}$ one then has

$$\widehat{\mathcal{X}}={\mathcal{X}}_{\mathrm{PPN}}+\sum _{i=\mathrm{p},\mathrm{c}}{A}_i{ϵ}_i+\sum _{i=\mathrm{p},\mathrm{c}}\sum _{j=\mathrm{p},\mathrm{c}}{A}_{ij}{ϵ}_i{ϵ}_j+\cdots ,$$

(40)

where $A_i$, $A_{ij}$, etc. are expected to depend on the specifics of the theory.

Table 1. Limits for strong-field PPN parameters from pulsar experiments. A “<” means that the limit applies to the absolute value. Some of these limits are significantly (orders of magnitude) tighter than limits on the corresponding (weak-field) PPN parameter obtained in the Solar System (see [3] for a comparison).

Parameter	Limit (95% C.L.)	Pulsar Experiment
$\widehat{\xi }$	<$3.9\times {10}^{-9}$	PSRs B1937+21 and J1744−1134 [135]
${\widehat{\alpha }}_1$	$-{0.4}_{-3.1}^{+3.7}\times {10}^{-5}$	PSR J1738+0333 [140]
${\widehat{\alpha }}_2$	<$1.6\times {10}^{-9}$	PSRs B1937+21 and J1744−1134 [141]
${\widehat{\alpha }}_3$	<$4\times {10}^{-20}$	PSR J1713+0747 [117]
${\widehat{\zeta }}_2$	<$1.3\times {10}^{-5}$	Combination of several pulsars [142]

lists the current pulsar limits on various strong-field PPN parameters. In the absence of a cancellation between the PPN term and the ${ϵ}_i$ dependent terms in Equation (40), one can consider the limits for the strong field generalizations also as a limit for the weak field PPN parameter. Furthermore, compared to weak-field tests, binary pulsar experiments are also sensitive to a violation of these gravitational symmetries that are limited to the strong self-fields of NSs. Combined with sufficiently accurate weak-field tests, pulsar limits can therefore be used to test for such strong-field specific deviations [104].

8. Future Outlook

In the previous sections, we have demonstrated that pulsars allow us to perform high-precision tests of aspects of strong-field gravity that are either unique or complement efforts with other experiments or observations. Even though we have mostly concentrated on binary pulsars, we note that also selected isolated pulsars help with this undertaking (see e.g., PSRs B1937+21 and J1744−1134 as listed in Table 1). Overall, it is therefore not surprising that gravity tests with pulsars remain a most important key science driver in pulsar astronomy. There are indeed many developments in the field that promise important progress, quantitatively in terms of improved measurements, but also qualitatively in testing new effects. In this section we will present some of those expected short and mid-term science goals. The list, naturally, has to be incomplete given the richness of that field of research. Especially since the history of pulsar astronomy demonstrates that some of the most important results came from discoveries that had not at all been anticipated.

8.1. Time and Sensitivity

As indicated in Section 2, the overall precision of pulsar gravity tests scales with the precision of the TOAs, which in turn scales with the sensitivity of our observations. At the same time, cadence, orbital phase coverage and, in particular, the length of the timing experiments—which in some cases spans several decades!—are crucial. For instance, our ability to measure the orbital period decay due to GW damping improves rapidly with the time span of our experiment, $T_{\mathrm{obs}}$. Given the same cadence in the observations, $\delta {\dot{P}}_{\mathrm{b}}\propto T_{\mathrm{obs}}^{-2.5}$ [13,14]. Hence, continued timing observations of the known systems, some of them introduced in the previous sections, will continue to improve the measurement precision.

The improvement due to longer time baselines is greatly supported by continuous advances in radio astronomical techniques. Especially, in the last decade, the new availability of fast digitizers, powerful real-time processing units, as well as affordable commodity computing power has allowed us to develop new broad band radio receivers. With significantly larger bandwidths, one does not only improve the sensitivity and hence timing precision, but one can also measure and mitigate more precisely the effects of the interstellar medium on the signal propagation described in Section 2.1.

New telescopes with bigger collecting area will also give an immediate improvement in the TOA precision. The distance based on the timing parallax, for instance, is one of those measurements particularly benefiting from such a development [143]. About two years ago the 500-m FAST radio telescope in the northern hemisphere (China) and the 64 × 13.9-m dishes of the MeerKAT interferometer in the southern hemisphere (South Africa) have started with their regular observations, both of which already have proven to provide superb timing data with promising early science results [144,145]. Over the next few years these two telescopes are expected to greatly benefit the field of gravity tests with pulsars. Of great importance here is the fact that the MeerKAT interferometer will soon be extended to MeerKAT+ in a collaboration between the South African Radio Astronomy (SARAO) and the Max-Planck Society (MPG)8 and eventually to SKA1. How this will improve our gravity tests has recently be demonstrated in [146] on the basis of extensive mock data simulations for the Double Pulsar. Many of the timing parameters of the Double Pulsar, including the timing parallax, are expected to improve significantly over the next 10 years. At some point (≳2030) they are in the reach of GW-related 3.5PN corrections in the equations of motion [147]. By then the Double Pulsar should allow, presumably for the first time in a binary pulsar, for a (moderately precise) test of the Lense-Thirring effect, via its contribution to the advance of periastron [31,148].9 In terms of improving GR tests, newly discovered, more relativistic systems, like PSRs J1946+2052 [151] (a more compact version of the Double Pulsar, but unfortunately not equally highly inclined) and J1757−1854 [152] (a more compact version of the Hulse-Taylor pulsar) are promising as well.

8.2. Complementary Information

Pulsar tests also often benefit from significant progress in other areas of astronomy providing important crucial information. For instance, at some point, for many of the relativistic systems the error in the distance measurement and the uncertainty in the Galactic potential will become the limiting factor in measuring the GW damping (cf. Section 4.1). However, better understanding of the Galactic gravitational potential is expected from new models based on GAIA data [153].

When it comes to dipolar radiation tests with pulsar-WD systems (Section 5), then in some of the cases, like PSRs J1738+0333 and J0348+0432, the limiting factor will be the uncertainty in the WD models restricting our capability of deriving the WD masses from high resolution spectroscopy. Improvements in our understanding of WDs are therefore also highly beneficial. However, pulsars like PSR J2222−0137 will not be limited by this, as there are additional PK parameters measured in that system, like $\dot{\omega }$ and the Shapiro parameters [102], whose errors will just get smaller with time and better telescopes. A promising pulsar in that context is PSR J1913+1102, which was discovered only a few years ago [154]. Although the companion is a NS, the high asymmetry in mass [155] also means an asymmetry in compactness and therefore in many alternatives to GR a significant gravitational dipole.

8.3. New Laboratories & New Avenues

New telescopes do not only promise an improvement in timing precision, they also come with new superior survey capabilities. This warrants high hopes for the discovery of new pulsars suitable for gravity tests. Apart from the discovery of more relativistic systems of the known kind, there is the potential of discovering systems which allow for qualitatively new tests. In the first place, there is the discovery of a pulsar-black hole (BH) system. Such a system could be used for precision tests of BH physics, complementary to other tests [156,157,158,159]. The ultimate system in this sense would be the discovery of a pulsar in a sufficiently relativistic orbit ($P_{\mathrm{b}}\lesssim 1\mathrm{yr}$) around Sgr A${}^{*}$, the supermassive BH in the Galactic center [160,161,162]. Gravity tests with such a pulsar could nicely complement the tests with the S-stars and the EHT [163]. Finding binary pulsars in the central region of our Galaxy, would already be of great interest for gravity tests, for instance to improve tests related to DM [124,164].

The most interesting binary systems to find are those that are most compact and that exhibit the largest accelerations of our pulsar clock. In those cases, search techniques developed over the last 50 years need to probe a highly dimensional parameter space that not only covers spin period and orbital parameters but also pulse width and the a priori unknown dispersion of the signal in the interstellar medium. Apart from benefiting from the steady increase in available computer power, astronomers have developed special techniques (“acceleration” and “jerk searches”, see [10]), have adopted algorithms successfully developed for GW detector experiments (see e.g., [165]) or have leveraged the power of Global Volunteer Computing like “Einstein@Home” [166]. Using these and other techniques, each search observation can result in a huge number of candidates, most of which will be due to random statistical fluctuations or the impact of man-made radio interference signals. Pulsar astronomers have therefore also developed techniques to sift these candidates for the real pulsar signals in an efficient manner. Without the introduction of techniques related to machine learning and artificial intelligence (e.g., [167]), it would nowadays be impossible to cope with the large data volumes that will increase even further in the era of SKA. As an example, the TRAPUM survey to discover pulsars and transients with MeerKAT [168] produces about 3 PB per observing night, which cannot be stored but needs to be searched in (quasi-) realtime. Another benefit of the larger collecting area of MeerKAT, FAST or finally SKA is that deeper sensitivity is often achieved with much smaller integration times, $T_{\mathrm{obs}}$. Hence, shorter parts of the orbits are covered during one pointing, which reduces the amount of acceleration space that needs to be searched enormously, since the number of acceleration steps scales as $T_{\mathrm{obs}}^3$. This alone is an important factor and gives reasons to believe that many more, and even more compact systems will be discovered in the near future.

Last but not least, pulsars can also be used in the detection and eventually study of ultra-low frequency GWs. Rather than serving as the GW emitters, they can form a Galaxy-sized GW detector in the ongoing effort of so-called pulsar timing array (PTA) experiments [169]. The most promising sources are binary supermassive BHs. Although detection is still pending, with new telescopes like the SKA one might eventually even be able to study the properties of these light-years long GWs in details, like their polarization modes and their dispersion relation, complementary to similar tests at high frequencies with ground-based GW observatories [170,171]. Prospects and details of PTA experiments (see e.g., [169]) are beyond the scope of this review.

9. Summary

In the last few decades since the discovery of the first binary pulsar, pulsar astronomy has provided some of the best experiments for gravitational physics, with GR having passed them all with flying colors. Binary pulsars provided the first verification of the existence of gravitational waves, having reached by now a precision of well below 0.1% in confirming the quadrupole formula of GR (Section 4). Quite generally, pulsars with NS companions could be utilized to test various aspects in the gravitational interaction of two strongly self-gravitating (material) objects, this includes the spin precession of pulsars due to relativistic spin-orbit coupling (Section 4.3) and the propagation of electromagnetic signals in the curved spacetime of a NS (Section 4.2). Pulsar WD systems have been used to constrain the existence of dipolar GWs, a prediction by many alternatives to GR (Section 5). Various pulsar experiments lead to tight constraints on potential violations of different building blocks of the SEP. There is the confirmation of the UFF of a strongly self-gravitating body within a hierarchical stellar triple system (Section 6.1), a test of the UFF towards DM (Section 6.2), and the verification of the constancy of the gravitational constant (Section 7.1), complementary to the weak-field limits from the Solar System. There are various pulsar experiments that lead to tight constraints on the strong-field generalizations of different PPN parameters, linked to violations of the gravitational LPI and LLI (two of the columns of the SEP) and the conservation of momentum (Section 7). In many ways, pulsar experiments are complementary to other gravity experiments, like the ones coming from the Solar System or from LIGO/Virgo merger events. The reason is multi-layered. In some cases it is because pulsars reach an unmatched precision, in other cases because they test specific effects only accessible in pulsar experiments. Pulsars certainly test a quite specific part of the parameter space of gravity (see Figure 7).

In the future we expect significant improvements in gravity tests with pulsars coming from two major directions, mostly driven by new facilities, like FAST, MeerKAT and the upcoming SKA, but also driven by new developments in field of big data and high-performance computing (Section 8). On the one hand there is the improvement in timing precision and the continuous extension of the timing baseline for known pulsars. On the other hand, there is the discovery of new, more relativistic or qualitatively different pulsar systems. Here, the discovery of a pulsar with a BH companion is on the top of the list of many pulsar astronomers, ideally a pulsar in orbit around Sgr A${}^{*}$, the supermassive BH in the Galactic center.