TIMING MEASUREMENTS OF THE RELATIVISTIC BINARY PULSAR PSR B1913+16

Pulsar PSR B1913+16 (PSR J1915+1606) was the first binary pulsar to be discovered (Hulse & Taylor 1975). Its 7.75 hr binary period and 300 km s⁻¹ projected orbital velocity suggested that it would exhibit a rich set of potentially measurable relativistic effects, and that has turned out to be the case. Subsequent studies have shown that the pulsar's companion must also be a neutron star. With tidal and rotational stellar distortions eliminated as complicating factors for any dynamical analysis, the system is an essentially clean laboratory for testing relativistic gravity. In this respect, the largest remaining complications depend on reference-frame accelerations related to the structure and dynamics of our Galaxy.

Over the past 35 years, we and our colleagues have timed PSR B1913+16 at Arecibo Observatory with steadily improving equipment and analysis techniques. Among the best known results are measurement of the general relativistic advance of periastron at a rate some 35,000 times that of Mercury in the solar system (Taylor et al. 1976); the effect of gravitational radiation damping, causing a measurable rate of orbital decay (Taylor et al. 1979); and the detection of changes in the pulse shape, resulting from geodetic spin precession (Weisberg et al. 1989). Our results continue to be fully consistent with general relativity and have placed strong constraints on alternative, previously viable, relativistic theories of gravity (Weisberg & Taylor 1981, 2002; Taylor & Weisberg 1982, 1989; Taylor et al. 1992; Clifton & Weisberg 2008). The purpose of this paper is to provide our latest analysis of the pulsar and its orbit, based on timing observations from 1974 through 2006.

Our highest quality data set consists of 7650 five-minute integrations of the pulse profile at wavelengths near 20 cm, each subsequently reduced to a topocentric time-of-arrival (TOA) measurement at Arecibo Observatory. These data span the years 1981 through 2006; earlier measurements are consistent with these but (owing to their much larger uncertainties) have little effect on the final results. In 2004 we made a special effort to observe regularly throughout the year, so as to improve the quality of the astrometric results. Parameters of the observing equipment and associated measurement uncertainties are presented in Table 1. Further details of the Princeton Mark I observing system are given by Taylor & Weisberg (1982), while the Mark II and III are described in Taylor & Weisberg (1989). The latest system, the Wideband Arecibo Pulsar Processor (WAPP), consists of four spectrometers optimized for pulsar work, each having 512 frequency channels across a 100 MHz bandwidth and synchronously accumulating the pulsar signal into 1024 phase bins across the full pulsar period.

Data were analyzed with the program TEMPO, using the JPL DE405 solar system ephemeris and the relativistic timing model of Damour & Deruelle (1986). The Damour–Deruelle approach is particularly desirable because it leaves open the question of the correct relativistic theory of gravity: one solves for phenomenological parameters whose precise relation to a gravitational theory can be investigated later. Further discussion of the Damour–Deruelle model and its application to binary pulsars in TEMPO can be found in Taylor & Weisberg (1989). Important quantities determined in this way include those related to the pulsar's celestial position, spin, and orbit. The first category includes right ascension α, declination δ, proper motions μ_α and μ_δ, epoch t₀, pulse repetition frequency f, spin-down rate $\dot{f}$, and a glitch epoch and frequency discontinuity Δf. Fitted orbital parameters include five Keplerian quantities: the projected semimajor axis of the pulsar orbit x ≡ a₁sin i, eccentricity e, epoch of periastron passage T₀, period P_b, and longitude of periastron ω₀; and the following relativistic or "post-Keplerian" parameters: average rate of periastron advance $\langle \dot{\omega }\rangle$, variations in gravitational redshift and time dilation γ, and orbital period derivative $\dot{P_b}$. As with most other pulsars that have been timed carefully over several decades, a number of "nuisance parameters" must also be measured to account for unmodeled long-term timing irregularities. These extra fitted terms have no clear physical interpretation beyond being somehow related to the poorly understood structure and dynamics of the spinning neutron star. Their mathematical form is somewhat arbitrary; as described further below, we have experimented with a number of higher-order time derivatives of the pulsar rotation frequency and one additional frequency discontinuity.

3.1. Astrometric and Spin Parameters

The astrometric and basic spin parameters determined from the full data set are listed in Table 2. A significant result for proper motion has been obtained for the first time. While the pulsar position is determined from the annual variation of TOAs as the Earth moves about its orbit, proper motion is determined from a secular variation of this annual signal. Timing noise or other systematic effects can contaminate such measurements. We now believe that the previously reported value of PSR B1913+16's proper motion (Taylor & Weisberg 1989) was biased by timing noise, a problem which was exacerbated by the concentrated but ∼biennially spaced observing campaigns frequently employed for this source. To circumvent these problems, we observed PSR B1913+16 several times over the course of the calendar year 2004 to achieve thorough data coverage around Earth's orbit. By merging the 2004 data with observations made in 1985–1988, which also had good coverage throughout those years, we have now obtained a robust measurement of the pulsar's proper motion. See Section 5 for further analyses of the proper motion result.

As in other pulsars, the measured spin-down rate $\dot{f}$ is assumed to be the result of an electromagnetic braking torque on the spinning, strongly magnetized neutron star. Its value is deterministic and (at the level of accuracy quoted in Table 2) independent of the span of data over which it is fitted. In addition, the pulsar experienced a well-defined "classical" glitch in 2003 May (see Figures 1(a) and 1(b)). The magnitude Δf/f = 3.7 × 10⁻¹¹ of the event is smaller than almost all seen in the population of normal (not recycled) pulsars, with only some of the glitches in PSR B0355+54 having a comparably small magnitude (Melatos et al. 2008; Chukwude & Urama 2010). Ours is only the second glitch to be detected in a recycled pulsar, with the other (the smallest known among all pulsars) occurring in the millisecond pulsar PSR B1821−24 in globular cluster M28 (Cognard & Backer 2004).

Table 2. Astrometric and Spin Parameters

Parameter	Valuea
t0 (MJD)b	52984.0
α (J2000)	19h15m2799928(9)
δ (J2000)	16°06'273871(13)
μα (mas yr−1)	−1.43(13)
μδ (mas yr−1)	−0.70(13)
f (s−1)	16.94053778563(15)
$\dot{f}$ (s−2)	−2.4761(9) ×10−15
Glitch epoch (MJD)	52770(20)
Δf (s−1)	6.2(2) ×10−10

Notes. ^aFigures in parentheses represent estimated uncertainties in the last quoted digit. The estimated uncertainties range from (3–10)× the formal fitted uncertainties, in order to also reflect variations resulting from different assumptions regarding timing noise, etc. ^bThis quantity is the epoch of the next six measurements tabulated here.

Download table as:  ASCIITypeset image

The remaining timing parameters—higher-order derivatives $\ddot{f}$, $\ddot{f}$, ..., and (in some fits) another small frequency discontinuity in mid-1992—were introduced to the timing solution in an ad hoc manner, in order to "whiten" the remaining post-fit residuals. Although we offer no clear or unique physical interpretation for these parameters, their combined effects are almost certainly a consequence of stochastic timing-noise processes in the neutron star interior (Cordes & Downs 1985; Arzoumanian et al. 1994; Urama et al. 2006). Unlike the case of the well-fitted 2003 May glitch, the values of these fitted parameters are not independent of the data span analyzed and cannot be expected to extrapolate the timing behavior accurately to future epochs. Identifying the mid-1992 behavior as a discrete event is highly uncertain, in part because of coarse sampling around that time. Our data can be fit almost as well by introducing several additional frequency derivatives instead of a second discrete event (see Figure 1(c) for post fit residuals). Similar results were obtained by fitting multiple harmonically related sinusoids to the timing noise with the routine "FITWAVES" (Hobbs et al. 2004) of program TEMPO2 (Hobbs et al. 2006; Edwards et al. 2006), instead of the multiple frequency derivatives described above. Since we are unable to converge on a unique glitch parameter solution for the mid-1992 behavior, we do not include this event in Table 2.

Zoom In Zoom Out Reset image size

3.2. Orbital Parameters

Fitted orbital parameters for PSR B1913+16 are listed in Table 3. As we have emphasized before (Taylor & Weisberg 1989), values for each of the Damour–Deruelle post-Keplerian parameters expected in general relativity can be expressed in terms of the Keplerian parameters and the initially unknown masses of the pulsar and its companion, m₁ and m₂. The appropriate expressions for $\langle \dot{\omega }\rangle$ and γ are

$$\begin{eqnarray} \langle \dot{\omega }\rangle &=& 3 G^{2/3} c^{-2} (P_b 2 \pi)^{-5/3} (1-e^2)^{-1} (m_1 + m_2)^{2/3} \nonumber \\ &=& 2.113323(2)\ \left[\frac{(m_1 + m_2)}{M_\odot }\right]^{2/3} \;{\rm deg\;yr}^{-1}, \\ \gamma &=& G^{2/3} c^{-2} e (P_b / 2 \pi)^{1/3} m_2 (m_1 + 2 m_2) (m_1 + m_2)^{-4/3} \nonumber \\ &=& 0.002936679(2) \!\left[\frac{m_2 (m_1 + 2 m_2) (m_1 + m_2)^{-4/3}}{M_\odot ^{2/3}}\!\right] \;{\rm s}.\nonumber\\ \end{eqnarray} \tag{ 1 }$$

In the second line of each equation we have substituted values for P_b and e from Table 3, and used the constants GM_☉/c³ = 4.925490947 × 10⁻⁶ s and 1 Julian yr =86400 × 365.25 s. The figures in parentheses represent uncertainties in the last quoted digit, determined by propagating the uncertainties listed in Table 3. In each case, the uncertainties are dominated by the experimental uncertainty in orbital eccentricity, e.

Table 3. Orbital Parameters

Parameter	Valuea
T0 (MJD)	52144.90097841(4)
x ≡ a1sin i (s)	2.341782(3)
e	0.6171334(5)
Pb (d)	0.322997448911(4)
ω0 (deg)	292.54472(6)
$\langle \dot{\omega }\rangle$ (deg yr−1)	4.226598(5)
γ (ms)	4.2992(8)
$\dot{P}_b$	−2.423(1) ×10−12

Note. ^aFigures in parentheses represent estimated uncertainties in the last quoted digit. The estimated uncertainties range from (2–6)× the formal fitted uncertainties, in order to reflect also the variations resulting from different assumptions regarding timing noise, etc.

Download table as:  ASCIITypeset image

Equation (1) may be solved for the total mass of the PSR B1913+16 system, yielding M = m₁ + m₂ = 2.828378 ± 0.000007 M_☉. The additional constraint provided by Equation (2) permits a solution for each star's mass individually, m₁ = 1.4398 ± 0.0002 M_☉ and m₂ = 1.3886 ± 0.0002 M_☉. As far as we know, these are the most accurately determined stellar masses outside the solar system. It is interesting to note that since the value of Newton's constant G is known to a fractional accuracy of only 1 × 10⁻⁴, M can be expressed more accurately in solar masses than in grams.

3.3. Gravitational Radiation Damping

According to general relativity a binary star system should radiate energy in the form of gravitational waves. Peters & Matthews (1963) showed that the resulting rate of change in orbital period should be

$$\begin{eqnarray} \dot{P}_b^{\rm GR} &=&-\!\frac{192\,\pi \,G^{5/3} }{5\,c^5} \left(\frac{P_b}{2\pi }\right)^{-5/3} \left(1 + \frac{73}{24} e^2 + \frac{37}{96} e^4\right) \nonumber\\ &&\times(1-e^2)^{-7/2} \ m_1\,m_2\,(m_1+m_2)^{-1/3} \nonumber \\ &=& -\!1.699451(8)\times 10^{-12} \left[\frac{m_1 m_2 (m_1+m_2)^{-1/3}}{M_\odot ^{5/3}}\right].\nonumber\\ \end{eqnarray} \tag{ 3 }$$

Inserting values obtained for m₁ and m₂ and propagating uncertainties appropriately, we obtain the general relativistic predicted value

$$\begin{equation} \dot{P}_b^{\rm GR} = -2.402531\pm 0.000014 \times 10^{-12}. \end{equation} \tag{ 4 }$$

Equations (3) and (4) apply in the orbiting system's reference frame. Relative acceleration of that frame with respect to the solar system barycenter will cause a small additional contribution to the observed $\dot{P}_b$. Damour & Taylor (1991) presented a detailed discussion of this effect and other possible contributions to $\dot{P}_b$. Recent progress in determining the galactic-structure parameters allows us to update the relevant quantities and compute a new value for the kinematic correction to $\dot{P}_b$. Using R₀ = 8.4 ± 0.6 kpc for the distance to the galactic center and Θ₀ = 254 ± 16 km s⁻¹ for the circular velocity of the local standard of rest (Ghez et al. 2008; Gillessen et al. 2009; Reid et al. 2009), and d = 9.9 ± 3.1 kpc for the pulsar distance (Weisberg et al. 2008), we obtain the kinematic contribution, $\Delta \dot{P}_{\rm b, gal}$:

$$\begin{equation} \Delta \dot{P}_{\rm b, gal}=-0.027\pm 0.005 \times 10^{-12}. \end{equation} \tag{ 5 }$$

Thus, we find the ratio of the observed-to-predicted rate of orbital period decay to be

$$\begin{equation} \frac{\dot{P}_b - \Delta \dot{P}_{\rm b, gal}}{\dot{P}_b^{\rm GR}} = 0.997\pm 0.002. \end{equation} \tag{ 6 }$$

Agreement between the observed orbital decay and the general relativistic prediction is illustrated in Figure 2, which shows how excess orbital phase (relative to an unchanging orbit) has accumulated since the pulsar's discovery in 1974. We note that the overall experimental uncertainty embodied in Equation (6) is now dominated by uncertainties in the galactic parameters and pulsar distance, not the pulsar timing measurements. Even better agreement between the observed and expected values of $\dot{P}_b$ would be obtained if the true value of R₀ or d were slightly smaller, or Θ₀ slightly larger. For example, observed and expected values agree if d = 6.9 kpc, which is within the Weisberg et al. (2008) error envelope. It will be interesting to see whether improved future estimates of these quantities will show one or more of these conditions to be true.

Zoom In Zoom Out Reset image size

Two other relativistic phenomena are potentially measurable in the PSR B1913+16 system: geodetic precession and gravitational propagation delay. Spin–orbit coupling should cause the pulsar's spin axis to precess (Damour & Ruffini 1974; Barker & O'Connell 1975a, 1975b), which should lead to observable pulse shape changes. Weisberg et al. (1989) first detected such changes, which were observed and modeled further by Kramer (1998). Weisberg & Taylor (2002) and Clifton & Weisberg (2008) found that the pulsar beam is elongated in the latitude direction and becomes wider in longitude with increasing distance from the beam axis in latitude. These models suggest that in the next decade or so, precession may move the beam far enough that the pulsar will become unobservable from Earth for some decades, before eventually returning to view.

The excess propagation delay (Shapiro 1964) caused by the passage of pulsar signals through the curved spacetime of the companion is largest at the pulsar's superior conjunction. The maximum amplitude varies with time because the impact parameter at superior conjunction strongly depends on the current value of ω. In this respect, the orbital geometry was particularly unfavorable in the mid-1990s (see Damour & Taylor 1992), but in coming years the propagation delay should start to become observable. Damour & Deruelle (1986) characterize the measurable quantities as range r = (Gm₂/c³) and shape s ≡ sin i of the Shapiro delay, where i is the orbital inclination. As orbital precession carries our line of sight deeper into the companion's gravitational well, future observations should permit the robust measurement of these two parameters, and hence two additional tests of relativistic theories of gravity (Damour 2009; Esposito-Farese 2009).

Our pulsar's proper motion measurement (Section 3.1), combined with the distance estimate discussed in Section 3.3, corresponds to a transverse velocity (with respect to the solar system barycenter) of 75 km s⁻¹ with a galactic position angle of 306°, i.e., directed 36° above the galactic plane. The ∼30% distance uncertainty places similar limits on velocity accuracies.

We can now estimate two components of the pulsar systemic velocity in its own standard of rest by combining the measured pulsar transverse velocity and distance, the solar motion with respect to our local standard of rest (Schönrich et al. 2010), and galactic quantities R₀ and Θ₀. The third component of motion, which is inaccessible via proper motion measurements, lies close to the direction of Galactic rotation at the pulsar's position.

The pulsar's galactic planar and polar velocity components relative to its standard of rest are 247 km s⁻¹ almost directly away from the galactic center and 51 km s⁻¹ toward the galactic North Pole, respectively. (This is significantly larger than the measured velocity in the solar system barycenter frame because the pulsar's standard of rest velocity fortuitously cancels much of the pulsar's peculiar velocity with respect to it.) The systemic velocity of B1913+16 is significantly larger than other well-measured double neutron star binary system velocities, including the J0737−3037 (transverse velocity 10 km s⁻¹; Stairs et al. 2006), J1518+4905 (transverse velocity 25 km s⁻¹; Janssen et al. 2008), and B1534+12 (transverse velocity 122 km s⁻¹; Thorsett et al. 2005) systems.

We have analyzed the full set of Arecibo timing data on pulsar B1913+16 to derive the best values of all measurable quantities. A significant proper motion has finally been determined. A small glitch was observed in the pulsar's timing behavior, the second known glitch in the population of recycled pulsars. The measured rate of orbital period decay continues to be almost precisely the value predicted by general relativity, providing conclusive evidence for the existence of gravitational radiation. Uncertainties in galactic accelerations now dominate the error budget in $\dot{P}_b$ and are likely to do so until the pulsar distance can be measured more accurately. We expect that the Shapiro gravitational propagation delay will yield additional tests of relativistic gravity within a few more years.

The three authors gratefully acknowledge financial support from the US National Science Foundation. Arecibo Observatory is operated by Cornell University under cooperative agreement with the NSF. We thank Joseph Swiggum for assistance with analyses of glitches in the pulsar population; and C. M. Ewers,⁴ A. de la Fuente, J. T. Green, and Z. Pei for assistance with observations.