Abstract

PSR J1903+0327, a millisecond pulsar in an eccentric (e= 0.44) 95-d orbit with an ∼1 M companion poses a challenge to our understanding of stellar evolution in binary and multiple-star systems. Here we describe optical and radio observations which rule out most of the scenarios proposed to explain formation of this system. Radio timing measurements of three post-Keplerian effects yield the most precise measurement of the mass of a millisecond pulsar to date: 1.667 ± 0.021 solar masses (99.7 per cent confidence limit). This rules out some equations of state for superdense matter; furthermore, it is consistent with the spin-up of the pulsar by mass accretion, as suggested by its short spin period and low magnetic field. Optical spectroscopy of a proposed main-sequence counterpart shows that its orbital motion mirrors the pulsar’s 95-d orbit; being therefore its binary companion. This finding rules out a previously suggested scenario which proposes that the system is presently a hierarchical triple. Conventional binary evolution scenarios predict that, after recycling a neutron star into a millisecond pulsar, the binary companion should become a white dwarf and its orbit should be nearly circular. This suggests that if PSR J1903+0327 was recycled, its present companion was not responsible for it. The optical detection also provides a measurement of the systemic radial velocity of the binary; this and the proper motion measured from pulsar timing allow the determination of the systemic 3D velocity in the Galaxy. We find that the system is always within 270 pc of the plane of the Galaxy, but always more than 3 kpc away from the Galactic Centre. Thus an exchange interaction in a dense stellar environment (like a globular cluster or the Galactic Centre) is not likely to be the origin of this system. We suggest that after the supernova that formed it, the neutron star was in a tight orbit with a main-sequence star and the present companion was a tertiary farther out. The neutron star then accreted matter from its evolving inner companion, forming a millisecond pulsar. The inner companion then disappeared, either due to a chaotic three-body interaction with the outer star (caused by the expansion of the inner orbit that necessarily results from mass transfer), or in the case of a very compact inner system, due to ablation/accretion by the newly formed millisecond pulsar. We discuss in detail the possible evolution of such a system before the supernova.

1 INTRODUCTION

PSR J1903+0327 was the first millisecond pulsar (MSP1) discovered in the ongoing Arecibo L-band Feed Array (ALFA) pulsar survey (Cordes et al. 2006). In the discovery paper Champion et al. (2008), we presented the results of phase-coherent radio timing of this pulsar carried out with the Green Bank and Arecibo radio telescopes. These quickly revealed that the pulsar was in a 95-d orbit around a 1 solar mass (M) companion. This object is remarkable for being the first (and thus far, the only) disc MSP known to have an eccentric (e = 0.44) orbit. In globular clusters (GCs) there are several binary MSPs with eccentric orbits; but those are thought to be caused by perturbations of the binary systems by occasional close interactions with other stars.

Coincident with the pulsar position derived from the timing, a star was found whose near-infrared magnitudes were consistent with a 1-M main-sequence star at the distance and reddening estimated for PSR J1903+0327. It was not known then whether this was just an unlikely (∼2.6 per cent) chance alignment or whether the star is genuinely associated with PSR J1903+0327, and if so whether it is the binary companion responsible for the 95-d orbit of the pulsar. Such a finding would be surprising, as the conventional understanding of MSP evolution posits the scenario that such a neutron star (NS) is spun up to high spin frequencies by accretion of matter and angular momentum from a companion star while the companion passes through a giant phase (Bhattacharya & van den Heuvel 1991). This circularizes the system, and a recycled MSP is left orbiting a low-mass white dwarf (the remnant core of the donor) in a low-eccentricity orbit (e < 10−3; Phinney 1992). Until the discovery of PSR J1903+0327, all known MSPs in the Galactic disc had such low-eccentricity orbits. For reviews, see Phinney & Kulkarni (1994), Stairs (2003), Tauris & van den Heuvel (2006).

For these reasons, Champion et al. (2008) proposed that PSR J1903+0327 may be part of a triple system where the 95-d orbit of the pulsar is caused by a massive unseen WD and the third member is the star detected in the near-infrared. The latter is in a long-period orbit and drives the eccentricity of the inner pair through the Kozai mechanism (Kozai 1962). An alternative possibility, also discussed in Champion et al. (2008), is that the companion to PSR J1903+0327 in the 95-d orbit is the star detected in the near-infrared, but that this eccentric, unusual system originated in an exchange interaction in a dense stellar environment, like a globular cluster.

In this paper, we present new optical measurements and further radio timing of PSR J1903+0327 obtained with the aim of testing these scenarios. The plan for the rest of this paper is as follows. The optical and radio observations are described in Section 2. The immediate results from these observations are described in Section 3. In Section 4 we discuss the implications of these results regarding the formation and evolution of this system. In Section 5 we discuss how this system might have formed. We summarize our main conclusions in Section 6.

2 OBSERVATIONS

2.1 Optical observations

The long-slit spectroscopy of the suspected counterpart to PSR J1903+0327 was obtained with FORS2 (Appenzeller et al. 1998), the low-dispersion spectrograph of ESO’s Very Large Telescope. Four spectra were obtained on 2008 June 21, three on 2008 August 23 and one a day later, on August 24. All spectra had exposure times of 46 min., and used a 1-arcsec slit combined with the 1028Z holographic grism, providing a wavelength coverage of 7830–9570 Å. The detectors were read out with 2 × 2 binning, yielding a resolution of 3.4 Å, sampled at 0.86 Å pixel−1. The slit was placed such that both the pulsar companion and a bright nearby star were centred on the slit. The observations were made during clear and photometric nights, with the seeing between 0.48 and 0.72 arcsec. The spectral observations were corrected for bias and were flat-fielded using lamp flats.

Spectral extraction is complicated by the bright star, henceforth star A, located 2.3 arcsec from the pulsar counterpart (see Fig. 1). The star is brighter by about 5.6 mag in the I band and as a result about 20 per cent of the detected counts at the spatial position of the pulsar counterpart are from the wings of the brighter star. The regular optimal extraction algorithm as described by Horne (1986) is not suited to extract blended spectra as it makes no assumptions about the spatial profile. Instead, we use a variation on the algorithm by Hynes (2002), which is based on the optimal extraction method of Horne (1986).

Figure 1

This 32 × 32 arcsec2 subsection of a 5-min SDSS i-band image taken with GMOS at Gemini North on September 20th, 2007 shows the location of the 1 arcsec slit. Besides the companion to PSR J1903+0327 denoted by ‘PSR’, spectra were also extracted of stars A, B and C.

Figure 1

This 32 × 32 arcsec2 subsection of a 5-min SDSS i-band image taken with GMOS at Gemini North on September 20th, 2007 shows the location of the 1 arcsec slit. Besides the companion to PSR J1903+0327 denoted by ‘PSR’, spectra were also extracted of stars A, B and C.

Hynes (2002) uses an analytic function to describe the spatial profile as a function of wavelength. This profile is fitted to the spectrum of an isolated template source and then used to simultaneously extract the spectra of the blended sources. We use a variation on this method. First, instead of using a Voigt function to describe the spatial profile, we use a Moffat (1969) function, essentially a modified Lorentzian with a variable exponent. Secondly, instead of removing the sky contribution before extraction, we include it in the fit, representing the sky as a first-order polynomial added to the Moffat profiles for each object in the blend. Finally, the absence of an isolated template source forced us to use star A as a template. We used an iterative scheme to converge the properties of the profile as a function of wavelength from star A while removing the contribution of the pulsar counterpart and another faint star (star B) that was also part of the blend. In addition to these three objects, a fourth star (star C) is located 28 arcsec north-west of PSR J1903+0327 on the slit. This star was not blended and extracted normally.

Arc lamp exposures obtained during daytime with the telescope pointing towards the zenith were used for wavelength calibration. However, comparison of wavelengths of the night sky emission lines from the science exposures showed wavelength offsets of up to 0.8 Å between different exposures, most likely due to flexure caused by the telescope pointing away from zenith. To correct for these offsets the wavelength calibration of June 21st was used to create a secondary line list of some 60 night sky emission lines in the first spectrum, which was subsequently used to calibrate the remaining seven spectra. Typical rms residuals of these fits were less than 0.05 Å.

Fig. 2 shows the averaged spectrum of the counterpart to PSR J1903+0327. The spectrum has a signal-to-noise ratio (S/N) of 15–20 in the wavelength range of 8000–9000 Å. The spectrum is generally featureless except for clear absorption lines around 8540 and 8660 Å. At these wavelengths the most common spectral lines are those belonging to the Ca ii triplet, and the Paschen series of hydrogen (see Cenarro et al. 2001). Though the Paschen series has lines on these wavelengths, the absence of the other lines in the series argues for their being due to the calcium triplet.

Figure 2

The spectrum of the counterpart of PSR J1903+0327. Shown with error bars are the average data from the eight individual spectra, normalized, shifted to zero velocity and binned by a factor of 2. The bestfitting normalized spectrum from the library of Munari et al. (2005) is shown as the solid black line, shifted downwards by 0.5 units. The average night-sky emission spectrum is shown by the solid grey line, showing the location of emission lines as regions of reduced sensitivity. The vertical scale is arbitrary. The location of stellar absorption lines belonging to the Ca ii triplet and the Paschen series of hydrogen are denoted by grey shade above the spectrum.

Figure 2

The spectrum of the counterpart of PSR J1903+0327. Shown with error bars are the average data from the eight individual spectra, normalized, shifted to zero velocity and binned by a factor of 2. The bestfitting normalized spectrum from the library of Munari et al. (2005) is shown as the solid black line, shifted downwards by 0.5 units. The average night-sky emission spectrum is shown by the solid grey line, showing the location of emission lines as regions of reduced sensitivity. The vertical scale is arbitrary. The location of stellar absorption lines belonging to the Ca ii triplet and the Paschen series of hydrogen are denoted by grey shade above the spectrum.

2.2 Radio timing observations

We made pulse time of arrival (TOA) measurements of PSR J1903+0327 using the 305-m Arecibo radio telescope and the 105-m Green Bank Telescope from 2006 December through 2010 January; the observation dates and radio frequencies are summarized in Fig. 3. The TOAs obtained before 2008 January were described and used by Champion et al. (2008). We now describe the observational set-ups used at the two observatories.

Figure 3

Top: radio frequencies of TOAs versus time (in years). Bottom: post-fit TOA residuals versus time (in years). The different colours indicate different observing systems. Light blue: Green Bank SPIGOT data, centred at 2100 MHz, Pink: Green Bank SPIGOT data, centred at 1800 MHz. This was not previously taken into account. Dark blue: Arecibo WAPP data taken with the ‘S-high’ receiver. Light green: Arecibo WAPP data taken with the ‘S-wide’ receiver. Red: Arecibo WAPP data taken with the ‘L-wide’ receiver. Black: Arecibo ASP data (also taken with the ‘S-wide’ receiver).

Figure 3

Top: radio frequencies of TOAs versus time (in years). Bottom: post-fit TOA residuals versus time (in years). The different colours indicate different observing systems. Light blue: Green Bank SPIGOT data, centred at 2100 MHz, Pink: Green Bank SPIGOT data, centred at 1800 MHz. This was not previously taken into account. Dark blue: Arecibo WAPP data taken with the ‘S-high’ receiver. Light green: Arecibo WAPP data taken with the ‘S-wide’ receiver. Red: Arecibo WAPP data taken with the ‘L-wide’ receiver. Black: Arecibo ASP data (also taken with the ‘S-wide’ receiver).

In this re-analysis we use all the TOAs obtained from processing of data taken with the Green Bank S-band receiver (with frequency coverage is from 1650 to 2250 MHz); these cover most of the year 2007. The data were acquired using the ‘Spigot’ pulsar back-end – a three-level autocorrelation spectrometer which samples observing bands of up to 800 MHz (Kaplan et al. 2005). Autocorrelation functions (ACFs) of length 2048 lags were accumulated with three-level precision and written to disc every 81.92 μs. The ACFs were subsequently Fourier transformed to synthesize 2048-channel power spectra. We separately analysed the lower and upper halves of the band, dedispersing each and deriving TOAs using the methods described in the Supplement to Champion et al. (2008). The TOA uncertainties are 6 μs in the upper half and 10 μs in the lower half of the band. In Champion et al. (2008), only data from the upper half of the band were used.

The Arecibo data were obtained with four Wideband Arecibo Pulsar Processors (WAPPs; Dowd, Sisk & Hagen 2000). With the ‘L-wide’ receiver, the WAPPs had a bandwidth of 50 MHz and were centred at 1520, 1570, 1620 and 1670 MHz; each produced 512-lag ACFs sampled with three-level precision, accumulated and integrated every 64 μs. With the ‘S-wide’ receiver, the WAPPs had a bandwidth of 50 MHz and were centred at 2125, 2175, 2225 and 2275 MHz, though the last one was sometimes disconnected in order to use the Arecibo Signal Processor (ASP, described below) in its place. At these frequencies, the WAPPs produced 256-lag ACFs which were integrated for 32 μs. With the ‘S-high’ receiver, the WAPPs were configured in dual-board mode, allowing 100 MHz of bandwidth in each of eight correlator boards, centred at 100 MHz intervals from 3150 to 3850 MHz. Each board produced 128-lag ACFs integrated for 32 μs.

As for the Spigot data, the ACFs from the WAPPs are Fourier transformed to produce power spectra. These data were then dedispersed and folded modulo the pulsar’s period using a routine written by one of us (IHS) specifically for this purpose; this folds the spectra using polynomial coefficients calculated specifically for the radio frequencies used in the observations. The folding results in 128-bin pulse profiles every 500 s for each WAPP.

During some of the S-wide observations we collected data with the Arecibo Signal Processor (ASP, Demorest 2007) in parallel with three WAPPs. The ASP band had a centre frequency of 2350 MHz. The ASP coherently dedisperses a maximum of 16 4-MHz bands, for a maximum bandwidth of 64 MHz. These data were then immediately folded at the pulsar’s rotational period, and 512-bin pulse profiles were stored for each of these bands every 120 s. However, because of the large dispersion measure (DM) of PSR J1903+0327 we can only coherently dedisperse 12 of these 4-MHz bands at any given time. Finally, we add all the bands and four 120-s integrations to produce pulse profiles with a good S/N.

All pulse profiles are cross-correlated in the Fourier domain (Taylor 1992) with a low-noise template derived from the sum of the pulse profiles obtained with the same spectrometer and at the same frequency. From this, we derive a total of 1872 topocentric TOAs. Both the high-resolution WAPP and ASP data provide TOA measurements with an average uncertainty of ∼ 1 μs.

In the next step we use the tempo2 software package (Hobbs, Edwards & Manchester 2006) for the TOA analysis. This applies the clock corrections intrinsic to the observatory, the Earth rotation data and the observatory coordinates to convert the TOAs to terrestrial time (TT), as maintained by the Bureau International des Poids et Mesures. To convert TT TOAs to coordinated barycentric time (TCB) TOAs, we used the DE/LE 421 Solar system ephemeris (Folkner, Williams & Boggs 2008). The program then minimizes the squares of TOA residuals — the difference between the observed and predicted TOAs. We used the ‘DD’ orbital model (Damour & Deruelle 1985, 1986); this is optimal for describing eccentric orbits in a theory-independent manner. The best-fitting parameters are listed in Table 1; their uncertainties are the 1σ values estimated by tempo2 (throughout this paper, the quoted uncertainties are 1σ, except when stated otherwise). The TOA residuals are displayed at the bottom plot of Fig. 3.

Table 1

Timing and derived parameters for PSR J1903+0327 using tempo and tempo2 as described in the text. All parameters are as measured at the Solar system barycentre. In parentheses we present the 1σ uncertainties, except in a few cases (b) where we present 99.7 per cent confidence limits.

Timing programme  tempo  tempo2 
Timing parameters 
Solar system ephemeris … … … …  DE 405/LE 405  DE 421/LE 421 
Reference time-scale … … … …  TDB  TCB 
Orbital model … … … …  Modified DD  DD 
Reference time (MJD) … … … …  55000  55000 
Right ascension, α (J2000) … … … …  19h03m05s793 296 (2)  19h03m05s793 213 (10) 
Declination, δ (J2000) … … … …  03°27′19″.21053 (6)  03°27′19″.20911 (6) 
Proper motion in α, μα (mas yr−1) … … … …  −2.01 (7)  −2.06 (7) 
Proper motion in δ, μδ (mas yr−1) … … … …  −5.20 (12)  −5.21 (12) 
Spin frequency, ν (Hz) … … … …  465.135 245 551 237 (10)  465.135 238 339 217 (9) 
First derivative of ν, forumla (10−15 Hz s−1) … … … …  −4.0719 (2)  −4.0719 (2) 
Dispersion measure, DM (cm−3 pc) … … … …  297.5244 (6)  297.5245 (6) 
First derivative of DM (cm−3 pc yr−1) … … … …  −0.0083 (6)  −0.0084 (6) 
Orbital period Pb (d) … … … …  95.174 117 277 (14)  95.174 118 753 (14) 
Projected semimajor axis, x (lt s) … … … …  105.593 4628 (5)  105.593 4643 (5) 
Eccentricity, e… … … …  0.436 678 410 (3)  0.436 678 409 (3) 
Longitude of periastron, ω (°) … … … …  141.653 1042 (2)  141.652 4786 (6) 
Time of passage through periastron, T0 (MJD) … … … …  550 15.581 404 51 (4)  550 15.581 588 59 (4) 
Derivative of x, forumla (10−15lt sforumla) … … … …  +20 (3)  +21 (3) 
Apsidal motion, forumla (forumla) … … … …  0.000 2400 (2)  0.000 2400 (2) 
‘Range’ parameter of the Shapiro delay, r/T (M) … … … …  1.03 (3) 
‘Shape’ parameter of the Shapiro delay, s… … … …  0.9760 (15) 
Orthometric amplitude of the Shapiro delay, h3 (μs) … … … …  2.602 (25) 
Orthometric ratio of the Shapiro delay, ς… … … …  0.803 (6) 
Limits (not fitted with other timing parameters) 
Second derivative of ν, forumla (10−26 Hz s−2) … … … …  4.8 (1.7)  6.1 (1.6) 
First derivative of e, forumla (10−16 s−1) … … … …  1.4 (6) 
First derivative of Pb, forumla (forumla) … … … …  −53 (33)  −64 (31) 
Derived parameters 
Spin period, P (ms) … … … …  2.149 912 364 349 21 (7) 
First derivative of spin period, forumla (forumla) … … … …  1.882 03 (11) 
Characteristic age, τc (109 yr) … … … …  1.9 
Dipolar magnetic flux density at the poles, B0 (108 G) ...  2.0 
Galactic longitude, l… … … …  37°.3363 
Galactic latitude, b… … … …  −1°.0136 
Distance, D (kpc) … … … …  6.4 ± 1.0 
Total proper motion, μ (mas yr−1) … … … …  5.60 (11) 
Galactic position angle of proper motion, Θμ… … … …  263°.9 ± 0.8° 
Velocities, solar system barycentre reference frame 
Transverse velocity, VT (km s−1) … … … …  168 ± 24 
Radial velocity, γ (km s−1) … … … …  42.0 ± 4.4 
Total 3D velocity, V (km s−1) … … … …  174 ± 25 
Velocities, Galactocentric reference frame (cylindrical) 
V ρ (km s−1) … … … …  17 ± 13 
V ϕ (km s−1) … … … …  −189 ± 5 
Vz (km s−1) … … … …  −11 ± 4 
Total 3D velocity, V (km s−1) … … … …  190 ± 5 
Velocities, relative to pulsar standard of rest (cylindrical) 
V ϕ (km s−1) … … … …  28 ± 5 
Total 3D velocity, V (km s−1) … … … …  37 ± 9 
Derived masses and orbital orientation 
Mass function, f (M) … … … …  0.139 558 441 (2) 
Orbital inclination, i (°) … … … …  77°.47 (15) or 102°.53 (15) 
Total mass of binary, Mt (M) … … … …  2.697 (29)b 
Companion mass, mc (M) … … … …  1.029 (8)b 
Pulsar mass, mp (M) … … … …  1.667 (21)b 
Timing programme  tempo  tempo2 
Timing parameters 
Solar system ephemeris … … … …  DE 405/LE 405  DE 421/LE 421 
Reference time-scale … … … …  TDB  TCB 
Orbital model … … … …  Modified DD  DD 
Reference time (MJD) … … … …  55000  55000 
Right ascension, α (J2000) … … … …  19h03m05s793 296 (2)  19h03m05s793 213 (10) 
Declination, δ (J2000) … … … …  03°27′19″.21053 (6)  03°27′19″.20911 (6) 
Proper motion in α, μα (mas yr−1) … … … …  −2.01 (7)  −2.06 (7) 
Proper motion in δ, μδ (mas yr−1) … … … …  −5.20 (12)  −5.21 (12) 
Spin frequency, ν (Hz) … … … …  465.135 245 551 237 (10)  465.135 238 339 217 (9) 
First derivative of ν, forumla (10−15 Hz s−1) … … … …  −4.0719 (2)  −4.0719 (2) 
Dispersion measure, DM (cm−3 pc) … … … …  297.5244 (6)  297.5245 (6) 
First derivative of DM (cm−3 pc yr−1) … … … …  −0.0083 (6)  −0.0084 (6) 
Orbital period Pb (d) … … … …  95.174 117 277 (14)  95.174 118 753 (14) 
Projected semimajor axis, x (lt s) … … … …  105.593 4628 (5)  105.593 4643 (5) 
Eccentricity, e… … … …  0.436 678 410 (3)  0.436 678 409 (3) 
Longitude of periastron, ω (°) … … … …  141.653 1042 (2)  141.652 4786 (6) 
Time of passage through periastron, T0 (MJD) … … … …  550 15.581 404 51 (4)  550 15.581 588 59 (4) 
Derivative of x, forumla (10−15lt sforumla) … … … …  +20 (3)  +21 (3) 
Apsidal motion, forumla (forumla) … … … …  0.000 2400 (2)  0.000 2400 (2) 
‘Range’ parameter of the Shapiro delay, r/T (M) … … … …  1.03 (3) 
‘Shape’ parameter of the Shapiro delay, s… … … …  0.9760 (15) 
Orthometric amplitude of the Shapiro delay, h3 (μs) … … … …  2.602 (25) 
Orthometric ratio of the Shapiro delay, ς… … … …  0.803 (6) 
Limits (not fitted with other timing parameters) 
Second derivative of ν, forumla (10−26 Hz s−2) … … … …  4.8 (1.7)  6.1 (1.6) 
First derivative of e, forumla (10−16 s−1) … … … …  1.4 (6) 
First derivative of Pb, forumla (forumla) … … … …  −53 (33)  −64 (31) 
Derived parameters 
Spin period, P (ms) … … … …  2.149 912 364 349 21 (7) 
First derivative of spin period, forumla (forumla) … … … …  1.882 03 (11) 
Characteristic age, τc (109 yr) … … … …  1.9 
Dipolar magnetic flux density at the poles, B0 (108 G) ...  2.0 
Galactic longitude, l… … … …  37°.3363 
Galactic latitude, b… … … …  −1°.0136 
Distance, D (kpc) … … … …  6.4 ± 1.0 
Total proper motion, μ (mas yr−1) … … … …  5.60 (11) 
Galactic position angle of proper motion, Θμ… … … …  263°.9 ± 0.8° 
Velocities, solar system barycentre reference frame 
Transverse velocity, VT (km s−1) … … … …  168 ± 24 
Radial velocity, γ (km s−1) … … … …  42.0 ± 4.4 
Total 3D velocity, V (km s−1) … … … …  174 ± 25 
Velocities, Galactocentric reference frame (cylindrical) 
V ρ (km s−1) … … … …  17 ± 13 
V ϕ (km s−1) … … … …  −189 ± 5 
Vz (km s−1) … … … …  −11 ± 4 
Total 3D velocity, V (km s−1) … … … …  190 ± 5 
Velocities, relative to pulsar standard of rest (cylindrical) 
V ϕ (km s−1) … … … …  28 ± 5 
Total 3D velocity, V (km s−1) … … … …  37 ± 9 
Derived masses and orbital orientation 
Mass function, f (M) … … … …  0.139 558 441 (2) 
Orbital inclination, i (°) … … … …  77°.47 (15) or 102°.53 (15) 
Total mass of binary, Mt (M) … … … …  2.697 (29)b 
Companion mass, mc (M) … … … …  1.029 (8)b 
Pulsar mass, mp (M) … … … …  1.667 (21)b 

We have performed a parallel analysis using the tempo software package.2 The resulting timing parameters and their uncertainties (as estimated by tempo) are also presented in Table 1. Several timing parameters are different from those estimated by tempo2. All time-like quantities (spin frequency, orbital period, projected semimajor axis, time of passage through periastron and ω, which is strongly covariant with the latter) are different because tempo converts the TOAs to barycentric dynamic time (TDB), a time-scale that is different from TCB. The right ascension and declination are different because the DE/LE 405 Solar system ephemeris (Standish 1998) uses an earlier celestial reference frame. Furthermore, we included in tempo a modified version of the DD orbital model that uses the orthometric parametrization of the Shapiro delay, as described in (Freire & Wex 2010). The ‘orthometric amplitude’ (h3) and the ‘orthometric ratio’ (ς) provide an improved description of the areas of the (mc, sin i) plane where the system is likely to be compared to the ‘range’ (r) and ‘shape’ (s) parameters used in the normal parametrization. In cases where we can measure other post-Keplerian (PK) parameters, the new parametrization yields an improved test of general relativity.

Using either tempo or tempo2 there are no apparent trends in the residuals and the normalized χ2 is 1.3, which is similar to what we observe in other MSPs timed with these observing systems. This means that the timing model provides a complete description of the observed TOAs, with no significant, unmodelled effects present. This also means that the 1σ uncertainties reported by tempo and tempo2 are reasonably accurate estimates of the real uncertainties. The non-unity value of the normalized χ2 likely reflects a slight underestimation of the TOA uncertainties by the cross-correlation analysis, but it could also be due to red noise in the pulsar rotation, as suggested by the marginal (almost 3σ) detection of forumla. To account for this, the timing parameters in Table 1 were obtained after the multiplication of the TOA uncertainty estimates by a scalefactor; this is calculated separately for each data set so that its normalized χ2 is 1. This scalefactor is 1.1 for all the Arecibo data and 1.3 for the GBT data.

3 RESULTS

3.1 On the nature of the companion

Radial velocities and spectral parameters were determined by comparing the observed spectra, both individual and averages at each epoch, with synthetic spectral templates from the spectral library of Munari et al. (2005). These models range in effective temperature Teff, surface gravity g, rotational velocity Vrot and metallicity [M/H] and are sampled at a resolution of λ/Δλ = 20 000. The wavelength range between 8400 and 8800 Å was used. To remove the effects of the finite resolution of the instrument, the synthetic spectra were convolved with a Gaussian profile with a width equal to the seeing, and truncated at the width of the slit. For each object, an iterative procedure was used to first determine the radial velocities using a starting template, then use those velocities to shift all spectra to zero velocity and create an average which was then used to obtain a better matching template for determining more accurate radial velocities.

The comparison of the average spectrum with synthetic spectra from the library of Munari et al. (2005) yields a temperature of Teff = 5825 ± 200 K (1σ) and sets a 2σ lower limit on the surface gravity of forumla, firmly excluding the possibility of a giant or subgiant star (see Fig. 4). Based on the shape of the absorption lines we estimate a 3σ upper limit on the rotational broadening of vrotsin i* < 140 km s−1, where i* is the inclination of the star. The observed spectrum does not constrain the metallicity. Stellar evolution models (Girardi et al. 2000) show that a main-sequence (MS) star with the observed mass and effective temperature of the companion of PSR J1903+0327 is consistent with the observed infrared colours for distances between 6 and 8 kpc, assuming the estimated reddening as a function of distance from 2MASS stars (Champion et al. 2008). This is in agreement with the 6.4-kpc distance estimated from the pulsar’s DM (see Table 1) using the Cordes & Lazio (2002) model of the Galactic electron distribution. The ages predicted from these models range from 4 to 7 Gyr.

Figure 4

Confidence contours for the effective temperature (Teff) and surface gravity [log (g cm s−2)] for the companion of PSR J1903+0327.

Figure 4

Confidence contours for the effective temperature (Teff) and surface gravity [log (g cm s−2)] for the companion of PSR J1903+0327.

3.2 Radial velocities

Table 2 shows the radial velocities of the pulsar counterpart relative to the Solar system barycentre. They are V1 = 92.4 ± 8.3 km s−1 during the first epoch, and V2 = 20.2 ± 5.2 km s−1 during the second. The velocity difference between the two epochs of V1V2 = 72.2 ± 9.8 km s−1 deviates from constant velocity at the 7.4σ level. Between the two epochs, the averaged radial velocities of stars A and C differ by V1V2=−14.2 ± 0.3 km s−1 for star A and −8.4 ± 0.4 km s−1 for star C. These differences are caused by slight differences in centring of the source on the slit (Bassa et al. 2006). Based on the location of star A with respect to the centre of the slit before and after each exposure, we estimate an average offset in radial velocity of about −8.6 ± 3.4 km s−1 for the first epoch and 3.7 ± 2.9 km s−1 for the second epoch. This is comparable to the weighted average in V1V2=−12.1 ± 0.3 km s−1 for stars A and C.

Table 2

The celestial position, JHKs magnitudes and averaged radial velocities (relative to the Solar system barycentre) at both epochs of the four stars located on the slit. The position and magnitudes are obtained from the data described in Champion et al. (2008).

ID  αJ2000  δJ2000  J  H  K s  V 1 (km s−1 V 2 (km s−1
forumla   +03°27′21″.31  14.85 ± 0.09  14.44 ± 0.10  14.23 ± 0.19  42.0 ± 0.2  56.3 ± 0.2 
forumla   +03°27′23″.83  19.49 ± 0.09  18.52 ± 0.10  18.08 ± 0.09  49.6 ± 6.1  48.6 ± 2.7 
forumla   +03°27′47″.50  15.07 ± 0.09  14.52 ± 0.10  14.27 ± 0.09  45.9 ± 0.2  54.3 ± 0.4 
PSR  forumla   +03°27′19″.18  19.22 ± 0.09  18.41 ± 0.10  18.03 ± 0.09  92.4 ± 8.3  20.2 ± 5.2 
ID  αJ2000  δJ2000  J  H  K s  V 1 (km s−1 V 2 (km s−1
forumla   +03°27′21″.31  14.85 ± 0.09  14.44 ± 0.10  14.23 ± 0.19  42.0 ± 0.2  56.3 ± 0.2 
forumla   +03°27′23″.83  19.49 ± 0.09  18.52 ± 0.10  18.08 ± 0.09  49.6 ± 6.1  48.6 ± 2.7 
forumla   +03°27′47″.50  15.07 ± 0.09  14.52 ± 0.10  14.27 ± 0.09  45.9 ± 0.2  54.3 ± 0.4 
PSR  forumla   +03°27′19″.18  19.22 ± 0.09  18.41 ± 0.10  18.03 ± 0.09  92.4 ± 8.3  20.2 ± 5.2 

Subtracting these radial velocity offsets from the measured velocities of the counterpart, we obtain V1 = 101.0 ± 8.7 km s−1 and V2 = 16.5 ± 5.9 km s−1, and a velocity difference of V1V2 = 85 ± 11 km s−1. This is consistent with the V1V2 = 88.72 km s−1 difference in radial velocity predicted from the pulsar timing ephemeris. This confirms that the optical counterpart to PSR J1903+0327 is the object in the 95-d orbit around the pulsar. Fig. 5 shows the radial velocity measurements of the companion and the radial velocity predictions based on the orbital parameters determined from pulsar timing. This provides an independent estimate of the mass ratio (R = 1.55 ± 0.20) and the systemic radial velocity of the binary, γ = 44.3 ± 4.9 km s−1, which cannot be derived from the radio timing.

Figure 5

The predicted radial velocities of the pulsar (solid line) and companion (dashed line) based on the pulsar timing ephemeris. These are vertically offset by the systemic radial velocity of the system relative to the Solar system barycentre γ, which cannot be determined from the radio measurements. The two optical radial velocity measurements are shown with their error bars. Their difference is consistent with the predictions from the pulsar timing ephemeris, and their absolute values allow the estimation of γ = 44.3 ± 4.9 km s−1.

Figure 5

The predicted radial velocities of the pulsar (solid line) and companion (dashed line) based on the pulsar timing ephemeris. These are vertically offset by the systemic radial velocity of the system relative to the Solar system barycentre γ, which cannot be determined from the radio measurements. The two optical radial velocity measurements are shown with their error bars. Their difference is consistent with the predictions from the pulsar timing ephemeris, and their absolute values allow the estimation of γ = 44.3 ± 4.9 km s−1.

3.3 Shapiro delay

Since the companion is a non-degenerate star, it could in principle have a strong stellar wind, which should be detectable as a variation of DM as a function of the orbital phase. Such a signal would produce a distortion in our measurement of the Shapiro delay. In Fig. 6 we display the DM averaged over 36 bins of the pulsar’s mean anomaly, after correction of the long-term DM variations. For some of these intervals we have smaller amounts of data or a small range of frequencies. If these result in DM determinations with uncertainties greater than forumla, they are not depicted. We detect no DM variations greater than forumla. This means that the mass estimates derived from the Shapiro delay are accurate, within their uncertainties.

Figure 6

Top: radio frequency of TOAs as a function of orbital phase. The orbital phases with the larger number and frequency spread of observations are those for which the DM measurements should be better. These are best at an orbital phase of 0.95, where superior conjunction happens. Bottom: DM as a function of orbital phase. All TOAs were divided into 10° orbital phase bins and a DM was derived for them, keeping all other parameters fixed.

Figure 6

Top: radio frequency of TOAs as a function of orbital phase. The orbital phases with the larger number and frequency spread of observations are those for which the DM measurements should be better. These are best at an orbital phase of 0.95, where superior conjunction happens. Bottom: DM as a function of orbital phase. All TOAs were divided into 10° orbital phase bins and a DM was derived for them, keeping all other parameters fixed.

The better timing precision, larger number of TOAs, longer data span and optimized orbital coverage of our new data set not only improve previous measurements of relativistic effects (namely the observed apsidal motion forumla and the s parameter of the Shapiro delay; see Champion et al. 2008) but also allow a precise measurement of the r parameter of the Shapiro delay, in addition to a measurement of the proper motion μ and the variation it causes on the apparent size of the orbit forumla (see Table 1). We now discuss the relevance of these parameters.

3.3.1 Shapiro delay and component masses

With the previously available timing data it was not possible to measure the companion mass mc from the Shapiro delay alone. The precise (and frequent) Arecibo timing has completely changed this situation, providing a very clear Shapiro delay signal (see Fig. 7). Assuming that general relativity (GR) is the correct theory of gravity, we obtain  

1
formula
where T=G M/c3 = 4.925 490 947 μs is the solar mass M in time units. As usual, G is Newton’s gravitational constant and c is the velocity of light. The orbital inclination derived from ς is i = 77°.4 ± 0°.4 or 102°.6 ± 0°.4. Given the mass function f, these values imply a pulsar mass mp = (1.67 ± 0.08) M, a total binary mass Mt = (2.70 ± 0.11) M and a mass ratio Rs = 1.62 ± 0.03.

Figure 7

Top: post-fit TOAs versus orbital phase. The Shapiro delay was not taken into account, but all the Keplerian parameters were fitted. Bottom: post-fit TOAs versus orbital phase, with Shapiro delay taken into account. The different colours indicate different observing systems, as in Fig. 3.

Figure 7

Top: post-fit TOAs versus orbital phase. The Shapiro delay was not taken into account, but all the Keplerian parameters were fitted. Bottom: post-fit TOAs versus orbital phase, with Shapiro delay taken into account. The different colours indicate different observing systems, as in Fig. 3.

3.3.2 Estimating the uncertainties

The uncertainties quoted above were estimated using the Bayesian technique described by Splaver et al. (2002). We assume that mc and cos i have constant a priori probability. For each point in a grid of mc, cos i values, we calculate (trivially) the r and s parameters to describe the Shapiro delay, assuming that general relativity is the correct theory of gravitation.

We then fit a timing solution similar to that of Table 1 to the TOAs but keep (r, s) fixed. We record the resultant χ2 and calculate the probability density for that (mc, cos i) space using  

2
formula
where χ2min is the minimum value of χ2 in the whole map. The contour levels that include the 68.3, 95.4 and 99.7 per cent of all probability are displayed in Fig. 8 as the thin contours. We then project the 2D probability density function (PDF) into the mc and cos i axes to calculate lateralized 1D PDFs.

Figure 8

Companion mass as a function of cos i and mp. The thick contour levels are derived from a 3D χ2 map of the cos i, mc and Ω space (where Ω is the position angle of the line of nodes) and then collapsed on the planes represented in the figure (see text for details). The thin contour levels represent a 2D χ2 map of the cos imc space calculated taking only the Shapiro delay into account. The lines represent the constraints derived from the spectroscopic mass ratio (R, black dashed), the apsidal motion (forumla, solid orange – here with increased uncertainty due to the proper motion, solid pink), the harmonic amplitude (h3) and harmonic ratio (ς) of the Shapiro delay (in purple) and finally an upper limit on the inclination given by forumla (pink solid line). The grey area in the mass–mass diagram is excluded by sin i≤ 1. In the marginal plots we can see that the 1D probability distribution functions for the pulsar and companion masses are much narrower when the apsidal motion (even with uncertainty caused by the proper motion) is taken into account (thick lines), but this assumes that there are no significant classical contributions to forumla. The latter must be <2.3 arcsec century−1 (1σ) given the agreement between the h3, ς and forumla bands.

Figure 8

Companion mass as a function of cos i and mp. The thick contour levels are derived from a 3D χ2 map of the cos i, mc and Ω space (where Ω is the position angle of the line of nodes) and then collapsed on the planes represented in the figure (see text for details). The thin contour levels represent a 2D χ2 map of the cos imc space calculated taking only the Shapiro delay into account. The lines represent the constraints derived from the spectroscopic mass ratio (R, black dashed), the apsidal motion (forumla, solid orange – here with increased uncertainty due to the proper motion, solid pink), the harmonic amplitude (h3) and harmonic ratio (ς) of the Shapiro delay (in purple) and finally an upper limit on the inclination given by forumla (pink solid line). The grey area in the mass–mass diagram is excluded by sin i≤ 1. In the marginal plots we can see that the 1D probability distribution functions for the pulsar and companion masses are much narrower when the apsidal motion (even with uncertainty caused by the proper motion) is taken into account (thick lines), but this assumes that there are no significant classical contributions to forumla. The latter must be <2.3 arcsec century−1 (1σ) given the agreement between the h3, ς and forumla bands.

We translate the 2D PDF of the (cos i, mc) space to a 2D PDF of the (mp, mc) space using the mass function f (Lorimer & Kramer 2005). This is then lateralized on to the mp axis, resulting in a 1D PDF for the pulsar mass. It is from the 1D PDFs that we derive the uncertainties for cos i, mp and mc.

3.4 Orientation of the orbit and component masses

3.4.1 Kinematic effects: orbital orientation

The additional, previously undetermined parameters in Table 1 are the proper motion (with total magnitude μ and position angle Θμ) and the apparent variation of the size of the orbit, forumla; this is given by  

3
formula
where i is the orbital inclination; Ω is the position angle of the line of nodes; forumla is the variation of x due to an intrinsic variation of i; forumla is a variation of the semimajor axis a; D is the distance to the pulsar; and al is the difference of Galactic accelerations of the binary and the Solar system, projected along the line of sight.

The first term is by far the largest; it is an effect of the proper motion, which constantly changes the viewing geometry. If this causes a change of i as viewed from the Earth then there will be a secular change of the projected size of the orbit, x=a sin i/c (Arzoumanian et al. 1996; Kopeikin 1996).

The second term describes the changing Doppler shift of the binary relative to the Solar system barycentre. For a nominal DM distance of 6.4 kpc, this amounts to forumla. This is 3 orders of magnitude smaller than the measurement error; this term is thus ignored in the following discussion. The third term results from any real changes in the orbital plane of the binary system. The only likely contribution to this is from spin–orbit coupling (Smarr & Blandford 1976; Lai, Bildsten & Kaspi 1995); the resulting forumla should be 1 order of magnitude smaller than the measurement error (see detailed discussion in Section 3.4.4). Finally, the fourth term results from a real change in the size of the orbit. The prediction for forumla caused by gravitational wave emission is of the order of forumla; many orders of magnitude too small. Therefore, only the first term is likely to make a significant contribution to the observed forumla.

In this case, for each of the two possible values of i discussed above, there are two possible (Θμ−Ω) solutions; these are depicted in Fig. 9 by the intersection of the forumla and ς lines. However, because of the uncertainty in the measurements of forumla and ς, these two pairs of solutions merge into two wide areas centred at (Θμ−Ω) ∼ 90° for i∼ 102°.6 and (Θμ−Ω) ∼ 270° for i∼ 77°.4.

Figure 9

Confidence contours for cos i and mc as a function of the orientation of the line of nodes relative to the direction of the proper motion Θμ−Ω. The contour levels include 66.3, 95.4 and 99.7 per cent of the total probability, which is derived from a 3D χ2 map of the (cos i, mc, Θμ−Ω) space and then collapsed on to the planes represented in the figure. The dashed contour lines indicate 90° < i < 180°. The vertical purple dashed lines are the constraints introduced by the harmonic ratio of the Shapiro delay (ς), the pink solid curves are the constraints derived from the measurement of forumla. In the marginal plots are displayed the 1D probability distribution functions. We can see that there is still a large uncertainty in (Θμ−Ω). This introduces a large uncertainty in the mass estimates (of which mc is shown in the plot).

Figure 9

Confidence contours for cos i and mc as a function of the orientation of the line of nodes relative to the direction of the proper motion Θμ−Ω. The contour levels include 66.3, 95.4 and 99.7 per cent of the total probability, which is derived from a 3D χ2 map of the (cos i, mc, Θμ−Ω) space and then collapsed on to the planes represented in the figure. The dashed contour lines indicate 90° < i < 180°. The vertical purple dashed lines are the constraints introduced by the harmonic ratio of the Shapiro delay (ς), the pink solid curves are the constraints derived from the measurement of forumla. In the marginal plots are displayed the 1D probability distribution functions. We can see that there is still a large uncertainty in (Θμ−Ω). This introduces a large uncertainty in the mass estimates (of which mc is shown in the plot).

3.4.2 Kinematic effects: apsidal motion and masses

The observed apsidal motion forumla arcsec century−1 is, in principle, a combination of the relativistic apsidal advance forumla; a kinematic contribution forumla due to the changing viewing geometry of the system; and a contribution caused by spin–orbit coupling forumla:  

4
formula
We will now consider each of these terms in detail.

The forumla is caused by the change of viewing geometry due to the proper motion μ (Kopeikin 1996):  

5
formula
where μ′=μ/sin i = 0.56 arcsec century−1. Given the aforementioned uncertainty of Θμ−Ω, it is impossible at the moment to know the exact value of forumla, but it has a very sharp upper limit given by |cos (Θμ−Ω)| ≤ 1. This means that uncertainty in forumla caused by the lack of precise knowledge of forumla is at most equal to ±μ′ = 0.56 arcsec century−1. This is more than six times larger than the current uncertainty in the measurement of forumla.

This affects the precision of our estimate of the total mass of the system. If forumla, then the total system mass could be determined directly using  

6
formula
where nb = 2π/Pb is the orbital frequency and Pb is the orbital period (e.g. Lorimer & Kramer 2005); the result would be Mt = (2.697 ± 0.004) M. Because of the uncertainty of forumla and forumla, we obtain instead Mt = (2.697 ± 0.029) M (99.7 per cent confidence limit, see Section 3.4.5).

3.4.3 Further uncertainty estimates

To estimate these uncertainties rigorously, we have made a second χ2 map that uses the principles described above, but includes a third dimension, the difference in the position angles of the proper motion and line of nodes Θμ−Ω. This map makes use of the extra information provided by the measurements of forumla and forumla.

As before, from the values of mc and cos i at each point we calculate r and s. But now we also use the mass function f to calculate mp. From the total mass Mt=mp+mc we calculate forumla using equation 6. From the values of cos i and (Θμ−Ω) at each point we calculate forumla using equation 3 and forumla using equation 5. This is made using only the best value for the proper motion μ, which has a relatively small uncertainty.

We then fit the resulting timing solution to the TOAs, but keep the four computed parameters fixed (r, s, forumla and forumla). We record the resultant χ2 and calculate the probability density for that point of the (mc, cos i, Θμ−Ω) space using an expression similar to equation 2. This 3D ‘cube’ of probability density is then projected on to its faces: the (cos i, mc) space (see thick contours in Fig. 8) and the (cos i, Θμ−Ω) and (mc, Θμ−Ω) spaces (Fig. 9). We then project this 3D PDF on to the three axes, calculating lateralized 1D PDFs for mc, cos i and Θμ−Ω. It is from these that we derive the uncertainties of these parameters.

In the (mc, Θμ−Ω) panel of Fig. 9 we can see that if Θμ−Ω were more precisely known, we would have a much better estimate of mc. This is the graphical demonstration that the kinematic effects are now the main source of uncertainty in the measurement of mc and mp, not the uncertainty intrinsic to the measurement of forumla. We can also see that as the precision of forumla improves, we will have a total of four possible values for Θμ−Ω, two for each possibility of cos i. The two solutions closer to the centre of Fig. 9μ−Ω = 180°) will have a large value of mc (and mp), since, according to equation 5, for this orbital orientation forumla, therefore forumla. The opposite will be true for the two solutions closer to Θμ−Ω = 0°. Unless we can independently determine cos i or Θμ−Ω, we will have two degenerate values for mc and mp. In principle, this degeneracy can be lifted if we can measure the orbital annual parallax. This should be possible with a fourfold improvement in timing precision.

3.4.4 Contribution from spin–orbit coupling

Thus far we have assumed that the spin–orbit contributions to forumla and forumla are negligible. Comparing equation (7.54) and the following in Will (1993) with the more general equation in Wex (1998), we find for forumla:  

7
formula
where a is the semimajor axis of the orbit of the companion as seen from the pulsar, Rc is the companion radius (which we assume to be one solar radius), Pb and e are as listed in Table 1, J2 is the quadrupolar moment of the star, θ is the angle between the orbital and companion’s angular momentum and Φ0 is the longitude of the ascending node in a reference frame defined by the total angular momentum vector. For PSR J1903+0327, we obtain a maximum value of forumla arcsec century−1.

As we have shown in Section 3.1, the companion star has a mass, age and temperature similar to that of the Sun. Such stars should rotate slowly like the Sun (Irwin et al. 2009); the latter has a rotational velocity of vrot, ⊙ = 2 km s−1. This should result in a quadrupolar moment similar to that of the Sun, J2, ⊙∼ 1.7 × 10−7. For this value of J2 we get forumla arcsec century−1, about seven times smaller than the measurement uncertainty of forumla.

However, because the companion’s rotation is not well constrained by the optical measurements (which indicate vrotsin i* < 140 km s−1, 3σ; see Section 2.1), we can use the agreement among h3, ς and forumla to constrain forumla assuming that GR is the correct theory of gravity. The minimum total mass compatible (at the 1σ level) with h3 and ς is 2.59 M. The corresponding minimum forumla (84.1 arcsec century−1) can then be used to derive a 1σ limit of forumla arcsec century−1 (here we are assuming a median expected value of 0 for forumla). Future optical/near-infrared measurements might be able to better constrain vrot. If this is small then the agreement of the three different PK parameters can be used directly as a test of general relativity.

To calculate forumla, we use equation (81) of (Wex 1998) to derive  

8
formula
Using the parameters above, we obtain forumla, which is 1 order of magnitude below the measurement error. Therefore, an improved measurement of forumla will likely improve the determination of the orbital orientation.

3.4.5 Component masses and orbital inclination

From the 3D map, we obtain the following parameters: i = 77°.47 ± 0°.15 or 102°.53 ± 0°.15 (1σ), Mt = (2.697 ± 0.029) M, mc = (1.029 ± 0.008) M, mp = (1.667 ± 0.021) M and R = 1.620 ± 0.008. The underlying probability distribution functions for these parameters have a shape that is very different from a Gaussian curve (see Fig. 8), therefore it is not very meaningful to refer to σ limits. For all the previous values we indicated instead 99.7 per cent confidence limits. These values and their uncertainties can be understood qualitatively as the result of the intersection of the 1σ bands of h3 and forumla. Their intersection results in measurement of i that is more precise than the value derived from s (or ς). Such a result cannot be understood using the regular r, s parametrization of the Shapiro delay.

These values are entirely consistent with those derived solely from the Shapiro delay in Section 3.3, but much more precise, despite the uncertainty introduced by forumla.

4 IMPLICATIONS

4.1 Neutron star mass

The recent mass measurement of PSR J1614−2230 [mp = (1.97 ± 0.4) M, 1σ; see Demorest et al. 2010] rules out the presence of hyperons, bosons and free quarks at densities comparable to the nuclear saturation density, and demonstrates that NSs can be stable at masses well above the Chandrasekhar mass. The mass measurement of PSR J1903+0327 (mp = (1.667 ± 0.021) M, see Section 3.4) supports the latter conclusion, being also inconsistent with the softest proposed equations of state for superdense matter. This plus the mass measurements of PSR J0437−4715 (Verbiest et al. 2008) and PSR J0621+1002 (Splaver et al. 2002) confirm beyond doubt the results of a previous statistical analysis of masses of several binary MSPs in globular clusters (Freire et al. 2008a), which showed that MSPs have a much wider mass distribution than observed among the members of double-neutron star systems. The likely cause is the accretion episode that spun up these NSs (Bhattacharya & van den Heuvel 1991; Lorimer 2008).

There are several reasons why PSR J1903+0327 appears to be a normal MSP (i.e., spun up by accretion of mass from a companion): (a) its spin period is typical of that of MSPs and ∼8 times shorter than that of the fastest young pulsar, PSR J0537−6910, near the N157B supernova remnant in the Large Magellanic Cloud (Marshall et al. 1998); (b) its magnetic field (B0) is similar to those observed among MSPs and 2 orders of magnitude smaller than the smallest B0 thus far observed for a young pulsar, 3 × 1010 G for PSR J1852+0040, a 105-ms X-ray pulsar near the centre of the supernova remnant Kesteven 79 (Halpern & Gotthelf 2010); (c) its mass suggests it has accreted mass from a former donor star. Precise mass measurements for components of double-neutron star systems range from 1.25 M (Kramer et al. 2006) to 1.44 M (Weisberg, Nice & Taylor 2010); this shows that a large majority of NSs form with masses near the Chandrasekhar limit (Lorimer 2008). If this was also the case for PSR J1903+0327 then it must have gained ΔM∼ 0.22–0.42 M during accretion, which is well within the range expected of low-mass and intermediate-mass binary pulsars (Pfahl, Rappaport & Podsiadlowski 2002). All of this is inconsistent with scenarios where PSR J1903+0327 forms directly from the collapse of a single star. Even alternative MSP formation scenarios (like merger of massive WDs or accretion induced collapse) require two stars to form the MSP.

4.2 PSR J1903+0327 is not a triple system

The discovery that the companion to PSR J1903+0327 is the MS star previously suspected of being associated with the system immediately rules out the detailed triple system scenario proposed in Champion et al. (2008). Our measurement of the variation of eccentricity had already ruled out the possibility that the Kozai mechanism is the origin of the observed eccentricity (Gopakumar, Bagchi & Ray 2009).

Despite this, it is clear that the present companion cannot be the star that recycled the pulsar. If it were, then it must have at some point filled the region where matter will remain bound to it (its ‘Roche lobe’) in order to transfer matter to the companion. The problem is that at the closest approach between the stars the companion is ∼23 times smaller than its present Roche lobe. Being an unevolved MS star it was even smaller in the past3 and even less likely to have supplied the material that recycled this pulsar.

Could the former donor have exchanged orbits with the main-sequence companion and now orbit the J1903+0327 binary system at a large distance? If there were a third component in the system, the J1903+0327 binary system would be accelerating towards it, with a line-of-sight component given by al. This should produce a variation of its orbital period given by forumla. Variations in al due to the relative motion of the J1903+0327 binary system and the outer component should produce a variation in the first derivative of the spin frequency (forumla), forumla (e.g. Backer, Foster & Sallmen 1993). Our measurements of forumla and forumla are not statistically significant (see Table 1), i.e. we detect no third component in this system. If the former donor still exists, it is no longer bound to this system.

4.3 Motion of PSR J1903+0327 in the Galaxy

A possible explanation for the absence of the former donor is that it might have been exchanged by the present companion. Before the discovery of PSR J1903+0327, this was the only mechanism known to produce MSP binaries with eccentric orbits. This can only occur in environments with high stellar densities such as the cores of globular clusters (Phinney 1992) or, hypothetically, the Galactic Centre. This is the reason why all other eccentric binary MSPs are observed in globular clusters (Freire et al. 2004; Ransom et al. 2004, 2005; Freire et al. 2008b). For this reason it was proposed in Champion et al. (2008) that PSR J1903+0327 might have formed in such an environment and subsequently have been ejected by the recoil produced by the presumed exchange interaction.

To investigate this possibility, we used the radial velocity γ and the pulsar’s Galactic coordinates, DM distance and proper motion to derive all three coordinates of its position in space and all three components of its velocity vector, as in Wex, Kalogera & Kramer (2000) and Lazaridis et al. (2009); the results are displayed in Table 1. We can thus calculate the past trajectory of the binary in the Galactic potential (Kuijken & Gilmore 1989). We implement this process in a Monte Carlo simulation of 10 000 orbits in the model Galactic potential. We take the uncertainties of μ, γ and distance4 into account by integrating many orbits with random initial conditions having a distribution of parameters consistent with the observed values and uncertainties. We also compare our results with those obtained using an alternative model for the Galactic potential (Paczyński 1990).

Our results are presented in Figs 10 and 11. The starting orbits have a range of proper motions and systemic velocities consistent with the uncertainties in Table 1. We also assume a 1σ relative uncertainty of 15 per cent in the distance. For each point we integrate the pulsar orbit back in time for about 300 Myr and record the maximum (Rmax) and minimum (Rmin) distances to the centre of the Galaxy. The diagonal black line shown in the top plot of Fig. 11 indicates circular orbits and the grey line eccentricities of 0.2. Rmax cannot be smaller than the minimum possible distance to the Galactic Centre at present (blue line); this is given by R0sin l, where R0 is the Sun’s distance to the Galactic Centre [7.7 kpc in the model used for this simulation (Kuijken & Gilmore 1989)] and l is the pulsar’s Galactic longitude, 37°.33. In this plot we can see that the minimum possible distance from the Galactic Centre is always larger than 3 kpc. This excludes formation in the Galactic Centre and in the bulge globular clusters.

Figure 10

The present Galactic positions (dots) and past orbits, integrated for 1 Gyr (solid curves) of the four MSPs with known radial velocities: PSR J1903+0327 (Red), PSR J1012+5307 (dark blue; Lazaridis et al. 2009), PSR J1023+0038 (light blue; Thorstensen & Armstrong 2005; Archibald et al. 2009) and PSR J1738+0333 (olive green). The Sun’s position is indicated by the yellow circle. In the top plot we look at the Galaxy from the North Galactic pole, in the bottom plot from the direction in the Galactic plane perpendicular to the line from the Sun to the Galactic Centre. We can see that PSR J1903+0327 is always near the Galactic plane, and that its orbital velocity about the Galaxy has a relatively small eccentricity.

Figure 10

The present Galactic positions (dots) and past orbits, integrated for 1 Gyr (solid curves) of the four MSPs with known radial velocities: PSR J1903+0327 (Red), PSR J1012+5307 (dark blue; Lazaridis et al. 2009), PSR J1023+0038 (light blue; Thorstensen & Armstrong 2005; Archibald et al. 2009) and PSR J1738+0333 (olive green). The Sun’s position is indicated by the yellow circle. In the top plot we look at the Galaxy from the North Galactic pole, in the bottom plot from the direction in the Galactic plane perpendicular to the line from the Sun to the Galactic Centre. We can see that PSR J1903+0327 is always near the Galactic plane, and that its orbital velocity about the Galaxy has a relatively small eccentricity.

Figure 11

Monte Carlo simulation of 10 000 pulsar orbits. Top: minimum and maximum distances from the centre of the Galaxy. We can see here that the pulsar is never within 3 kpc from the centre of the Galaxy Bottom: starting velocities in the previous simulation relative to each simulated pulsar’s average of local velocities (i.e., the system’s peculiar velocity).

Figure 11

Monte Carlo simulation of 10 000 pulsar orbits. Top: minimum and maximum distances from the centre of the Galaxy. We can see here that the pulsar is never within 3 kpc from the centre of the Galaxy Bottom: starting velocities in the previous simulation relative to each simulated pulsar’s average of local velocities (i.e., the system’s peculiar velocity).

It is also apparent from our simulations that the distance from the Galactic plane never exceeds 270 pc (99.7 per cent confidence limit), which is unusual for MSPs. The globular clusters observed outside the bulge have a wide distribution of Galactic heights (Harris 1996), which implies that most of them orbit the Galaxy well outside the plane. The formation of PSR J1903+0327 in such a cluster would require that the system is ejected exactly at the time that the cluster is crossing the Galactic plane, and that the resulting velocity is very nearly in the Galactic plane and close to the average local Galactic rotation (the system’s current peculiar velocity is 37 ± 9 km s−1; see the bottom plot of Fig. 11). We find this possibility highly unlikely. The probability of formation in an exchange interaction appears therefore to be extremely small.

5 FORMATION OF PSR J1903+0327

To summarize, the measured mass of PSR J1903+0327 (plus its small spin period and magnetic field) make any scenario where it forms as it is unlikely; all pulsar characteristics are consistent with it having been recycled (Section 4.1). Our limits on the variation of the eccentricity and our optical detection exclude the triple system scenario proposed in Champion et al. (2008; Section 4.2). In particular, we find that the present companion is not the donor star that recycled the pulsar; if this former donor still exists it is no longer bound to the system. Furthermore, our calculations of the 3D systemic motion (Section 4.3) show that it is unlikely that the system originated in an environment with high stellar density. Therefore, it is unlikely that the former donor star was exchanged by the present companion. So how did PSR J1903+0327 form? And what happened to the former donor star?

5.1 Multiple star spiral-in

We suggest that the system started as a triple (indicated below in square brackets), where the outer companion (the present MS star companion) was initially on a much wider orbit than at present. The inner binary system (indicated in parentheses) had a much shorter orbital period and consisted of two somewhat more massive MS stars:  

formula
(the numbers are conjectural, we show them for illustrative purposes only). The more massive star (the progenitor of the pulsar) evolves and becomes a red supergiant (RSG):  
formula
At this point, it starts transferring mass to the other component in the inner binary. In such a situation, mass transfer is generally unstable, and the envelope of the RSG is expected to engulf the companion and lead to a common envelope (CE) phase, where the helium core of the RSG and the inner companion are embedded in the envelope of the RSG. Because of friction with the envelope, the orbit of this embedded binary will decay, releasing orbital energy in the process. If this energy is sufficient to eject the envelope, this will leave a much closer binary consisting of a helium star and the more-or-less unaffected inner companion (see, e.g. Bhattacharya & van den Heuvel 1991):  
formula
If the outer component of the triple system is close enough, it may also be affected by this CE phase, since the envelope will greatly expand because of the orbital energy that is being deposited within it (to radii of ∼5–10 au; see, e.g. Podsiadlowski 2001). It will also be engulfed by the expanding envelope and will also experience a spiral-in phase (Eggleton & Verbunt 1986). This will produce a much closer triple system. Without realistic hydrodynamical simulations, it is difficult to determine the post-CE parameters. Here, we will consider two potential outcomes, which we will discuss in more detail later:  
formula
or  
formula
Eventually, the He core will explode in a supernova and produce an NS. The sudden mass-loss associated with this event will make the orbits of the other two components wider and somewhat eccentric:  
formula
or  
formula
This event could be either a normal Fe-core collapse supernova (SN) or an electron-capture (e-capture) supernova (see Podsiadlowski et al. 2004, and further references therein). Given the fact that the system remained bound after the explosion, it is not likely that the SN produced a major kick, nor large fractional mass-loss, particularly under scenario ‘a’, where the system is wider and acquires a small eccentricity. This is consistent with the observed peculiar velocity for the system, which is smaller than for the other systems where we can make a complete determination of this quantity (see Fig. 10). We note too that such SN characteristics are expected of e-capture supernovae.

5.2 Post-supernova evolution

At this stage the inner binary resembles the progenitor of a typical low- or intermediate-mass X-ray binary (L/IMXB). Once the normal stellar component of the inner binary has evolved to fill its Roche lobe, mass transfer to the NS will start to increase its mass and spin it up to millisecond periods (see, e.g. Pfahl et al. 2002). Initially, when the separation of the inner binary is sufficiently small compared to the periastron of the outer component, the evolution of the L/IMXB will not be significantly affected by the presence of the outer component.

If the mass-transfer rate is very high, mass-loss from the inner binary will cause the orbit of the outer component to widen to some degree. However, even for conservative mass transfer the inner binary is generally expected to widen more drastically, in particular in the 2.5-d scenario (scenario ‘a’ in Section 5.1; see Podsiadlowski, Rappaport & Pfahl 2002),  

formula
until its separation reaches a critical fraction of the periastron distance of the outer companion. This leads to a chaotic three-body interaction (see, e.g. Hut 1984; Phillips 1993). One of the most probable outcomes is the ejection of the least massive component – in this case most likely the mass donor in the inner binary, leaving a bound, eccentric binary system consisting of the MSP and the formerly outer companion, in an orbit tighter than before:  
formula
If the starting orbital period is small (PB < 8 h, scenario b in Section 5.1), then gravitational radiation might lead to the destruction of the remnant of the former donor by the newly formed MSP, very much in the same way as PSR 1957+20, the black-widow pulsar (Fruchter, Stinebring & Taylor 1988) is evaporating its companion (although in the latter case the process might take several Hubble times). This is a possible formation channel for isolated MSPs. Since these represent ∼20 per cent of the known MSPs in the Galactic disc, this should not be too unlikely a scenario. The end result would also resemble the J1903+0327 binary system.

5.3 Alternative scenarios

One can imagine alternative triple scenarios to the one outlined above. For example, it is possible that the progenitor of the MSP was initially the outer component of a stable triple system. When this massive star evolved to become a red supergiant, it could have engulfed the inner binary completely, leading to a spiral-in of that binary inside the common envelope. As it spirals in, its orbit could be disrupted by the tide induced by the more massive component; both of its components then start orbiting the core of the more massive He star independently, but at different distances given their different orbital velocities prior to disruption. This would produce a result quite similar to the double star spiral-in. Another possibility is that this binary survived the CE phase, but it was later disrupted when its orbital radius increased due to a mass transfer stage.

6 CONCLUSIONS

Our Arecibo and Green Bank radio timing observations include a full determination of the relativistic Shapiro delay, a very precise measurement of the apsidal motion and new constraints of the orbital orientation of the system. Through a detailed analysis of all of these, we derive new constraints on the masses of the pulsar and companion. The mass of the pulsar (1.667 ± 0.021 M) adds to a growing population of NSs with masses which significantly exceed 1.4 M, and rules out some of the proposed equations of state of superdense matter. This high mass also suggests that the pulsar was recycled.

Our observations with the Very Large Telescope reveal shifts in the spectral lines of the suspected companion to PSR J1903+0327, which show that it has the 95-d orbital motion expected for the pulsar’s binary companion. This star has characteristics very similar to the Sun; making this system unique among all known binary MSPs. A detailed discussion of the results of the radio timing and optical spectroscopy rules out most proposed formation mechanisms. We are then forced to accept the only remaining possibility, i.e. that the system originated from a triple system with a compact inner component and the presently observed MS companion as the outer component. The inner compact binary then formed the observed MSP, with the elimination of the former donor star (i.e., evolution towards an ‘isolated’ MSP), or its ejection from the system.

Detailed stellar evolution simulations of the formation of this system might provide new insights on how MSPs form. As an example, the mass-loss and kick velocities associated with the supernova event that formed this NS must have been small enough to allow the whole system to remain bound. Such a gentle formation would produce a small recoil velocity for the whole system; this is consistent with the small peculiar velocity observed at present. This would generally favour alternative small-kick formation mechanisms, like electron capture supernovae. This means that a dynamical study of this system could lead to a conclusive demonstration of alternative formation channels for NSs in general. Some of these exotic mechanisms (for example, white dwarf mergers) could be especially important for explaining the existence of isolated MSPs.

1
We define a millisecond pulsar as a pulsar with spin period P≤ 20 ms and with a low surface magnetic field, B∼ 108 − 9 G.
3
If this star had ever filled its Roche lobe, the orbit would have been circularized. This is inconsistent with the present large and non-varying eccentricity.
4
We used a 1σ relative uncertainty in the distance of 15 per cent such that it is similar to the range of optically derived distances.

This work was supported by the NSF through a cooperative agreement with Cornell University to operate the Arecibo Observatory. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Pulsar research at UBC is supported by an NSERC Discovery Grant. VMK receives support via an NSERC Discovery Grant, CIFAR, FQRNT, and from the Lorne Trottier Chair and a Canada Research Chair. JWTH is a Veni Fellow of The Netherlands Organisation for Scientific Research (NWO). MAM and DRL are supported by the West Virginia EPSCoR programme and the Research Corporation for Scientific Advancement. JPWV is supported by the European Union under Marie Curie Intra-European Fellowship 236394. DJN was supported by NSF grant AST 0647820 to Bryn Mawr College. We also thank the referee, Peter P. Eggleton, for the careful review and many constructive suggestions.

REFERENCES

Appenzeller
I.
et al.,
1998
,
The Messenger
 ,
94
,
1
Archibald
A. M.
et al.,
2009
,
Sci
 ,
324
,
1411
Arzoumanian
Z.
Joshi
K.
Rasio
F.
Thorsett
S. E.
,
1996
, in
Johnston
S.
Walker
M. A.
Bailes
M.
, eds,
IAU Colloq. 160, Pulsars: Problems and Progress
 .
Astron. Soc. Pac.
, San Francisco, p.
525
Backer
D. C.
Foster
R. F.
Sallmen
S.
,
1993
,
Nat
 ,
365
,
817
Bassa
C. G.
van Kerkwijk
M. H.
Koester
D.
Verbunt
F.
,
2006
,
A&A
 ,
456
,
295
Bhattacharya
D.
van den Heuvel
E. P. J.
,
1991
,
Phys. Rep.
 ,
203
,
1
Cenarro
A. J.
Cardiel
N.
Gorgas
J.
Peletier
R. F.
Vazdekis
A.
Prada
F.
,
2001
,
MNRAS
 ,
326
,
959
Champion
D. J.
et al.,
2008
,
Sci
 ,
320
,
1309
Cordes
J. M.
Lazio
T. J. W.
,
2002
, preprint (0207156)
Cordes
J. M.
et al.,
2006
,
ApJ
 ,
637
,
446
Damour
T.
Deruelle
N.
,
1985
,
Ann. Inst. H. Poincaré (Physique Théorique)
 ,
43
,
107
Damour
T.
Deruelle
N.
,
1986
,
Ann. Inst. H. Poincaré (Physique Théorique)
 ,
44
,
263
Demorest
P. B.
,
2007
, PhD thesis, Univ. California, Berkeley
Demorest
P.
Pennucci
T.
Ransom
S.
Roberts
M.
Hessels
J. W. T.
,
2010
,
Nat
 ,
467
,
1081
Dowd
A.
Sisk
W.
Hagen
J.
,
2000
, in
Kramer
M.
Wex
N.
Wielebinski
R.
, eds,
IAU Colloq. 177, Pulsar Astronomy – 2000 and Beyond
 .
Astron. Soc. Pac.
, San Francisco, p.
275
Eggleton
P. P.
Verbunt
F.
,
1986
,
MNRAS
 ,
220
,
13p
Folkner
W. M.
Williams
J. G.
Boggs
D. H.
,
2008
.
JPL Planetary and Lunar Ephemeris, de421
 .
Technical Report IOM 343R-08-003, NASA Jet Propulsion Laboratory
Freire
P. C. C.
Wex
N.
,
2010
,
MNRAS
 ,
409
,
199
Freire
P. C.
Gupta
Y.
Ransom
S. M.
Ishwara-Chandra
C. H.
,
2004
,
ApJ
 ,
606
,
L53
Freire
P. C. C.
Wolszczan
A.
van den Berg
M.
Hessels
J. W. T.
,
2008a
,
ApJ
 ,
679
,
1433
Freire
P. C. C.
Ransom
S. M.
Bégin
S.
Stairs
I. H.
Hessels
J. W. T.
Frey
L. H.
Camilo
F.
,
2008b
,
ApJ
 ,
675
,
670
Fruchter
A. S.
Stinebring
D. R.
Taylor
J. H.
,
1988
,
Nat
 ,
333
,
237
Girardi
L.
Bressan
A.
Bertelli
G.
Chiosi
C.
,
2000
,
A&A
 ,
141
,
371
Gopakumar
A.
Bagchi
M.
Ray
A.
,
2009
,
MNRAS
 ,
399
,
L123
Halpern
J. P.
Gotthelf
E. V.
,
2010
,
ApJ
 ,
709
,
436
Harris
W. E.
,
1996
,
AJ
 ,
112
,
1487
. Updated version at http://www.physics.mcmaster.ca/resources/globular.html
Hobbs
G. B.
Edwards
R. T.
Manchester
R. N.
,
2006
,
MNRAS
 ,
369
,
655
Horne
K.
,
1986
,
PASP
 ,
98
,
609
Hut
P.
,
1984
,
ApJS
 ,
55
,
301
Hynes
R. I.
,
2002
,
A&A
 ,
382
,
752
Irwin
J.
Aigrain
S.
Bouvier
J.
Hebb
L.
Hodgkin
S.
Irwin
M.
Moraux
E.
,
2009
,
MNRAS
 ,
392
,
1456
Kaplan
D. L.
et al.,
2005
,
PASP
 ,
117
,
643
Kopeikin
S. M.
,
1996
,
ApJ
 ,
467
,
L93
Kozai
Y.
,
1962
,
AJ
 ,
67
,
591
Kramer
M.
et al.,
2006
,
Sci
 ,
314
,
97
Kuijken
K.
Gilmore
G.
,
1989
,
MNRAS
 ,
239
,
571
Lai
D.
Bildsten
L.
Kaspi
V. M.
,
1995
,
ApJ
 ,
452
,
819
Lazaridis
K.
et al.,
2009
,
MNRAS
 ,
400
,
805
Lorimer
D. R.
Kramer
M.
,
2005
,
Handbook of Pulsar Astronomy
 .
Cambridge Univ. Press
, Cambridge
Lorimer
D. R.
,
2008
,
Living Rev. Relativity
 . http://relativity.livingreviews.org/Articles/lrr-2008-8/
Marshall
F. E.
Gotthelf
E. V.
Zhang
W.
Middleditch
J.
Wang
Q. D.
,
1998
,
ApJ
 ,
499
,
L179
Moffat
A. F. J.
,
1969
,
A&A
 ,
3
,
455
Munari
U.
Sordo
R.
Castelli
F.
Zwitter
T.
,
2005
,
A&A
 ,
442
,
1127
Paczyński
B.
,
1990
,
ApJ
 ,
348
,
485
Pfahl
E.
Rappaport
S.
Podsiadlowski
P.
,
2002
,
ApJ
 ,
573
,
283
Phillips
J. A.
,
1993
, in
Phillips
J. A.
Thorsett
S. E.
Kulkarni
S. R.
, eds,
ASP Conf. Ser. Vol. 36
.
Astron. Soc. Pac.
, San Francisco, p.
321
Phinney
E. S.
,
1992
,
Philos. Trans. R. Soc. Lond. A
 ,
341
,
39
Phinney
E. S.
Kulkarni
S. R.
,
1994
,
Ann. Rev. Astr. Ap.
 ,
32
,
591
Podsiadlowski
P.
,
2001
, in
Podsiadlowski
P.
Rappaport
S.
King
A. R.
D’Antona
F.
Burderi
L.
, eds, ASP Conf. Ser. Vol. 229,
Evolution of Binary and Multiple Star Systems
 .
Astron. Soc. Pac.
, San Francisco, p.
239
Podsiadlowski
P.
Rappaport
S.
Pfahl
E. D.
,
2002
,
ApJ
 ,
565
,
1107
Podsiadlowski
P.
Langer
N.
Poelarends
A. J. T.
Rappaport
S.
Heger
A.
Pfahl
E.
,
2004
,
ApJ
 ,
612
,
1044
Ransom
S. M.
Stairs
I. H.
Backer
D. C.
Greenhill
L. J.
Bassa
C. G.
Hessels
J. W. T.
Kaspi
V. M.
,
2004
,
ApJ
 ,
604
,
328
Ransom
S. M.
Hessels
J. W. T.
Stairs
I. H.
Freire
P. C. C.
Camilo
F.
Kaspi
V. M.
Kaplan
D. L.
,
2005
,
Sci
 ,
307
,
892
Smarr
L. L.
Blandford
R.
,
1976
,
ApJ
 ,
207
,
574
Splaver
E. M.
Nice
D. J.
Arzoumanian
Z.
Camilo
F.
Lyne
A. G.
Stairs
I. H.
,
2002
,
ApJ
 ,
581
,
509
Stairs
I. H.
,
2003
,
Living Rev. Relativ.
 ,
6
,
5
Standish
E. M.
,
1998
,
JPL Planetary and Lunar Ephemerides, DE405/LE405, Memo IOM 312.F-98-048, (Pasadena: JPL)
 . http://ssd.jpl.nasa.gov/iau-comm4/de405iom/de405iom.pdf
Tauris
T. M.
van den Heuvel
E. P. J.
,
2006
.
Formation and Evolution of Compact Stellar X-ray Sources
 . p.
623
Taylor
J. H.
,
1992
,
Philos. Trans. R. Soc. Lond. A
 ,
341
,
117
Thorstensen
J. R.
Armstrong
E.
,
2005
,
AJ
 ,
130
,
759
Verbiest
J. P. W.
et al.,
2008
,
ApJ
 ,
679
,
675
Weisberg
J. M.
Nice
D. J.
Taylor
J. H.
,
2010
,
ApJ
 ,
722
,
1030
Wex
N.
,
1998
,
MNRAS
 ,
298
,
67
Wex
N.
Kalogera
V.
Kramer
M.
,
2000
,
ApJ
 ,
528
,
401
Will
C. M.
,
1993
,
Theory and Experiment in Gravitational Physics
 .
Cambridge Univ. Press
, Cambridge