# Studying the solar system with the International Pulsar Timing Array

R.N. Caballero, Y.J. Guo, K.J. Lee, P. Lazarus, D.J. Champion, G. Desvignes, M. Kramer, K. Plant, Z. Arzoumanian, M. Bailes, et al.

► To cite this version:

R.N. Caballero, Y.J. Guo, K.J. Lee, P. Lazarus, D.J. Champion, et al. Studying the solar system with the International Pulsar Timing Array. Monthly Notices of the Royal Astronomical Society, 2018, 481 (4), pp.5501-5516. 10.1093/mnras/sty2632. hal-01929067

HAL Id: hal-01929067

https://hal.science/hal-01929067v1

Submitted on 8 Mar 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Studying the Solar system with the International Pulsar Timing Array

R. N. Caballero,¹★ Y. J. Guo,² K. J. Lee,¹,²★ P. Lazarus,¹ D. J. Champion ◦,¹

G. Desvignes,¹ M. Kramer,¹,³ K. Plant,⁴ Z. Arzoumanian,⁵ M. Bailes,⁶ C. G. Bassa,⁷

N. D. R. Bhat,⁸ A. Brazier,⁹,¹⁰ M. Burgay ◦,¹¹ S. Burke-Spolaor,¹²,¹³ S. J. Chamberlin,¹⁴

S. Chatterjee,¹⁰,¹⁵ I. Cognard,¹⁶,¹⁷ J. M. Cordes,¹⁰,¹⁵ S. Dai,¹⁸ P. Demorest,¹⁹ T. Dolch,²⁰

R. D. Ferdman,²¹ E. Fonseca,²² J. R. Gair,²³ N. Garver-Daniels,¹²,¹³ P. Gentile,¹²,¹³

M. E. Gonzalez,²⁴,²⁵ E. Graikou,¹ L. Guillemot,¹⁶,¹⁷ G. Hobbs,¹⁸ G. H. Janssen,⁷,²⁶

R. Karuppusamy,¹ M. J. Keith,³ M. Kerr ◦,²⁷ M. T. Lam,¹²,¹³ P. D. Lasky,²⁸

T. J. W. Lazio,²⁹ L. Levin,³,¹² K. Liu,¹ A. N. Lommen,³⁰ D. R. Lorimer,¹²,¹³

R. S. Lynch,³¹ D. R. Madison,³¹ R. N. Manchester,¹⁸ J. W. McKee ◦,¹,³

M. A. McLaughlin,¹²,¹³ S. T. McWilliams,¹²,¹³ C. M. F. Mingarelli,¹,³² D. J. Nice,³³

S. Osłowski ◦,¹,⁶,³⁴ N. T. Palliyaguru,³⁵ T. T. Pennucci ◦,³⁶ B. B. P. Perera ◦,³

D. Perrodin,¹¹ A. Possenti,¹¹ S. M. Ransom,³¹ D. J. Reardon,¹⁸,²⁸ S. A. Sanidas,³,³⁷

A. Sesana,³⁸ G. Shaifullah ◦,⁷ R. M. Shannon,⁶,³⁹ X. Siemens,⁴⁰ J. Simon,²⁹

R. Spiewak,⁶,⁴⁰ I. Stairs,²⁵ B. Stappers,³ D. R. Stinebring,⁴¹ K. Stovall,¹⁹

J. K. Swiggum,⁴⁰ S. R. Taylor,⁴² G. Theureau,¹⁶,¹⁷,⁴³ C. Tiburzi ◦,¹,³⁴ L. Toomey,¹⁸

R. van Haasteren,²⁹ W. van Straten ◦,⁴⁴ J. P. W. Verbiest,¹,³⁴ J. B. Wang,⁴⁵

X. J. Zhu ◦²⁸,³⁹ and W. W. Zhu¹,⁴⁶

*Affiliations are listed at the end of the paper*

Accepted 2018 September 24. Received 2018 September 5; in original form 2018 May 25

## ABSTRACT

Pulsar-timing analyses are sensitive to errors in the Solar-system ephemerides (SSEs) that timing models utilize to estimate the location of the Solar-system barycentre, the quasi-inertial reference frame to which all recorded pulse times-of-arrival are referred. Any error in the SSE will affect all pulsars, therefore pulsar timing arrays (PTAs) are a suitable tool to search for such errors and impose independent constraints on relevant physical parameters. We employ the first data release of the International Pulsar Timing Array to constrain the masses of the planet–moons systems and to search for possible unmodelled objects (UMOs) in the Solar system. We employ 10 SSEs from two independent research groups, derive and compare mass constraints of planetary systems, and derive the first PTA mass constraints on asteroid-belt objects. Constraints on planetary-system masses have been improved by factors of up to 20 from the previous relevant study using the same assumptions, with the mass of the Jovian system measured at 9.5479189(3) × 10⁻⁴ M☉. The mass of the dwarf planet Ceres is measured at 4.7(4) × 10⁻¹⁰ M☉. We also present the first sensitivity curves using real data that place generic limits on the masses of UMOs, which can also be used as upper limits on the mass of putative exotic objects. For example, upper limits on dark-matter clumps are comparable to published limits using independent methods. While the constraints on planetary masses derived with all employed SSEs are consistent, we note and discuss differences in the associated timing residuals and UMO sensitivity curves.

**Key words:** pulsars: general – methods: data analysis – methods: statistical – ephemerides.

<footnote>* E-mail: caballero.astro@gmail.com (RNC); kjlee@pku.edu.cn (KJL)</footnote>
# 1 INTRODUCTION

Millisecond pulsars (MSPs) are the most stable rotators known to date in the observable Universe. Pulsar timing (see e.g. Lorimer & Kramer 2005) is a powerful technique through which we record the times-of-arrival (TOAs) of the pulses and use a sophisticated model to convert the topocentric TOA, or site arrival time, to the pulse time-of-emission in the pulsar's co-moving reference frame. The success of the model's fit is assessed from the timing residuals – the difference between the observed and the model-predicted TOAs – which capture all the information unaccounted for in the timing model. Timing residuals from contemporary high-precision timing of the brightest and most stable MSPs observed are at levels of a few hundreds of nanoseconds (see e.g. Verbiest et al. 2016).

The high precision with which MSPs can be timed has made them the primary targets for studies of gravity in the (quasi-stationary) strong field regime, primarily through the studies of their orbital behaviour and especially when in tight orbits with other neutron stars (see e.g. Damour 2009, for a review). Additionally, MSPs can be used as reference clocks to study interesting phenomena that affect their TOAs but are extrinsic to their rotational and orbital behaviour. It is self-evident that it is to our advantage to use multiple MSPs to observe such extrinsic phenomena when possible, especially when trying to measure an effect which is expected to affect TOAs from all MSPs and depends on the pulsar's sky position. We refer to such an ensemble of regularly observed MSPs as a pulsar timing array (PTA; Foster & Backer 1990). The primary scientific goal of PTA researchers is the direct detection of low-frequency gravitational waves (GWs), at nHz frequencies, including stochastic GW backgrounds (GWBs). Three collaborations are currently leading these efforts, namely the European Pulsar Timing Array (EPTA; Desvignes et al. 2016), the North-American Nanohertz Observatory for Gravitational Waves (NANOGray; Arzoumanian et al. 2015), and the Parkes Pulsar Timing Array (PPTA; Reardon et al. 2016). These collaborations work together under the International Pulsar Timing Array consortium (IPTA; Verbiest et al. 2016) in an effort to combine data, resources, and expertise to maximize their scientific output.

While the timing model includes the pulsar's rotational, astrometric, and orbital parameters, and accounts for the time-delay effects of the interstellar medium on the pulsed-signal propagation, it is in fact the transformation of the observation site arrival time to the arrival time at the Solar-system barycentre (SSB) that may introduce correlated signals in the TOAs most likely to interfere with the GWB searches. Such correlated signals may arise from errors in the terrestrial time standards and the Solar-system ephemeris (SSE) used to predict the position of the SSB at any given time of interest. The correlated signals from these two types of errors result in monopolar and dipolar spatial correlations, respectively (see Tiburzi et al. 2016), leading to cross-correlations in the timing residuals of pulsar pairs that may resemble those from a GWB, which have their basis on the quadrupole angular correlation pattern caused by GWs (Hellings & Downs 1983). The presence of signals from clock and SSE errors increases the false-alarm probability of GW detection with PTAs (Tiburzi et al. 2016). In principle, these signals are distinguishable from each other if the data are sufficiently informative, and to manage this, it is especially important to increase the number of MSPs contributing to the analysis and to the sampling of the cross-correlation curve (Siemens et al. 2013; Taylor et al. 2016). While we examine methods to minimize these errors and mitigate their effects when searching for GWs in the PTA data, one can also use the data to extract scientific information on topics other than

GWs. In particular, PTA data have been employed to develop an independent pulsar-based time standard (Hobbs et al. 2012) and to constrain the masses of the Solar-system planetary systems (SSPS) (Champion et al. 2010, henceforth CHM10).

The SSEs that we use for pulsar timing are constructed via numerical integrations of the equations of motion for the known Solar-system bodies. These integrations are subject to a wealth of observational data from telescopes, radio and laser ranging, and spacecrafts orbiting the planets and their moons, when available. Such input data also include estimates of the masses of the planets and other important Solar-system bodies. Observationally, it is not the mass, *M*, but the gravitational parameter of the bodies that is determined, i.e. *GM* where *G* is the universal constant of gravitation. This parameter can be determined with much higher precision than the gravitational constant (see e.g. Petit & Luzum 2010), a fact that limits the precision of measurements of *M* in SI units. For this reason, the masses of Solar-system bodies such as the planets are expressed as the ratios of their gravitational parameters to the solar gravitational parameter (heliocentric gravitational constant), *GM*⊙.

New data are added over time, so that newer SSEs are subject to data of better accuracy and observational sampling. Many of the involved parameters are fitted and adjusted while creating the final SSE. As noted in CHM10, while this process gives accurate predictions for the positions of the planetary systems with respect to the Earth–Moon system, they do not manage to constrain the masses much better than the measurements used as initial values. This is reflected by the fact that typically the ratios of the gravitational parameters of the planetary systems with respect to the solar parameter are held fixed during numerical integrations. What changes between SSE versions with respect to the reference planetary masses is either the initial mass values of the planetary systems, for example after new mass estimates by spacecraft fly-bys, and/or the estimate of the solar gravitational parameter, which can be a fitted parameter in the SSE. Therefore, the input planetary masses in principle differ between the various SSE versions. With this in mind, CHM10 search for errors in the masses of the planetary systems, as the most possible errors that pulsar-timing data could identify.

In this paper we focus on extending the work of CHM10 using the first IPTA data release (IPTA DR 1; Verbiest et al. 2016). In addition to improving the constraints on the SSPS masses, we provide the first PTA constraints on the most massive asteroid-belt objects (ABO) and employ a recently published algorithm (Guo, Lee & Caballero 2018, henceforth GLC18) to search for possible unmodelled objects (UMOs) in Keplerian orbits in the Solar system. We also make a quantitative comparison of SSEs provided by two independent groups, namely the Institut de Mécanique Céleste et de Calcul des Éphémérides (IMCCE) and the Jet Propulsion Laboratory (JPL).

The rest of the paper is organized as follows. In Section 2 we briefly overview the IPTA DR 1 and list the MSPs used in the present study and their basic observational properties. In Section 3 we describe the analysis methods, which includes the details of single-pulsar noise analysis and analyses for constraining the mass of the eight planetary systems and ABOs, as well as the masses of UMOs. The results of the analyses are presented in Section 4. We finally discuss scientific implications of our results and our conclusions in Section 5.

# 2 DATA SET: THE IPTA DR 1

IPTA DR 1 is described in Verbiest et al. (2016), and we only give a brief overview in this section. The full data set consists of TOAs from 49 MSPs. Data were collected by the three regional
PTAs over a total time-span of up to 27.1 yr using seven telescopes across the world, namely the Effelsberg Radio Telescope, the Lovell Telescope, the Nançay Radio Telescope, and the Westerbork Synthesis Radio Telescope by the EPTA, the Arecibo Observatory and the Green Bank Telescope by NANOGrav, and the Parkes Radio Telescope by the PPTA. The IPTA DR 1 data set was constructed by combining data that were published in Kaspi, Taylor & Ryba (1994), Demorest et al. (2013), Manchester et al. (2013), Zhu et al. (2015), and Desvignes et al. (2016).

It is important to note that the TOAs from the different telescopes and different studies were calculated with various methods. Although all TOA calculations were based on template-matching methods (Taylor 1992), where each observed profile is cross-correlated with a profile template of arbitrary phase, there are technical differences regarding issues such as the methods to create the pulse-profile templates, and algorithms for optimal template matching. There are also different approaches with regards to the way that the recorded information is used. For example, in some cases the total intensity profiles were used, which are created by summing the flux of all polarization modes, frequency bands, and sub-integrations, while in others cases, one TOA was calculated per frequency band. These choices reflect differences in the sensitivity of instruments over time and analysis methods which have developed to address them. For example, data from a receiver with limited total bandwidth would use total intensity profiles to achieve useful signal-to-noise ratio, while a more modern broadband receiver can achieve sufficient signal-to-noise ratio with sub-bands of the total bandwidth. In this case, one may opt to produce TOAs per sub-band as a way to mitigate, for example, effects of possible evolution of the pulse profile over the observing frequency (see e.g. Xilouris et al. 1996; Kramer et al. 1999), or possible noise that is limited in certain sub-bands (Cordes, Shannon & Stinebring 2016; Lentati et al. 2016). One may also opt to not integrate profiles over time for short-period pulsars in order to better sample the orbit (e.g. Desvignes et al. 2016).

For MSPs for which data from more than one PTA were available, the IPTA data combination increased the time-span of the MSP data, as well as their cadence and observing-frequency coverage. Increased time-span and cadence allows improved sampling of orbits at longer and shorter periods, respectively. They also lead to better characterization of low- and high-frequency noise properties. Noise mitigation is further aided by improved observing-frequency coverage which is particularly crucial in measuring chromatic noise processes related with the turbulent ionized interstellar medium. The combination of data from multiple telescopes, when available, also offers the chance to use individual data sets in the same observing-frequency bands to search for noise due to systematics (Lentati et al. 2016).

The IPTA DR 1 served as a first testing ground for the use of pulsar-timing noise models that were more complex compared to previous studies such as Arzoumanian et al. (2015), Caballero et al. (2016), or Reardon et al. (2016), which only used data from individual PTAs. It was exactly the aforementioned properties of the IPTA DR 1 that motivated the inclusion of additional noise components in the noise analysis presented in Lentati et al. (2016). The analysis was made in particular to attempt to distinguish between noise specific to each pulsar and noise due to systematics in the data of a given observing system, or noise that is associated with a specific observing frequency band. The intent of introducing the latter noise term is to probe chromatic noise that does not follow the dispersive law of cold homogeneous plasma, associated with

<table><caption>**Table 1.** General characteristics of the IPTA DR 1.0 data (Verbiest et al. 2016) for the MSPs used in this study (note that PSR J0437-4715 was not used to derive mass limits of Solar-system bodies; see Section 4.1). For each pulsar we note the total time-span, *T*, the average cadence, the number of telescopes contributing data, and the weighted RMS of the timing residuals (after subtracting the waveform of the DM variations) The residual RMS was derived using the planetary ephemeris DE421.</caption><thead><tr><th>PSR Name (J2000)</th><th><i>T</i> (yr)</th><th>Average cadence (d)</th><th>Number of telescopes</th><th>Residual weighted RMS (μs)</th></tr></thead><tbody><tr><td>J0437-4715</td><td>14.9</td><td>5.1</td><td>1</td><td>0.3</td></tr><tr><td>J0613-0200</td><td>13.7</td><td>4.3</td><td>6</td><td>1.2</td></tr><tr><td>J1012+5307</td><td>14.4</td><td>6.3</td><td>5</td><td>1.7</td></tr><tr><td>J1713+0747</td><td>17.7</td><td>5.1</td><td>7</td><td>0.3</td></tr><tr><td>J1744-1134</td><td>17.0</td><td>8.4</td><td>6</td><td>1.1</td></tr><tr><td>J1909-3744</td><td>10.8</td><td>4.4</td><td>3</td><td>0.2</td></tr></tbody></table>

temporal dispersion measure (DM) variations (see e.g. Keith et al. 2013; Lee et al. 2014).

In this paper we study the timing data from six MSPs in total and employed data from five of these to constrain the masses of Solar-system bodies. The MSPs were selected based on the contribution of each MSP to the overall results as discussed in Section 4.1. The key observational properties of the data for each of these pulsars are presented in Table 1. By comparison to the IPTA DR 1 data, the one change we have made is related to PSR J1713+0747. The large number of TOAs (19 972) would make the current analysis significantly computationally expensive. This large number of TOAs primarily stems from NANOGrav data, which are not averaged over the observing frequency band, resulting in one TOA per frequency channel. To reduce the computational cost for PSR J1713+0747 we employed the TEMPO2 routine AVERAGEDATA and produced an average TOA for each epoch per observing frequency band by summing up all channels across the frequency band.

## 3 ANALYSIS METHODS

We have implemented two methods to study the Solar system with the IPTA DR 1. Both methods rely on searching for residuals induced by the periodic oscillation of the SSB due to the presence of a mass in orbit that is not accounted for by the pulsar timing models. This mass can either be a difference from the real mass of a Solar-system body to that assumed by the SSE, for which we employed the method discussed in Section 3.2, or the mass of a UMO not included in the SSE, for which we employed the method discussed in Section 3.3.

We clarify here that in this study we are only modelling possible errors in the SSE reference masses. We do not examine the effects of positional errors. Under this model, possible small errors in orbital elements could be absorbed in the mass-error parameter. As we noted in Section 1, we expect that mass errors are more likely to be detected first with pulsar timing analysis, however, sensitivity to errors in orbital elements are not excluded. The GLC18 method can also be focused on applying upper limits on orbital parameter of UMOs, but this is beyond the scope of this study. We discuss further work in pulsar-timing research that attempts to extend PTA studies to orbital elements of planets in Section 5.

Prior to discussing the SSE analysis, we first give an overview of the single-pulsar timing and noise analysis.
## 3.1 Single-pulsar timing and noise analysis

As with other applications of PTAs, constraining the masses of known or unknown bodies in orbit around the SSB requires good characterization of the noise in individual pulsar data (see Cordes 2013; Verbiest & Shaifullah 2018, for reviews on sources of noise in pulsar timing), as noise components may have significant power at frequencies related to a planetary orbit. Insufficient accounting of the noise can lead to significant bias on the measured values of the timing parameters and their uncertainties (Coles et al. 2011; van Haasteren & Levin 2013). CHM10 pointed out these effects in the context of constraining planets' masses and specifically did not include one of the four pulsars they used, PSR J0437-4715, when estimating the mass error of Mars. Specifically, CHM10 argued that its noise model was not sufficient to account for spectral features close to the orbital frequency of Mars, and including this pulsar would thus bias the solution for the specific planet.

For the work presented in this paper, for each pulsar we created different timing and noise models for each SSE. The initial phase-coherent timing models were obtained using the timing software TEMPO2 (Hobbs, Edwards & Manchester 2006). TEMPO2 uses a previously derived timing model (which could be as simple as the pulsar discovery position and rotational frequency) and iteratively performs a least-squares fit of the model to the TOAs until the reduced chi-squared of the residuals is minimized. TEMPO2 applies a linearized approximation to calculate the small, linear offsets of model parameters from the pre-fit value (see also Edwards, Hobbs & Manchester 2006). The least-squares fit can be unweighted or weighted according to the TOA uncertainties. Throughout this work, our timing solutions use weighted fits. These initial individual pulsar-timing models do not include parameters related to errors in the SSE or any noise components. We then employed TEMPONEST (Lentati et al. 2014) to perform a Bayesian (simultaneous) timing and noise analysis, with the same noise modelling used, for example, in Caballero et al. (2016). TEMPONEST samples the joint parameter space of the timing and noise parameters using MULTINEST (Feroz, Hobson & Bridges 2009), a Bayesian inference algorithm based on nested sampling (see Skilling 2004), while evaluating the timing model at each point of the parameter space using the TEMPO2 algorithms.

Before proceeding to the correlated-signal analysis, we produced the final noise models employing the analysis package that we use to make the search for errors in the SSE in order to have a consistent mathematical noise-model parametrization. During this last stage, we performed a Bayesian noise analysis while analytically marginalizing over the timing parameters, also using MULTINEST as the sampler. In brief, the noise model consists of the following components:

(i) Uncorrelated noise terms, modelled with a pair of corrections to the TOA uncertainties per observing system (white noise). The TEMPONEST analysis includes an EFAC (for Error FACTor, a multiplicative factor) and an EQUAD (for Error in QUADrature, a factor added in quadrature). The application of these terms attempts to create a timing solution where appropriate relative weights between the different observing systems are given, since TOA uncertainties calculated via template-matching methods do not always fully account for the TOA scatter. EFACs are used to correct underestimation of the uncertainty, for example due to low signal-to-noise ratio of the observed pulse profile, differences in the pulse profile and the template or presence of noise other than white radiometer noise in the profile, at significant levels. EQUADs are primarily used to account for additional scatter in the TOAs due to physical processes such

as pulse phase jitter (e.g. Liu et al. 2012; Shannon et al. 2014). The corrected TOA uncertainty, $\hat{\sigma}$, and initial uncertainty, $\sigma$, are then related as

$$ \hat{\sigma}^2 = (\sigma \cdot \text{EFAC})^2 + \text{EQUAD}^2. \quad (1) $$

During the final noise analysis, we applied a 'global' EFAC per pulsar, to regularize the white-noise level against the other noise components.

(ii) An achromatic (observing-frequency independent) low-frequency stochastic component (red noise) per pulsar, modelled as a wide-sense stationary stochastic process with a one-sided power-law spectrum of the form

$$ S(f) = \frac{A^2}{f} \left( \frac{f}{f_c} \right)^{2\hat{\alpha}}, \quad (2) $$

where $f$ is the Fourier frequency, $f_c = 1 \text{ yr}^{-1}$ is a reference frequency, $A$ is the amplitude in units of time, and $\hat{\alpha}$ is the spectral index. This noise component is added to model primarily physical noise from irregularities in the rotation of the pulsar, often referred to as 'spin noise' (e.g. Kramer et al. 2006; Shannon & Cordes 2010). In the absence of other dedicated model components (see Lentati et al. 2016), this component will also include noise due to possible systematics in the data.

(iii) A chromatic (observing-frequency dependent) low-frequency stochastic component (DM noise). It has the same spectral properties as the red-noise component, but with the restriction that the induced residuals reflect TOA delays that follow the dispersive law of cold homogeneous plasma (e.g. Landau & Lifshitz 1960), i.e. the time delay of a signal at two observing frequencies, $v_1$ and $v_2$, along a line of sight with DM value, $D_l$, is

$$ \Delta T_{D_l} = \kappa \frac{D_l}{\text{pc cm}^{-3}} \left[ \left( \frac{v_1}{\text{GHz}} \right)^{-2} - \left( \frac{v_2}{\text{GHz}} \right)^{-2} \right], \quad (3) $$

where $\kappa = 4.15 \times 10^{-3}$ s.

The power-law power spectra used to describe the stochastic noise components have sharp cut-offs at $f = 1/T$, with $T$ the data span. This cut-off reflects the fact that power at frequencies below $1/T$ is fitted out by the timing model, as discussed in previous works (van Haasteren et al. 2009; Lee et al. 2014). In particular, we fit for the rotational period and period derivative to remove the low-frequency power of the red noise, and the DM first and second derivatives to remove the low-frequency power of the DM-variations noise. As such, a linear and a quadratic term for DM-variations are always implemented in the (deterministic) timing models used in this work.

Finally, we note that the timing model also needs to take into account the dispersive delays from the plasma of the solar wind (You et al. 2007). Our timing models implement the standard TEMPO2 solar-wind model (Edwards et al. 2006), that assumes a spherical distribution of free electrons with a nominal density of $4 \text{ cm}^{-3}$ at 1 au. Deviations of the electron density distribution from this value (e.g. due to solar activity or deviations from the assumed electron density distribution and/or density at 1 au) will induce additional delay signals that become significant when the line of sight to the pulsar is close to the solar disc. The result is then induced residuals with annual signatures, which peak at epochs where the pulsar is at its smallest elongation. Pulsars with low ecliptic latitudes are more susceptible to such effects (only PSR J1744-1134 falls into that category from the pulsars used in this study). Unmitigated solar-wind signals have complex power spectra and show spatial correlations similar to that caused by SSE errors (Tiburzi et al. 2016) so that they could interfere with the sensitivity of PTA data
at high frequencies. More careful modelling of the solar wind is planned for future work. Data from new observing campaigns at lower frequencies (see e.g. Tiburzi & Verbiest 2018) can provide valuable input for better modelling and mitigation of dispersive delays from the solar wind.

## 3.2 Analysis method for known Solar-system bodies

We first discuss our approach in searching for coherent waveforms in the MSPs from possible errors in the SSPS masses assumed in the SSE. We employ a frequentist analysis using a code that implements the method described in CHM10.

The method considers small errors, $\delta m$, in SSPS masses, $m$, so that $\delta m \ll m$. Such errors will induce residuals due to periodic linear shifts in the SSB position with the period of the planetary orbit. In such a small mass-error case, we can neglect higher order effects on the residuals due to the SSB motion. CHM10 examined the extent of secondary effects using a modified version of the DE421 SSE where the mass of Jupiter deviated the real value by $7 \times 10^{-11} M_\odot$ (an amount compatible to the precision that current PTAs can probe the Jovian system mass, as one can see from the results in the next section) and concluded that such effects were negligible in the case of Jupiter after fitting for the timing model. We further investigated these secondary effects with methods similar to the work in CHM10 and reached similar conclusions. The cases of the inner planets, Mercury and Venus, show additional complexity because of the effects on the orbit of the Earth-Moon system that errors in these planetary masses would cause. The induced residuals from such effects, however, fall into different frequencies to the orbital frequencies of the inner planets. Consequently, although a fully dynamical model could make use of such signals as additional information in constraining the planetary masses, we have verified that these signals, if present in the data, do not affect the results and conclusions from the narrow-frequency signal search employed in this work.

In the first-order CHM10 approximation, the induced residuals from the erroneous mass are then only associated with the (Solar-system related) Rømer delay, the geometric vacuum delay of the TOA at the observatory and at the SSB. The induced residuals will reflect the shift in the position of the SSB along the barycentric position vector of the SSPS, $b$, associated with an error in the pulse time-of-emission. For the multipulsar and multi-SSPS case, this error is calculated for each time epoch as

$$ \tau_b^{n,k} \approx \frac{1}{c M_T} \sum_{i,j}^{n,k} \delta m_i (b_i \cdot \hat{R}_j), \quad (4) $$

where indices $i$ and $j$ refer to the $i$-th (out of $n$) SSPS and the $j$-th (out of $k$) pulsar, respectively, $b$ is the barycentric position vector of the SSPS, $\hat{R}_j$ is the unit barycentric position vector of the pulsar (or pulsar binary) barycentre, $c$ is the speed of light, and $M_T$ is the total mass of the Solar system, which was approximated by $M_T \approx M_\odot$.

Since equation (4) is linear, $\tau_b$ can be directly added to the linear timing model of TEMPO2. Although a single pulsar can provide measurements of the $\delta m$ parameters, equation (4) shows how the measurement precision is dependent on the pulsar's sky position and, therefore, better and less biased measurements can be made by fitting for these parameters simultaneously with many pulsars. In CHM10, this was performed using TEMPO2, which was appropriately modified to allow the $\delta m$ parameters to be fitted in a 'global' timing analysis. Such an analysis is a simultaneous fit of the timing models of various pulsars, where a subset of the parameters, which

we call global, are common for all pulsars. In this example, the $\delta m$ parameters are the global parameters.

The covariance matrix for each pulsar, $C$, is constructed using the maximum-likelihood values of the posterior distribution of the Bayesian noise analysis. It is defined as

$$ C = C_w + C_r + C_{dm}, \quad (5) $$

where the constituent matrices are the white-, red- and DM-noise covariance matrices. $C_w$ is a diagonal matrix with the main diagonal populated with the variances of the TOAs (after application of EFACs and EQUADs). The red- and DM-noise covariance matrices are populated by elements defined, respectively, as (Lee et al. 2014)

$$ C_{r,ij} = \int_{1/T}^{\infty} S_r(f) \cos(2\pi f t_{ij}) df, \quad \text{and} \quad (6) $$

$$ C_{dm,ij} = \frac{\kappa^2 \int_{1/T}^{\infty} S_d(f) \cos(2\pi f t_{ij}) df}{v_r^2 v_j^2}. \quad (7) $$

In the above equations, the $i$ and $j$ indices refer to observing epochs, $f$ is the Fourier frequency, $v$ denotes the observing frequency and $t_{ij}$ is the time lag between the corresponding time epochs.

We can now proceed to search for coherent waveforms as predicted by equation (4) via a global timing analysis. During our analysis, apart from the global parameters, for each pulsar we only fitted for a limited number of timing parameters to ensure that the condition numbers of the design matrices (discussed below) were small and matrix inversions are computationally stable. The timing parameters fitted for are the rotational frequency and its derivative, the DM and derivatives (first and second included in timing models), the pulsar position and parallax. The rotational frequency, DM, and their derivatives correlate with low-frequency noise parameters and $\delta m$ parameters related to the planets with the longest periods. Pulsar position and parallax are also significantly affected by changes to the SSEs (see also fig. 1 in Caballero 2018). We have done so after confirming that the timing models were not influenced by this practice.

Using standard linear-algebra methods we fitted the timing parameters denoted with the column matrix, $\epsilon$, as

$$ \epsilon = (A_r^T A_r)^{-1} A_r^T A_q^T C^{-\frac{1}{2}} t, \quad (8) $$

and the corresponding variances are given by

$$ \sigma_{\epsilon}^2 = \operatorname{diag}(\mathbf{J}^{-1}), \quad \mathbf{J} = \mathbf{D}^T \mathbf{C}^{-1} \mathbf{D}. \quad (9) $$

In these equations, $t$ is the column matrix of the timing residuals, $\mathbf{D}$ is the design matrix (calculated with TEMPO2 during the individual pulsar timing analysis), $\mathbf{C}$ is the covariance matrix, and $\mathbf{J}$ is the Fisher-information matrix. $A_q$ and $A_r$ are the Q and R compositions of matrix $A = \mathbf{C}^{-\frac{1}{2}} \mathbf{D}$, respectively. The $T$, $-1$, and $-\frac{1}{2}$ superscripts denote the transpose, inverse, and inverse of the square root of a given matrix, respectively. All matrices are the total matrices, for all pulsars; $t$ is formed by appending all pulsar timing residuals, $\mathbf{C}$ is the block-diagonal matrix of all pulsar covariance matrices, and $\mathbf{D}$ is formed by appending the SSPS-waveform column matrices to the block-diagonal matrix of all pulsar design matrices. In this way, the SSPS-waveforms act as global parameters to the fit.

The columns with the global SSPS $\delta m$ waveforms are calculated using equation (4). The position of the pulsar is known from the timing model. The position vector of the SSPS for a given observing epoch is calculated based on the information for the SSPS
orbits provided by the used SSE. IMCCE and JPL provide libraries that contain modules and functions that read in the data from the ephemerides and calculate the positions and velocities of the SSPSs for given times. IMCCE and JPL provide the CALCEPH<fnref n="1" /> (Gastineau et al. 2015) and SPICE<fnref n="2" /> libraries, respectively. Having confirmed that both libraries give completely consistent results, we used the CALCEPH in all related work, except the calculations regarding mass errors of ABOs, as discussed in Section 4.2.1.

### 3.3 Analysis method for unknown Solar-system bodies

The approximation used in CHM10 can also be employed in the case where instead of errors in the SSE's reference mass of the SSPS, we consider the mass of UMOs, for which we then also need to model the dynamics of their motion. Such an analysis is beneficial for different reasons. First, it gives the potential to PTAs to probe the masses and dynamics of any object in orbit around the SSB and to impose constraints on physical parameters of proposed or hypothetical objects (see Section 4.3), such as Planet Nine (Brown & Batygin 2016) or dark matter in the Solar system with specified mass distributions (Loeb & Zaldarriaga 2005; Pitjev & Pitjeva 2013; Pitjeva & Pitjev 2013). In this study we focus on a simple model which assumes small bodies in Keplerian orbits around the SSB in order to probe to first order the sensitivity of the real PTA data set to orbiting masses in the Solar system. While not specifically applied in order to constrain the parameter space of specific proposed objects, the analysis assumes orbits that approximate those of most Solar-system bodies (excluding perturbations) and the results can serve as a confirmation of our mass constraints on known bodies and as a means to compare the different SSEs at first order.

For this analysis, we implement the algorithm presented in GLC18, which searches for coherent waveforms from bodies in Keplerian orbits around the SSB, in the TOAs of all pulsars. The details of the approach to search for UMOs, including the mathematical framework, the choice of prior distributions and the analysis algorithm, can be found in GLC18. The algorithm solves the dynamical problem of bodies in Keplerian orbits. By neglecting higher-order effects due to the SSB motion as in CHM10 and any perturbations on the UMO from any object except the Sun, the algorithm is currently restricted to searches of small objects and that are not in orbit around a major planet. The dynamical model contains seven unknown parameters, i.e. the mass of the UMO, *m*, and the six Keplerian orbital parameters, i.e. the semimajor axis, *a*, the eccentricity, *e*, the longitude of the ascending node, $\Omega$, the inclination of the orbit, *i*, the argument of perihelion, $\omega$, and the reference phase, $\phi_0$. For a set of values for these parameters, the model determines the barycentric position vector of the UMO, *b*, and uses equation (4) to calculate the induced signal in the TOAs $S(\xi)$, where we use $\xi$ to denote the seven unknown parameters.

The UMO-induced waveform is now an unknown waveform in the data, and no longer part of the timing parameters. The analysis now uses the reduced likelihood (van Haasteren et al. 2009), which is used when solving the problem while analytically marginalizing over the parameters that are not of interest (often referred to as nuisance parameters). In this case, these are all the timing parameters, $\epsilon$. For the multipulsar case, where we search for a coherent waveform *S* in all pulsars, the reduced likelihood function can be written

as

$$ \Lambda \propto \frac{1}{\sqrt{|\mathbf{CC}'|}} \times \exp \left( -\frac{1}{2} \sum_{i,j,I,J} (\mathbf{t}_{I,i} - S(\xi)_{I,i})^{\mathrm{T}} \mathbf{C}'_{I,J,i,j} (\mathbf{t}_{J,j} - S(\xi)_{J,j}) \right), \quad (10) $$

where the *I*, *J* indices denote pairs of pulsars, the *i*, *j* indices denote pairs of time epochs, and $\mathbf{C}' = \mathbf{C}^{-1} - \mathbf{C}^{-1}\mathbf{D}(\mathbf{D})$. We note that one can use the alternative formulation of the likelihood introduced in Lentati et al. (2013). By applying Bayes's theorem, one can proceed to perform Bayesian parameter estimation as

$$ P(\xi|X) \propto P(\xi)\Lambda. \quad (11) $$

In this compact notation, *X* is the data and $\xi$ are all the model parameter we want to sample, that is, the Keplerian orbital parameters of the UMO and the pulsar-noise parameters we opt to fit simultaneously. Therefore, $P(\xi|X)$ is the posterior probability distribution of the parameter(s) of interest, and $P(\xi)$ is the prior probability distribution of the parameter(s). The parameter space is explored using MULTINEST.

The analysis algorithm for UMOs offers flexibility in the analysis, allowing analytical marginalization over the timing parameters and limiting the prior range of orbital parameters or fixing them to a given value. For the work presented in this paper, we analytically marginalize over the timing parameters and simultaneously search over the UMO orbital parameters and pulsar-noise parameters.

Following the same procedure as in GLC18, we first ran an analysis using the least informative priors for the parameters in order to get the posterior distributions from which we can determine whether we have a possible detection of a UMO. These priors are uniform in the log-space for the parameters with dimension and uniform for dimensionless parameters. In the non-detection case, as is the case in all our IPTA DR 1 analyses, we proceeded to a follow-up, upper-limit analysis to determine the data's sensitivity to any given UMO mass at any semimajor axis value. For the upper-limit analysis, we changed the mass priors to uniform in linear space and performed Bayesian inference for a grid of fixed semimajor axis values. We will refer to these upper limits of the mass as a function of the semimajor axis as the mass sensitivity curves.

## 4 ANALYSES AND RESULTS

The analysis with our implementation of the CHM10 method used 10 different SSEs, five from IMCCE (designation 'INPOP') and five from JPL (designation 'DE'). An overview of the SSEs we employed can be found in Table 2. Before proceeding to searches for correlated SSE-error signals across pulsars, we performed some preliminary searches for errors in masses of SSPSs using single-pulsar data to check the effects of the noise model we select and to compare the performance of our implementations of the CHM10 method with that of TEMPO2. We also made a first-order comparison of the effects on pulsar timing from choosing a different SSE during the analysis.

We tested whether using the noise model described in Section 3.1 produced significantly different results than when using the more complex models published in Lentati et al. (2016). In that work, the SSE DE421 was used, so we used this SSE for a proper comparison. We used single-pulsar constraints on $\delta m$ of the SSPSs using TEMPO2, which can use both types of noise models for single-pulsar

<footnote n="1">¹http://www.imcce.fr/fr/presentation/equipes/ASD/inpop/calceph/</footnote>

<footnote n="2">²https://naif.jpl.nasa.gov/naif/toolkit.html</footnote>
<table><caption>**Table 2.** List of SSEs used in the analyses.</caption><thead><tr><th>IMCCE ephemerides</th><th>Reference</th></tr></thead><tbody><tr><td>INPOP06C</td><td>Fienga et al. (2008)</td></tr><tr><td>INPOP08</td><td>Fienga et al. (2009)</td></tr><tr><td>INPOP10E</td><td>Fienga et al. (2013)</td></tr><tr><td>INPOP13C</td><td>Fienga et al. (2014)</td></tr><tr><td>INPOP17A</td><td>Viswanathan et al. (2017)</td></tr><tr><td>JPL ephemerides</td><td>Reference</td></tr><tr><td>DE405</td><td>Standish (1998)</td></tr><tr><td>DE418</td><td>Folkner et al. (2007)</td></tr><tr><td>DE421</td><td>Folkner, Williams & Boggs (2009)</td></tr><tr><td>DE430</td><td>Folkner et al. (2014)</td></tr><tr><td>DE435</td><td>Folkner, Park & Jacobson (2016)</td></tr></tbody></table>

cases to constrain the mass error. This test was useful for investigating whether any of the pulsars had such noise properties that using a simpler noise model would create a significant bias in the multipulsar, correlated search for errors in the SSE input masses of Solar-system bodies. We did not find any statistically significant differences between the $\delta m$ measurements using the different noise models. We then proceeded to compare the single-pulsar results using TEMPO2 and the algorithm described here, implementing the noise model used in this work. We found the $\delta m$ measurements to be consistent using the two different codes.

We carried out a first-order examination of the effects of our choice of SSE during the timing analysis. As the timing residuals are the primary metric of the completeness of the timing model, we compared the residuals' weighted root mean square (RMS) for each pulsar when using different SSEs. The results for six MSPs (see next section for the selection of pulsars) are summarized in Fig. 1. If we assume that the residual RMS will be minimal for the best-performing SSE, the SSE ranking varies for different pulsars, suggesting that the SSE performance is dependent on the sky position. It is known that the differences between the pairs of SSEs have various sky patterns, an effect that can be illustrated using simulated data (see Caballero 2018).

It is important to keep in mind that SSE related residuals can be fitted out by a number of timing parameters if they have power at those frequencies (Blandford, Narayan & Romani 1984). We are aware that this happens with parameters such as the annual term of the position of the pulsar and other astrometric parameters (see e.g. Madison, Chatterjee & Cordes 2013; Wang et al. 2017). Residual signals due to possible SSE imperfections may also be covariant with pulsar noise parameters. As a result, in the absence of independent constraints on pulsar timing parameters, the SSE ranking based on the timing residuals RMS does not necessarily mean overall better accuracy on the data used to construct the SSE. Madison et al. (2013) also demonstrated that the ability of the noise models included in the timing analysis to prevent leakage of residuals in astrometric parameters depends on the total timespan of the pulsar data set. Therefore, a given SSE may perform differently in terms of the residual RMS for pulsars with different time-spans, even when their true noise properties are similar, since de-correlating pulsar noise, SSE residuals and astrometric parameters requires sufficient data length. Additionally, a given SSE may be over- or underperforming by comparison to another SSE for different Solar-system bodies when used in pulsar timing, so the data-span can further influence the overall performance of an SSE.

While at present the differences in the RMS values of the residuals using different SSEs are within the noise-fluctuation levels, it is

clear that without a full account of such effects in the timing model, cross-checking our results using various SSEs makes studies such as the one presented in this paper more meticulous and robust. A direct consequence of the issues discussed above is that a result regarding the constraints on planetary masses becomes more reliable when using pulsars at as many sky positions as possible and with comparable timing precision and overall data quality, when possible.

## 4.1 Selection of pulsars for analysis

The last point to consider before proceeding to the analysis is which pulsars to use. Searching correlated signals with many pulsars is a computationally intensive task. It has thus been common practice to attempt a ranking of the pulsars available, in order to choose those expected to contribute the most to the analysis. The type of signal sought, the noise characteristics of each pulsar as well as the details of each pulsar's data quality (cadence, time-span, observing frequencies, etc.) play crucial roles in the ranking.

We made a single ranking of the pulsars that we used for both analysis methods described in Section 3 so that we are able to directly compare the results of the analysis for modelled and un-modelled Solar-system objects. Our approach was to use the GLC18 Bayesian code described in Section 3.3 to determine the sensitivity curves of the single-pulsar data to UMO masses. For this, we used the SSE DE421. Table 3 shows a breakdown of the average sensitivity in four intervals of semimajor axis, chosen to be equal in logarithmic space. One can see that the relative sensitivity between pulsars can change over the semimajor axis or equivalently over the period of the Keplerian orbit. Given that the time-spans of our pulsar data sets are between the orbital periods of Jupiter and Saturn while the cadence for all pulsars is much shorter than the period of Mercury, we anticipate that the sensitivity of the pulsars in the third semimajor axis interval ($5 < a/au < 17$) is the most impactful to our results. We therefore made a priority list according to the average sensitivity in that interval. Beyond the top six pulsars, the average sensitivity drops significantly and we therefore decided not to use more MSPs for this work.

The set we used to derive mass limits eventually consisted of five pulsars (highlighted in Table 3). Despite the fact that PSR J0437-4715 is fourth in the ranking, we decided to not include it in this analysis. This is because examination of the posterior distribution of the orbital parameters from the single-pulsar Bayesian analysis for UMOs using PSR J0437-4715 revealed possible systematics in the high-frequency regime (i.e. for small values of the semimajor axis) which resulted in the calculated sensitivity curve violating the analytic sensitivity curve (see GLC18 for details on the analytic sensitivity curve). Our analysis revealed systematics in the range of $1 < a/au < 5$ which complicated the upper limit analysis, and worsen the UMO mass upper limits when including this pulsar in the multipulsar analysis, in contrast to the expectation from the pulsar's noise properties and analytical sensitivity curve. To avoid the potential effects and complications due to these systematics, which we reproduced using multiple SSEs, we did not include PSR J0437-4715 in constraining masses of Solar-system bodies. This pulsar is very bright and as such has very small TOA uncertainties, but it is known to suffer from multiple sources of time-correlated noise (see e.g. Lentati et al. 2016), which gave significant effects on our analysis exactly because of the low TOA uncertainties. We remind the reader that, as discussed in Section 3.1, CHM10 assumed their noise model for PSR J0437-4715 was not precise enough around the orbit of Mars (~1.5 au). We will focus on the results
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The RMS of the timing residuals of the six MSPs listed in Table 1, using the 10 SSEs listed in Table 2. The dashed red, dotted green, and dash-dotted blue lines represent RMS values which are 10, 20, and 30 ns larger than the smallest RMS achieved for the MSP in question. Note that PSR J0437−4715 was not included when calculating mass constraints for Solar-system bodies (see Section 4.1).</figcaption></figure>

<table><caption>**Table 3.** Average sensitivity to mass of UMOs in Keplerian orbits in four ranges of the semimajor axis, *a*, for single-pulsar cases. The table reports the sensitivity as the logarithms of the average of the $1\sigma$ upper limits on the mass of UMOs within each semimajor axis range. The MSPs are listed in order of sensitivity (best to worse) in the interval $a \in [5, 17]$. Given the data set's cadence and time-span, this is the interval where the analysis performance is expected to impact mostly on our results. MSPs in boldface were selected to derive the mass constraints of SSPS, ABOs, and UMOs (see discussion in main text).</caption><thead><tr><th rowspan="2">PSR<br>Name<br>(J2000)</th><th colspan="4">log (M/M<sub>☉</sub>)</th></tr><tr><th>a (au)<br>∈[0.4, 1.4]</th><th>a (au)<br>∈[1.4, 5]</th><th>a (au)<br>∈[5, 17]</th><th>a (au)<br>∈[17, 60]</th></tr></thead><tbody><tr><td>J1713+0747</td><td>-9.921</td><td>-9.933</td><td>-8.514</td><td>-5.824</td></tr><tr><td>J1909-3744</td><td>-10.040</td><td>-10.436</td><td>-8.317</td><td>-6.077</td></tr><tr><td>J1744-1134</td><td>-9.337</td><td>-9.520</td><td>-8.200</td><td>-5.720</td></tr><tr><td>J0437-4715</td><td>-9.737</td><td>-9.244</td><td>-8.091</td><td>-5.861</td></tr><tr><td>J1012+5307</td><td>-9.113</td><td>-9.372</td><td>-7.764</td><td>-5.323</td></tr><tr><td>J0613-0200</td><td>-9.323</td><td>-9.600</td><td>-7.645</td><td>-5.135</td></tr></tbody></table>

without PSR J0437−4715, to directly compare the results of the analysis for modelled and unmodelled Solar-system objects. Examining the exact origins of the systematics is beyond the scope of this paper and is left for future work.

## 4.2 Constraints on masses of known Solar-system bodies

We used our implementation of the CHM10 method, as described in Section 3.2 on the five-pulsar subset of IPTA DR 1, using the 10 SSEs noted in Table 2. As explained in Section 3.2, our analysis seeks possible errors in the input masses, assuming that the mass error is small such that only geometric delays of the pulse propagation due to errors in the estimated position of the SSB are significantly affecting the timing residuals. The SSE input values

were taken directly from the header information of the SSEs using the CALCEPH_INSPECTOR tool of the CALCEPH library.

Fig. 2 shows the results of the analysis for all 10 SSEs. The analysis included all planetary systems (excluding the Earth–Moon system). For planets with moons we refer to the position and mass of the system’s barycentre. The results from the various SSEs are statistically consistent. We also observe that despite the fact that most $\delta m$ measurements are consistent with zero near the $1\sigma$ level, for each planet the central values from the various SSEs are not randomly distributed around zero but have consistent, systematic biases, i.e. are either positive or negative. The only exceptions are INPOP17A for Jupiter, and INPOP08 and DE405 for Mars, although this can be compensated by the very small values with respect to the uncertainties. The most likely reason for these systematic biases is that at this given level of timing precision the results are almost completely constrained by the data, rather than from differences between SSEs within the limits of the random noise from the measurements they use as input data.

Fig. 3 shows the distribution of the significance (central value divided over the $1\sigma$ uncertainty) of the measurements. We see that 38.57 per cent of the cases (27 out of 70) show a measurement with a significance above $1\sigma$. This distribution of errors is very close to a Gaussian distribution (where the corresponding percentage would be at most 31.73). The small difference from the expected error distribution can be due to correlations of long-orbital $\delta m$ signals and low-frequency noise, together with the fact that the analysis assumes symmetric uncertainties. That is because when a $\delta m$ signal correlates with noise parameters, its probability distribution may in fact be asymmetric and the uncertainty would be larger on one side of the median value than the other. Full Monte Carlo sampling of the SSE and noise parameters could be implemented in future work to have a better understanding of these correlations.

Results on Saturn and the ice giants are largely inconclusive. The orbital periods of Uranus and Neptune (84 and 165 yr, respectively)
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The derived central values and $1\sigma$ uncertainties for errors on the masses of the planetary systems with respect to each SSE's input values, for analyses using the 10 SSEs listed in Table 2. The figure on the left-hand-side includes the ice giants to emphasize the much larger uncertainties on their derived masses. The figure on the right-hand-side excludes the ice giants for clarity. The SSE superscript is used to denote that each result is tied to the values of the Sun's gravitational mass of the given SSE (see main text, Section 4.2 for details).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The normalized histogram of the distribution of the significance (central value divided over the $1\sigma$ uncertainty) of the SSPS mass measurements. 38.57 per cent of the cases show a significance over unity, compared to the 31.73 per cent expected for a Gaussian distribution. The vertical, black, dashed line indicates the significance of 1. The red, solid line shows the corresponding cumulative distribution.</figcaption></figure>

are more than six times longer than the data time-span, while their masses are more than five time smaller than the mass of Saturn. It is therefore expected, a priori, that our data set will be completely insensitive to any possible small errors in their masses. We include them nevertheless in our analysis, since the uncertainties of $\delta m$ for these planets are a good indication of the goodness of the uncertainties calculated in the presence of time-correlated noise in the pulsar data and the sufficiency of the pulsar noise models we use. In the presence of low-frequency noise, if the models underestimate the

noise levels, one would expect to see significant detections of $\delta m$ for planets with periods longer than the data set's time-span. The results are as expected, with the uncertainties on $\delta m$ of the ice giants being orders of magnitude larger than the rest of the planets.

Using the results of the analysis, we derived the mass constraints of the planetary systems in the Solar system using the IPTA DR 1 and the 10 SSEs employed in this study. The results are summarized in Table 4. Since the solar gravitational parameter is a fitted quantity in the SSEs, and is therefore different in each case, we express all results as ratios of the planetary gravitational parameters as derived using a specific SSE (superscript SSE) with respect to the nominal solar gravitational parameter, $GM_\odot^N = 1.3271244 \times 10^{20} \text{ m}^3 \text{ s}^{-2}$, in compliance with the guidelines from the 2015 IAU Resolution B3³ (Mamajek et al. 2015). We follow this approach in all mass constraints results we present. While at the precision of the current data set this does not cause any differences in the results within the uncertainties, we nevertheless adopt this approach to allow correct comparisons with future results. A first observation is the consistency in the uncertainties, despite fluctuations in the central values of $\delta m$. The IPTA DR 1 data set is sensitive to mass differences of a few times $10^{-11} \text{ M}_\odot$ for systems up to the Jovian, which constitutes a significant improvement from the approximately $10^{-10} \text{ M}_\odot$ reported in CHM10. We note that for the Saturnian system, this sensitivity is approximately $3 \times 10^{-10} \text{ M}_\odot$ while for the ice giants, the sensitivity drops significantly to approximately $10^{-8} \text{ M}_\odot$.

To evaluate our results, we compare them with the results from CHM10 and with the current best estimates<fnref n="4" /> (CBEs) adopted by the

<footnote n="3">³Available at: https://www.iau.org/static/resolutions/IAU2015_English.pdf</footnote>

<footnote n="4">⁴Up-to-date information at: http://maia.usno.navy.mil/NSFA/NSFA_cbe.html</footnote>
<table><caption>**Table 4.** The mass constraints on the planetary systems derived with the IPTA DR 1, using 10 different SSEs, expressed as ratios of their gravitational masses to that of the nominal solar gravitational mass (see main text, Section 4.2 for details). Numbers in brackets indicate the uncertainty in the last digit quoted. All results are consistent at the 1σ level.</caption><thead><tr><th rowspan="2">Solar-system<br>ephemeris</th><th colspan="7">(GM)<sub>IPTA1</sub><sup>SSE</sup>/(GM)<sub>☉</sub><sup>N</sup></th></tr><tr><th>Mercury</th><th>Venus</th><th>Mars</th><th>Jupiter</th><th>Saturn</th><th>Uranus</th><th>Neptune</th></tr></thead><tbody><tr><td>DE405</td><td>1.6600(3) × 10<sup>-7</sup></td><td>2.44782(2) × 10<sup>-6</sup></td><td>3.2271(1) × 10<sup>-7</sup></td><td>9.5479196(4) × 10<sup>-4</sup></td><td>2.858863(3) × 10<sup>-4</sup></td><td>4.368(1) × 10<sup>-5</sup></td><td>5.144(6) × 10<sup>-5</sup></td></tr><tr><td>DE418</td><td>1.6599(3) × 10<sup>-7</sup></td><td>2.44783(2) × 10<sup>-6</sup></td><td>3.2273(1) × 10<sup>-7</sup></td><td>9.5479195(3) × 10<sup>-4</sup></td><td>2.858859(2) × 10<sup>-4</sup></td><td>4.368(1) × 10<sup>-5</sup></td><td>5.143(6) × 10<sup>-5</sup></td></tr><tr><td>DE421</td><td>1.6599(3) × 10<sup>-7</sup></td><td>2.44782(2) × 10<sup>-6</sup></td><td>3.2273(1) × 10<sup>-7</sup></td><td>9.5479195(3) × 10<sup>-4</sup></td><td>2.858860(2) × 10<sup>-4</sup></td><td>4.368(1) × 10<sup>-5</sup></td><td>5.144(6) × 10<sup>-5</sup></td></tr><tr><td>DE430</td><td>1.6599(3) × 10<sup>-7</sup></td><td>2.44782(2) × 10<sup>-6</sup></td><td>3.2272(1) × 10<sup>-7</sup></td><td>9.5479193(4) × 10<sup>-4</sup></td><td>2.858861(3) × 10<sup>-4</sup></td><td>4.367(1) × 10<sup>-5</sup></td><td>5.145(6) × 10<sup>-5</sup></td></tr><tr><td>DE435</td><td>1.6598(3) × 10<sup>-7</sup></td><td>2.44782(2) × 10<sup>-6</sup></td><td>3.2272(1) × 10<sup>-7</sup></td><td>9.5479193(3) × 10<sup>-4</sup></td><td>2.858860(2) × 10<sup>-4</sup></td><td>4.367(1) × 10<sup>-5</sup></td><td>5.147(6) × 10<sup>-5</sup></td></tr><tr><td>INPOP06C</td><td>1.6599(3) × 10<sup>-7</sup></td><td>2.44782(2) × 10<sup>-6</sup></td><td>3.2273(1) × 10<sup>-7</sup></td><td>9.5479194(4) × 10<sup>-4</sup></td><td>2.858863(3) × 10<sup>-4</sup></td><td>4.367(1) × 10<sup>-5</sup></td><td>5.145(7) × 10<sup>-5</sup></td></tr><tr><td>INPOP08</td><td>1.6600(3) × 10<sup>-7</sup></td><td>2.44782(2) × 10<sup>-6</sup></td><td>3.2271(1) × 10<sup>-7</sup></td><td>9.5479193(5) × 10<sup>-4</sup></td><td>2.858863(3) × 10<sup>-4</sup></td><td>4.368(1) × 10<sup>-5</sup></td><td>5.144(7) × 10<sup>-5</sup></td></tr><tr><td>INPOP10E</td><td>1.6599(3) × 10<sup>-7</sup></td><td>2.44782(2) × 10<sup>-6</sup></td><td>3.2272(1) × 10<sup>-7</sup></td><td>9.5479193(4) × 10<sup>-4</sup></td><td>2.858860(3) × 10<sup>-4</sup></td><td>4.367(1) × 10<sup>-5</sup></td><td>5.146(6) × 10<sup>-5</sup></td></tr><tr><td>INPOP13C</td><td>1.6597(3) × 10<sup>-7</sup></td><td>2.44782(2) × 10<sup>-6</sup></td><td>3.2273(1) × 10<sup>-7</sup></td><td>9.5479193(5) × 10<sup>-4</sup></td><td>2.858861(3) × 10<sup>-4</sup></td><td>4.367(1) × 10<sup>-5</sup></td><td>5.149(7) × 10<sup>-5</sup></td></tr><tr><td>INPOP17A</td><td>1.6599(3) × 10<sup>-7</sup></td><td>2.44783(2) × 10<sup>-6</sup></td><td>3.2273(1) × 10<sup>-7</sup></td><td>9.5479189(5) × 10<sup>-4</sup></td><td>2.858862(3) × 10<sup>-4</sup></td><td>4.367(1) × 10<sup>-5</sup></td><td>5.149(7) × 10<sup>-5</sup></td></tr></tbody></table>

International Astronomical Union (IAU) for the planet-moons systems. The CBEs, denoted with the CBE superscript, are selected from the literature and are derived directly from spacecraft data. For the comparison, we also expressed the CBE results with respect to the nominal solar gravitational parameter. We note again, that such an approach does not change our results within the uncertainties due to the current data precision, but we follow this practice to allow better comparisons with future results and follow the recommended best practices by the IAU. Compared to CHM10, the mass constraints have improved by factors of 5.7, 8.5, 20, 6.7, and 4 for the planetary systems of Mercury, Venus, Mars, Jupiter, and Saturn, respectively. In Table 5, we also compare our results with the IAU CBEs. The precision of the mass constraints derived in this study for planetary systems is lower by factors that range from of ~3 for the case of Jupiter, up to ~10³ for Mercury. In the case of Mercury, the large difference reflects the very significant improvement in the planet's gravity field measurements by the MESSENGER spacecraft. The CBEs for Mercury's gravitational mass (Mazarico et al. 2014) are about a factor 10³ more precise than the previous CBEs (Anderson et al. 1987).

### 4.2.1 Asteroid-belt objects

The main asteroid belt hosts small bodies with masses that reach up to order $10^{-10}$ M☉. With the IPTA DR 1 having sensitivity to mass errors of the order $10^{-11}$ to $10^{-10}$ M☉ between the orbits of Mars and Jupiter (see also next section), it is logical to attempt constraining the masses of the largest bodies of the main belt. This is the first time that PTA data are used to derive mass constraints on ABOs. As our data are only beginning to be sensitive to ABO masses, in this work we perform a pilot study and use only one SSE. Future work with more sensitive data can focus more on comparisons between the pulsar-timing constraints on ABO mass using different SSEs. We employed the SSE DE435 together with additional, high-precision positional data for the ABOs from the New Horizons SPICE Data Archive,<fnref n="5" /> which were used for the New Horizons spacecraft mission. These auxiliary data are provided by JPL in the SPICE kernel format and for this reason, for this application we use the SPICE library and tools.

The δm measurements are shown in Fig. 4 and Table 6 presents the mass constraints derived. We produced IPTA mass constraints on the

three ABOs included in the IAU body constants, namely the dwarf planet Ceres and the asteroids Pallas and Vesta and additionally for another two large asteroids, Juno and Hygiea. For both Ceres and Pallas, the IPTA mass constraint is only slightly over an order of magnitude larger than the IAU CBEs. On the other hand, the IPTA precision on the mass of Vesta is five orders of magnitude worse. This is because of a very precise new determination of the asteroid's mass, orbital, and orientation parameters by Konopliv et al. (2014), which increased the precision of the mass measurement by a factor $10^5$ from the previous best estimate. This was achieved by measurements made with radiometric tracking and optical data from the Dawn spacecraft (Russell & Raymond 2011). The Dawn space mission was specifically designed to send the spacecraft in orbit around Ceres and Vesta for detailed studies. We note that although not yet adopted by the IAU, a publication has recently appeared presenting results for Ceres by the Dawn mission, which has also improved the precision of its mass measurement by a factor of 100 (Konopliv et al. 2018). For the asteroids Juno and Hygiea the uncertainty is equal or higher than the mass constraint and therefore we can only assume upper limits of $9 \times 10^{-11}$ and $6 \times 10^{-11}$ M☉ on their masses, respectively, at the 68 per cent confidence level.

## 4.3 Constraints on masses of UMOs

We used the same five-pulsar list as in the analysis for the SSPs and ABOs in the previous section and employed the method outlined in Section 3.3 to conduct the Bayesian analysis to search for UMOs. This is the first time that such an analysis has been conducted using real PTA data. Given the very high consistency in the results produced using the 10 SSEs in the previous section, we opted to focus on three SSEs, namely DE421, DE435, and INPOP17A. The first was chosen for comparison reasons, since it is the SSE used in CHM10 and the IPTA DR 1 data release and noise-analysis papers (Verbiest et al. 2016; Lentati et al. 2016), while the other two were chosen because they are the latest from each SSE family among those used in this study. Table 7 gives an overview of the types of prior probability distributions used for the sampled parameters, as well as the ranges of their values.

We performed a blind orbital analysis, i.e. we fully searched over the UMO mass and orbital parameters. Our analysis was restricted to circular and eccentric orbits. For all three SSE cases we derived a non-detection result, and produced the mass sensitivity curves, which we present in Fig. 5. In Fig. 6 we overplot the results from the three cases for direct visual comparison. The results show how

<footnote n="5">⁵https://ssd.jpl.nasa.gov/x/spk.html</footnote>
<table><caption>**Table 5.** Comparison between the mass constraints on the planetary systems from this work (IPTA1), the CHM10 results and the CBEs adopted by the IAU. Numbers in brackets indicate the uncertainty in the last digit quoted, the different results are expressed in terms of the nominal solar gravitational mass (see main text, Section 4.2 for details). The sensitivity of the methods can be compared via the ratio of their uncertainties ($\sigma$). For IPTA values, we used the case with the highest uncertainty for each planetary system. Where multiple SSE cases gave the same uncertainty, we note the mass constraint derived with the most recent SSE. The IPTA and IAU have the most comparable mass uncertainties in the case of the Jovian system. The largest difference in the case of Mercury.</caption><thead><tr><th>Planetary system</th><th>(GM)<sup>SSE</sup><sub>IPTA1</sub>/(GM)<sup>N</sup><sub>☉</sub></th><th>σ<sub>CHM10</sub>/σ<sub>IPTA1</sub></th><th>(GM)<sup>CBE</sup><sub>IAU</sub>/(GM)<sup>N</sup><sub>☉</sub></th><th>σ<sub>IPTA1</sub>/σ<sub>IAU</sub></th></tr></thead><tbody><tr><td>Mercury</td><td>1.6599(3) × 10<sup>-7</sup></td><td>5.5</td><td>1.66012099(6) × 10<sup>-7</sup></td><td>5.3 × 10<sup>3</sup></td></tr><tr><td>Venus</td><td>2.44783(2) × 10<sup>-6</sup></td><td>8.5</td><td>2.44783824(4) × 10<sup>-6</sup></td><td>50.0</td></tr><tr><td>Mars</td><td>3.2273(1) × 10<sup>-7</sup></td><td>20</td><td>3.2271560(2) × 10<sup>-7</sup></td><td>500.0</td></tr><tr><td>Jupiter</td><td>9.5479189(5) × 10<sup>-4</sup></td><td>6.7</td><td>9.54791898(16) × 10<sup>-4</sup></td><td>3.13</td></tr><tr><td>Saturn</td><td>2.858863(3) × 10<sup>-4</sup></td><td>4.0</td><td>2.85885670(8) × 10<sup>-4</sup></td><td>37.5</td></tr><tr><td>Uranus</td><td>4.367(1) × 10<sup>-5</sup></td><td>n.a.</td><td>4.366249(3) × 10<sup>-5</sup></td><td>333.3</td></tr><tr><td>Neptune</td><td>5.149(7) × 10<sup>-5</sup></td><td>n.a.</td><td>5.151383(8) × 10<sup>-5</sup></td><td>875.0</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** The derived central values and $1\sigma$ uncertainties for errors on the masses of five massive ABOs with respect to the SSE's input values, for an analysis using the DE435 SSE and updated high-precision positional data from the *New Horizons SPICE Data Archive*. The SSE superscript is used to denote that each result is tied to the values of the Sun's gravitational mass of the given SSE (see main text, Section 4.2 for details).</figcaption></figure>

the relative sensitivity of the data at various distances from the SSB changes when using different SSEs. While for semimajor axis values up to the orbit of Mars the three SSEs are in very close agreement, for wider orbits the results are less consistent, with DE421 showing overall higher sensitivity, i.e. giving the lowest upper limits. DE435 and INPOP17A show their biggest differences in the semimajor axis range 4–8 au, i.e. in the asteroid belt, and around the orbit of Jupiter, and become fully consistent for distances beyond 20 au. Table 8 presents the upper limits on the mass of UMOs at selected values of the semimajor axis, for the three SSEs used.

Direct comparison to the results for known SSPs at the same semimajor axis values is only approximate, since this analysis assumes unperturbed Keplerian orbits, in contrast to the analysis in the previous section which follows the exact orbits based on observations. It is nevertheless useful to make the comparison as a

cross-check, since the much larger degrees of freedom in the search for UMOs should always result in worse sensitivity by comparison to that of known bodies for the same semimajor axis values. This is indeed the case in our analysis, with the upper limits from the blind search being ~2–14 times higher. One could also extend the upper-limit analysis to wider orbits in order to retrieve, for example, an upper limit on the mass of Planet Nine. In GLC18, the results using simulated data show that the precision of the IPTA DR 1 is not sufficient to give informative constraints on the mass of Planet Nine. We therefore did not attempt this, but reserve such effort for future work.

As discussed in GLC18, the results from this type of analysis directly provide upper limits on the presence of any type of massive objects in orbit around the SSB. As such, our results are also applicable to more exotic objects such as dark matter clumps (Loeb & Zaldarriaga 2005) or cosmic strings (Blanco-Pillado, Olum & Shlaer 2014). For distances above 2 au from the SSB (where the sensitivity is maximum) we can exclude (with a 68 per cent confidence level) the presence of dark matter clumps (in eccentric, Keplerian orbits) with masses up to $1.2 \times 10^{-11}$ M☉. For distances up to Saturn's orbit (≈9.6 au), the upper limits range between $4 \times 10^{-10}$ and $2 \times 10^{-9}$ M☉ (depending on the used SSE). For comparison, we note that Pitjev & Pitjeva (2013) and Pitjeva & Pitjev (2013) present upper limit of $1.7 \times 10^{-10}$ M☉ for the dark matter mass in the sphere within Saturn's orbit, using independent data and methodology. Their approach searches for perturbations on the orbital motion of planets due to the acceleration by an assumed dark matter distribution in the interplanetary space. This comparison is only indicative, since that work assumes that dark matter has a continuous distribution that is spherically symmetric relative to the Sun, with a fixed central density and exponential drop with increased distance Pitjev & Pitjeva (2013) and Pitjeva & Pitjev (2013). Other density distributions are also discussed, but none of those models assumes clumps as we did in this study. We note that Pitjev & Pitjeva (2013) employ an SSE independent of the ones used in this study. Specifically, they use the EPM2011 (Pitjeva 2013) which is published by the Institute of Applied Astronomy of the Russian Academy of Sciences.

# 5 DISCUSSION AND CONCLUSIONS

In the work described in this paper we have employed previously published methods on a subset of the first IPTA data release in order to constrain the masses of Solar-system bodies using ten different
<table><caption>**Table 6.** Comparisons of the mass constraints for the five most massive ABOs derived in this work with the IAU CBEs. When IAU CBEs are unavailable (noted with * superscript), we use the values from Carry (2012). The IPTA masses were derived using the SSE DE435 and updated high-precision positional data from the *New Horizons SPICE Data Archive*. For the comparison, the different results are expressed in terms of the nominal solar gravitational mass (see main text, Section 4.2 for details). Numbers in brackets indicate the uncertainty in the last digit quoted.</caption><thead><tr><th>Name</th><th>Minor planet category</th><th>(GM)<sup>SSE</sup><sub>IPTA1</sub>/(GM)<sup>N</sup><sub>☉</sub></th><th>(GM)<sup>CBE</sup><sub>IAU</sub>/(GM)<sup>N</sup><sub>☉</sub></th><th>σ<sub>IPTA1</sub>/σ<sub>IAU</sub></th></tr></thead><tbody><tr><td>1 Ceres</td><td>Dwarf planet</td><td>4.8(4) × 10<sup>-10</sup></td><td>4.72(3) × 10<sup>-10</sup></td><td>13.3</td></tr><tr><td>2 Pallas</td><td>Asteroid</td><td>1.4(4) × 10<sup>-10</sup></td><td>1.03(3) × 10<sup>-10</sup></td><td>13.3</td></tr><tr><td>3 Juno*</td><td>Asteroid</td><td>4(5) × 10<sup>-11</sup></td><td>1.37(1) × 10<sup>-11</sup></td><td>500</td></tr><tr><td>4 Vesta</td><td>Asteroid</td><td>1.1(1) × 10<sup>-10</sup></td><td>1.3026846(9) × 10<sup>-10</sup></td><td>1.1 × 10<sup>5</sup></td></tr><tr><td>10 Hygiea*</td><td>Asteroid</td><td>3(3) × 10<sup>-11</sup></td><td>4.3(3) × 10<sup>-11</sup></td><td>100</td></tr></tbody></table>

<table><caption>**Table 7.** Prior types and ranges for the Bayesian analysis to constrain the masses of UMOs. Two sets of priors are shown, one for the blind search of UMOs and one for the mass upper limit analysis (see discussion on priors in Section 3.3).</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Prior range</th></tr><tr><th>Blind search</th><th>Upper-limit analysis</th></tr></thead><tbody><tr><td>m (M<sub>☉</sub>)</td><td>log-uniform in [10<sup>-25</sup>, 10<sup>-5</sup>]</td><td>Uniform in [0, 10<sup>-5</sup>]</td></tr><tr><td>a (au)</td><td>log-uniform in [0.1, 10]</td><td>Fixed in [0.4, 60]</td></tr><tr><td>e</td><td>Uniform in [0, 0.99]</td><td>Uniform in [0, 0.99]</td></tr><tr><td>Ω</td><td>Uniform in [0, 2π]</td><td>Uniform in [0, 2π]</td></tr><tr><td>i</td><td>Uniform in [0, π]</td><td>Uniform in [0, π]</td></tr><tr><td>ω</td><td>Uniform in [0, 2π]</td><td>Uniform in [0, 2π]</td></tr><tr><td>φ<sub>0</sub></td><td>Uniform in [0, 2π]</td><td>Uniform in [0, 2π]</td></tr></tbody></table>

SSEs, five from IMCCE and five from JPL. Using a new computational implementation of the method first described in CHM10, we have derived new mass constraints for the SSPSs, which were found to be statistically consistent using all 10 SSEs. While the biases from the SSE reference values appear consistent for each SSPS, the results appear to be dominated by data noise. Within the uncertainties, our results are in agreement with the CBEs from the IAU which overall have significantly lower uncertainties. For the first time, PTA data were also used to significantly constrain the masses of the most massive ABOs. A Bayesian method from GLC18 was also employed for the first time on real data to provide generic sensitivity limits on the mass of UMOs in the Solar system using pulsar timing.

The new mass constraints on all planetary systems show improvements of factors 4–20 from the last work that used the same method, namely CHM10, emphasizing the fact that increasing the precision, cadence, frequency coverage and time-span of the pulsar-timing data allows for constant improvements of PTA sensitivity to potential errors in SSEs. As such, the IPTA greatly serves this research since the combination of independent data sets improves the data overall in all these aspects. As noted in Caballero (2018), the use of the IPTA combined data improved the sensitivity to planetary masses by factors up to ~4 by comparison to only using EPTA data, when using the same pulsars in both cases. Additionally, the IPTA combined data set also allowed more options with regards to choosing MSPs for the analyses and this study has benefited from using a larger and different sample of pulsars than CHM10. We note that constraining planetary masses with pulsar-timing data helps us cross-check the data quality and pulsar noise models using information on physical properties that are measured completely

independently. Large deviations from the SSE's reference masses or unexplained signals present only in one pulsar's data, which are not detected with multipulsar searches for correlated signals, can indicate insufficiencies of the noise models or possible systematics in the data of a given pulsar.

In this paper we have additionally demonstrated with real data the ability of algorithms that search for UMOs of any type in the Solar system, to provide generic mass sensitivity curves using pulsar timing. While with certain limitations, the GLC18 code applied in this paper highlights differences between SSEs. As we saw in Section 4.3, the main differences in UMO-mass sensitivity curves between DE421, DE435, and INPOP17A appear in the asteroid belt and around Jupiter. While the details of the differences between SSEs are beyond the scope of this study, we note that these results may be due to changes in the way that ABO masses and their perturbations on each other and on Mars are estimated, as well as recent updates in the positional data of Mars, Jupiter, and Saturn that JPL and IMCCE have been implementing (e.g. Folkner et al. 2014; Viswanathan et al. 2017). For example, when estimating perturbations of the ABOs on the orbit of Mars, both DE421 and DE435 used data for the 343 ABOs identified to be dominant. However, while in DE421 only eleven ABO masses were individually calculated (for the rest either the initial values were kept fixed or values were fixed to approximate values derived densities assumed per taxonomic class), for DE435 (also the case for DE430) the individual masses were calculated for all 343 ABOs. Future work with more precise data sets could focus more on the effects of such difference on pulsar timing and applications.

As noted in CHM10, since PTAs are sensitive to the total mass of the SSPSs, if PTAs in the future measure differences in the masses with statistical significance, those differences may reflect differences in the masses or total number of moons taken into account when estimating the position of the planet–moons barycentre and total mass estimations. Although the GLC18 algorithm is not directly applicable to bodies in orbit around major planets, the differences in the sensitivity curves around semimajor axis values close to planetary orbits may still be associated with such errors. The mass of UMOs at these semimajor axis values may also reflect differences in the SSEs regarding the positional data of the planets and moons, since with the applied methodology such effects could potentially be absorbed by the UMO mass parameter. These results underline the potential of pulsar timing and PTA research to also provide feedback and independent checks to groups developing SSEs, and add information in the future for SSE development.

While at the precision that the IPTA DR 1 can probe the masses of SSPSs we have confirmed that the SSEs give consistent results, we found that for any given MSP the timing residuals resulting from us-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Sensitivity curves for unmodelled masses in Keplerian orbits for the IPTA DR 1, using three different SSEs. From top to bottom, the figures show the $1\sigma$ (solid, black), $2\sigma$ (dashed, blue), and $3\sigma$ (dotted, red) upper limits (corresponding to the 68, 95, and 99.7 per cent credible intervals of the posterior distribution) for the mass for the DE421, DE435, and INPOP17A SSEs, respectively. The dot-dashed lines show the expected amplitude of the residuals induced by a given mass in Keplerian orbit at any given semimajor axis value. The grey shaded region shows the position of the asteroid belt. The cyan stars indicate the official IAU masses for the planetary systems. The corresponding IAU uncertainties and uncertainties from our analysis of known SSPs on the planetary masses are plotted as yellow triangles and green circles, respectively, for comparison.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** A comparison of the $1\sigma$ upper limits for the three used SSEs presented in Fig. 5.</figcaption></figure>

<table><caption>**Table 8.** Derived upper limits for the mass of UMOs in Keplerian orbits around the SSB. The upper limits quoted correspond to the 95 per cent credible intervals of the posterior distributions.</caption><thead><tr><th rowspan="2">Semimajor<br>axis<br>(au)</th><th colspan="3">(Gm)<sup>SSE</sup><sub>IPTA1</sub>/(G<sub>M</sub>)<sup>N</sup><sub>☉</sub></th></tr><tr><th>INPOP17A</th><th>DE435</th><th>DE421</th></tr></thead><tbody><tr><td>0.5</td><td>1.14815362 × 10<sup>-10</sup></td><td>1.00000000 × 10<sup>-10</sup></td><td>1.07151931 × 10<sup>-10</sup></td></tr><tr><td>1.4</td><td>7.24435960 × 10<sup>-11</sup></td><td>5.88843655 × 10<sup>-11</sup></td><td>6.02559586 × 10<sup>-11</sup></td></tr><tr><td>5.0</td><td>2.23872114 × 10<sup>-10</sup></td><td>1.28824955 × 10<sup>-10</sup></td><td>1.20226443 × 10<sup>-10</sup></td></tr><tr><td>10</td><td>2.95120923 × 10<sup>-09</sup></td><td>2.39883292 × 10<sup>-09</sup></td><td>1.86208714 × 10<sup>-09</sup></td></tr><tr><td>17</td><td>9.12010839 × 10<sup>-08</sup></td><td>5.62341325 × 10<sup>-08</sup></td><td>3.31131121 × 10<sup>-08</sup></td></tr><tr><td>60</td><td>2.57039578 × 10<sup>-06</sup></td><td>2.81838293 × 10<sup>-06</sup></td><td>1.86208714 × 10<sup>-06</sup></td></tr></tbody></table>

ing different SSEs can vary at different levels. For individual pulsars we have noted that the differences between the RMS deviations of residuals formed using various SSEs (see Fig. 1) were up to ≈65 ns (which corresponds to relative differences of up to 22 per cent). The consistency in the SSPS masses from the various SSEs in the presence of the timing-residual differences means that, to a large extent, the levels of noise are such that the $\delta m$ uncertainties compensate for these residual differences. These results, however, motivate further research into the role of SSEs in pulsar timing models, beyond the effects on the mass constraints of Solar-system bodies.

It is worth noting that SSEs are regularly being updated with new data [DE436 is also available and the reader can see it applied to PTA data in Arzoumanian et al. (2018)] and the IAU regularly evaluates new data and updated the published CBEs. As we have seen in Section 4.2.1 for the cases of Mercury and Vesta, new data from space missions can indeed give at times dramatic improvements in the measurements of masses and other physical properties of Solar-system bodies. As such we will regularly have to check the impact of SSE updates on pulsar-timing and PTA applications and compare our results to updated CBEs. Pulsar-timing results, on the other hand, are also able to show strong improvements over time. As was noted already in CHM10, pulsar timing has the benefit of being able to improve on mass constraints with more accumulation of data, even if the data quality remains constant. The example of Saturn is important to highlight, since within the next 10 yr our data sets will be long enough to fully sample its orbit, which has a period of 29.5 yr. This will allow us to fully de-correlate signals from Saturn from those of timing and noise parameters and the uncertainties of the Saturnian $\delta m$ will be very significantly reduced. General
predictions on expected improvements in probing parameters of Solar-system bodies with future instruments can be made using generic sensitivity curves for UMOs. In GLC18, it was shown that regular observations of 20–40 MSPs with future radio telescopes, which may achieve timing residuals of average RMS levels of order ~100 ns for 20 yr, can potentially improve our mass constraints on the Jovian system by another two orders of magnitude, at levels below the current constraints by space missions. Therefore, despite the anticipated improvements in the CBE mass values of the Jovian system as a result of the analysis of data by the JUNO and JUICE space missions, our predictions are indicative of the potential for interesting results that can be produced in the future with respect to Solar-system studies using PTA data.

Finally, we comment on the relation between research into SSEs, Solar-system studies and GW searches with PTAs. Results of studies such as the present can give hints on which GW frequencies one can expect most of the differences in the limits by PTAs when using different SSEs. For example, the results shown in Fig. 6 suggests that searches for GWs using DE421 and DE435 would be mostly affected by SSE choice around GW frequencies ~2–8 nHz. The sensitivity of PTAs to the dimensionless strain of stochastic GWBs is currently of the order of $10^{-15}$ at reference frequency $1 \text{ yr}^{-1}$ (e.g. Verbiest et al. 2016). As an example, let us consider the case of a GWB formed by the superposition of GWs for a large number of GW-driven supermassive black hole binaries, which have a dimensionless strain that scales with the GWB frequency as $f_{\text{gwb}}^{-2/3}$, (e.g. Sesana 2013). The RMS of such a GWB signal would then be at levels $\lesssim 200$ ns. Depending on the pulsar position, differences between SSEs are shown to vary between ~15 and 450 ns (Caballero 2018). It is therefore rational to anticipate that once other sources of noise, such as IISM related chromatic noise, are mitigated, the GWB searches and limits will begin to depend more clearly on the choice of SSE. This is indeed the case and already GWB upper limits and detection statistics are being affected by the choice of SSEs, which leads to the need of introducing SSE-related parameters in the model in order to mitigate such effects (Arzoumanian et al. 2018). Vice versa, one expects the presence of a GWB to influence our results when trying to constrain planetary masses with PTA data. As discussed in Section 1, with more MSPs of high data precision, a GWB and a $\delta m$ signal should be distinguishable on the basis of their different angular correlation. Further work is underway to understand the impact of SSE selection to the timing and noise models and to bridge the systematics of each SSE to a given analysis.

Early attempts to do this included non-physical, generic error vectors of the position of the SSB (Deng et al. 2013; Tiburzi et al. 2016), which was implemented with real EPTA data in Lentati et al. (2015). The latter study only used one SSE, since there was no evidence of the GWB-strain upper limits being influenced by adding this SSE component in the model. In this case, the effects of the ephemeris error added in the Bayesian model and the relevant parameters were simultaneously sampled with GWB parameters. Further work in this direction was demonstrated in Taylor et al. (2017). Using simulated data, they recovered the signal induced by an error in the SSE's input Jupiter mass, correctly estimated the value of the error in the mass, and were able to distinguish the mass-error signal from a GWB signal. More recently, Arzoumanian et al. (2018) implemented a physical SSE-perturbation model that allows a combination of coordinate-frame drifts, gas-giant mass perturbations (as in this paper), and Jupiter orbital-element perturbations. Their findings indicated that upper limits and signal-versus-noise odds ratios for a GWB can vary significantly depending on the choice of SSE. The new model led to identical SSE-marginalized

GWB statistics, regardless of the initial SSE model (both JPL and IMCEE models were used). Both this work and Arzoumanian et al. (2018) the models were limited in the use of parameters describing linear mass-perturbation effects on the TOAs, but further components will be employed in the future using the upcoming new IPTA data releases.

## ACKNOWLEDGEMENTS

Part of this work is based on observations with the 100-m telescope of the Max-Planck-Institut für Radioastronomie (MPIfR) at Effelsberg. The Nançay Radio Observatory is operated by the Paris Observatory, associated with the French Centre National de la Recherche Scientifique (CNRS). Pulsar research at the Jodrell Bank Centre for Astrophysics and the observations using the Lovell Telescope are supported by a consolidated grant from the STFC in the UK. The Westerbork Synthesis Radio Telescope is operated by the Netherlands Institute for Radio Astronomy (ASTRON) with support from The Netherlands Foundation for Scientific Research NWO. The Green Bank Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Arecibo Observatory is operated by SRI International under a cooperative agreement with the NSF (AST-1100968), and in alliance with Ana G. Méndez-Universidad Metropolitana and the Universities Space Research Association. The Parkes radio telescope is part of the Australia Telescope National Facility which is funded by the Australian Government for operation as a National Facility managed by CSIRO. This work was supported by the MPG funding for the Max-Planck Partner Group. Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA). Work at NRL is supported by NASA. We acknowledge financial support from ‘Programme National de Cosmologie and Galaxies’ (PNCG) and from ‘Programme National Gravitation, Références, Astronomie, Métrologie’ (PNGRAM) funded by CNRS/INSU-IN2P3-INP, CEA, and CNES, France. We acknowledge financial support from the project ‘Opening a new era in pulsars and compact objects science with MeerKat’ in the context of INAF grant ‘SKA-CTA 2016’. Part of this research was funded by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), CE170100004. The Flatiron Institute is supported by the Simons Foundation. The computation was performed using the HERCULES cluster at the Max Planck Computing and Data Facility, Garching, the DIRAC cluster in KIAA and the TIANHE II supercomputer at Guangzhou.

RNC acknowledges the support of the International Max Planck Research School Bonn/Cologne and the Bonn-Cologne Graduate School for part of this work. KJL is supported by the National Natural Science Foundation of China (Grant No. 11373011). KJL is supported by XDB23010200, National Basic Research Program of China, 973 Program, 2015CB857101 and NSFC U15311243, 11690024. We are also supported by the MPG funding for the Max-Planck Partner Group. PL acknowledges the support of the International Max Planck Research School Bonn/Cologne. GD and KL acknowledge financial support by the European Research Council (ERC) for the ERC Synergy Grant BlackHoleCam under contract no. 610058. ZA, AB, SB-S, SC, JMC, PD, TD, EF, NG-D, PG, MTL, TJWL, ANL, DRL, RSL, DRM, MAM, STWM, CMFM, DJN, NTP, TTP, SMR, XS, JS, IS, DRS, KS, JKS, SRT, and RvH are members of the NANOGrav Physics Frontiers Center, which is supported by NSF award 1430284. GJ and GS acknowledge support from the Netherlands Organisation for Scientific Research NWO
(TOP2.614.001.602). JW is supported by Qing Cu Hui of Chinese Academy of Sciences (CAS). AS is supported by a University Research Fellowship of the Royal Society. SAS acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 694745; PI B. W. Stappers). SO was supported by the Alexander von Humboldt Foundation and acknowledges Australian Research Council grant Laureate Fellowship FL150100148. RvH acknowledges support by NASA through Einstein Fellowship grant PF3-140116. WWZ is supported by the Chinese Academy of Science Pioneer Hundred Talents Program and the Strategic Priority Research Program of the Chinese Academy of Sciences Grant No. XDB230000000. We thank Alessandro Ridolfi for a careful review of this paper and useful comments, as well as Agnès Fienga and Mickael Gastineau for useful pointers regrding the use of CALCEPH and the anonymous referee for helpful comments and suggestions.

## REFERENCES

Anderson J. D., Colombo G., Espsitio P. B., Lau E. L., Trager G. B., 1987, *Icarus*, 71, 337

Arzoumanian Z. et al., 2015, *ApJ*, 813, 65

Arzoumanian Z. et al., 2018, *ApJ*, 859, 47

Blanco-Pillado J. J., Olum K. D., Shlaer B., 2014, *Phys. Rev. D*, 89, 023512

Blandford R., Narayan R., Romani R. W., 1984, *J. Astrophys. Astron.*, 5, 369

Brown M. E., Batygin K., 2016, *ApJ*, 824, L23

Caballero R. N., 2018, in Weltevrede P., Perera B. B. P., Preston L. L., Sanidas S., eds, Proc. IAU Symp. 337, Pulsar Astrophysics the Next Fifty Years, Cambridge University Press, p. 154

Caballero R. N. et al., 2016, *MNRAS*, 457, 4421

Carry B., 2012, *Planet. Space Sci.*, 73, 98

Champion D. J. et al., 2010, *ApJ*, 720, L201 (CHM10)

Coles W., Hobbs G., Champion D. J., Manchester R. N., Verbiest J. P. W., 2011, *MNRAS*, 418, 561

Cordes J. M., 2013, *Class. Quantum Gravity*, 30, 224002

Cordes J. M., Shannon R. M., Stinebring D. R., 2016, *ApJ*, 817, 16

Damour T., 2009, in Colpi M., Casella P., Gorini V., Moschella U., Possenti A., eds, Astrophysics and Space Science Library, Vol. 359, Physics of Relativistic Objects in Compact Binaries: From Birth to Coalescence, Springer, Netherlands, p. 1

Demorest P. B. et al., 2013, *ApJ*, 762, 94

Deng X. P. et al., 2013, *Adv. Space Res.*, 52, 1602

Desvignes G. et al., 2016, *MNRAS*, 458, 3341

Edwards R. T., Hobbs G. B., Manchester R. N., 2006, *MNRAS*, 372, 1549

Feroz F., Hobson M. P., Bridges M., 2009, *MNRAS*, 398, 1601

Fienga A., Manche H., Laskar J., Gastineau M., 2008, *A&A*, 477, 315

Fienga A. et al., 2009, *A&A*, 507, 1675

Fienga A., Manche H., Laskar J., Gastineau M., Verma A., 2013, preprint (arXiv:1301.1510)

Fienga A., Manche H., Laskar J., Gastineau M., Verma A., 2014, preprint (arXiv:1405.0484v2)

Folkner W. M., Standish E. M., Williams J. G., Boggs D. H., 2007, IOM 343R-07-005

Folkner W. M., Williams J. G., Boggs D. H., 2009, Interplanetary Network Progress Report, 178, 1

Folkner W. M., Williams J. G., Boggs D. H., Park R. S., Kuchynka P., 2014, Interplanetary Network Progress Report, 196, 1

Folkner W. M., Park R. S., Jacobson R. A., 2016, IOM 392R-16-003.

Foster R. S., Backer D. C., 1990, *ApJ*, 361, 300

Gastineau M., Laskar J., Manche H., Fienga A., 2015, Astrophysics Source Code Library, record ascl:1505.001

Guo Y. J., Lee K. J., Caballero R. N., 2018, *MNRAS*, 475, 3644 (GLC18)

Hellings R. W., Downs G. S., 1983, *ApJ*, 265, L39

Hobbs G. B., Edwards R. T., Manchester R. N., 2006, *MNRAS*, 369, 655

Hobbs G. et al., 2012, *MNRAS*, 427, 2780

Kaspi V. M., Taylor J. H., Ryba M. F., 1994, *ApJ*, 428, 713

Keith M. J. et al., 2013, *MNRAS*, 429, 2161

Konopliv A. S. et al., 2014, *Icarus*, 240, 103

Konopliv A. S. et al., 2018, *Icarus*, 299, 411

Kramer M., Lange C., Lorimer D. R., Backer D. C., Xilouris K. M., Jessner A., Wielebinski R., 1999, *ApJ*, 526, 957

Kramer M., Lyne A. G., O'Brien J. T., Jordan C. A., Lorimer D. R., 2006, Science, 312, 549

Landau L. D., Lifshitz E. M., 1960, Electrodynamics of Continuous Media. Pergamon Press, Oxford

Lee K. J. et al., 2014, *MNRAS*, 441, 2831

Lentati L., Alexander P., Hobson M. P., Taylor S., Gair J., Balan S. T., van Haasteren R., 2013, *Phys. Rev. D*, 87, 104021

Lentati L., Alexander P., Hobson M. P., Feroz F., van Haasteren R., Lee K. J., Shannon R. M., 2014, *MNRAS*, 437, 3004

Lentati L. et al., 2015, *MNRAS*, 453, 2576

Lentati L. et al., 2016, *MNRAS*, 458, 2161

Liu K., Keane E. F., Lee K. J., Kramer M., Cordes J. M., Purver M. B., 2012, *MNRAS*, 420, 361

Loeb A., Zaldarriaga M., 2005, *Phys. Rev. D*, 71, 103520

Lorimer D. R., Kramer M., 2005, Handbook of Pulsar Astronomy. Cambridge Univ. Press, Cambridge

Madison D. R., Chatterjee S., Cordes J. M., 2013, *ApJ*, 777, 104

Mamajek E. E. et al., 2015, preprint (arXiv:1510.07674)

Manchester R. N. et al., 2013, *PASA*, 30, e017

Mazarico E., Genova A., Goossens S., Lemoine F. G., Neumann G. A., Zuber M. T., Smith D. E., Solomon S. C., 2014, *J. Geophys. Res. (Planets)*, 119, 2417

Petit G., Luzum B., 2010, IERS Technical Note. p. 36

Pitjev N. P., Pitjeva E. V., 2013, *Astron. Lett.*, 39, 141

Pitjeva E. V., 2013, *Solar Syst. Res.*, 47, 386

Pitjeva E. V., Pitjev N. P., 2013, *MNRAS*, 432, 3431

Reardon D. J. et al., 2016, *MNRAS*, 455, 1751

Russell C. T., Raymond C. A., 2011, *Space Sci. Rev.*, 163, 3

Sesana A., 2013, *MNRAS*, 433, L1

Shannon R. M., Cordes J. M., 2010, *ApJ*, 725, 1607

Shannon R. M. et al., 2014, *MNRAS*, 443, 1463

Siemens X., Ellis J., Jenet F., Romano J. D., 2013, *Class. Quantum Gravity*, 30, 224015

Skilling J., 2004, in Fischer R., Preuss R., von Toussaint U., eds, AIP Conf. Proc. Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering. Am. Inst. Phys., New York, p. 395

Standish E. M., 1998, IOM 312F-98-048

Taylor J. H., 1992, *Phil. Trans. R. Soc. A*, 341, 117

Taylor S. R., Vallisneri M., Ellis J. A., Mingarelli C. M. F., Lazio T. J. W., van Haasteren R., 2016, *ApJ*, 819, L6

Taylor S. R., Lentati L., Babak S., Brem P., Gair J. R., Sesana A., Vecchio A., 2017, *Phys. Rev. D*, 95, 042002

Tiburzi C., Verbiest J., 2018, in Weltevrede P., Perera B. B. P., Preston L. L., Sanidas S., eds, Proc. IAU Symp. 337, Pulsar Astrophysics the Next Fifty Years, Cambridge University Press, p. 279

Tiburzi C. et al., 2016, *MNRAS*, 455, 4339

van Haasteren R., Levin Y., 2013, *MNRAS*, 428, 1147

van Haasteren R., Levin Y., McDonald P., Lu T., 2009, *MNRAS*, 395, 1005

Verbiest J. P. W., Shaifullah G. M., 2018, *Class. Quantum Gravity*, 35, 133001

Verbiest J. P. W. et al., 2016, *MNRAS*, 458, 1267

Viswanathan V., Fienga A., Gastineau M., Laskar J., 2017, Notes scientifiques et techniques de l'Institut de Méchanique et Céleste, ISSN: 1621-3823, ISBN: 2-910015-79-3

Wang J. B. et al., 2017, *MNRAS*, 469, 425

Xilouris K. M., Kramer M., Jessner A., Wielebinski R., Timofeev M., 1996, *A&A*, 309, 481

You X. P., Hobbs G. B., Coles W. A., Manchester R. N., Han J. L., 2007, *ApJ*, 671, 907

Zhu W. W. et al., 2015, *ApJ*, 809, 41
<footnote n="1">¹Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="2">²Kavli institute for Astronomy and Astrophysics, Peking University, Beijing 100871, P. R. China</footnote>

<footnote n="3">³Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="4">⁴Cahill Center for Astrophysics, California Institute of Technology, MC 249-17, 1200 E. California Boulevard, Pasadena, CA 91125, USA</footnote>

<footnote n="5">⁵X-Ray Astrophysics Laboratory, NASA Goddard Space Flight Center, Code 662, Greenbelt, MD 20771, USA</footnote>

<footnote n="6">⁶Centre for Astrophysics and Supercomputing, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="7">⁷ASTRON, the Netherlands Institute for Radio Astronomy, Postbus 2, NL-7990 AA, Dwingeloo, the Netherlands</footnote>

<footnote n="8">⁸International Centre for Radio Astronomy Research, Curtin University, Bentley, WA 6102, Australia</footnote>

<footnote n="9">⁹Cornell Center for Advanced Computing, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="10">¹⁰Cornell Center for Astrophysics and Planetary Science, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="11">¹¹INAF – Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047 Selargius (CA), Italy</footnote>

<footnote n="12">¹²Department of Physics and Astronomy, West Virginia University, Morgantown, WV 26506, USA</footnote>

<footnote n="13">¹³Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA</footnote>

<footnote n="14">¹⁴Institute for Gravitation and the Cosmos, Department of Physics, The Pennsylvania State University, University Park, PA 16802, USA</footnote>

<footnote n="15">¹⁵Astronomy Department, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="16">¹⁶Laboratoire de Physique et Chimie de l'Environnement et de l'Espace LPC2E UMR7328, Université d'Orléans, CNRS, F-45071 Orléans, France</footnote>

<footnote n="17">¹⁷Station de Radioastronomie de Nançay, Observatoire de Paris, PSL University, CNRS, Université d'Orléans, F-18330 Nançay, France</footnote>

<footnote n="18">¹⁸CSIRO Astronomy and Space Science, Australia Telescope National Facility, Box 76, Epping, NSW 1710, Australia</footnote>

<footnote n="19">¹⁹National Radio Astronomy Observatory, PO Box O, Socorro, NM 87801, USA</footnote>

<footnote n="20">²⁰Department of Physics, Hillsdale College, 33 E. College Street, Hillsdale, MI 49242, USA</footnote>

<footnote n="21">²¹Faculty of Science, University of East Anglia, Norwich Research Park, Norwich NR4 7TJ, UK</footnote>

<footnote n="22">²²Department of Physics, McGill University, 3600 University St., Montreal, QC H3A 2T8, Canada</footnote>

<footnote n="23">²³School of Mathematics, James Clerk Maxwell Building, Peter Guthrie Tait Road, University of Edinburgh, Edinburgh EH9 3FD, UK</footnote>

<footnote n="24">²⁴Vancouver Coastal Health, Department of Nuclear Medicine, 899 W 12th Ave, Vancouver, BC V5Z 1M9, Canada</footnote>

<footnote n="25">²⁵Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada</footnote>

<footnote n="26">²⁶Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, NL-6500 GL Nijmegen, the Netherlands</footnote>

<footnote n="27">²⁷Space Science Division, Naval Research Laboratory, Washington, DC 20375-5352, USA</footnote>

<footnote n="28">²⁸Monash Centre for Astrophysics (MoCA), School of Physics and Astronomy, Monash University, VIC 3800, Australia</footnote>

<footnote n="29">²⁹Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr, M/S 67-201, Pasadena, CA 91109, USA</footnote>

<footnote n="30">³⁰Department of Physics and Astronomy, Haverford College, 370 Lancaster Ave, Haverford, PA 19041, USA</footnote>

<footnote n="31">³¹National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA</footnote>

<footnote n="32">³²Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA</footnote>

<footnote n="33">³³Physics Department, Lafayette College, Easton, PA 18042, USA</footnote>

<footnote n="34">³⁴Fakultät für Physik, Universität Bielefeld, Postfach 100131, D-33501 Bielefeld, Germany</footnote>

<footnote n="35">³⁵Physics Department, Texas Tech University, Box 41051, Lubbock, TX 79409, USA</footnote>

<footnote n="36">³⁶Hungarian Academy of Sciences MTA-ELTE 'Extragalactic Astrophysics' Research Group, Institute of Physics, Eötvös Loránd University, Pázmány P. s. I/A, Budapest 1117, Hungary</footnote>

<footnote n="37">³⁷Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, NL-1098 XH Amsterdam, the Netherlands</footnote>

<footnote n="38">³⁸School of Physics and Astronomy, The University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="39">³⁹ARC Centre of Excellence for Gravitational Wave Discovery (OzGrav), Swinburne University of Technology, John St, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="40">⁴⁰Center for Gravitation, Cosmology and Astrophysics, Department of Physics, University of Wisconsin-Milwaukee, PO Box 413, Milwaukee, WI 53201, USA</footnote>

<footnote n="41">⁴¹Physics and Astronomy Department, Oberlin College, Oberlin, OH 44074, USA</footnote>

<footnote n="42">⁴²TAPIR Group, MC 350-17, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="43">⁴³Laboratoire Univers et Théories (LUTh), Observatoire de Paris, PSL University, CNRS, Université Paris-Diderot, F-92190 Meudon, France</footnote>

<footnote n="44">⁴⁴Institute for Radio Astronomy & Space Research, Auckland University of Technology, Private Bag 92006, Auckland 1142, New Zealand</footnote>

<footnote n="45">⁴⁵Xinjiang Astronomical Observatory, Chinese Academy of Sciences, 150 Science 1-Street, Urumqi, Xinjiang 830011, China</footnote>

<footnote n="46">⁴⁶National Astronomical Observatories, Chinese Academy of Science, 20A Datun Road, Chaoyang District, Beijing 100012, China</footnote>

<footnote>This paper has been typeset from a TEX/LATEX file prepared by the author.</footnote>