The NANOGrav 11 yr Data Set: Constraints on Planetary Masses Around 45 Millisecond Pulsars

In 1992, the first two confirmed exoplanets were found around the isolated pulsar B1257+12 using pulsar timing techniques (Wolszczan & Frail 1992). This 6.2 ms pulsar was initially found to host two planets, each being a few earth masses with periods of 2–3 months. Two years after this initial discovery, another periodic signal was discovered and later confirmed to be a third exoplanet around the same pulsar, and to this day, this planet is still the least-massive confirmed exoplanet (Wolszczan 1994), with a mass of only about two moon masses, or about 10²³ kg. In 2000, a circumbinary planet was discovered around PSR B1620−26 and its white dwarf companion, and this planet was also announced to be the oldest planet discovered with an age of about 12.6 Gyr (Ford et al. 2000).

Methods of pulsar-planet detection are described in Kramer (2018), and early work in the field of pulsar planets was reviewed in Phillips & Thorsett (1994). There have been few recent systematic searches for pulsar planets, particularly around millisecond pulsars. Kerr et al. (2015) searched for planets around a large number of pulsars, but they were non-millisecond pulsars. Caballero et al. (2018) used millisecond pulsar data to search for unknown planets in our own solar system, but their method, cross-correlating data across pulsars, precluded detection of planets around individual pulsars.

In total, six planetary-mass bodies have been confirmed to be orbiting pulsars (Bailes et al. 2011; Spiewak et al. 2018), but the mechanism of forming pulsar–planet systems is unclear. If these planets formed in the protoplanetary disk phase, we would not expect them to survive the subsequent supernova explosion, making the discovery of pulsar planets a somewhat unexpected result. One theory proposes that when the supernovae that produce these pulsars explode, some of the matter that is ejected is then gravitationally recaptured by the pulsar and forms a rotating disk, functioning similarly to the protoplanetary disks with which we are familiar (Lin et al. 1991).

Another theory that has been proposed for pulsars in binary orbits involves the planet matter actually coming from the pulsar's companion. Many millisecond pulsars (MSPs) that are thought to have been spun up to their current rotational period by accreting matter from a low-mass donor star and receiving an influx of angular momentum that eventually vaporizes their companions, earning the nickname "black widow" pulsars (Bhattacharya & van den Heuvel 1991; Arzoumanian et al. 1999; Roberts 2011). It has been suggested that the high-energy relativistic wind from these pulsars could ablate their low-mass companions, but the matter being chipped away is then recaptured by the pulsar and forms a rotating disk (Perets 2010). Like the previous model, planetesimals would then form from this disk similarly to how they form in protoplanetary disks around newly formed stars. It is even possible that mechanisms in the theories described above could lead to the formation of asteroid belts around pulsars (Shannon et al. 2013).

In this Letter, we discuss and analyze our search for planets around MSPs in the 11 yr data set (Arzoumanian et al. 2018) from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav). NANOGrav is a collaboration of scientists from the United States and Canada whose goal is to use pulsars as a galactic-scale detector for long-period gravitational waves (Ransom et al. 2019). This data set contains approximately 11 yr of time of arrival (TOA) data for 45 MSPs and their timing ephemerides. In the pursuit of gravitational wave detection, NANOGrav has compiled a data set of some of the most highly stable MSPs in the sky with timing precision of under 1–2 μs, which makes this data set an excellent resource for not only gravitational wave studies, but for other astrophysics as well.

This Letter describes our use of the NANOGrav 11 yr data set to search for exoplanets around pulsars as well as to set lower limits on the masses of exoplanets that could be detected using this data set. In Section 2, we expand upon the qualities of the NANOGrav data set. In Section 3 we discuss how we expect to see planetary signals in pulsar timing residuals (the difference between observed and predicted TOAs) and how we optimize the data set for exoplanet searches. In Section 4 we introduce our detection scheme as well as our methods for determining mass constraints. In Section 5 we present the results of our search and our efforts to define detection limits, and finally in Section 6 we summarize our work and propose some future work in this field.

The NANOGrav 11 yr data set contains high-precision TOA measurements as well as timing models for 45 MSPs taken over the span of roughly 11 yr, though observations of some pulsars span longer or shorter periods of time. Observations of these pulsars were taken with a roughly monthly, though sometimes weekly, cadence on the 305 m William E. Gordon Telescope at the Arecibo Observatory as well as on the 100 m Robert C. Byrd Green Bank Telescope at the Green Bank Observatory. Of these 45 pulsars, 31 are in binary systems, the orbits of which we have modeled with both Keplerian and post-Keplerian parameters to account for apparent orbital deviations from Keplarian motion (Damour & Taylor 1992).

In this Letter, we primarily consider the pulsar timing residuals of the NANOGrav data to analyze the variation between the data and our models. These models were created using a weighted, linear least-squares fit generated using TEMPO (Nice et al. 2015), a software package commonly used in pulsar timing. For most of the pulsars included in this data set, the rms of the timing residuals is below 2 μs; the only exceptions include those pulsars for which the collected data exhibit strong red noise, or low-frequency noise (Arzoumanian et al. 2018).

To determine if exoplanets are orbiting the MSPs in this data set, we analyze the timing residuals to search for a characteristic variation that we would expect to see from our models as a result of a planetary presence. Due to the reflex motion of a pulsar in response to a planetary orbit, we expect to see fluctuations in pulsar TOAs as the pulsar orbits the center of mass (COM) of the pulsar–planet system. The pattern of these fluctuations would repeat periodically and, assuming a circular orbit, would therefore be seen as a sinusoidal structure in the pulsar's residuals. By extracting the amplitude and frequency of these sinusoids, we can gather information about the pulsar's distance from the COM and the frequency of its orbits, which is also equal to the planet's orbital frequency. We can input these parameters into the equations for Kepler's third law

$$\begin{eqnarray}&&{P}^{2}=\displaystyle \frac{4{\pi }^{2}{\left(R+r\right)}^{3}}{G\left(M+m\right)}\end{eqnarray} \tag{ 1 }$$

and the definition of the COM

$$\begin{eqnarray}&&{MR}={mr},\end{eqnarray} \tag{ 2 }$$

where P is the planet's orbital period; R and r are the distances to the system's COM from the pulsar (and its companion, if present) and the planet, respectively; M is the central mass of the pulsar system, either the pulsar's mass M₁ by itself in the case of an isolated pulsar or the combined mass of the pulsar and the mass of its companion M₁ + M₂; and m is the mass of the planet. We use these equations to ultimately calculate the mass of the planet associated with the aforementioned sinusoidal signal and find the relationship m ∝ P^−2/3.

To determine some characteristics of the planetary signals that we might see in our data, we carefully consider the orbital frequencies that would be detectable and would make the most physical sense given the parameters of our systems. As a result of the, at most, weekly sampling of the NANOGrav data, we determine that we would be unlikely to detect any planets with orbital periods of less than one week, so we define 7 days as our generic minimum orbital period to be probed. For the pulsars in binary orbits, additional analysis is done to account for the various possible planetary configurations given the characteristics of our binaries as described below.

We consider two possible planetary configurations for the binary systems in the NANOGrav 11 yr release, both of which are described in detail in Chambers et al. (2002). The first configuration that we consider is a planet in a satellite S-type orbit, in which the planet orbits only the pulsar, and the companion orbits both the planet and the pulsar around the outskirts of the system. Holman & Wiegert (1999) defined a critical semimajor axis a_c for planets in this configuration in which planets with semimajor axis a ≤ a_c will be stable against perturbations caused by the binary companion. For S-type orbits, the semimajor axes follow the relationship

$$\begin{eqnarray}&&\begin{array}{l}a\leqslant [(0.464\pm 0.006)+(-0.380\pm 0.010)\mu \\ \ +\left(-0.631\pm 0.034\right)e+\left(0.586\pm 0.061\right)\mu e\\ \ +\left(0.150\pm 0.041\right){e}^{2}+\left(-0.198\pm 0.074\right)\mu {e}^{2}]\ {a}_{b},\end{array}\end{eqnarray} \tag{ 3 }$$

where e is the eccentricity of the binary orbit; μ is the mass ratio M₂/(M₁ + M₂); and a_b is the semimajor axis of the pulsar's orbit. For the majority of the binary systems in this data set, e ≈ 0, M₁ ≈ 1.5 M_⊙, and M₂ ≈ 0.2 M_⊙ (μ ≈ 1/9). Given these parameters, Equation (3) reduces to

$$\begin{eqnarray}&&a\lessapprox 0.4{a}_{b}.\end{eqnarray} \tag{ 4 }$$

PSR J1903+0327 is the one notable exception to this approximation due to its eccentric orbit (e ≈ 0.44) with a main-sequence star (Fonseca et al. 2016).

However, because of the relatively compact orbits of the binaries in the NANOGrav 11 yr data set, the periods (listed in Table 1) corresponding to critical semimajor axes of planets in this configuration around a NANOGrav pulsar would be extremely short and therefore below the one-week limit set by NANOGrav's sampling frequency. We determine that a planet in an S-type configuration would not be detectable in the NANOGrav data set.

Table 1.  Orbital Characteristics of Planets around Binary MSPs and Mass Limits for Planets with P = 100 Days

Pulsar	Ppl,max, S-type	${P}_{\mathrm{pl},\min }$, P-type	${M}_{\mathrm{pl}}$, P = 100 days
	(days)	(days)	(Mmoon)
J0023+0923	2.58 × 10−5	4.70 × 10−4	1.00
J0030+0451	⋯	⋯	0.47
J0340+4130	⋯	⋯	1.90
J0613−0200	4.53 × 10−3	0.08	0.28
J0636+5128	⋯	⋯	1.14
J0645+5158	⋯	⋯	0.91
J0740+6620	⋯	⋯	1.45
J0931−1902	⋯	⋯	2.68
J1012+5307	1.76 × 10−3	0.03	0.87
J1024−0719	⋯	⋯	1.08
J1125+7819	⋯	⋯	3.64
J1453+1902	⋯	⋯	8.39
J1455−3330	0.73	13.30	2.19
J1600−3053	0.10	1.89	0.22
J1614−2230	0.15	2.74	0.38
J1640+2224	1.63	29.77	0.67
J1643−1224	0.50	9.08	0.26
J1713+0747	0.73	13.29	0.16
J1738+0333	7.99 × 10−4	0.01	1.09
J1741+1351	0.14	2.64	0.44
J1744−1134	⋯	⋯	0.56
J1747−4036	⋯	⋯	1.44
J1832−0836	⋯	⋯	1.00
J1853+1303	1.03	18.80	0.97
B1855+09	0.11	2.03	0.21
J1903+0327	1.51	137.84	0.32
J1909−3744	0.01	0.19	0.13
J1910+1256	0.39	7.02	0.53
J1911+1347	⋯	⋯	1.10
J1918−0642	0.10	1.74	0.72
J1923+2515	⋯	⋯	1.33
B1937+21	⋯	⋯	0.01
J1944+0907	⋯	⋯	0.93
B1953+29	0.70	12.72	0.88
J2010−1323	⋯	⋯	0.85
J2017+0603	0.01	0.23	0.54
J2033+1734	⋯	⋯	1.45
J2043+1711	0.01	0.15	0.27
J2145−0750	0.13	2.34	1.53
J2214+3000	5.70 × 10−5	1.04 × 10−3	1.98
J2229+2643	⋯	⋯	1.84
J2234+0611	⋯	⋯	1.30
J2234+0944	⋯	⋯	1.51
J2302+4442	1.46	26.67	1.90
J2317+1439	0.01	0.25	0.45

Note. ${P}_{\mathrm{pl},\max }$ is the maximum orbital period possible for planets in an S-type configuration. ${P}_{\mathrm{pl},\min }$ is the minimum orbital period of planets in P-type configurations. ${M}_{\mathrm{pl}}$ is the minimum planetary mass detectable in a circular orbit using the 11 yr data set given a 90% confidence level and adjusted by the factor of 2.5 discussed in Section 5.

Download table as:  ASCIITypeset image

We then consider bodies in planetary P-type orbits, where the planet orbits both the pulsar and its companion as if they are a single mass. Holman & Wiegert (1999) defined a lower limit for the semimajor axes of planets in this configuration, and thus possible values follow the relationship

$$\begin{eqnarray}\begin{array}{rcl}a & \geqslant & \left[(1.60\pm 0.04)+(5.10\pm 0.05)e\right.\\ & & \quad +(-2.22\pm 0.11){e}^{2}+(4.12\pm 0.09)\mu \\ & & \quad +(-4.27\pm 0.17)e\mu +(-5.09\pm 0.11){\mu }^{2}\\ & & \quad +\left.(4.61\pm 0.36){e}^{2}{\mu }^{2}\right]\ {a}_{b}.\end{array}\end{eqnarray} \tag{ 5 }$$

As with Equation (3), using our assumptions that e ≈ 0 and μ ≈ 1/9, we can reduce Equation (5) to

$$\begin{eqnarray}&&a\gtrapprox 2{a}_{b}.\end{eqnarray} \tag{ 6 }$$

Due to the close orbits of the NANOGrav binaries, the majority of the minimum planetary periods in P-type configurations (also listed in Table 1) are smaller than the 7 day limit defined by our observation frequency. For these systems, we adopt 7 days as the lowest orbital period to probe; for systems with minimum periods that are greater than 7 days, we use the constraints set by their binary orbits to define the minimum orbital period to investigate.

We also investigate period constraints set by possible interactions with the timing models at lower frequencies. Certain patterns in pulsar timing data recur periodically as a result of Earth's rotation around the Sun. In order to remove these systematics to better model the behavior of the pulsar, the timing models fit for the pulsar's position and proper motion, which ultimately erases power from any periodic signals with a frequency of 1 yr⁻¹, and for parallax, which will remove power at frequencies of 2 yr⁻¹. Though fitting for parallax does not result in a significant chance that signals with frequencies of 2 yr⁻¹ will go undetected, we find that it would be nearly impossible to detect signals of about 1 yr⁻¹ due to TEMPO's fitting procedure for position and proper motion. We also define a lower limit for potential orbital frequencies as 1/T, where T is the total span of the data.

The NANOGrav data set timing models use a combination of white and red noise in the model of the timing of each pulsar.

In the NANOGrav data set timing models, the white noise is modeled by three factors: EFAC, a multiplicative factor applied to measured TOA uncertainties; EQUAD, a factor added in quadrature to measured TOA uncertainties; and ECORR, a factor that accounts for noise correlated across frequencies in data collected across a wide band. Such white noise models could hide the presence of signals from pulsar planets. Therefore, for this Letter, we excluded the white noise model (setting EFAC to 1 and EQUAD and ECORR to 0). Theoretically, this reduction in TOA uncertainties could lower our detection threshold and produce spurious detections of planet signals. In practice, as described below, we had no detections of planets, so we believe the reduction in TOA uncertainties does not contaminate our results.

NANOGrav also uses parameters to account for red noise. However, we only expect red noise to be significant on timescales as long as or longer than the data set, and we only search for planets with periods up to 2000 days, so the strong red noise exhibited by some of the pulsars in this data set is unlikely to affect our analysis. We therefore adopt the red noise parameters found by Arzoumanian et al. (2018) and use "unwhitened" residuals, or residuals off of which the red-noise model has not been subtracted. We then rerun TEMPO on all the pulsars in the 11 yr data set to generate new timing models given these noise parameters.

As discussed above, we would expect to see sinusoidal variation in a pulsar's timing residuals if a planet is present in the system. To detect these sinusoids, we use Astropy's LombScargle²⁶ module (named after Lomb 1976 and Scargle 1982) to perform a periodogram analysis for detecting sinusoids in our data that may not be visible by eye (VanderPlas et al. 2012; Astropy Collaboration et al. 2013; VanderPlas & Ivezić 2015). The Lomb–Scargle periodogram is a tool specifically designed to be used with unevenly sampled data with unequal errors bars to fit sinusoids of various frequencies to a data set. The quality of a sinusoid's fit to the data is determined by how dramatically the reduced-χ² of the data improves if a fitted sinusoidal signal were removed. The Astropy module produces dimensionless, normalized power spectra that reflect these determinations with high spikes in power at the frequencies associated with the possible sinusoids identified in the data.

To use this module to detect potential planetary signals, we define a detection threshold based on the maximum power that could be achieved by residuals without a planetary signal. We simulate 10,000 realizations of white-noise residuals for each pulsar with each data point being taken from a Gaussian distribution centered at 0 that has a standard deviation equal to the error on that TOA. We also set the error on this simulated TOA equal to the error on the observed TOA. We set our detection threshold as the highest power value generated from the periodograms run on these simulated residuals. We then run the real residual data for each pulsar through the Astropy LombScargle module and define a potential planetary detection as the presence of a power spike greater than the threshold value defined above.

We also use the LombScargle module to set lower limits on the planetary masses that can be detected using pulsar timing by injecting planetary signals into simulated residuals. Because we cannot be sure no undetected planetary signal exists in the NANOGrav 11 yr data, we again simulate residuals that mirror the noise levels of each pulsar but are guaranteed to contain only white noise. We define a frequency grid of 1000 linearly spaced frequencies with the maximum and minimum defined as in Section 3 above. We then inject planetary signals with these frequencies and known amplitudes into the data and attempt to detect them using Lomb–Scargle periodograms and the detection technique described above.

At each chosen frequency, we lower the amplitude of the injected sinusoid until it can no longer be detected more than 90% of the time over 1000 simulations of noisy residuals. Figure 1 shows an example of simulated residuals injected with a planetary signal with an amplitude at the lower limit of our detectable range. If a periodogram spike is detected with a value greater than our noise-only threshold, we also identify whether the periodogram has recovered the correct period of the injected signal by determining whether the number of planetary cycles of the recovered signal over the span of the data set matches that of the injected signal. After using the amplitude and frequency data to determine the lowest detectable mass for each orbital period, we perform a linear, least-squares fit in logarithmic space to obtain a best-fit line for each pulsar's mass-period relationship (see Figure 3).

Zoom In Zoom Out Reset image size

To ensure that our method for determining lower limits is realistic, we also run TEMPO on pulsar residuals that include an injected planetary signal, to see if covariances with other timing parameters will affect planet detectability. Due to the extremely time-consuming nature of running thousands of TEMPO iterations, we choose to test the validity of our detection scheme with a single simulated pulsar under the assumption that the results will be reflective of all the pulsars in this data set. We simulate pulsar TOAs spanning about 4600 days that reflect NANOGrav's uneven sampling rate and have uniform 1 μs errors, making it a fairly typical NANOGrav MSP, and we add a sinusoid to the TOAs that is representative of a planetary signal. We then run TEMPO again to obtain a model for the simulated pulsar with its injected signal and proceed to attempt to detect the signal as described above.

No exoplanets are discovered around the pulsars in this data set. Figure 2 shows two examples of power spectra from the Lomb–Scargle periodograms. As shown, no significant peaks stand out in the spectra, and all power values are well below the thresholds determined for each pulsar. We therefore conclude that no planets orbiting with periods within our range of sensitivity exist around the pulsars in the NANOGrav data set.

Zoom In Zoom Out Reset image size

The results of the simulated planetary signal injections reinforce our conclusion regarding the presence of exoplanets around these pulsars. The injections demonstrate the incredible sensitivity of NANOGrav pulsar timing data to planetary perturbations, which is shown in Figure 3, with the mass constraints for planets with periods of 100 days listed in Table 1. Within the range of orbital periods investigated here, all of the NANOGrav pulsars demonstrate a sensitivity to planetary masses well below an earth mass and even as low as a fraction of a moon mass. These limits also demonstrate a −2/3 slope as a function of orbital period, which results from Equation (1) as well as NANOGrav's roughly equal sensitivity to planetary perturbations across frequencies.

Zoom In Zoom Out Reset image size

The detections made after running TEMPO on simulated residuals with an injected planetary signal are also shown in Figure 3. There is a significant loss of sensitivity at 1 yr due to fitting for pulsar position and proper motion. Our ability to detect planetary signals wanes further at periods of about 2000 days or greater when the injected signal's period becomes equal to or longer than the range of the data set.

We also find that running TEMPO on data that already contains a planetary signal, yet with no planets in the timing model, will result in a small loss of sensitivity. TEMPO's model fitting procedure with the incorrect timing model absorbs some of that planetary signal, causing us to become less sensitive to planets. We used synthetic TOAs from an isolated pulsar with a constant dispersion measure to calibrate the mass-detection penalty and find it to be approximately a factor of 2.5 over the planetary orbital periods of interest (i.e., when searching TEMPO-derived residuals for an injected signal in the pulse arrival times, the planets we are able to detect are 2.5 times more massive than those whose signals are injected directly into the residuals themselves). As a result, we increase the lower limits found using our original method by this factor of 2.5 in order to more accurately reflect the lowest planetary masses that would be detected after using TEMPO to fit and model the NANOGrav pulsars. Additionally, because no red noise parameters are used to model the noise of the simulated pulsar, the post-TEMPO simulated data limit appears slightly more sensitive than the constraints for the real pulsars, whose models account for red noise, after they are adjusted by a factor of 2.5.

Figure 3 also indicates how the sensitivity of the NANOGrav pulsars compares with that of other exoplanet detection methods. Using data from the NASA Exoplanet Archive, the least-massive 10% of exoplanets found using the listed methods are plotted over the mass limits we find using timing of the NANOGrav pulsars. In general, the NANOGrav pulsars demonstrate that pulsar timing is far more sensitive to planetary detections than other methods, sometimes by orders of magnitude.

As a result of the exquisite sensitivity stemming from the high timing precision and long baselines of the NANOGrav data set, we determine that pulsar timing is capable of detecting planetary masses orders of magnitude smaller than any other detection method. Our results are primarily limited by the sampling frequency of the NANOGrav data. In Section 3, our investigations show that there is a physical possibility of planets existing with orbital periods shorter than those that can be detected given this data set's sampling frequency. More closely spaced timing observations, taken several times a day over the span of weeks or months, would allow us to probe higher-frequency orbits to which we are not currently sensitive. It is worth noting, however, that if such short-period planets existed in our data with fairly massive planetary companions, the white-noise parameters from our standard NANOGrav fits would indicate substantial excess noise in our data (i.e., via large EFACs, for instance), which is something that we do not presently see.

We are of course also limited by the frequency with which planets exist around pulsars, which is not currently known. As only a handful of planetary-mass bodies have to date been discovered around three or four pulsars, it seems likely that their existence around millisecond pulsars is quite rare. Further investigations are needed to understand the demographics of this population of planets.

The NANOGrav project receives support from National Science Foundation (NSF) Physics Frontiers Center award number 1430284. The National Radio Astronomy Observatory and the Green Bank Observatory are facilities of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The Arecibo Observatory is a facility of the National Science Foundation operated under cooperative agreement (#AST-1744119) by the University of Central Florida in alliance with Universidad Ana G. Méndez (UAGM) and Yang Enterprises (YEI), Inc. Portions of this work performed at NRL are supported by the Chief of Naval Research. S.M.R. is a CIFAR Fellow. The NASA Exoplanet Archive is operated by the California Institute of Technology under contract with the National Aeronautics and Space Administration.