MEASURING THE MASS OF SOLAR SYSTEM PLANETS USING PULSAR TIMING

The technique of pulsar timing can provide precise measurements of the rotational, astrometric, and orbital parameters of a pulsar by modeling the observed pulse times of arrival (TOAs). The basic timing analysis provides a fittable parametric model of delays associated with variations in the Euclidean distance between the pulsar and the Earth (resulting from Earth's orbital motion, the proper motion of the pulsar, and its binary motion), dispersive delays in the interstellar medium, and general relativistic time dilation of clocks in the observatory and pulsar frames and along the propagation path (see, e.g., Edwards et al. 2006). The largest variable delay term is the so-called Roemer delay: the modulation caused by the orbital motion of the Earth relative to the solar system barycenter (SSB). The amplitude of this delay is up to ∼500 s, while pulse TOAs for many pulsars are measurable with an uncertainty of much less than 1 μs. This delay is compensated for by using a numerical solar system ephemeris (e.g., Standish 1998). However, the solar system ephemerides cannot be perfect and, at some level, will introduce systematic effects into the timing process. In addition to their use in pulsar timing, these ephemerides are used to provide guidance information for space missions (in fact, this was the original motivation for their development), and hence there is considerable interest in improving their accuracy.

The measured TOAs, t_i, are related to the rotational phase, ϕ_i, of the pulsar at the time of emission as follows:

$$\begin{equation} \phi _i = \nu {\textsc t}_i + \frac{\dot{\nu }{\textsc t}_i^2}{2} + \cdots, \end{equation} \tag{ 1 }$$

where

$$\begin{equation} {\textsc t}_i = t_i + \frac{(\bm {s_i}+\bm {r_i})\cdot \bm {R}}{c} - \Delta _i. \end{equation} \tag{ 2 }$$

Here, ${\textsc t}_i$ is the time of pulse emission, s_i and r_i are, respectively, the vectors from the SSB to the geocenter and from the geocenter to the observatory at time t_i, and R is a unit vector from the SSB toward the pulsar. Δ_i accounts for numerous other delays not relevant to the present discussion (see, e.g., Edwards et al. 2006). Equation (1) expresses the rotational behavior of the pulsar as a Taylor series, which for most millisecond pulsars receives significant non-stochastic contributions from only the two terms shown. Equation (2) relates the times of emission and reception, explicitly including the variations in light-travel time resulting from the motion of the observatory with respect to the SSB. If the parameters of the timing model are perfect, then ϕ_i is always an integer. The differences between the observed phase and that predicted by the timing model are referred to as the "timing residuals," usually expressed in time units through division by ν. The best-fit timing model is generally that which minimizes the weighted sum of the squared residuals, where the weights are the reciprocals of the squared measurement uncertainties in t_i.

The most commonly used solar system ephemerides for pulsar timing are from NASA's Jet Propulsion Laboratory (JPL). They are constructed by numerical integration of the equations of motion and adjustment of the model parameters to fit data from optical astrometry, astrolabe measurements, observations of transits and occultations of the planets and their rings, radar ranging of the planets, radio astrometry of the planets using very long baseline interferometry, radio ranging and Doppler tracking of spacecraft, and laser ranging of the Moon (Standish 1998). These observations constrain the motion of solar system bodies with respect to the Earth, however they do not tightly constrain the planetary masses. This is reflected in the fact that the planetary/solar mass ratios are normally held fixed in the fit.

If the vector between the observatory and SSB is not correctly determined, then systematic timing residuals will be induced. For instance, if the mass of the Jovian system is in error, then sinusoidal timing residuals with a period equal to Jupiter's orbital period will be induced. The identification of such residuals therefore provides a method to limit or detect planetary mass errors in the solar system ephemeris.

In this Letter, we use data taken as part of the international effort to detect gravitational waves (Manchester 2008; Janssen et al. 2008; Jenet et al. 2009; Hobbs et al. 2010) using an array of pulsars to constrain the masses of the solar system planetary systems. These data sets are described in Section 2 and their analysis is discussed in Section 3. In Section 4 the results are presented, and in Section 5 we discuss the potential of future observations and the constraints on unknown solar system bodies.

The pulsars used in this analysis (listed in Table 1) were chosen from the sample observed as part of the International Pulsar Timing Array project (Hobbs et al. 2010). The four pulsars were selected based upon the precision of their measured TOAs, the magnitude of timing irregularities, and on the length of the data set. The data sets for PSRs J0437−4715, J1744−1134, and J1909−3744 are those published by Verbiest et al. (2009) except for a reweighting as described below.

For each observation of a pulsar, typically of 1 hr duration, the data are folded at the rotation period of the pulsar and summed to produce a single pulse profile of relatively high signal-to-noise ratio. The TOA for each profile was obtained by cross-correlating the profile with a high signal-to-noise ratio template and adjusting the start time of the observation for the phase offset between the template and observed profiles. The psrchive (Hotan et al. 2004) and tempo2 (Hobbs et al. 2006) software packages were used to process the data and to obtain timing solutions.

The data set for PSR J1857+0943 is a combination of the previously published TOA data from the Arecibo telescope (Kaspi et al. 1994) in addition to new data from Arecibo, Parkes, and Effelsberg. This combined data set is over 22 years long. Even though this data set is nearly 10 years longer than the other data sets in our sample, accurate TOAs were not obtained for just over three years during the upgrade of the Arecibo telescope. The lack of useful data means that an arbitrary offset has to be included between the pre- and post-upgrade data sets. This arbitrary offset absorbs low-frequency power in the residuals which reduces the sensitivity of the fit to low-frequency terms.

The uncertainties for the parameters produced by the standard weighted least-squares fit implemented into tempo2 assume that the reduced χ² of the fit is unity. In most pulsar data sets the reduced χ² of the fit is significantly larger than one. There are a number of possible reasons for this, including: radio frequency interference causing subtle shape changes in the profile, variations in the interstellar propagation path, intrinsic variations in the pulse profile or the pulsar rotation rate, instrumental artifacts, errors in the clocks used to timestamp the data, or gravitational waves. Many of these effects have a steep-spectrum or "red" character and manifest approximately as low-order polynomials in the timing residuals. In order to improve the estimate of the TOA uncertainty (and therefore the uncertainty of the parameters in the fit), it is common practice to introduce a multiplier that is applied to the TOA uncertainties at fitting. This is usually determined by fitting a polynomial to "whiten" (i.e., flatten) the residuals and then calculating the multiplier required to bring the reduced χ² to unity. Because of the ad hoc nature of this process and because we are searching for long-period signals (i.e., signals with periods similar to the length of our data sets), we use an improved technique to whiten the data and obtain accurate timing model parameters in the presence of red noise and with poorly known TOA uncertainties. This technique is called "Cholesky whitening" and is summarized briefly in Section 3, but will be described more fully in an upcoming paper.

The position of the SSB in a Euclidean frame can be written as a sum over all solar system bodies (including the Sun), where M_j is the mass of the body and b_j the vector position of the body (where the i subscripts used in Equations (1) and (2) have been dropped for clarity):

$$\begin{equation} \bm {b_B} = \sum _j \bm {b_j}\frac{M_j}{M_T}, \end{equation} \tag{ 3 }$$

where M_T = ∑_jM_j. An erroneous set of masses M'_j = M_j − δ_j leads to an erroneous estimate of the barycenter vector

$$\begin{equation} \bm {b^{\prime }_B} - \bm {b_B} \approx -\sum _j \bm {b_j}\frac{\delta _j}{M_T}, \end{equation} \tag{ 4 }$$

where it has been assumed that δ_j ≪ M_j and, consequently, that ∑_jM'_j ≈ ∑_jM_j. If we take the origin of the reference frame to be at the SSB, then b_B = 0 and

$$\begin{equation} \bm {s^{\prime }} - \bm {s} = -\bm {b^{\prime }_B} \approx \sum _j \bm {b_j}\frac{\delta _j}{M_T}. \end{equation} \tag{ 5 }$$

The error in the model time of emission is then

$$\begin{eqnarray} {\textsc t}^{\prime } - {\textsc t} &=& \frac{(\bm {s}^{\prime } - \bm {s}) \cdot \bm {R}}{c}\qquad\qquad\quad \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &\approx & \frac{1}{cM_\odot } \sum _j \delta _j (\bm {b_j} \cdot \bm {R}), \end{eqnarray} \tag{ 7 }$$

that is, we approximate the effect of a change in the mass of a planet as a relocation of the SSB along the vector from the original SSB to that planet. The tempo2 software package has been modified to include the right side of Equation (7) as additional terms in the model. The modified timing model obtains b_j from a specified version of the JPL series of solar system ephemerides (in this work, DE421, Folkner et al. 2009). The model parameters δ_j measure the difference between the best-fit masses and the values assumed by the chosen solar system ephemeris. Indices j of 1–9 refer to the planets (and Pluto) in ascending order of mean distance from the Sun (note that b₃ ≡ s). Examples of the induced timing residuals resulting from an increase in the Jovian system mass of 5 × 10⁻¹⁰ M_☉ are given in Figure 1 for PSRs J0437−4715 and J1857+0943 (ignoring any effects caused by the fitting procedures).

Zoom In Zoom Out Reset image size

In order to deal with all our data sets, we adapted tempo2 to fit multiple pulsars simultaneously. Fitting for the pulsar specific parameters is based solely on the TOAs of that pulsar, whereas pulsar-independent parameters, such as a planetary mass, are fitted globally over all data sets. This procedure reduces the impact of timing noise in individual pulsars on the derived values for the global parameters.

For data sets whose post-fit residuals have a reduced χ² value close to 1.0, it is possible simply to fit for the planetary system mass. To determine realistic uncertainties for the TOAs, we selected short sections of data (∼30 days long depending upon the sampling) for each pulsar, observatory, and back-end instrument used. Weighted fits to each of these short sections of data gave independent estimates of the correction factors which were subsequently averaged and applied to the data set for that combination. This procedure avoided contamination of the correction factors by non-white noise in the data sets. In the presence of non-white noise, standard fitting procedures lead to biased parameter estimates and underestimated uncertainties (see, e.g., Verbiest et al. 2008).

The red noise in each data set was modeled following the method outlined by Verbiest et al. (2008). For PSR J0437−4715 the model developed by Verbiest et al. (2008) was used, while for PSRs J1744−1134, J1857+0943, and J1909−3744, the red noise was fitted by a power-law function, P ∝ f^−α, with exponents of α = 0.9, 0.7, and 0.55, respectively. These noise models provide two methods to determine the parameter uncertainties. First, conservative estimates of the parameter uncertainties were obtained using a Monte Carlo simulation as described by Verbiest et al. (2008). Second, we implemented a new technique that both whitens the residuals and modifies the function being fitted before obtaining the parameter values and uncertainties using a Cholesky factorization of the data covariance matrix. This procedure, known as "Cholesky whitening," will be fully described in a forthcoming paper.

To test our analysis technique, a new ephemeris was created that had identical parameters to the DE421 ephemeris, except for a small decrease in the mass of Jupiter by 7 × 10⁻¹¹ M_☉. The effect of this change was investigated by simulating TOAs that are predicted exactly by a given timing model and the DE421 ephemeris. These simulated data were then analyzed using tempo2 with the modified ephemeris. The resulting pre-fit residuals show the expected sinusoid at Jupiter's orbital period together with an annual term of about half the amplitude of the Jupiter term. Changing the mass of Jupiter has many secondary effects in the modified ephemeris. These include a slight variation in the Astronomical Unit which leads to the annual sinusoid. This small effect is likely to be undetectable in real data and in any case would be absorbed as an offset in the position of the pulsar by ∼0.1 mas. The tempo2 fitting correctly recovered the simulated offset in Jupiter's mass.

Using the DE421 ephemeris, we have obtained timing residuals for the four pulsars listed in Table 1 and fitted for a mass difference for each of the planetary systems from Mercury to Saturn. The resulting mass measurements are listed in Table 2, where the 1σ uncertainties given in parentheses are in the last quoted digit. All results from this work are consistent with the best current measurements; the number of standard deviations between the masses derived in this work and the best-known masses are given in last column.

The mass measurement for Mars was determined without the use of the PSR J0437−4715 data. A spectral analysis of the data set shows significant power in a broad feature around the period of the Martian orbit which could contaminate a fit for the narrow feature that would indicate an error in the mass of Mars. The simple red-noise model used to calculate the correct uncertainties is not detailed enough to account for this feature and so this data set was not used.

Our current data sets are sensitive to mass differences of approximately 10⁻¹⁰ M_☉, independent of the planet. Consequently, our most precise fractional mass determination is for the Jovian system. We therefore check our result by comparing the Jovian system mass obtained using different subsets of our data. In Figure 2, we show the fitted mass difference compared with the value used for the DE421 ephemeris, 9.5479191563 × 10⁻⁴ M_☉, for each pulsar separately and the weighted mean. For comparison, we also show the best Jovian system mass from the Pioneer and Voyager (Campbell & Synnott 1985) and the Galileo (Jacobson et al. 2000) spacecraft. The results obtained by fitting to individual pulsar data sets show a small scatter around the DE421 mass value with no pulsar showing more than a 2σ deviation. The weighted mean deviates from the best-known measurement by only 1.1σ and has considerably smaller uncertainties than the mass determination derived from Pioneer and Voyager (Campbell & Synnott 1985).

Zoom In Zoom Out Reset image size

While the result presented here for the Jovian system is more precise than the best measurement derived from the Pioneer and Voyager spacecraft by a factor of ∼4, the result from the Galileo spacecraft is still better by a factor of ∼20. For a pulsar timing array of 20 pulsars, regularly sampled every two weeks, with white data sets and an rms timing residual of 100 ns, the uncertainty of the mass estimate decreases with increasing data span such that the mass uncertainty of the Galileo measurement for Jupiter would be reached in ∼7 years of observations; see the solid line in Figure 3. Note that this curve does not follow a simple power-law function because of the fitting procedures that are undertaken when dealing with pulsar data sets. Figure 3 also shows that, with ∼13 years of data, the uncertainty of the current Cassini measurement for Saturn is reached.

Zoom In Zoom Out Reset image size

These predictions rely on pulsar data sets remaining "white" over timescales of a decade or more at very high levels of timing precision. While this has yet to be demonstrated, the indications from recent decade-long data sets (Verbiest et al. 2009) are encouraging.

Analysis of data from current and future spacecraft will produce improved measurements of planetary masses. For example, NASA's New Frontiers Mission to Jupiter, Juno, is expected to reach Jupiter in 2016 and orbit it for more than a year. A major scientific objective of this mission is to probe Jupiter's gravitational field in detail, thereby providing a very precise mass for Jupiter.

While spacecraft measurements are likely to continue to provide the most precise mass measurements for most of the planets, at least for the next decade, it should be noted that the pulsar measurements are independent with different assumptions and sources of uncertainty. Independent methods are particularly important for high-precision measurements where sources of systematic error may not be well understood. Furthermore, while spacecraft such as Juno are very sensitive to the mass of the individual body being orbited, they tell us very little about the satellite masses of that body. Only five of the Jovian and nine of the Saturnian satellites are included in the system mass assumed for DE421 (R. A. Jacobson 2010, private communication). Since the pulsar technique is sensitive to the mass of the entire planetary system, it can provide a measure of the mass undetermined by spacecraft observations.

By combining the pulsar and satellite measurements, it will be possible to test the inverse-square relation of gravity and distance for Jupiter masses and distances between 0.1 and 5 AU. Although no deviations apart from known general relativistic effects are expected, it is important to place limits on such effects where possible.

The pulsar timing technique is also sensitive to other solar system objects such as asteroids and currently unknown bodies, e.g., trans-Neptunian objects (TNOs). Measurements of anomalous period derivatives and binary period derivatives for a number of millisecond pulsars have already been used to place limits on the acceleration of the solar system toward nearby stars or undetected massive planets (e.g., Zakamska & Tremaine 2005; Verbiest et al. 2008). Pulsar timing array experiments with a wide distribution of pulsars on the sky will be sensitive to the dipolar spatial dependence resulting from any error in the solar system ephemeris, including currently unknown TNOs. Any ephemeris error will be distinguishable from the effects of gravitational waves passing over the Earth as the latter have a quadrupolar spatial signature. Limits for unknown masses have also been placed by spacecraft using deviations from their predicted trajectories. Doppler tracking data from the two Pioneer spacecraft were searched for accelerations due to an unknown planet. The anomalous acceleration detected in these data, a_P = (8.7 ± 1.3) × 10⁻¹⁰ m s⁻¹, is attributed to non-gravitational sources (Anderson et al. 2002) and is not detected in planetary measurements (Folkner 2010).

An exciting possibility for the future is the creation of a solar system ephemeris that includes pulsar timing data in the overall fit. Such a fit would be able to determine the masses of the planetary systems while simultaneously fitting for orbital parameters.

We have used the four longest and most precise data sets taken for pulsar timing array projects to constrain the masses of the solar system planetary systems from Mercury to Saturn. In all cases, these measurements are consistent with the best-known measurements. For the Jovian system, our measurement improves on the Pioneer and Voyager spacecraft measurement and is consistent with the mass derived from observations of the Galileo spacecraft as it orbited the planet between 1995 and 2003. Pulsar timing has the potential to make the most accurate measurements of planetary system masses and to detect currently unknown solar system objects such as TNOs. In the future, pulsar timing data can be included in the global solutions used to derive solar system ephemerides, thereby improving their precision.

We thank R. A. Jacobson for useful information about the planetary mass determinations, and F. A. Jenet, S. Oslowski, E. Splaver, K. Xilouris, and X. P. You for assistance with the observations and analysis. We also thank the referee for constructive comments. This work is undertaken as part of the Parkes Pulsar Timing Array project which was supported by R.N.M.'s Australian Research Council Federation Fellowship (FF0348478). G.H. is the recipient of an Australian Research Council QEII Fellowship (DP0878388), D.N. is the recipient of an NSF grant (AST-0647820) and A.L. is the recipient of a National Science Foundation CAREER award (AST-0748580). Pulsar research at UBC is funded by an NSERC Discovery Grant. The Parkes radio telescope is part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a National Facility managed by the Commonwealth Scientific and Industrial Research Organisation.