Relativistic Measurements from Timing the Binary Pulsar PSR B1913+16
AAccepted by APJ 2016 Jun 1
Preprint typeset using L A TEX style emulateapj v. 5/2/11
RELATIVISTIC MEASUREMENTS FROM TIMING THE BINARY PULSARPSR B1913+16
J. M. Weisberg and Y. Huang
Department of Physics and Astronomy, Carleton College, Northfield, MN 55057
Draft version June 10, 2016
ABSTRACTWe present relativistic analyses of 9257 measurements of times-of-arrival from the first binary pul-sar, PSR B1913+16, acquired over the last thirty-five years. The determination of the “Keplerian”orbital elements plus two relativistic terms completely characterizes the binary system, aside froman unknown rotation about the line of sight; leading to a determination of the masses of the pulsarand its companion: 1 . ± .
001 M (cid:12) and 1 . ± .
001 M (cid:12) , respectively. In addition, the completesystem characterization allows the creation of tests of relativistic gravitation by comparing measuredand predicted sizes of various relativistic phenomena. We find that the ratio of observed orbital perioddecrease due to gravitational wave damping (corrected by a kinematic term) to the general relativisticprediction, is 0 . ± . δ θ , although its intrinsicvalue is obscured by currently unquantified pulsar emission beam aberration. We have also marginallymeasured the time derivative of the projected semimajor axis, which, when improved in combinationwith beam aberration modelling from geodetic precession observations, should ultimately constrainthe pulsar’s moment of inertia. Subject headings: binaries: close — gravitation — gravitational waves — pulsars: individual (PSRB1913+16) INTRODUCTION
Pulsar B1913+16 was the first binary pulsar discov-ered (Hulse & Taylor 1975). The system consists oftwo neutron stars (one an observed pulsar) orbiting ina very tight, highly eccentric orbit, and it remains oneof the best for studying relativistic gravitation (Weis-berg & Taylor 1981; Taylor & Weisberg 1982, 1989; andWeisberg, Nice, & Taylor 2010, hereafter WNT). In thispaper, we update WNT with the addition of post-2006data and with further relativistic timing analyses. Theaddition of significant quantities of data acquired withmodern data-acquisition devices has enabled us to mea-sure several additional relativistic phenomena for thefirst time in this system, while also refining previously-measured ones. Among the parameters newly measuredwith various degrees of accuracy are the Shapiro gravita-tional propagation delay, a relativistic correction to theelliptical orbital shape, and the time derivative of pro-jected pulsar semimajor axis. All of the the data usedin this study are published with this paper, while ouranalysis software is published in an online repository.We describe the observations used in this work in § § §
4, and their applications fortests of relativistic gravitation are discussed in §
5. Weconclude in § [email protected], orcid https://orcid.org/0000-0001-9096-6543 DATA
The data for our analyses consists of 9257 pulse timesof arrival (TOAs) derived from five-minute integrationsof the pulsar signal at frequencies near 1400 MHz mea-sured at Arecibo Observatory from 1981 - 2012. The pa-rameters of the various observing systems and the num-ber of TOAs from each through epoch 2006 are tabu-lated in WNT; WAPP spectrometer observations sincethen have added another 1652 TOAs to the total, eachacquired by three WAPPs deployed simultaneously at ap-proximately contiguous 100-MHz bands near 1400 MHz.Geodetic precession of the pulsar spin axis has inducedpulse profile changes (Weisberg et al. 1989; Kramer 1998;Weisberg & Taylor 2002; Clifton & Weisberg 2008) thathave lately grown increasingly larger, presumably as ourline of sight approaches the edge of the pulsar beam.Nevertheless, for purposes of uniformity, we use only asingle profile template while finding TOAs for all WAPPdata. This procedure induces time offsets into our TOAdataset between different sessions and frequencies, whichhave grown to a level where they should be compensatedfor. To do so, we adopted the following process. First,we formed a pulse profile at each frequency band for eachtwo-week session. Each resulting “session - band” stan-dard profile has much greater
S/N than does a singlefive - minute integration, while still being short enoughto avoid the secular changes we are trying to measure.Next, we measured the offset of the midpoint of this ses-sion - band standard profile with respect to the grandstandard profile. (The midpoint is assumed to corre-spond to a fixed longitude on the pulsar regardless of a r X i v : . [ a s t r o - ph . H E ] J un J. M. Weisberg and Y. Huangprofile shifts).Then we fitted out a “primary” linear model of theprofile offsets as a function of time at each band. In thisfashion, we provided an empirically determined, modelindependent, first-order TOA correction that accountsfor the effects of profile changes, thereby minimizing thelong-term effects of profile shifts that might be mistakenas the signature of other phenomena. We next fittedthe timing model to all such “primarily-offsetted” TOAsand chose the resulting dispersion measure as our nom-inal value. Finally, secondary offsets were determinedthrough single session fits, where the dispersion measurewas fixed at its nominal value and the residual offset ofeach band was fitted for and then removed. This processensures that the infinite-frequency TOAs calculated fromeach band and session are self-consistently de-dispersed and offsetted.To verify that the above profile-variation-correctionprocess does not contaminate our parameter measure-ments, we also employed an alternate approach, fittingfor an offset for each band in each session simultaneouslywith all other parameters (Demorest et al. 2013). Thisprocedure yields parameters that agree with our methodto within 1 σ for all parameters, suggesting that our mea-surements are robust with respect to the methods usedto remove profile-shift-induced timing offsets. RELATIVISTIC ANALYSIS OF TOAS
Using an augmented version of the TEMPO soft-ware program, we fitted the relativistic timing modelof Damour & Deruelle (1986, hereafter DD); or, in cer-tain cases, the DD model augmented by the Freire &Wex (2010, hereafter FW) Shapiro parametrization (see § to our TOAs. In these models, the pulsar signalencounters several distinct types of delays on its journeyfrom the orbiting pulsar to the solar system barycenter,such that the infinite-frequency pulse arrival time at thesolar system barycenter, t ssbc , is given by: t ssbc = D − [ T + ∆ Roemer ( T ) + ∆ Einstein ( T )+∆ Shapiro ( T ) + ∆ Aberration ( T )]; (1)where each delay is a function of the pulsar proper time ofpulse emission, T , and whose details depend on a numberof physical parameters that are fitted for. (The Dopplerfactor D accounts for the relative motion of the solar sys-tem and binary system barycenters.) The various termsin Eq. 1 are detailed in DD and Damour & Taylor (1992,hereafter DT92), and we will comment further on the lasttwo terms of Eq. 1 in the following two sections.Among the fitted parameters, we determined improvedvalues of the pulsar spin and orbital parameters thatwere published in WNT, plus a number of new ones. Forthe first time, we have successfully fitted for the Shapiro(1964) gravitational propagation delay while also placingconstraints on two additional ones: a relativistic correc-tion to the quasi-elliptical shape of the orbit and the This procedure also absorbs TOA variations induced by DMfluctuations at the levels and timescales expected from studies ofmillisecond pulsars (You et al. 2007). See http://sourceforge.net/projects/TEMPO/ for our aug-mented version of TEMPO, which contains our fitting routine forthe FW Shapiro parameters in high-eccentricity binaries. shrinkage rate of its projected semimajor axis, as de-scribed further below.The improved, previously fitted parameters include thepulsar spin frequency and derivative(s) f , ˙ f , . . . ; five“Keplerian” orbital elements defined as projected pulsarsemimajor axis x ≡ a sin i where i is the orbital incli-nation, orbital period P b , eccentricity e , reference epoch T , the reference epoch’s longitude of periastron ω ; rel-ativistic “post-Keplerian” parameters defined as meanrate of periastron advance (cid:104) ˙ ω (cid:105) , gravitational redshift -time dilation term γ , and orbital period derivative ˙ P b .The newly fitted post-Keplerian parameters includethe following: (i) two Shapiro delay terms called shape s and range r in the DD parametrization, or two (different)quantities ς and h in the alternate FW parametrization;(ii) the orbital elliptical shape correction parameter δ obs θ (to our knowledge never previously fitted for in any bi-nary system), which appears in DT92’s full expression forEq. 1’s Roemer term ; and (iii) ˙ x obs and ˙ e obs , the ob-served time derivatives of x and e . All of the new param-eters are discussed in greater detail in § § §
4, and relativistic tests resulting from these measure-ments are discussed in §
5. A set of TEMPO input filescontaining input parameters and the TOAs is availableas an online companion to this paper, on arXiv, and atzenodo: http://dx.doi.org/10.5281/zenodo.54764.
Fitting for the Shapiro Gravitational PropagationDelay via two Different Parametrizations
Until this work, we have been unable to explicitlymeasure the two Shapiro gravitational propagation de-lay terms that characterize Eq. 1’s ∆
Shapiro ( T ), owingto their relatively small timing signature and to their co-variance. The two terms are identified as s and r in DD;while FW recently developed an alternate parametriza-tion of the phenomenon, wherein their two fitted parame-ters, ς and h , are orthonormal. Our software implemen-tation of their parametrization for high-eccentricity pul-sars is available in our augmented version of TEMPO. The measurements of either pair of Shapiro param-eters, ( s, r ) or ( ς, h ), can be utilized in either of twodifferent manners, as described below.First, if general relativity is assumed to be the correcttheory of gravitation, then either pair of Shapiro mea-surables can be utilized as independent constraints onthe orbital inclination and binary companion mass. Wesummarize the theory here, and then apply it in § s and r translate directly intosin i and m (the companion mass), respectively:sin i = s, (2)and m = (cid:18) c G (cid:19) r = (cid:18) r T (cid:12) (cid:19) M (cid:12) , (3)with c the speed of light, G the Newtonian gravitationalconstant, T (cid:12) = G M (cid:12) /c = 4 .
925 490 947 µs . Before it can be utilized for tests of relativity, the δ obs θ param-eter must be corrected for a comparable aberration term which iscurrently undeterminable (see § iming Relativistic Binary Pulsar PSR B1913+16 3The alternate FW parametrization of the Shapiro delaygives sin i = 2 ςς + 1 , (4)while m is a combination of the two measurables ( h , ς ): m = (cid:18) c (cid:19) h ς = (cid:18) h /ς T (cid:12) (cid:19) M (cid:12) . (5)Alternatively, each measured parameter of the Shapiropair can be considered to be an independent test of rel-ativistic gravitation. We apply this procedure in § Determination of the Relativistic Orbital ShapeCorrection δ θ in the Presence of the AberrationDelay In order to successfully measure the intrinsic value of δ θ , which nominally quantifies a relativistic correctionto the shape of the approximately elliptical orbit in Eq.1’s Roemer delay expression, one must compensate theobserved value for the influence of a phenomenon thatcomparably affects TOAs, namely the orbital-phase de-pendent aberration of the pulsar beam, as described byDD and DT92. Those authors provide a prescription forcalculating and eliminating the confounding aberrationsignature from the observed value of δ θ , if the aberrationgeometry is known. In principle, the necessary infor-mation can be gleaned from studies of profile changesresulting from geodetic precession of the pulsar spin axis(Weisberg & Taylor 2002; Clifton & Weisberg 2008). Inthis section, we summarize the theoretical expressions re-quired to quantify δ θ and aberration, and we will applyour obervations to these results in § Aberration in Eq. 1, resulting fromaberration of the rotating pulsar beam, is dependent onthe time-variable transverse component of the pulsar’sorbital velocity. DD and DT92 parametrize the instan-taneous delay via the aberration parameters A ( t ) and B ( t ):∆ Aberration = A ( t ) { sin( ω + A e ( u )) + e sin ω } + B ( t ) { cos( ω + A e ( u )) + e cos ω } , (6)where A e ( u ) is a true-anomaly-like quantity, and A ( t )and B ( t ) are dependent on the precessing spin-axis ge-ometry: A ( t ) = − f − P b x sin i (1 − e ) / sin η sin λ , (7) B ( t ) = − f − P b x sin i (1 − e ) / cos i cos η sin λ , (8)with λ and η the geodetically-precessing polar angles ofthe pulsar spin axis with respect to the line of sight andline of nodes, respectively. (DD and DT92 suggestedthe substitution of a single fixed parameter, A , for thetwo parameters A ( t ) and B ( t ), because observations atthe time suggested that the spin and orbital angular mo-menta are aligned. However, subsequent observationsof pulse profile changes have shown that this is not thecase.)While the above equations provide a complete de-scription of the calculation of ∆ Aberration at any proper emission time, DD and DT92 also provide an alternateapproach that focuses on aberration parameters thatchange slowly (on precession timescales) as a result ofspin-axis precession. This procedure, detailed below, ismore closely tailored to parameters determinable fromTOA analyses.DT92 show that aberration will bias the observed valueof the relativistic orbital shape parameter δ obs θ with re-spect to its intrinsic value δ intr θ : δ obs θ = δ intr θ − (cid:15) A , (9)where the small parameter (cid:15) A is defined as (cid:15) A ≡ A ( t ) x . (10)Hence the observational bias can be removed, givenmeasurements of the aberration parameter A ( t ) and theKeplerian quantity x . The corrected value δ intr θ couldthen serve as an additional test of gravitation theory. Other Parameters Affected by the Aberration Delay
In addition to affecting δ θ , DT92 show that aberrationalso affects the observed x and e values. However, thefractional corrections to x and e are tiny. More inter-esting is the effect of geodetic spin-axis precession on thetime-derivatives of these parameters, because DT92 showthat they are potentially measurable. The precessionalmotion will cause the aberration geometry to change, re-sulting in secular changes to A ( t ) and hence to (cid:15) A onprecession timescales: d(cid:15) A dt = 1 x dA ( t ) dt = − f − P b i (1 − e ) / Ω geodetic1 sin λ × (sin i cos λ sin 2 η + cos i sin λ cos η ) , (11)where 2 π/ Ω geodetic1 is the geodetic precession period ofthe pulsar spin axis. The observed, normalized timederivative of e results from the sum of two phenomena: (cid:18) ˙ ee (cid:19) obs = d(cid:15) A dt + (cid:18) ˙ ee (cid:19) GW , (12)where “GW” designates effects due to gravitationalwaves; while the observed, normalized time derivativeof x stems from a combination of five terms: (cid:18) ˙ xx (cid:19) obs = d(cid:15) A dt + (cid:18) ˙ a a (cid:19) GW + (cid:18) cot i didt (cid:19) SO − µ cot i sin(Θ µ − Ω) − ˙ DD , (13)where “SO” refers to spin-orbit coupling. (See Lorimer &Kramer (2004) for an expression that includes additionalterms needed for some other binary pulsars.)The quantity di/dt in the third term of Eq. 13, result-ing from pulsar spin-orbit coupling, is developed here J. M. Weisberg and Y. Huangfrom expressions in DT92: didt = G c (cid:26) m m (cid:27) S a (1 − e ) / sin λ cos η = 1 c (cid:26) m m (cid:27) I (2 πf )(2 π/P b ) ( m + m )(1 − e ) / sin λ cos η, (14)where m is the pulsar mass, a R is the semimajor axisof the relative orbit, S = I (2 πf ) is the magnitude ofthe pulsar spin angular momentum, and I is its mo-ment of inertia. The fourth term of Eq. 13 results fromthe changing projection of the line of sight onto the or-bital plane due to proper motion, with µ and Θ µ respec-tively the amplitude and position angle of proper motionand Ω the position angle of the line of nodes (Kopeikin1996). The final term of Eq. 13, involving changes in theDoppler factor D of Eq. 1, is caused by the relative line-of-sight galactic accelerations of the solar system and thebinary system.The above equations demonstrate that measurementsof ˙ e or ˙ x , along with experimental or theoretical determi-nations of some of the other quantities appearing therein,can usefully constrain yet others. RESULTS OF THE FITS
We fitted the parameters discussed above to the fullset of TOAs, using the TEMPO software, as modifiedby us. See Tables 1 and 2 for our results and their es-timated uncertainties. The uncertainties quoted thereinrepresent the standard errors from the TEMPO fit (ex-cept as noted). This convention differs from our previouspractice, wherein many uncertainties were instead esti-mated from fitted parameter variations across multiplereasonable fits. While the old procedure facilitated theincorporation of some systematic uncertainties into theerror budget; the more stable recent instrumental config-urations appear to minimize such effects.Some of the fitted parameters shifted by several σ withrespect to the values reported in Weisberg et al. (2010).The shifts can all be attributed to the new incorporationof a frequency and time offset for each WAPP observ-ing session and center frequency in order to account forgeodetic-precession-induced profile changes (see § x . The latter procedure also led toa significantly larger uncertainty in the fitted value of γ and in quantities derived therefrom.The astrometric and spin solutions are listed in Table1. These are quite similar to those given in Weisberget al. (2010), except that our longer post-glitch base-line made it clear that the previously discovered glitch atMJD ≈ derivative, ∆ ˙ f . There remainsonly one known glitch having a significantly smaller valueof ∆ f /f [in globular cluster millisecond PSR B1821-24(Mandal et al. 2009)], although several of magnitude sim-ilar to the one tabulated here are now known. [See theonline Jodrell Bank Pulsar Glitch Catalogue (Espinozaet al 2011)]. Note that, as with Weisberg et al. (2010),ten higher-order spin derivatives were also fitted for to Table 1
Astrometric and Spin ParametersParameter Value a t (MJD) b . . . . . . . . . . . . 52984.0 α (J2000) . . . . . . . . . . . . . 19 h m . s δ (J2000) . . . . . . . . . . . . . 16 ◦ (cid:48) . (cid:48)(cid:48) µ α (mas yr − ) . . . . . . . . − µ δ (mas yr − ) . . . . . . . . − f (s − ) . . . . . . . . . . . . . . . 16.940537785677(3)˙ f (s − ) . . . . . . . . . . . . . . . − × − Glitch ParametersGlitch epoch (MJD). . . 52777(2)∆ f (s − ) . . . . . . . . . . . . . 5.49(3) × − ∆ ˙ f (s − ) . . . . . . . . . . . . . . − × − Figures in parentheses represent formal TEMPOstandard errors in the last quoted digit, except forthe glitch parameters. The stated uncertainty inglitch epoch results from empirically varying theglitch epoch until ∆ χ corresponds to the 68%confidence level; the quoted uncertainties in theother glitch parameters were derived from theirvariations as the glitch epoch was varied over thechosen range. b This quantity is the epoch of the next six mea-surements tabulated here.
Table 2
Orbital ParametersParameter Value a T (MJD) . . . . . . . 52144.90097849(3) x ≡ a sin i (s). . . 2.341776(2) e . . . . . . . . . . . . . . . 0.6171340(4) P b (d). . . . . . . . . . . 0.322997448918(3) ω (deg) . . . . . . . . 292.54450(8) (cid:104) ˙ ω (cid:105) (deg / yr) . . . 4.226585(4) γ (ms) . . . . . . . . . . 0.004307(4)˙ P obs b . . . . . . . . . . . . − × − δ obs θ . . . . . . . . . . . . . 4.0(25) × − ˙ x obs . . . . . . . . . . . . . − × − ˙ e obs (s − ) . . . . . . . 0.0006(7) × − Shapiro Gravitational Propagation Delay ParametersDamour & Deruelle (1986) Parametrization s . . . . . . . . . . . . . . . 0 . +0 . − . r ( µ s) . . . . . . . . . . 9 . +2 . − . Freire & Wex (2010) Parametrization ς . . . . . . . . . . . . . . . 0.38(4) h . . . . . . . . . . . . . . 0.6(1) × − Figures in parentheses represent formal TEMPO stan-dard errors in the last quoted digit. The DD Shapiro pa-rameters s and r , which are highly covariant, and theiruncertainties, were refined through a process illustratedin Fig. 2. eliminate the effects of timing noise. Their values arenot shown in the Table as they do not correspond tomeaningful physical parameters.Table 2 displays the results of our fit to orbital param-eters, including the eight final entries, which are fittedhere for the first time in this system. Note that the firstiming Relativistic Binary Pulsar PSR B1913+16 5 Figure 1.
The Shapiro gravitational propagation delay variationaround the orbit at three epochs. The curve represents the ex-pected delay based on a general relativistic calculation, while thepoints and their error bars result from combining all residuals toa special fit (see text) near the given epoch into one of 20 orbitaltime bins. Time is reckoned with respect to T Conj , the epoch ofthe pulsar’s superior conjunction with the companion. Each curvepeaks at that epoch, when the pulsar’s earthbound signals plungemost deeply into the companion’s gravitational well. The ampli-tude and shape of the curves evolve, due to relativistic precessionof the orbital ellipse, as quantified by the advancing longitude ofperiastron ω . two of these eight new parameters, namely δ obs θ and ˙ x obs ,are measured at the marginal 1 . σ level, while the third,˙ e obs , is only an upper limit. All others in this Table, in-cluding the new Shapiro terms, are measured with highconfidence. In the next sections, we discuss importantorbital measurables, including corrections that must bemade to some of the observed quantities in order to de-termine their intrinsic values. The Observed and Intrinsic Orbital PeriodDerivative
The observed orbital period derivative ˙ P obsb , must becorrected by a term ˙ P galb , resulting from the relativegalactic accelerations of the solar system and the bi-nary system (Damour & Taylor [1991, hereafter DT91],DT92), in order to yield the intrinsic derivative, ˙ P intrb :˙ P intrb = ˙ P obsb − ˙ P galb . (15)Using galactic parameters of R = 8 . ± .
16 kpc andΘ = 240 ± P galb = − (0 . ± . × − . Inserting also ˙ P obsb from Table 2 into Eq. 15, we calculate that ˙ P intrb = − (2 . ± . × − . The uncertainty in this resultis dominated by the error in ˙ P galb , which in turn is setprincipally by the pulsar distance uncertainty. A VLBAparallax campaign on the pulsar, currently in progress,will hopefully improve these uncertainties. First Successful Measurement of the ShapiroGravitational Propagation Delay Parameters inPSR B1913+16
Because of relativistic precession of the elliptical orbit,the Shapiro delay has recently grown to an amplitudeof ∼ µs around the orbit, rendering it relatively eas-ily measurable. Fig. 1 illustrates the enhancement ofthe Shapiro delay signal around the orbit over the lastdozen years, during which time the WAPP receivers alsocame into use, thereby increasing our observing band-width ten-fold. The curves in the Figure illustrate ageneral relativistic calculation of the expected Shapirodelay variation around the orbit, while the data pointsare residuals to a TEMPO fit freezing all parameters attheir best-fit values, except for the Shapiro parameter r (corresponding in General Relativity to the companionmass). The latter quantity was artificially set to zero tosimulate the absence of the Shapiro delay in the fit. Thepattern of residuals systematically matches theoreticalexpectations for the Shapiro delay.We have now successfully determined the two Shapiroterms in both the DD and FW parametrizations (see Ta-ble 2). While the Shapiro delay has been observed in sev-eral other binary pulsar systems, these results mark thefirst successful detection of a pair of Shapiro measurablesin the PSR B1913+16 system. [The two DD parameters, s and r , had been jointly constrained by Taylor & Weis-berg (1989). For the ensuing decade, unfavorable orbitalgeometry rendered its amplitude unmeasurably small;while the last decade has seen both improving orbital ge-ometry and advances in observing instrumentation.] Weaccounted for the significant nonlinear covariance of thethe s and r parameters in estimating their values anduncertainties in Table 2, using a process delineated bySplaver et al. (2002). The procedure is illustrated graph-ically in Fig. 2, which also shows the best-fitting FWShapiro parameters. (Tighter constraints on the inclina-tion and companion mass can be derived indirectly fromother measurements. See § Best Determination of Component Masses andOrbital Inclination
The measurement of the first seven quantities in Ta-ble 2 enables the precise general relativistic determina-tion of the component masses and orbital inclination.Specifically, our measurements of (cid:104) ˙ ω (cid:105) and γ , along withthe Keplerian elements, leave only the two unknowns m (cid:104) ˙ ω (cid:105) ,γ and m (cid:104) ˙ ω (cid:105) ,γ in the following two general rel-ativistic equations: (cid:104) ˙ ω (cid:105) =3 G / c − ( P b / π ) − / (1 − e ) − ( m (cid:104) ˙ ω (cid:105) ,γ + m (cid:104) ˙ ω (cid:105) ,γ ) / = 3 T / (cid:12) ( P b / π ) − / (1 − e ) − (cid:18) m (cid:104) ˙ ω (cid:105) ,γ + m (cid:104) ˙ ω (cid:105) ,γ M (cid:12) (cid:19) / , (16) J. M. Weisberg and Y. Huang Figure 2. (a) Measured constraints on | cos i | and m result-ing from TEMPO fits for two different parametrizations of theShapiro gravitational propagation delay within the context of gen-eral relativity. The Damour & Deruelle (1986) s and r param-eters map directly onto the displayed | cos i | [= + √ − s ] and m [= ( r/ T (cid:12) )M (cid:12) ] axes, respectively; the black contours show joint1, 2, and 3 σ (∆ χ = 2 . , . , .
8) confidence limits on those quan-tities, derived from a set of TEMPO fits to a large grid of (fixed)( | cos i | , m ). The alternate Freire & Wex (2010) best-fit parame-ter constraints and their ± σ limits are shown in green. Their fit-ted parameter ς transforms directly into the displayed | cos i | axis(see Eq. 4), whereas their h parameter does not map uniquelyonto either of the axes (see Eq. 5). The marginal distributionsin (b) and (c) result from collapsing the resulting two-dimensionalDD probability distribution onto the | cos i | and m axes respec-tively, in which the mean (solid black) and the 1 σ bounds (greyregion) are displayed, yielding | cos i | = 0 .
73 (+0 . , − .
11) and m = 1 .
95 (+0 . , − .
71) M (cid:12) (68.3% confidence). and γ = G / c − e ( P b / π ) / m (cid:104) ˙ ω (cid:105) ,γ ( m (cid:104) ˙ ω (cid:105) ,γ + 2 m (cid:104) ˙ ω (cid:105) ,γ ) × ( m (cid:104) ˙ ω (cid:105) ,γ + m (cid:104) ˙ ω (cid:105) ,γ ) − / = T / (cid:12) e ( P b / π ) / m (cid:104) ˙ ω (cid:105) ,γ M (cid:12) × (cid:18) m (cid:104) ˙ ω (cid:105) ,γ + 2 m (cid:104) ˙ ω (cid:105) ,γ M (cid:12) (cid:19) (cid:18) ( m (cid:104) ˙ ω (cid:105) ,γ + m (cid:104) ˙ ω (cid:105) ,γ M (cid:12) (cid:19) − / . (17)Solving simultaneously for the two component masses,we find that m (cid:104) ˙ ω (cid:105) ,γ = 1 . ± .
001 M (cid:12) and m (cid:104) ˙ ω (cid:105) ,γ =1 . ± .
001 M (cid:12) . These values agree with WNTto within 2 σ , while our precision is poorer due toa less-precisely determined γ (see § x and P b measure- ments: (sin i ) (cid:104) ˙ ω (cid:105) ,γ = 0 . ± . | cos i | (cid:104) ˙ ω (cid:105) ,γ = 0 . ± . m and sin i or | cos i | valuesdetermined directly from the Shapiro propagation delaymeasurements of § Toward the First Published Measurement of theRelativistic Orbital Shape Correction, δ θ , in anySystem We have successfully measured the apparent post-Keplerian orbital shape correction term, δ obs θ (see Table2). As noted in § (cid:15) A (see Eq. 9).Geodetic spin-precession modeling of this system shouldin principle determine the necessary aberration parame-ters by specifying the spin axis orientation (specifically,its polar angles η and λ ) over time (see Eq. 7). However,we find that the currently available pulseshape variationfits (Kramer 1998; Weisberg & Taylor 2002; Clifton &Weisberg 2008) yield inconsistent solutions for these pa-rameters. Consequently, although we have successfullymeasured δ obs θ , it is not yet possible to determine δ intr θ nor hence to use its measured value as an additional testof relativistic gravitation.Despite our current inability to measure δ intr θ , we candetermine its expected value δ GR θ via a general relativisticcalculation (DD; DT92): δ GR θ = G / c (cid:18) P b π (cid:19) − / ( m + m ) − / (cid:26) m + 6 m m + 2 m (cid:27) =T (cid:12) / (cid:18) P b π (cid:19) − / (cid:18) m + m M (cid:12) (cid:19) − / × (cid:40) (cid:18) m M (cid:12) (cid:19) + 6 (cid:18) m m M (cid:12) (cid:19) + 2 (cid:18) m M (cid:12) (cid:19) (cid:41) , (18)yielding δ GR θ = (6 . ± . × − . Eq. 9 can thenbe inverted to give the aberration parameter (cid:15) A = (2 . ± . × − . This timing-derived value of (cid:15) A will providea modest consistency check on future geodetic precessionmodeling. Implications of Fits for the First Time-Derivativesof e and x As noted above, Eqs. 12 and 13 demonstrate that thesuccessful measurement of ˙ e obs and ˙ x obs would lead toconstraints on other quantities of interest. It is thereforeuseful to further investigate the various terms composingthese equations.The first (aberration) term of both equations, d(cid:15) A /dt ,was defined in Eqs. 10 and 11. However, as discussed in § (cid:15) A , additional progress in understanding the pul-sar’s spin axis orientation is needed before d(cid:15) A /dt can beconfidently determined. Nevertheless, our current geode-tic precession modeling suggests that its value is in the ∼ − s − range, and varying with spin-precessionalphase. Constraints from ˙ e obs iming Relativistic Binary Pulsar PSR B1913+16 7The second and final term of Eq. 12 involves the timeevolution of e induced by gravitational wave (“GW”)emission. For the Eq. 12 second term, Peters (1964)shows that this term, (cid:18) ˙ ee (cid:19) GW = − c a µ ( m + m ) × (1 − e ) − / (cid:18) e (cid:19) = − (cid:12) / m M (cid:12) m M (cid:12) ( m + m M (cid:12) ) / (cid:18) P b π (cid:19) − / × (1 − e ) − / (cid:18) e (cid:19) = − . × − s − , (19)with µ the reduced mass. This term is negligible com-pared to the expected value of d(cid:15) A /dt , except at fortu-itous precessional phases where the latter can drop tozero. Consequently, the d(cid:15) A /dt term dominates Eq. 12,so that a successful measurement of ( ˙ e/e ) obs would pro-vide a unique measurement of d(cid:15) A /dt . Unfortunately,our fitted value of ( ˙ e/e ) obs is currently not significantlydifferent from zero, although its upper limit is in the10 − s − range (see Table 2) expected for d(cid:15) A /dt . Constraints from ˙ x obs The second term in Eq. 13 delineates the gravitationalwave-induced orbital shrinkage rate, which can be eval-uated from measurables via (cid:18) ˙ a a (cid:19) GW = 23 (cid:32) ˙ P b P b (cid:33) GW = − . × − s − . (20)The third (spin-orbit) term of Eq. 13 varies approxi-mately sinusoidally on the geodetic precession timescale,with an amplitude of ∼ ± × − s − . Details await arobust determination of geodetic precession parameters(see Eq. 14).The fourth (Kopeikin 1996) term of Eq. 13 has a max-imum amplitude of ∼ . × − s − ; while the fifthand final term, − ˙ DD = + ˙ P b , gal P b = − × − s − (21)(see § P b , gal ).In summary, the first (aberration) and third (spin-orbit) terms dominate Eq. 13, so all others may be ig-nored. However, neither of these two terms is currentlyaccurately determinable. We do have a marginally signif-icant measurement of ˙ x obs (see Table 2). Consequently,if either of the two terms becomes well-determined inthe future along with an improved value of ˙ x obs , thenthe other term will also become accessible. For exam-ple, there are two possible paths toward determination of The exact value depends on the unknown alignment on the skyof the line of nodes. the aberrational term: First, additional geodetic preces-sion observations and modeling should better constrain d(cid:15) A /dt ; and second, additional observations could betterdetermine ˙ e obs , which, as noted in § d(cid:15) A /dt . At this point,the spin-orbit term would be calculable, leading to anexciting measurement of the pulsar’s moment of inertia,which has important implications for neutron star equa-tions of state (Lattimer & Schutz 2005). With the mea-surement precision of ˙ x and ˙ e improving with time t as t − / , another decade or so of observations are required.Unfortunately, geodetic spin axis precession may causethe pulsar to disappear before that time. TESTS OF RELATIVISTIC GRAVITATION
The determination of seven particular independentquantities suffices to fully determine the dynamics of abinary system within the context of a particular theoryof gravitation. For example, the most accurate deter-mination of component masses and orbital inclination in § δ GR θ in § Gravitational Radiation Emission and ˙ P b Gravitational radiation emission should cause the orbitto decay as orbital energy is radiated away. The quantity˙ P GRb is the resulting orbital period derivative expectedfrom a general relativistic calculation of this phenomenon(Peters and Mathews 1963):˙ P GR b = − π G / c (cid:18) P b π (cid:19) − / (cid:18) e + 3796 e (cid:19) (1 − e ) − / × m m ( m + m ) − / = − π T / (cid:12) (cid:18) P b π (cid:19) − / (cid:18) e + 3796 e (cid:19) (1 − e ) − / × (cid:18) m M (cid:12) (cid:19) (cid:18) m M (cid:12) (cid:19) (cid:18) m + m M (cid:12) (cid:19) − / (22)Inserting our measured and derived values and theiruncertainties into Eq. 22 , we find that ˙ P GR b =( − . ± . × − . To verify our estimate ofthe error in ˙ P GR b that was derived via propagation of un-certainty, we also employed a Monte-Carlo method withCholesky decomposition of the covariance matrix. In thisfashion, we simulated the joint normal distribution ofmeasured parameters { γ, ˙ ω, P b , e } , and then constructed This value may also be calculated directly from the first sevenorbital measurables in Table 2 alone, without the use of derivedquantities such as the masses (DT91).
J. M. Weisberg and Y. Huang
Figure 3.
The orbital decay of PSR B1913+16 as a function oftime. The curve represents the orbital phase shift expected fromgravitational wave emission according to General Relativity. Thepoints, with error bars too small to show, represent our measure-ments. a histogram of 1000000 derived ˙ P GR b and inferred theuncertainty therefrom.Consequently, we find that˙ P intr b ˙ P GR b = ( − . ± . × − ( − . ± . × − = 0 . ± . . (23)This result demonstrates that the system is losing en-ergy to gravitational radiation within ∼ σ of the ratepredicted by general relativity (see also Fig. 3 and thered curve in Fig. 4). The above number represents asignificant improvement over the value determined byWNT, 0 . ± . σ discrep-ancy between our measurements and general relativity.Interestingly, the new galactic parameters of Reid et al.(2014) are the principal reason for the improvement (viaa change in ˙ P galb ), while our measured values themselveschanged little.In addition to confirming general relativistic radiationdamping at this level, our result rules out large parameterspaces in plausible scalar-tensor theories of gravity. Inrecent years, however, other pulsars in neutron star-whitedwarf binary systems have overtaken PSR B1913+16 inconstraining these alternatives (Freire et al. 2012).DD92 point out that this test is a “mixed” strong-field probe in that it involves a combination of radiativeeffects (via ˙ P obs b ) and quasi-static phenomena (through (cid:104) ˙ ω (cid:105) and γ , whose values are needed in order to make aprediction of the expected ˙ P GR b ). Consequently, addi-tional tests, such as those described in the next section,that probe different aspects of strong-field gravitation,are also useful in constraining viable alternatives to gen-eral relativity. Shapiro Gravitational Propagation Delay
Each of the two newly measured Shapiro parametersrepresents another independent test of relativistic gravi-tation. As with the ˙ P b test of § (cid:104) ˙ ω (cid:105) and γ inorder to make a testable prediction for the value of theShapiro parameters. In this case, unlike the ˙ P b test, allof the post-Keplerian quantities probe quasi-static phe-nomena in strong fields. While Shapiro parameters havealready been measured in several other binary systems, itis especially useful to constrain theories via systems suchas this one and PSRs B1534+12 and J0737-3039A, whereat least three “excess” post-Keplerian parameters beyond (cid:104) ˙ ω (cid:105) and γ (one gravitational radiation parameter andtwo Shapiro quantities) are measurable. Although theprecision of the binary pulsar Shapiro parameter mea-surements is well below their measurement precision inthe weak solar gravitational field, it is this simultaneousdetermination of several parameters in strong-field con-ditions in each of these binary pulsar systems that leadsto the important constraints on relativistic theories ofgravitation.In the DD formulation of the Shapiro delay within gen-eral relativity, the Shapiro measurables s and r map di-rectly onto sin i and m , respectively. Hence we can testGeneral Relativity by comparing the Shapiro determina-tion of sin i ( ≡ s ) with that determined from (cid:104) ˙ ω (cid:105) and γ (called sin i (cid:104) ˙ ω (cid:105) ,γ ; see § s sin i (cid:104) ˙ ω (cid:105) ,γ = 0 . +0 . − . . ± . . +0 . − . . (24)Similarly, we can test General relativity by comparingthe Shapiro determination of m ( ≡ r M (cid:12) / T (cid:12) ) with thatdetermined from (cid:104) ˙ ω (cid:105) and γ : r M (cid:12) / T (cid:12) m (cid:104) ˙ ω (cid:105) ,γ = (1 . +0 . − . ) M (cid:12) (1 . ± . (cid:12) = 1 . +0 . − . . (25)The consistency (albeit with a rather low level of pre-cision) of these Shapiro determinations of sin i and m with those measured via the other post-Keplerian terms,and hence their confirmation of general relativity, is alsographically depicted in Fig. 4. The Shapiro terms havealso been measured in several other binary pulsar sys-tems with higher precision, and have also been shown tobe in agreement with general relativity. CONCLUSIONS
We report here on the measurements and relativisticanalyses of 9309 TOAs in over thirty years of high-qualityArecibo data on binary pulsar PSR B1913+16. We fit-ted for a number of previously unmeasurable parametersfor the first time in this system (and in one case, for thefirst time anywhere), which enabled us to significantlyadvance our relativistic analyses of the system. We pro-vide our newest measurents or derivations of all relevantphysical quantities of the binary system, with the excep-tion of Ω, the position angle of the line of nodes. Werigorously ascertained the uncertainties in the fitted andderived parameters. Having fully characterized the sys-tem, we proceeded to use it in several tests of generalrelativity in strong-field conditions.We have measured a gravitational radiation-inducediming Relativistic Binary Pulsar PSR B1913+16 9
Figure 4.
Constraints on the masses of the pulsar and the com-panion as a function of five post-Keplerian measurables, within thecontext of General Relativity. The width of each curve represents ± σ error bounds. The mutual near-intersection of all curves il-lustrates the agreement of our observations with general relativityin the strong-field conditions at the binary system. Table 3
Comparison of Gravitational Radiation-Induced Orbital Decaywith GR Prediction in Binary PulsarsPSR ˙ P intrb / ˙ P GRb
Ref.J0348+0432 1 . ± .
18 Antoniadis et al. (2013)J0737-3039 1 . ± .
014 Kramer et al. (2006)J1141-6545 1 . ± .
06 Bhat, Bailes, & Verbiest (2008)B1534+12 0 . ± .
06 Stairs et al. (2002)J1738+0333 0 . ± .
13 Freire et al. (2012)J1756-2251 1 . ± .
03 Ferdman et al. (2014)J1906+0746 1 . ± . a van Leeuwen et al. (2015)B1913+16 0 . ± . . ± .
03 Jacoby et al. (2006) a Assumes negligible proper motion. orbital period decrease whose rate agrees with the gen-eral relativistic expectation to within ∼ σ , which iscloser than found by WNT, largely as a result of an im-proved galactic correction resulting from more accurategalactic parameters (Reid et al. 2014).Similar orbital decay tests have now been performedwith several other binary pulsars (see Table 3 for pub-lished measurements of ˙ P intrb / ˙ P GRb ). The orbital de-cays of PSRs J0348+0432, J0737-3039, J1141-4565,J1738+0333, J1906+0746, B1913+16, and B2127+11Call exhibit agreement between observation and generalrelativity to within (or very close to) the authors’ stateduncertainties. PSR B1913+16 currently has the mostprecise such determination, and interferometric parallaxmeasurements, currently in progress, will hopefully fur-ther tighten the precision. Among the other two, PSR B1534+12 and PSR J1756-2251, various systematic effects such as an incorrect dis-tance in the galactic acceleration correction may explainthe small observed discrepancies, although it is possiblethat an incompleteness of general relativity or some un-known physical effect is responsible. See the work of Fer-dman et al. (2014) for an especially thorough descriptionof the most significant deviation of the orbital decay ratefrom the general relativistic prediction, found in PSRJ1756-2251.Our new (for this system) measurements of Shapirogravitational propagation delay parameters representtwo additional tests of relativistic gravitation, and arefully consistent with general relativity, although their rel-ative precision is currently far lower than the orbital de-cay test. This binary now joins several other systems, in-cluding PSRs J0737-3039A, B1534+12, and J1756-2251,in each providing at least three independent tests of rel-ativistic, strong-field gravitation.We have also marginally measured the orbital shapeparameter δ θ for the first time anywhere, but its intrin-sic value is corrupted by a comparable, undeterminedaberration delay. Future geodetic spin-orbit precessionmeasurements should lead to an accurate characteriza-tion of the aberration and then an additional relativisticgravitational test via the comparison of the aberration-corrected δ intr θ with δ GR θ .In addition, we fitted for the time derivative of orbitaleccentricity e and projected semimajor axis of the pulsarorbit, x , and we achieved an upper limit on the formerand a marginal detection of the latter. We discussedand quantified the various physical phenomena that cancontribute to these parameters. Unless the pulsar dis-appears in the next few years due to geodetic spin axisprecession, future timing observations should better de-fine these quantities, allowing for a determination of thepulsar’s moment of inertia I .We have placed online a subroutine and modificationsto the TEMPO TOA fitting software, which codes theFreire & Wex (2010) parametrization of the Shapiro de-lay for high-eccentricity binary pulsars such as the PSRB1913+16 system. An online companion of this paperprovides the TEMPO input files and TOAs upon whichthese analyses are based. The same data are on arXiv,and at zenodo: http://dx.doi.org/10.5281/zenodo.54764.Much of this experiment was pioneered by J. H. Tay-lor, to whom we owe our deepest thanks. D. J. Niceassisted with observing and analyses, and A. A. Chaelassisted in developing the FW analysis package. The au-thors gratefully acknowledge financial support from theUS National Science Foundation. The Arecibo Obser-vatory is operated by SRI International under a coop-erative agreement with the National Science Foundation(AST-1100968), and in alliance with Ana G. Mendez-Universidad Metropolitana, and the Universities SpaceResearch Association. Facility:
Arecibo
REFERENCESAntoniadis, J., Freire, P. C. C., Wex, N., et al. 2013, Science, 340,4480 J. M. Weisberg and Y. Huang