TESTING THEORIES OF GRAVITATION USING 21-YEAR TIMING OF PULSAR BINARY J1713+0747

Our data set consists of pulse timing observations of PSR J1713+0747 at the Arecibo Observatory, from 1992 through 2013, and at the Robert C. Byrd Green Bank Telescope (GBT), from 2006 through 2013, using several generations of data acquisition systems as described below (Table 1). The first 12 yr of these data (1992 August through 2004 May) were previously reported by Splaver et al. (2005) using the Mark III, Mark IV, and Arecibo-Berkeley Pulsar Processor (ABPP) systems. In this work, we incorporate an extension of the Mark IV data set (reduced using the same process as Splaver et al. 2005); data collected at Arecibo from two newer systems; and data collected at the GBT.

The earliest observations (Mark III) were made using a single receiver operating over a limited bandwidth, and hence intrinsic pulsar behavior cannot be separated from effects of the ISM. Subsequent observations at both telescopes were made using two receivers at widely spaced frequencies (1410 and 2380 MHz at Arecibo; 800 and 1410 MHz at the GBT) to allow for separation of these effects.

The earliest observations (1992–1994) used the Princeton Mark III (Stinebring et al. 1992), which collected dual-polarization data with a filter bank of 32 spectral channels each 1.25 MHz wide. Observations between 1998 and 2004 used the Princeton Mark IV (Stairs et al. 2000) instrument and the ABPP (Backer et al. 1997) system in parallel. The Mark IV system collected 10 MHz passband data using 2 bit sampling. The data were coherently dedispersed and folded at the pulse period offline. The ABPP system sampled voltages with 2 bit resolution and filtered the passband into 32 spectral channels (1.75 MHz per channel and 56 MHz in total for 1410 MHz band; 3.5 MHz per channel and 112 MHz in total for 2380 MHz band), and applied coherent dedispersion to each channel using 3 bit coefficients.

From 2004 to 2011/12, pulsar data were collected with the Astronomical Signal Processor (ASP; Demorest 2007) and its Green Bank counterpart GASP (Demorest 2007). The (G)ASP systems recorded 8 bit sampled ∼64 MHz bandwidth data and applyed real-time coherent dedispersion and pulse period folding. The resulting data contain 2048-bin full-Stokes pulse profiles integrated over 1–3 minutes. When observing with ASP we used 16 channels each 4 MHz wide between 1440 and 1360 MHz, and 16 channels between 2318 and 2382 MHz. When observing with GASP we used 12 channels between 1386 and 1434 MHz, and 16 channels between 822 and 886 MHz. The J1713+0747 ASP/GASP data were also reported in Demorest et al. (2013).

We started using the Green Bank Ultimate Pulsar Processing Instrument (GUPPI; DuPlain et al. 2008) for GBT observations in 2010 and its clone the Puerto-Rican Ultimate Pulsar Processing Instrument (PUPPI) for Arecibo observations in 2012. GUPPI and PUPPI use 8 bit sampling in real-time coherent dedispersion and pulse period folding mode and produce 2048-bin full-Stokes pulse profiles integrated over 10 second intervals. When observing with the 800 MHz receiver at the GBT, we use GUPPI to collect data from 62 spectral channels each 3.125 MHz wide, covering 724–918 MHz in frequency. With the L-band receiver, GUPPI uses 58 spectral channels each 12.5 MHz wide, covering the 1150–1880 MHz band. When observing with the L-band receiver at Arecibo, we use PUPPI to collect data from 1150–1765 MHz using 50 spectral channels each 12.5 MHz wide. When observing with the S band, PUPPI takes data from 1770–1880 MHz and 2050–2405 MHz using 38 spectral channels each 12.5 MHz wide. The spectral bandwidth and resolution provided by GUPPI and PUPPI are crucial for resolving the pulse profile evolution in frequency described in Section 3.4.

We combine the pulse times of arrivals (TOAs) used in Splaver et al. (2005) and those of the later observations, for a data span of 21 years, with a noticeable gap between 1994 and 1998, during the Arecibo upgrade. We used daily averaged TOAs from Splaver et al. (2005), because the original 190 s integration TOAs are not accessible. Data timestamps are derived from observatory masers and retrocorrected to universal Coordinated Time via GPS and then further corrected to the TT(BIPM) timescale using the 2012 version BIPM clock corrections with extrapolations to 2013. The TOAs are measured from the observational data through a series of steps. First the data are folded, as they were being taken, into pulse profiles using an ephemeris known to be good enough for predicting the pulse period for the duration of the observation. The folded profiles from different frequency channels and sub-integrations are often summed together to improve the signal-to-noise ratio of the profiles. Orthogonal polarizations are summed to produce a total-intensity profile. The summed profiles are then compared with a well-measured standard pulse profile from the appropriate frequency band. We employ frequency-domain cross-correlation techniques (Taylor 1992) to determine the phase of the pulse peak relative to the midpoint of the observation. The final TOA of a summed profile is then calculated by adding the mid-observation time and the product of pulse period and the measured peak phase. The flux density of the pulsar in these observations can also be measured by comparing the signal strength in the data with that of a calibration observation taken right before or after the pulsar observation in which a signal with known strength was injected. For the post-upgrade Arecibo data and all the GBT data, the flux density of the calibration signal is calibrated every month by comparing it with an astronomical object of known and constant flux density, in this case, the AGNs J1413+1509 and B1442+09.

The telescopes and data acquisition systems introduce varying degrees of computational and electronic delays into the measured TOAs. Further, in many cases, different standard pulse profiles were used to create TOAs from different data acquisition systems. As a result, the TOAs from different systems differ by small time offsets. When a transition is made from one system to another at a telescope, data are typically collected in parallel during a period of around a year, and those data are used to measure the offset between data taken with the instruments. Here we describe the measured offsets in more detail.

There was no overlap between data taken with Mark III (through 1994) and data taken with Mark IV and ABPP (beginning in 1998). The offset between these two systems was treated as a free parameter when fitting timing solutions to the full data set.

Mark IV and ABPP were used in parallel for J1713+0747 observations between 1998 and 2004. They collected data with different bandwidths (Table 1) and were computed using different profile templates. To align the ABPP ToAs with Mark IV, we fitted a phase offset between the 1410 MHz TOAs of the two instruments, and found that ABPP ToAs trail those of Mark IV by 0.46791 ± 0.00009 in pulse phase. In the timing modeling, we fix the phase offsets between Mark IV and ABPP to this value in both 1410 and 2300 MHz, and fit an extra time offset for 2300 MHz ABPP TOAs to account for any extra template misalignment.

Mark IV/ABPP and ASP were used in parallel for several epochs. We fit across the overlap data to determine the offset between the 1410 MHz TOAs of these data sets to measure an offset of 2.33 ± 0.10 μs (Mark IV trailing ASP). The offset between the 2300 MHz TOAs of Mark IV and ASP is treated as a free parameter in the full timing solution.

For the offsets between ASP and PUPPI at Arecibo, and GASP and GUPPI at the GBT, analysis of simultaneous observations of many pulsars, including both pulsar signals and radiometer noise, were used to measure very precise offsets between the instruments at each observatory; details will be given in a forthcoming paper (Arzoumanian et al. 2015). Further, the same standard pulse profiles were used in any given band. Thus these data sets form a continuous 9-year collection of TOAs with no arbitrary offsets.

Finally, the offset between the Arecibo 1410 MHz TOAs and GBT Bank 1410 MHz TOAs was treated as a free parameter when fitting timing solutions to the full data set. The offset between the 800 and 1410 MHz GBT TOAs and the offset between the 1410 and the 2300 MHz Arecibo TOAs are also fitted as free parameters.

The date span, number of observation epochs, and specifications of the systems are listed in Table 1.

The noise model used in this analysis is a parameterized model that is a function of several unknown quantities describing both correlated and uncorrelated noise sources. Uncorrelated noise is independent from one TOA to another, while the correlated noise is not. For instance, the template matching error mostly due to radiometer noise are uncorrelated in time, but the pulse jitter noise (Cordes & Shannon 2010), which affects the multi-frequency TOAs measured simultaneously across the band, is correlated in time. There is also time-correlated noise, such as red timing noise that is correlated from epoch to epoch. Among the various types of noise only the template matching error σ can be estimated when we compute the TOAs. Other sources of noise must be modeled separately. The solution to this problem has been discussed extensively in van Haasteren & Levin (2013), Ellis (2013), van Haasteren & Vallisneri (2014, 2015), Arzoumanian et al. (2014), and J. A. Ellis (PhD thesis 2014). In this section, we summarize the noise model used in this analysis (see, e.g., Arzoumanian et al. 2015, for more details).

To model uncorrelated noise, we use the standard EFAC and EQUAD parameters for each backend/receiver system (e.g., PUPPI backend with L-band receiver). These parameters simply re-scale the original TOA uncertainties

$$\begin{eqnarray}&&{\sigma }_{i,k}\to {E}_{k}{({\sigma }_{i,k}^{2}+{Q}_{k}^{2})}^{1/2},\end{eqnarray} \tag{ 1 }$$

where E_k and Q_k denote the EFAC and EQUAD parameters, respectively, and the subscript i is the TOA number and the subscript k denotes the backend/receiver system.

To model correlated noise we use the new ECORR parameter and a power-law red noise spectrum. The ECORR term describes short timescale noise that is completely correlated for all TOAs in a given observation but completely uncorrelated between observations. This term could be described as pulse phase jitter (Cordes & Shannon 2010) but could also have other components. The ECORR term manifests itself as a block diagonal term in the noise covariance matrix where the size of the blocks is equal to the number of TOAs in a given observation. The exact details of the implementation are described in Arzoumanian et al. (2015); however, the term essentially acts as a observation-to-observation variance. Finally, we model the red noise as a stationary Gaussian process that is parameterized by a power spectrum of the form

$$\begin{eqnarray}&&P(f)={A}_{\mathrm{red}}^{2}{(\displaystyle \frac{f}{{f}_{\mathrm{yr}}})}^{{\gamma }_{\mathrm{red}}},\end{eqnarray} \tag{ 2 }$$

where ${A}_{\mathrm{red}}$ is the amplitude of the red noise in μs ${\mathrm{year}}^{1/2}$ and ${\gamma }_{\mathrm{red}}$ is the spectral index.

These noise parameters are included in a joint likelihood that contains all timing model parameters. For the purposes of noise modeling, we analytically marginalize over the linear timing model parameters and explore the space of noise parameters via MCMC. We then use the MCMC results to determine the maximum likelihood noise parameters which are subsequently used as inputs to TEMPO/TEMPO2 GLS fitting routines. In our noise model we include EFAC, EQUAD, and ECORR parameters for data collected by different backend and receiver systems. However, we do not model EFAC and ECORR of the Mark III, Mark IV, and ABPP data because there were not enough TOAs in the legacy data set to constrain both EFAC and EQUAD. Furthermore, there was only one TOA per observation so we cannot constrain the observation-correlated noise modeled by ECORR. Instead, we set EFAC values to 1 and ECORR to 0 for these data sets, and use only EQUAD to model the white noise in these observations. The maximum likelihood values of the white noise parameters are presented in Table 4.

Shannon & Cordes (2012) studied the pulse arrival times from a single long exposure of PSR J1713+0747, and found that this pulsar's single pulses showed random jitter of $\simeq 26\;\mu {\rm{s}}$. A similar result of $\simeq 27\;\mu {\rm{s}}$ was found by Dolch et al. (2014) from a more recent study using a 24 hr continuous observation of PSR J1713+0747 conducted with major telescopes around the globe. Therefore, by averaging many pulses collected in the typical ∼20 minutes NANOGrav observation, one expects $\sim 27\;\mu {\rm{s}}/\sqrt{1200\nu }$ = 51 ns of jitter noise. Tables 2 and 3 show the best-fit timing parameters before and after we applied our noise model to the data. It is clear that the jitter-like observation-correlated noise affected the arrival time of the pulses, such that some timing parameters changed significantly after including the jitter model. The optimal jitter parameters (ECORR, as shown in Table 4) from our noise modeling are mostly consistent with the prediction from Shannon & Cordes (2012), with some of them being higher. This could be due to the covariance between the jitter parameters and the EQUAD parameters.

In Figure 1 we show the red noise realization based on our best noise model (Table 4) and compare it to the post-fit residuals of a TEMPO GLS fit. The bottom panels of Figure 1 show the one- and two-dimensional posterior probability plots of the red noise. This noise model describes the data well as we can see in Figure 2 in which the maximum likelihood realizations of red and jitter noise are subtracted out. We see from the figure that both the high and low frequency residuals (with red and jitter noise realizations subtracted) are white (described by our EFAC and EQUAD parameters) and the weighted residuals follow a zero-mean unit-variance Gaussian distribution. We do note that the normalized residuals do not seem to follow exactly the Gaussian distribution outside of the 3σ range, this affects only $\lesssim 0.3$% of the TOAs and will not significantly affect the results presented here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Red noise was previously reported for PSR J1713+0747 by Splaver et al. (2005), who modeled it using an eighth-order polynomial. The largest residual in the present data set appears in the time period 1999–2005 (Figure 1). This is when dual-frequency observations begin, and it is likely that the red noise model is absorbing unmodeled DM variations in singe-frequency data collected before 1999, thus making physical interpretations of red noise difficult.

We reanalyzed the Splaver et al. (2005) data sets (Mark III, Mark IV, and ABPP data) with the new red noise modeling technique, and show that the timing results from the new method are mostly consistent with those from Splaver et al. (2005), except for α, δ, ν, and $\dot{\nu }$, which probably changed due to the new red noise model, and Keplerian parameters, which probably changed due to the fact that we used TEMPO2's T2 model instead of TEMPO's DD model used by Splaver et al. (2005). We fit the Keplerian and post-Keplerian parameters simultaneously using the T2 model, whereas Splaver et al. (2005) fit only for the Keplerian parameters while holding the post-Keplerian parameters to their best-fit values using the DD model. Therefore, the uncertainties reported in Splaver et al. (2005) do not reflect the covariance between the two sets of parameters but ours do.

The red noise signal found by our noise model applied to the Splaver et al. (2005) data set is consistent with that found in the 21 year data set, both in terms of the noise parameters (Tables 2, 3, and 5), and in terms of the shapes of the red noise realization (Figure 1), while the red noise modeled by high-order frequency polynomials varies significantly depending on the order of the polynomial or and the observation time span.

Table 5.  The Timing Results from Splaver et al. (2005) and from a Re-analysis of the Splaver et al. (2005) Data Set Using New Red Noise Analysis Technique

Parameter	Splaver et al. (2005)	Red Noise Modela
R.A., α (J2000)	17:13:49.5305335(6)	17:13:49.5305321(6)
decl., δ (J2000)	7:47:37.52636(2)	7:47:37.52626(2)
Spin frequecy ν (s−1)	218.8118439157321(3)	218.811843915731(1)
Spin down rate $\dot{\nu }$ (s−2)	$-4.0835(2)\times {10}^{-16}$	$-4.0836(1)\times {10}^{-16}$
Proper motion in α, ${\mu }_{\alpha }=\dot{\alpha }\mathrm{cos}\delta$ (mas yr−1)	4.917(4)	4.917(4)
Proper motion in δ, ${\mu }_{\delta }=\dot{\delta }$ (mas yr−1)	−3.93(1)	−3.93(1)
Parallax, ϖ (mas)	0.89(8)	0.84(4)
Dispersion measure (pc cm−3)	15.9960	15.9940
Orbital period, ${P}_{{\rm{b}}}$ (day)	67.8251298718(5)b	67.825129921(4)
Change rate of ${P}_{{\rm{b}}}$, ${\dot{P}}_{{\rm{b}}}$ (10−12 s s−1)	0.0(6)	−0.2(7)
Eccentricity, e	0.000074940(1)b	0.000074940(1)
Time of periastron passage, T0 (MJD)	51997.5784(2)b	51997.5790(6)
Angle of periastron, ω (deg)	176.192(1)b	176.195(3)
Projected semimajor axis, x (lt-s)	32.34242099(2)b	32.3424218(3)
Cosine of inclination, $\mathrm{cos}i$	0.31(3)	0.32(2)
Companion mass, Mc (${M}_{\odot }$)	0.28(3)	0.30(3)
Position angle of ascending node, Ω (deg)	87(6)	89(4)
Solar system ephemeris	DE405	DE405
Reference epoch for α, δ, and ν (MJD)	52000	52000
Pulsar mass, ${M}_{\mathrm{PSR}}$ (${M}_{\odot }$)	1.3(2)	1.4(2)
Red noise amplitude (μs ${\mathrm{year}}^{1/2}$)	⋯	0.004
Red noise spectral index	⋯	5.14

Notes.

^aWe used TEMPO2's T2 binary model, which models the Keplerian (${P}_{{\rm{b}}}$, x, e, T₀, and ω) and post-Keplerian orbital elements ($\mathrm{cos}i$, Ω, and m₂ ) simultaneously. ^bSplaver et al. (2005) uses TEMPO's DD model and reports the uncertainties of the Keplerian parameters with the post-Keplerian ones fixed to their bestfit values.

Download table as:  ASCIITypeset image

The DM of a pulsar reflects the number of free electrons between the pulsar and the telescopes and it varies because our sight-line through the turbulent ISM and solar wind is changing as the pulsar, the Sun, the Earth, and the ISM all move with respect to each other. DM variation can affect the timing of high-precision pulsars significantly.

We fit simultaneously with other parameters the time-varying DM using the DMX model in TEMPO. This model fits independent DM values for TOA groups taken within 14 day intervals, except for the L-band-only Mark III TOAs. We grouped the Mark III TOAs together as a single group, because their frequency resolution and timing precision are not sufficient for measuring epoch-to-epoch DM changes.

We turned off the solar wind model for TEMPO and TEMPO2 by setting the solar wind electron density (at 1 AU from the Sun) parameter ${n}_{0}$ to 0 cm⁻³ (the default value is 10 cm⁻³), and used the DMX model to model all DM variations including contribution from solar wind.

Figure 3 shows the measured DM variation of PSR J1713+0747. The sudden dip and recovery of DM around 2008 (MJD 54800) is due to changes either in the ISM or in the solar wind. This DM dip is also observed independently by the Parkes observatory (Keith et al. 2013) and the European Pulsar Timing Array (G. Desvignes et. al. 2015, in preparation).

Zoom In Zoom Out Reset image size

Spectrum analysis of the time variation of flux, pulse arrival phase, and DM have been employed to study the turbulent nature of the ISM (e.g., Cordes et al. 1986; Rickett & Lyne 1990). It has been shown that DM variations of some pulsars are consistent with those expected from an ISM characterized by a Kolmogorov turbulence spectrum (Cordes et al. 1990; Rickett 1990; Kaspi et al. 1994; You et al. 2007; Keith et al. 2013; Fonseca et al. 2014). One can calculate the structure function of the varying DM:

$$\begin{eqnarray}&&{D}_{\phi }(\tau )=(\displaystyle \frac{2\pi K}{{f}^{2}})\langle {[\mathrm{DM}(t+\tau )-\mathrm{DM}(t)]}^{2}\rangle ,\end{eqnarray} \tag{ 3 }$$

where τ is a given time delay, $K=4.148\times {10}^{3}$ MHz² pc⁻¹ cm³ s, and f is the observing frequency in MHz. We expect, under the simplest assumptions, this function to follow a Kolmogorov power law ${D}_{\phi }{(\tau )=(\tau /{\tau }_{0})}^{\beta -2}$, where $\beta =11/3$ and ${\tau }_{0}$ is a characteristic timescale related to the inner scale of the turbulence. The pulsars with DM variations that fit this theory generally have large DM variations on timescale of years. However, PSR J1713+0747 does not show significant long-term DM variation (Figure 3). Conversely, it went through a steep drop and recovery around 2008. If such rapid DM changes are the result of variations in the ISM along the light of sight, such ISM variations do not fit the general characteristics of a Kolmogorov medium.

The term "timing noise" in pulsar timing generally refers to the non-white noise left in the timing residuals. An important contribution to timing noise is expected to come from the pulsar's spin irregularity, i.e., its long-term deviation from a simple linear slow down. Spin irregularity is often significant in younger pulsars, and may be modeled with high-order frequency polynomials (such as $\ddot{\nu }$, where ν is the pulsar's spin frequency). Potential causes of irregular spin behavior include unresolved micro-glitches, internal superfluid turbulence, magnetosphere variations, or external torques caused by matter surrounding the pulsar (Hobbs et al. 2010; Yu et al. 2013; Melatos & Link 2014). These mechanisms could lead to accumulative random perturbations in the pulsar's pulse phase, spin rate, or spin-down rate. Shannon & Cordes (2010) pointed out that one could model these types of timing noise using random walks. Random walks in phase (RW₀) would grow over time (T) proportionally to ${T}^{1/2}$, random walks in ν grow proportionally to ${T}^{3/2}$, random walks in $\dot{\nu }$ grow proportionally to ${T}^{5/2}$. Such spin noise would likely have a steep power spectrum with more power in the lower frequencies. This is considered as one of the main sources of "red" noise in pulsar timing.

The timing noise of radio pulsars has been studied by Cordes & Helfand (1980), Cordes & Downs (1985), Arzoumanian et al. (1994), D'Alessandro et al. (1995), Matsakis et al. (1997), and later by Hobbs et al. (2010) and Shannon & Cordes (2010) with large samples. Matsakis et al. (1997) adopted a generalized Allen Variance (traditionally used in measuring clock stability) to characterize the timing instability of pulsars:

$$\begin{eqnarray}&&{\sigma }_{z}(\tau )=\displaystyle \frac{{\tau }^{2}}{2\sqrt{5}}\langle {c}^{2}{\rangle }^{1/2},\end{eqnarray} \tag{ 4 }$$

where $\langle {c}^{2}\rangle$ denotes the sum of squares of the cubic terms fitted to segments of length τ. Hobbs et al. (2010) found a best-fit scaling model of ${\sigma }_{z}(10\;\mathrm{year})$ from a large sample of pulsars, including canonical pulsars (CPs) and MSPs:

$$\begin{eqnarray}&&{\mathrm{log}}_{10}[{\sigma }_{z}(10\mathrm{year})]=-1.37{\mathrm{log}}_{10}[{\nu }^{0.29}| \dot{\nu }{| }^{0.55}]+0.52,\end{eqnarray} \tag{ 5 }$$

where ν, $\dot{\nu }$ are the pulsar's spin and spin-down rate. We find that Hobbs et al. (2010)'s scaling model (${\sigma }_{z,10\mathrm{year}}^{\mathrm{model}}\simeq 1\times {10}^{-12}$) over-predicted ${\sigma }_{z,10\mathrm{year}}^{\mathrm{measured}}=5\times {10}^{-16}$ for PSR J1713+0747 by more than three orders of magnitude.

Cordes & Helfand (1980) defined a different timing noise characteristic ${\sigma }_{\mathrm{TN},2}^{2}$ based on the root mean square of residuals ${\sigma }_{{\mathcal{R}},2}^{2}$ from a timing fit that does not include any higher order spin parameters like $\ddot{\nu }$. The timing noise term is related to ${\sigma }_{{\mathcal{R}},2}^{2}$:

$$\begin{eqnarray}&&{\sigma }_{{\mathcal{R}},2}^{2}(T)={\sigma }_{\mathrm{TN},2}^{2}(T)+{\sigma }_{W}^{2},\end{eqnarray} \tag{ 6 }$$

where ${\sigma }_{W}^{2}$ is a time-independent term caused by white noise in the data. In this definition, timing noise ${\sigma }_{\mathrm{TN},2}^{2}(T)$ grows bigger over time while white noise stays constant.

Shannon & Cordes (2010) studied the ${\sigma }_{\mathrm{TN},2}^{2}$ from a large sample of CPs and MSPs. They found a scaling model:

$$\begin{eqnarray}&&\mathrm{ln}({\hat{\sigma }}_{\mathrm{TN},2})=1.6-1.4\mathrm{ln}(\nu )+1.1\mathrm{ln}| {\dot{\nu }}_{-15}| +2\mathrm{ln}({T}_{\mathrm{yr}}),\end{eqnarray} \tag{ 7 }$$

where ${\dot{\nu }}_{-15}$ is $\dot{\nu }$ in units of ${10}^{-15}$ s⁻², and ${T}_{\mathrm{yr}}$ is the observation time span in years. This scaling model predicts that, for 21-year timing of PSR J1713+0747, the residual rms without removing timing noise ${\sigma }_{\mathrm{TN},2}^{2}$ would be ∼400 ns. The measured rms of the red noise residual ${\sigma }_{{\mathcal{R}},\mathrm{RN}}^{2}=364$ ns, is consistent with the extrapolation from Shannon & Cordes (2010). The best-fit scaling law also indicates that the residuals of the sampled pulsars ${\hat{\sigma }}_{\mathrm{TN},2}$ seem to grow linearly with ${T}_{\mathrm{yr}}^{2}$. If the timing noise of the sampled pulsars is due to the accumulation of spin noise, and the spin noise is caused by the same physical processes, then this rms growth rate would imply that the spin noise of pulsars has a frequency power spectrum of power-law index ${\gamma }_{\mathrm{red}}\simeq -5$. This spectral index is consistent with the ${\gamma }_{\mathrm{red}}$ from our noise model. This can be seen in the bottom right plot of Figure 1 by noting that a spectral index of −5 is consistent with the posterior at the one-sigma level.

It is inconclusive whether or not the observed red noise can be interpreted as pulsar spin irregularity. Other sources of noise also could have contributed significantly. If we do assume that they are from spin irregularity, the estimated maximum likelihood red noise spectral index of ∼−3 favors that the pulsar spin irregularities come from random walks in either spin phase or spin rate, although other explanations cannot be ruled out due to the substantial uncertainty on the red noise spectral index (Figure 1).

Finally, Shannon & Cordes (2010) showed that the significance of timing noise coming from gravitational wave (GW) background could be estimated as ${\sigma }_{\mathrm{GW},2}\;\approx$ $1.3{A}_{0}{({T}_{\mathrm{yr}})}^{5/3}$ ns, where A₀ is the characteristic strain at f = 1 year⁻¹ and ${T}_{\mathrm{yr}}$ is the observational time span in years. The current best upper limit on GW characteristic strain is $2.4\times {10}^{-15}$ (Shannon et al. 2013), which predicts an upper limit on timing noise of ∼500 ns from GW background. Therefore, we cannot rule out the contribution of GWs in the timing noise.

After removing the dispersion that causes TOA delays proportional to ${f}^{-2}$, where f is the observing frequency, we still see small remaining frequency-dependent residuals from wide-band observations using different instruments and telescopes (Figure 4). It appears that the low-frequency (∼800 MHz) signals lead the high-frequency (1400 and 2300 MHz band) signals by microseconds. The cause of such TOA evolution is not clear. It could be a change in the pulsar's radiation pattern with frequency, or it could be the use of different standard profiles in different frequencies. Pulsar radiation of different frequencies may originate from different parts of the star's magnetosphere, and the radiation region of the pulsars' magnetosphere may be slightly distorted, leading to a frequency-dependent radiation pattern. Pennucci et al. (2014) and Liu et al. (2014) extensively discussed this phenomenon and developed TOA extraction techniques based on phase-frequency 2D pulse profile matching. This technique is not yet applied to our data set.

Zoom In Zoom Out Reset image size

Demorest et al. (2013) allowed an arbitrary offset between TOAs taken with different observing systems and at different frequencies in order to model profile evolution with frequency. However, the number of frequency channels has increased by a factor of ten with the modern wide-band instruments, making it much harder to mitigate profile-frequency evolution using frequency channel offsets. Instead, we used the FD model, a polynomial of the logarithm of frequency: ${\rm{\Delta }}{t}_{\mathrm{FD}}={\displaystyle \sum }_{i=1}^{n}{c}_{i}{(\mathrm{log}\;f)}^{i}$ (Arzoumanian et al. 2015; solid line in Figure 4) to fit for and remove the profile-frequency evolution, where ${\rm{\Delta }}{t}_{\mathrm{FD}}$ is the profile evolution term in units of seconds, f is the observing frequency in unit of GHz, and c_i are the FD model parameters. We employed an F-test with significance value of 0.0027 to determine how many FD parameters are needed to model profile frequency evolution. PSR J1713+0747 only requires n = 4 FD parameters.

The timing model of PSR J1713+0747 has been significantly improved by the 21-year timing effort. Most notably, the pulsar and the companion masses have been more precisely constrained (Tables 2 and 3, Figure 5) through Shapiro delay measurements. The companion's mass ${M}_{{\rm{c}}}=0.286\pm 0.012$ M_⊙ and the pulsar ${M}_{\mathrm{PSR}}=1.31\pm 0.11$ M_⊙ are in good agreement with the previously measured values (Splaver et al. 2005). Furthermore, we have carried out an MCMC run that uses the nonlinear timing model in order to map out any non-Gaussian correlations in parameter space. We find that the nonlinear model gives nearly identical results to the GLS method of TEMPO/TEMPO2. The covariance matrix used in our GLS fitting contains terms come from both the correlated and the uncorrelated noise; therefore, the timing parameter uncertainties we get have taken into account the contribution from the noise model. We note that, without the noise model, the derived pulsar masses would be substantially, and perhaps unrealistically, lower (with ${\rm{\Delta }}{M}_{\mathrm{PSR}}$ ranges from 0.3 to 0.4 ${M}_{\odot }$) than the values with the noise model (Tables 2 and 3), suggesting that correlated noise would significantly impact the accuracy of high precision timing analysis.

The pulsar's mass is compatible with the distribution of pulsar masses in other neutron star-white dwarf systems, and in good agreement with the distribution of pulsar masses found in recycled binaries (Özel et al. 2012; Kiziltan et al. 2013). The precise measurement of neutron star masses may eventually help us understand the properties of matter of extreme density (Demorest et al. 2010; Lattimer 2012; Antoniadis et al. 2013).

In the standard picture of binary evolution, an MSP with a low-mass white dwarf companion must have been spun up through accretion when the white dwarf was a giant star filling its Roche lobe. This should lead to a strong correlation between the binary period and the mass of the white dwarf companion (Rappaport et al. 1995; Tauris & Savonije 1999; Podsiadlowski et al. 2002). Indeed, this picture has been supported by the measurements of several pulsar binary systems (e.g., Kaspi et al. 1994; van Straten et al. 2001; Ransom et al. 2014; Tauris & van den Heuvel 2014). The orbital period and companion mass of PSR J1713+0747 fit this correlation very well, thus supporting the standard MSP evolution theory.

We have observed an apparent change in orbital period from PSR J1713+0747, ${\dot{P}}_{{\rm{b}}}=(0.36\pm 0.17)\times {10}^{-12}$ s s⁻¹ (Tables 2 and 3). This orbital period change is not intrinsic to the pulsar binary, but rather the result of the relative acceleration between the binary and the observer, i.e., the combination of differential acceleration in the Galactic gravitational potential (Damour & Taylor 1991) and the "Shklovskii" effect (Shklovskii 1970) which is caused by the transverse motion of the pulsar binary relative to Earth. We have good measurements of the distance and proper motion of the binary system, which allow us to remove these effects and study the system's intrinsic orbital decay.

The apparent change in orbital period due to differential acceleration in Galactic gravitational potential can be derived from

$$\begin{eqnarray}&&{\dot{P}}_{{\rm{b}}}^{\mathrm{Gal}}=\displaystyle \frac{{A}_{{\rm{G}}}}{c}{P}_{{\rm{b}}}=(-0.10\pm 0.02)\times {10}^{-12}\;{\rm{s}}\;{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 8 }$$

where ${A}_{{\rm{G}}}$ is the line of sight acceleration of the pulsar binary. ${A}_{{\rm{G}}}$ is obtained using Equation (5) in Nice & Taylor (1995), Equation (17) in Lazaridis et al. (2009), the local matter density of Galactic disk around solar system (Holmberg & Flynn 2004) and the Galactic potential model by Reid et al. (2014).

The Shklovskii effect causes ${P}_{{\rm{b}}}$ to change by

$$\begin{eqnarray}&&{\dot{P}}_{{\rm{b}}}^{\mathrm{Shk}}=({\mu }_{\alpha }^{2}+{\mu }_{\delta }^{2})\displaystyle \frac{d}{c}{P}_{{\rm{b}}}=(0.65\pm 0.02)\times {10}^{-12}\;{\rm{s}}\;{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 9 }$$

Therefore, the pulsar's intrinsic orbital decay is ${\dot{P}}_{{\rm{b}}}^{\mathrm{Int}}={\dot{P}}_{{\rm{b}}}^{\mathrm{Obs}}-{\dot{P}}_{{\rm{b}}}^{\mathrm{Shk}}-{\dot{P}}_{{\rm{b}}}^{\mathrm{Gal}}=(-0.20\pm 0.17)\times {10}^{-12}$ s s⁻¹, and is consistent with zero.

Due to the very long ∼68 day orbit, the binary's decay due to the emission of gravitational radiation is expected to be undetectable: ${\dot{P}}_{{\rm{b}}}^{\mathrm{GR}}=-6\times {10}^{-18}$ s s⁻¹ (Lorimer & Kramer 2005). Therefore, the insignificant intrinsic orbital decay rate is entirely consistent with the description of quadrupolar gravitational radiation within General Relativity (GR).

Other than the gravitational radiation, two classical effects could have played a role in ${\dot{P}}_{{\rm{b}}}^{\mathrm{Int}}$. One, ${\dot{P}}_{{\rm{b}}}^{\dot{M}}$, is caused by mass loss in the binary system, and the other, ${\dot{P}}_{{\rm{b}}}^{{\rm{T}}}$, is the contribution from tidal effects. The pulsar and the white dwarf both could lose mass due to their magnetic dipole radiation; the maximum mass loss rate due to this effect can be estimated from the star's rotational energy loss rate. In the case of the pulsar, ${\dot{M}}_{\mathrm{PSR}}=\dot{E}/{c}^{2}$, which is measurable through the spin down rate of the pulsar. The white dwarf generally loses mass at a much lower rate than the pulsar. Therefore, orbital change due to mass loss can be estimated as ${\dot{P}}_{{\rm{b}}}^{\dot{M}}\sim 1\times {10}^{-14}$ s s⁻¹ (Damour & Taylor 1991; Equations (9) and (10) of Freire et al. 2012b). This is an order of magnitude smaller than the measured uncertainties on ${\dot{P}}_{{\rm{b}}}^{\mathrm{Int}}$. The tidal effect in this binary system is expected to be ${\dot{P}}_{{\rm{b}}}^{{\rm{T}}}\ll 3\times {10}^{-14}$ s s⁻¹ based on the most extreme scenarios (the white dwarf spins at its break-up velocity and the tidal synchronizing timescale equals the characteristic age of the pulsar; see Equation (11) in Freire et al. 2012b and references therein). Both of these extra terms are much smaller than the observed uncertainties on ${\dot{P}}_{{\rm{b}}}^{\mathrm{Int}}$.

Based on the measurement of "excess" orbital period change ${\dot{P}}_{{\rm{b}}}^{\mathrm{exc}}={\dot{P}}_{{\rm{b}}}^{\mathrm{Int}}-{\dot{P}}_{{\rm{b}}}^{\dot{M}}-{\dot{P}}_{{\rm{b}}}^{{\rm{T}}}-{\dot{P}}_{{\rm{b}}}^{\mathrm{GR}}$, Damour et al. (1988) derived a phenomenological limit for $\dot{G}$ without considering the binding energy of the stars: $\dot{G}/G\simeq {-\dot{P}}_{{\rm{b}}}^{\mathrm{exc}}/(2{P}_{{\rm{b}}})=(1.0\pm 2.3)\times {10}^{-11}$ yr⁻¹ using the timing of binary PSR B1913+16. Since then ${\dot{P}}_{{\rm{b}}}^{\mathrm{exc}}$ of pulsar binaries, including PSR J1713+0747, have been used to constrain $\dot{G}/G$ (Kaspi et al. 1994; Nice et al. 2005; Deller et al. 2008; Lazaridis et al. 2009; Freire et al. 2012b). So far all pulsar observations show $\dot{G}/G$ consistent with being zero, with upper limits largely determined by the uncertainties in orbital period change rate, distance, and proper motions. PSR J1713+0747 has the smallest known ${\dot{P}}_{{\rm{b}}}^{\mathrm{exc}}/(2{P}_{{\rm{b}}})\simeq (-0.5\pm 0.9)\times {10}^{-12}$ yr⁻¹ (Section 4.2) and is particularly useful for constraining the time variability of the gravitational constant.

Nordtvedt (1990), Lazaridis et al. (2009), and Freire et al. (2012b) showed that a generic test of $\dot{G}/G$ can be achieved using pulsar binaries in a more rigorous fashion by incorporating the binding energy of the neutron stars. The binding energy of a compact star changes with the G, resulting in a changing mass, this will also affect the binary orbit. In a generic form, we could characterize this effect using a self-gravity "sensitivity" parameter s_p (Will 1993). The changing G will now change the orbital period of a pulsar binary system (Nordtvedt 1990; Lazaridis et al. 2009):

$$\begin{eqnarray}&&{\dot{P}}_{{\rm{b}}}^{\dot{G}}=-2\displaystyle \frac{\dot{G}}{G}[1-(1+\displaystyle \frac{{m}_{c}}{2M}){s}_{p}]{P}_{{\rm{b}}}.\end{eqnarray} \tag{ 10 }$$

This formalism is more generic in the sense that it incorporates the compactness of the neutron star, but it also assumes that non-perturbative effects are absent and higher-order contributions in the self-gravity sensitivity can be neglected.

Meanwhile, in the framework of an alternative gravitation theory that violates SEP, a binary system may emit dipole gravitational radiation (Will 1993, 2001; Lazaridis et al. 2009; Freire et al. 2012b and references therein). Such effects arise when the two bodies are very different in terms of their self-gravity, i.e., their compactness. Under the aforementioned assumptions of neglecting non-perturbative effects and higher order contributions of self-gravity sensitivity, this extra dipole radiation could lead to an extra orbital change term:

$$\begin{eqnarray}&&{\dot{P}}_{{\rm{b}}}^{{\rm{D}}}\simeq -4\pi \displaystyle \frac{{T}_{\odot }\mu }{{P}_{{\rm{b}}}}{\kappa }_{{\rm{D}}}{S}^{2},\end{eqnarray} \tag{ 11 }$$

(Will 1993; Lazaridis et al. 2009), where ${T}_{\odot }={{GM}}_{\odot }/{c}^{3}=4.925490947$ $\mu {\rm{s}}$, μ is the reduced mass (${m}_{p}{m}_{c}/M$) of the system, ${\kappa }_{{\rm{D}}}$ is a dipole gravitational radiation "coupling constant," and S is the difference between the self-gravity "sensitivities" of the two bodies ($S={s}_{p}-{s}_{c}$; ${s}_{p}\sim 0.1{m}_{{\rm{p}}}/{M}_{\odot }$ according to Damour & Esposito-Farèse 1992; and ${s}_{c}\ll {s}_{p}$). In Einstein's GR ${\kappa }_{{\rm{D}}}=0$—there is no self-gravity induced dipole gravitational radiation, but it is generally not the case in alternative theories that violate the SEP.

PSR J1713+0747 has a wider binary orbit than most other high-timing-precision pulsar binaries, making its ${\dot{P}}_{{\rm{b}}}^{{\rm{D}}}$ very small. Conversely, ${\dot{P}}_{{\rm{b}}}^{\dot{G}}$ is larger when ${P}_{{\rm{b}}}$ is large. This makes PSR J1713+0747 the best pulsar binary system for constraining the effect of the changing gravitational constant $\dot{G}$. Limits on both $\dot{G}$ and ${\kappa }_{{\rm{D}}}$ can be estimated in the same fashion as in Lazaridis et al. (2009): by solving $\dot{G}$ and ${\kappa }_{{\rm{D}}}$ simultaneously from the equation ${\dot{P}}_{{\rm{b}}}^{\mathrm{exc}}={\dot{P}}_{{\rm{b}}}^{{\rm{D}}}+{\dot{P}}_{{\rm{b}}}^{\dot{G}}$ (Equation (29) of Lazaridis et al. 2009) of different pulsars. We applied this method to four pulsars: PSR J0437−4715, PSR J1012+5307, PSR J1738+0333, and PSR J1713+0747 using timing parameters reported in Lazaridis et al. (2009), Freire et al. (2012b), and this work. The resulting confidence region of $\dot{G}$ and ${\kappa }_{{\rm{D}}}$ is shown in Figure 6. We found, at 95% confidence limit, $\dot{G}/G\;=\;$ $(0.6\pm 1.1)\times {10}^{-12}\;{\mathrm{yr}}^{-1}$; ${\kappa }_{{\rm{D}}}\;=\;$ $(-0.9\pm 3.3)\times {10}^{-4}$. This constraint on $\dot{G}$ is more stringent than previous pulsar-based constraints (Freire et al. 2012b), and close to one of the best constraints of this type ($\dot{G}/G=(-0.07\pm 0.76)\times {10}^{-12}$ yr⁻¹) from the Lunar Laser Ranging (LLR) experiment (Hofmann et al. 2010), which measured Earth-moon distance to $\sim {10}^{-11}$ precision through 39 years of laser ranging. Fienga et al. (2014) showed that $\dot{G}/G$ can be constrained to $(0.01\pm 0.18)\times {10}^{-12}$ yr⁻¹ through the analysis of solar system planetary ephemerides. The ${\kappa }_{{\rm{D}}}$ limit, which is not constrained by solar-system tests, is also slightly improved by using PSR J1713+0747. The pulsar-timing $\dot{G}$ and ${\kappa }_{{\rm{D}}}$ limits are particularly interesting in the testing of SEP-violating alternative theories, because they arise from a test using objects of strong self-gravitation. For example, in some classes of the scalar-tensor theories, the effect of $\dot{G}/G$ could be significantly enhanced in pulsar-white dwarf binaries (Wex 2014).

Zoom In Zoom Out Reset image size

The SEP states that the gravitational effect on a small test body is independent of its constitution, and in particular, that bodies of different self-gravitation should behave the same in the same gravitational experiments. This principle is violated in alternative theories of gravitation like the aforementioned Jordan–Fierz–Brans–Dicke scalar-tensor theory. The PSR J1713+0747 binary is an excellent laboratory for testing effects of SEP violation. If the SEP is violated, the neutron star and the white dwarf will be accelerated differently by the Galactic gravitational field, causing the binary orbit to be polarized toward the center of the Galaxy. The excess eccentricity is expected to be (Damour & Schäfer 1991):

$$\begin{eqnarray}&&| {{\boldsymbol{e}}}_{F}| =\displaystyle \frac{1}{2}\displaystyle \frac{{\rm{\Delta }}\;{g}_{\perp }{c}^{2}}{{\mathcal{F}}{\mathcal{G}}{({M}_{\mathrm{PSR}}+{M}_{{\rm{c}}})(2\pi /{P}_{{\rm{b}}})}^{2}},\end{eqnarray} \tag{ 12 }$$

where ${\mathcal{F}}$ is a factor accounts for potential changes in the periastron advance rate due to deviations from GR, ${\mathcal{G}}$ is the effective gravitational constant in the interaction between the pulsar and the white dwarf, and ${g}_{\perp }$ is the projection of Galactic acceleration on the orbital plane and Δ is the dimensionless factor that characterizes the significance of SEP violation. We assume ${\mathcal{F}}{\mathcal{G}}\approx G$ here and after. The Galactic acceleration of the pulsar system is derived from Holmberg & Flynn (2004), Reid et al. (2014).

GR predicts that there is no preferred reference frame in the universe but this may be different in alternative theories. The Parameterized Post-Newtonian (PPN) formalism parameterized possible deviations from GR into a set of parameters (see Will 2014 for the list of them), some of which are associated with the PFE. In this work, we test ${\alpha }_{3}$, one of the PFE-related parameters. If ${\alpha }_{3}\ne 0$, this would lead to both the presence of a PFE and the breaking of momentum conservation. A rotating body would be accelerated perpendicular to its spin axis and its absolute velocity in the preferred reference frame. In a pulsar binary, this effect would cause an excess in eccentricity, which can be estimated by (Damour & Esposito-Farèse 1992; Bell & Damour 1996):

$$\begin{eqnarray}&&| {{\boldsymbol{e}}}_{F}| ={\hat{\alpha }}_{3}\displaystyle \frac{{c}_{p}| {\boldsymbol{w}}| {P}_{{\rm{b}}}^{2}}{24\pi P}\displaystyle \frac{{c}^{2}}{{\mathcal{F}}{\mathcal{G}}({M}_{\mathrm{PSR}}+{M}_{{\rm{c}}})}\mathrm{sin}\beta ,\end{eqnarray} \tag{ 13 }$$

where ${\boldsymbol{w}}$ is the absolute velocity of the binary system relative to the preferred frame of reference, typically taken as that of the cosmic microwave background (CMB), P is the pulsar's spin period, β is the angle between ${\boldsymbol{w}}$ and the spin axis of the pulsar, and ${\hat{\alpha }}_{3}$ is the strong-field version of the PPN parameters ${\alpha }_{3}$. Here ${\boldsymbol{w}}={{\boldsymbol{w}}}_{\odot }+{{\boldsymbol{v}}}_{\mathrm{PSR}}$, where ${{\boldsymbol{w}}}_{\odot }=369\pm 0.9$ km s⁻¹ is the velocity of the solar system relative to the CMB (Hinshaw et al. 2009), and the term ${{\boldsymbol{v}}}_{\mathrm{PSR}}$ is the relative speed of the pulsar to our solar system. ${{\boldsymbol{v}}}_{\mathrm{PSR}}$ is only partially known because we can measure the pulsar system's proper motion on the sky but we cannot measure its line-of-sight velocity (${v}_{{\rm{r}}}$).

Fortunately, many variables in these equations are measurable in the case of the PSR J1713+0747 binary. It is possible to constrain Δ and ${\hat{\alpha }}_{3}$ using Bayesian techniques by assuming certain fiducial priors for ${v}_{{\rm{r}}}$ (Splaver et al. 2005; Stairs et al. 2005; Gonzalez et al. 2011). In our case, we assumed a Gaussian prior for ${v}_{{\rm{r}}}$ centred at zero with a width equal to the system's proper motion speed. Based on our 21-year timing of J1713+0747 alone, we find 95% confidence limits on the violations of SEP and Lorentz invariance ${\rm{\Delta }}\lt 0.01$ and ${\hat{\alpha }}_{3}\lt 2\times {10}^{-20}$, slightly improving the single pulsar limits from earlier data on this pulsar (Splaver et al. 2005; Stairs et al. 2005; Gonzalez et al. 2011). Stronger limits can be found by combining the results from multiple similar pulsar systems (Wex 2000; Stairs et al. 2005; Gonzalez et al. 2011).