The following article is Open access

Multimessenger Gravitational-wave Searches with Pulsar Timing Arrays: Application to 3C 66B Using the NANOGrav 11-year Data Set

, , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , and

Published 2020 September 7 © 2020. The Author(s). Published by the American Astronomical Society.
, , Citation Zaven Arzoumanian et al 2020 ApJ 900 102 DOI 10.3847/1538-4357/ababa1

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

0004-637X/900/2/102

Abstract

When galaxies merge, the supermassive black holes in their centers may form binaries and emit low-frequency gravitational radiation in the process. In this paper, we consider the galaxy 3C 66B, which was used as the target of the first multimessenger search for gravitational waves. Due to the observed periodicities present in the photometric and astrometric data of the source, it has been theorized to contain a supermassive black hole binary. Its apparent 1.05-year orbital period would place the gravitational-wave emission directly in the pulsar timing band. Since the first pulsar timing array study of 3C 66B, revised models of the source have been published, and timing array sensitivities and techniques have improved dramatically. With these advances, we further constrain the chirp mass of the potential supermassive black hole binary in 3C 66B to less than (1.65 ± 0.02) × 109 M using data from the NANOGrav 11-year data set. This upper limit provides a factor of 1.6 improvement over previous limits and a factor of 4.3 over the first search done. Nevertheless, the most recent orbital model for the source is still consistent with our limit from pulsar timing array data. In addition, we are able to quantify the improvement made by the inclusion of source properties gleaned from electromagnetic data over "blind" pulsar timing array searches. With these methods, it is apparent that it is not necessary to obtain exact a priori knowledge of the period of a binary to gain meaningful astrophysical inferences.

Export citation and abstract BibTeX RIS

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

1. Introduction

Continuous gravitational waves (GWs), defined by single-source cyclic GW emission, are expected to arise from the supermassive black hole binaries (SMBHBs) that form during a galaxy merger. When an SMBHB evolves such that it emits GWs in the microhertz-to-nanohertz GW band (orbital periods of weeks to several decades), a sufficiently massive and/or nearby SMBHB may be detectable by pulsar timing arrays (PTAs; e.g., Aggarwal et al. 2019, hereafter A19).

While GWs from individual sources in the PTA regime have been sought after in multiple works (Jenet et al. 2004; Arzoumanian et al. 2014; Feng et al. 2019) through a variety of methods, none have been detected. However, numerous advances have been made in the field of pulsar timing. As PTA experiments gain longer time baselines and higher cadences and the number of millisecond pulsars grows, sensitivity to GW sources increases. The notable ongoing PTA programs in the world include the European PTA, Parkes PTA, and the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) (e.g., Manchester et al. 2013; Desvignes et al. 2016; Arzoumanian et al. 2018b, respectively). Altogether, these PTAs time approximately 100 pulsars to high precision with the goal of GW detection, among other endeavors (e.g., Caballero et al. 2018; Hobbs et al. 2020).

In addition, sophisticated GW detection methods have been developed to detect quadrupolar continuous-wave (CW) signals in the data of coordinated timing arrays (e.g., Zhu et al. 2014; Babak et al. 2016; A19). However, past analyses that used the most up-to-date methods have used "blind" detection methods; that is, the software did not consider any binary model information gained from electromagnetic data to directly benefit the search. In comparison, most works that do target limits on specific sources using electromagnetic information have used smaller data sets consisting of a single pulsar with a periodogram approach (Jenet et al. 2004; Feng et al. 2019) or have been derived from the stochastic GW background (Zhu et al. 2019) rather than the full GW analysis pipeline. Here we have combined these methods in the first search of this type, where we used the entire NANOGrav array of pulsars and full GW search analysis, while incorporating electromagnetic data to conduct a more informed search for GWs from our test source, 3C 66B.

Since the report of a hypothesized orbital motion in the core of the galaxy 3C 66B by Sudou et al. (2003, hereafter S03), it has been an ideal test case for searches for GWs from SMBHBs. Using long-baseline interferometry, the authors found apparent elliptical motions in 3C 66B's radio core, modeling this motion as the gyration of the jet nozzle due to an orbit-induced precession of the smaller black hole's jet. They proposed a period and chirp mass for the binary of 1.05 ± 0.03 years and 1.3 × 1010 M, respectively. Given the relatively small redshift of the galaxy (z = 0.02126), a binary with those properties would be emitting gravitational radiation well within the sensitivity of PTAs.

As such, 3C 66B has long been a prime candidate for continuous GW detection. It was the first object targeted for CW detection, as reported by Jenet et al. (2004, hereafter J04), in which 7 years of Arecibo timing data from PSR B1855+09 (Kaspi et al. 1994) were used to search the Fourier domain timing residuals (commonly referred to as a Lomb–Scargle periodogram), using harmonic summing (Press et al. 1992), for a GW signal consistent with the binary period modeled by S03. With these methods, they did not see evidence of a significant signal and were able to place an upper limit of 7 × 109 M on the chirp mass of the system at a binary eccentricity of e = 0.

Since the study of J04, Iguchi et al. (2010, hereafter I10) reported a 93 day variability in the active galactic nucleus's millimeter light that was interpreted as likely due to Doppler boosting of a relativistic outflow that is modulated by orbital motion (its period differs due to geometric effects). The new model assumed the 1.05-year orbital period from S03 but predicted an updated chirp mass of 7.9 × 108 M, almost a full order of magnitude lower than the upper limit set by J04.

The galaxy 3C 66B was also one of the objects targeted by Zhu et al. (2019), who used a novel approach to test 3C 66B indirectly by using the source to predict the GW background strength implied by this source's existence. They concluded that the I10 model produced GW backgrounds that were larger than are currently probed by PTAs, implying that the source was not likely to be a binary with parameters as proposed by I10.

The work reported here presents a Bayesian cross-validation framework in which we use 3C 66B's binary parameter measurements as priors for our CW search. Our search has resulted in the most stringent direct GW-derived limit to date on the chirp mass of 3C 66B's candidate SMBHB. We also test, more generically within our search framework, what sensitivity improvements can be gained by knowing the GW frequency of a target to increasingly good precision.

Therefore, we have quantified the improvement made by searching for GWs from a specific source, including cases where the orbital period is only known with a large error or not known at all.

Note that because J04 used only one pulsar in their study, they were unable to perform a formal experiment to detect 3C 66B, as the use of one pulsar precludes the ability to demonstrate the quadrupolar signature that is unique to the influence of GWs. Thus, our study here is the first formal targeted detection experiment for 3C 66B using a PTA.

This paper is laid out as follows. In Section 2, we describe our data, mathematical model, and software pipeline. In Section 3 we report the detection Bayes factor and chirp mass upper limit for 3C 66B, as well as the results of new test methods. In Sections 4 and 5, we present our conclusions, as well as discuss implications for future detection prospects of this and other SMBHBs.

2. Analysis Methods

2.1. Pulsar Timing and Electromagnetic Data

We make use of the NANOGrav 11-year Data Set (Arzoumanian et al. 2018b), which provides high-precision timing of 45 ms pulsars. Only the 34 pulsars with baselines of at least 3 years are used for GW detection analyses (Arzoumanian et al. 2018a). We describe slight differences in the use of the data set in this work as compared to other papers in Section 2. However, the majority of the data are treated similarly to A19. Due to the 11-year timing baseline, the data set is most sensitive to binaries with orbital periods of less than a decade.

The electromagnetic data we incorporate into our models are mainly derived from S03 and I10, as well as the location from the NASA/IPAC Extragalactic Database (NED).40 These values are summarized in Table 1. The R.A., decl., and luminosity distance are taken as constants throughout the analysis, as the PTA sensitivity to sky location and distance is much lower than any associated errors. For consistency with earlier work, we take the luminosity distance of 3C 66B to be 85 Mpc, as in S03. Therefore, all calculations use ${H}_{0}=75\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$. Note that minor differences in the distance due to different reports of redshift or H0 cause only a small fractional variation in the results. If the fractional change in the luminosity distance is defined as

Equation (1)

any GW strain limit can be converted to the reader's preferred distance by multiplying the strain by d85 and ${ \mathcal M }$ limits by multiplying by ${d}_{85}^{3/5}$.

Table 1.  GW Model Values and Uncertainties

Parameter Value Reference
Chirp mass (${ \mathcal M }$) ${7.9}_{-4.5}^{+3.8}\times {10}^{8}\,{M}_{\odot }$ I10
GW frequency (${f}_{\mathrm{GW}})$ 60.4 ± 1.73 nHz S03
Redshift (z) 0.02126 Huchra et al. (1999)
R.A. 02h 23 m 11.4112 s Fey et al. (2004)
Decl. +42d 59 m 31.384s Fey et al. (2004)
GW strain (h) ${7.3}_{-5.8}^{+6.8}\times {10}^{-15}$ S03; I10

Download table as:  ASCIITypeset image

2.2. Signal Model

We use the methods presented in A19 for the generation of expected pulsar timing residuals influenced by a signal from a continuous GW from a circular SMBHB. While we will not present the full derivation, we will summarize below the relevant equations needed to follow our analysis of the NANOGrav data and refer the reader to A19 for more detail. Note that throughout this section, equations are written in natural units (where G = c = 1).

Pulsar timing residuals describe the deviation of an observed pulse arrival time from that predicted from a model based on spin, astrometric, interstellar delay, and, if needed, binary parameters of the pulsar. These are the basic data product that we use to search for GWs, which will not be included in the pulsar's timing model. A vector of timing residuals (δt) that is fit without a GW for each pulsar is modeled as

Equation (2)

where M is the design matrix, which describes the timing model, and epsilon is a vector of the linearized timing model parameter offsets from the best-fit solution. In other words, the timing model, which was originally derived without the presence of a GW, must now be adjusted. We write a vector describing the white noise in the data as nwhite, and the same for the red noise, nred, which is correlated over long timescales. The noise terms are described in more detail in Section 2.3.

The signal s can be derived as follows. For a GW source whose sky location is described by polar and azimuthal angles θ and ϕ, the strain induced by the emitted GWs is written in terms of two polarizations as

Equation (3)

where ${h}_{+,\times }$ are the polarization amplitudes and ${e}_{{ab}}^{+,\times }$ are the polarization tensors, which we write in the solar system barycenter (SSB) frame as

Equation (4)

Equation (5)

(Wahlquist 1987). In these equations, we define $\hat{{\rm{\Omega }}}$ as a unit vector pointing from the GW source to the SSB, written as

Equation (6)

We define the vectors $\hat{m}$ and $\hat{n}$ as

Equation (7)

Equation (8)

The pulsar's response to the GW source is described by the antenna pattern functions (Sesana & Vecchio 2010; Ellis et al. 2012; Taylor et al. 2016 and references therein),

Equation (9)

Equation (10)

where $\hat{p}$ is a unit vector pointing from the Earth to the pulsar.

Finally, we write the signal s induced by the GW, as seen in the pulsar's residuals, as

Equation (11)

Here ${\rm{\Delta }}{s}_{+,\times }$ represents the difference between the signal induced at the Earth (the Earth term) and that at the pulsar (the pulsar term) and can be written as

Equation (12)

where t is the time at which the GW passes the SSB, and tp is the time at which the GW passes the pulsar.41 These times can be related from geometry by

Equation (13)

where L is the distance to the pulsar.

For a circular binary at zeroth post-Newtonian order, ${s}_{+,\times }$ is given by (Wahlquist 1987; Corbin & Cornish 2010; Lee et al. 2011)

Equation (14)

where i is the inclination angle of the SMBHB, ψ is the GW polarization angle, dL is the luminosity distance to the source, and ${ \mathcal M }$ is the chirp mass, which is related to the two black hole masses as

Equation (15)

It is important to note that ${ \mathcal M }$ and ω, in this case, refer to the observed redshifted values.

For a circular binary, we relate the orbital angular frequency to the GW frequency with ${\omega }_{0}=\pi {f}_{\mathrm{GW}}$, where ${\omega }_{0}=\omega \left({t}_{0}\right)$. For this work, as in A19, we define t0 as the last MJD in the 11-year data set (MJD 57387). The orbital phase and frequency of the SMBHB are given by

Equation (16)

Equation (17)

where Φ0 and ω0 are the initial orbital phase and frequency. As in A19, we use the full expression for ω(t) to maintain consistency across runs, as this form is needed to model the signal at the higher frequencies sampled in some runs, as described in Section 2.4.3.

2.3. Software and Analyses

In this work, we make use of NANOGrav's GW detection package, enterprise,42 an open-source code written fully in Python that contains a built-in interface with the pulsar timing data and noise models required to perform Bayesian GW analysis (Arzoumanian et al. 2018a; limits and detection). The basic algorithms for Bayesian CW analysis are described in detail in a number of past works (e.g., Ellis 2013; Ellis & Cornish 2016).

Using enterprise, we can use a priori constraints on a binary system that come from electromagnetic observation (for instance, the period of 3C 66B) to set priors on GW parameters that are derived from the binary model. Within enterprise, we can easily add these priors to the timing and noise models to obtain a full model of the signal. We then perform the Markov Chain Monte Carlo (MCMC) methods implemented in PTMCMCSampler43 to find the posterior distribution for each of the free parameters. For "blind" CW searches, as in A19, we typically set uninformative priors, which are uniform across the allowed range of values, for the binary system's parameters, such as sky location, frequency, mass, and distance to the source. Thus, the methods here could be considered a "targeted" search by our use of informed priors.

For instance, in the simplest treatment of 3C 66B, a specific binary model has been hypothesized, with measurements and associated unknowns in the mass, mass ratio, and orbital frequency (e.g., S03; I10). We can use these electromagnetically constrained parameters, in addition to knowledge of the location of this object on the sky, to restrict our priors.44

Assuming an SMBHB with a circular orbit, a continuous GW signal can be characterized by eight of the following nine parameters,

Equation (18)

which represent the GW source's

  • 1.  
    and 2. position on the sky $(\theta ,\phi );$
  • 3.  
    GW frequency, related to the orbital frequency at some reference time $({f}_{\mathrm{GW}});$
  • 4.  
    orbital phase at some reference time (Φ0);
  • 5.  
    GW polarization angle (ψ);
  • 6.  
    orbital inclination (i);
  • 7.  
    chirp mass (${ \mathcal M }$);
  • 8.  
    luminosity distance (dL); and
  • 9.  
    strain amplitude (h0), which is related to the chirp mass, GW frequency, and luminosity distance.

The ninth parameter is redundant, as the strain amplitude h0 can be defined by

Equation (19)

(Sathyaprakash & Schutz 2009), where

and related to other physical parameters by

Equation (20)

Since the strain is entirely determined by ${ \mathcal M }$, fGW, and dL, a limit on h0 based on a PTA search can be translated into constraints on these source parameters. Since the uncertainties on θ, ϕ, and dL are much smaller than the PTA sky localization accuracy, by targeting a specific source with a known position and redshift, we can set these parameters as constant values and therefore reduce the number of search parameters to five.

In all runs, there is also a set of free parameters associated with each pulsar included in the PTA that are varied in the analysis. The first of these is the pulsar distance, which has a Gaussian prior in all cases. In pulsars where the distance is reported in Verbiest et al. (2012), the Gaussian is defined using the recognized distance and associated error. For the remaining pulsars, the Gaussian prior is set to a fiducial 1.0 ± 0.2 kpc, which is consistent with the distribution of distances and uncertainties obtained from Verbiest et al. (2012). Although this range does not necessarily encompass the actual distances to most of these pulsars, it works as a proxy value, and the choice of this value does not affect our results. As in A19, this assumption can be seen to hold in the posteriors for these pulsars, as the prior is returned in all cases, meaning that this analysis cannot inform on the distances for these pulsars. This is expected, as these pulsars are largely those with shorter observation baselines that are influencing the PTA to a smaller degree. The recovered pulsar distances also affect the GW frequency difference between the Earth and the pulsar, which therefore will be related to the chirp mass. When a wide range of chirp mass values are allowed by the data, the uncertainty in the pulsar distances is not significant to the final result of the search. Additionally, for small chirp masses, for even a large change in the distance to the pulsar, the change induced in the GW frequency at the pulsar is well below the resolution limit of the PTA (1/Tobs). This angular frequency at the pulsar can be calculated as

Equation (21)

where dp is the distance from the Earth to the pulsar.

Also included is the GW phase at the pulsar. While this quantity could be calculated geometrically from the other parameters, including it as a search parameter mitigates potential issues in sampling the complex parameter space that arise due to the large uncertainty on the distances to the pulsars compared to the GW wavelength.

As is standard for these types of analyses (e.g., Arzoumanian et al. 2018a; A19), the white noise of each pulsar (described as EFAC, EQUAD, and ECORR) is held fixed. The power spectral density of the pulsar intrinsic red noise is modeled as

Equation (22)

where Ared (the red-noise amplitude) and γ (the red-noise spectral index) are also allowed to vary in each pulsar in our MCMC simulation. Here fyr is 1/(1 year) in Hz. To assist the sampler, empirical distributions of the red-noise parameters were made from single pulsar noise run posteriors and used to create jump proposals. These determine how steps in the MCMC are taken through generating proposed samples and were added to significantly improve sampling and decrease burn-in time for our analyses. For a more detailed description, see Appendix A of A19.

Our treatment of the red noise in one pulsar, J0613−0200, required additional noise modeling. As described in A19, this pulsar possesses extra unmodeled noise processes that, in the 11-year CW search, presented as an increase in the strain upper limit at a frequency of 15 nHz. In this work, this manifested as poor sampling in the CW parameters, particularly in fGW. Because of this poor sampling, the fGW parameter would periodically get stuck near this frequency. Due to this pulsar's location relative to 3C 66B, which places it among the 10 pulsars with the highest antenna pattern response amplitudes, it is important to find a robust solution to these issues rather than remove the pulsar from the analysis. To mitigate this effect, we applied more sophisticated noise modeling techniques to allow the red noise to deviate from the typical power law, with corresponding jump proposals to assist sampling. The noise model that was chosen is a t-process spectrum, which allows for "fuzziness" in the typical power-law spectrum by scaling the power spectral density by a variable factor for each frequency. This model is created by generalizing the typical Gaussian process prior to a Student's t-distribution. This process will be discussed in more depth in J. Simon (2020, in preparation) and, due to increasingly complex data, will likely become more typical in future analyses.

Even with this model, poor sampling in the fGW parameter still occurred and can be attributed to unmodeled noise due to changes in the dispersion measure of pulsar J1713+0747 caused by variations in the interstellar medium along the line of sight (Lam et al. 2018; Hazboun et al. 2020). While this pulsar is NANOGrav's most sensitive in general, it is not particularly sensitive to 3C 66B, as shown in Figure 1; thus, excluding it did not significantly affect the upper limit on target 3C 66B. As such, this pulsar was removed from our search.

Figure 1.

Figure 1. Sky map depicting the antenna pattern response amplitude (${F}_{\times }^{2}+{F}_{+}^{2}$) due to a GW located at the sky position of 3C 66B. Also plotted are the locations of the 34 pulsars used in GW analyses of the NANOGrav 11-year data set, with the two pulsars in need of special attention noted with different colors.

Standard image High-resolution image

The above procedure is used for all enterprise runs, as described in detail in the next subsection.

2.4. Four Distinct Tests

We constructed several separate setups for enterprise for the purpose of testing distinct hypotheses. The purpose of each of these, and the difference in procedures within enterprise, is described below.

2.4.1. Detection

To determine if a CW from 3C 66B is detected, we conduct an enterprise search using a single frequency, with a value corresponding to the 1.05-year orbital period for a circular binary, making the final set of search parameters

Equation (23)

Due to the frequency resolution of the PTA, which is defined by the timing baseline, it is reasonable to set a parameter with errors of this magnitude (Table 1) to a constant value. However, we will explore the relaxation of this assumption in later sections. Note that the I10 and S03 models make assumptions about the electromagnetic data that may or may not be correct; our model simply tests the presence of an SMBHB in this system at a period of 1.05 years.

The detection prior on ${ \mathcal M }$ is log-uniform in the range 107 to 1010 M and sampled in log space. This prior is convenient for calculating Bayes factors as a measure of detection significance using the Savage–Dickey formula (Dickey 1971),

Equation (24)

Here ${{ \mathcal H }}_{1}$ is the model with a GW signal plus individual pulsar red noise, and ${{ \mathcal H }}_{0}$ is the model with individual pulsar red noise only. The prior and posterior volumes at h0 = 0 are $p\left({h}_{0}=0| {{ \mathcal H }}_{1}\right)$ and $p\left({h}_{0}=0| { \mathcal D },{{ \mathcal H }}_{1}\right)$, respectively. We are able to apply the Savage–Dickey formula because these models are nested (${{ \mathcal H }}_{0}$ is ${{ \mathcal H }}_{1}$ where h0 = 0), and $p\left({h}_{0}=0| { \mathcal D },{{ \mathcal H }}_{1}\right)$ is approximated as the fraction of quasi-independent samples in the lowest-amplitude bin of a histogram of h0. The error in the Bayes factor is computed as

Equation (25)

where n is the number of samples in the lowest-amplitude bin. This process is done once the samples in the GW strain are calculated from the directly sampled parameters. In the detection analyses, the red-noise amplitude is sampled with a matching prior (log-uniform in ${A}_{\mathrm{red}}$). All other GW parameters are searched with a uniform prior.

2.4.2. Upper Limits

To set an upper limit on the chirp mass of 3C 66B, we again conduct an enterprise search using a single frequency with a value corresponding to the 1.05-year orbital period, making the final parameter set as in the previous section (Equation (23)). However, in contrast with the case for detection, the upper-limit prior on ${ \mathcal M }$ is uniform (rather than log-uniform), meaning the prior set on ${\mathrm{log}}_{10}{ \mathcal M }$ exponentially increases over the range {7, 10}. This is done as an astrophysically reasonable prior, as we expect SMBHBs to lie anywhere in this mass range, while still allowing for efficient sampling. Additionally, this prior choice allows the derived upper limit to be as conservative as possible by allowing a higher proportion of high chirp mass samples and to be independent from the choice of lower prior bound. In the upper-limit analyses, the red-noise amplitude is sampled with a matching prior (uniform in Ared). Upper limits are taken to be the value of the 95th percentile of the posterior distribution. Following the approach of Arzoumanian et al. (2018a), we calculate the error on the upper-limit calculations as

Equation (26)

where x = 0.95, and Ns is the number of effective samples in the chain, which is estimated by dividing the total number of samples by the autocorrelation length of the chain.

2.4.3. Frequency Prior Testing

In addition to the tests described above of the S03 and I10 models, where the GW frequency is fixed to discrete values, as in other CW searches (A19; Arzoumanian et al. 2014), it is also crucial to test frequencies within the confidence region of these values. For this aim, we have developed methods to directly sample in fGW. These include specialized parameter groupings and jump proposals to help the sampler move through the more complex parameter space. Using these techniques, we are able to obtain an upper limit from the ${ \mathcal M }$ posteriors for a variety of frequency priors from various enterprise setups.

When searching over the GW frequency, a log-uniform chirp mass prior is used, and the samples are reweighted during upper-limit calculations to modify the prior choice from a uniform-in-log distribution of masses to a uniform-in-linear distribution of masses, the latter of which is more common in upper-limit analyses by virtue of insensitivity to the lower sampling boundary. This both assists with sampling and maintains a consistent prior on the GW strain, which is not directly sampled. To match the ${ \mathcal M }$ prior, a log-uniform prior is used on Ared. Since we are no longer fixing fGW to a single value, our final parameter set for these searches is

Equation (27)

In addition, we also choose to limit our GW frequency prior to a range of 1–100 nHz, rather than the 1–300 nHz used in A19. Besides the PTA's insensitivity at these high frequencies, we expect a source to remain in these frequency bins for very little time, with residence timescales as small as months, so their detection prospects are minimal (Burke-Spolaor et al. 2019; Hazboun et al. 2019b).

Using the three priors shown in Table 2, we are able to find reweighted upper limits for a variety of scenarios. These include

  • 1.  
    the GW frequency is known and set to a single value;
  • 2.  
    the GW frequency is known with large errors, and the error region is searched over; and
  • 3.  
    the GW frequency is not known or has significant uncertainty, and the entire PTA sensitivity band is searched over.

Then we examined the change in the reweighted chirp mass upper limit as a function of frequency prior width. In addition to allowing for possible errors in the orbital period measured by S03, these widened priors allow us to test the feasibility of this process on a less constrained source. Additionally, if there was any significant frequency evolution in the source, a signal would still have the chance to be detected in either of these setups. In addition to a single value and a uniform prior across the PTA sensitivity bandwidth, we also use 10 times the uncertainty on the predicted frequency as an example of a search with significant uncertainty. We also bin the samples of the widest fGW search to interpolate between these three individual prior widths. The results of this examination are described in Section 3 and summarized in Figure 5.

Table 2.  Frequency Prior Testing Weighted Upper Limits

Scenario fGW Prior Weighted ${ \mathcal M }$ Upper Limit (109 M)
1 Constant 1.57 ± 0.02
2 10σ 1.54 ± 0.01
3 Log-uniform 8.68 ± 0.07

Download table as:  ASCIITypeset image

2.4.4. Test of a Specific Binary Model

To directly test the consistency of the model presented in I10 with the NANOGrav data, we create priors for an enterprise run corresponding to the values presented (see the first line of Table 3). For fGW, we are able to use a Gaussian prior, where the error on the measured value from I10 directly corresponds to the standard deviation of the prior. However, ${ \mathcal M }$ has uneven error bars, so a more complicated prior is needed. Here we fit a skewed normal distribution to the reported value and error and construct a skewed normal prior based on this distribution; we also keep a log-uniform prior on Ared. Therefore, the final parameter set for this search was

Equation (28)

Table 3.  Model Testing Prior and Posterior Values

  log(Frequency) log(Chirp Mass)
Iguchi (prior) −7.219 ± 0.012 ${8.90}_{-0.24}^{+0.21}$
This work (posterior) $-{7.217}_{-0.013}^{+0.012}$ ${8.87}_{-0.24}^{+0.16}$

Download table as:  ASCIITypeset image

To analyze the amount of information gained between the prior and posterior models, we employed the Kullback–Leibler (K–L) divergence (Kullback & Leibler 1951). We calculate this information gain in bits between the posterior $p(x| d)$ and the prior p(x) as

Equation (29)

This is done for the distributions for both ${ \mathcal M }$ and fGW. To maintain consistency between forms of the posterior and the prior, we fit a skewed normal distribution to both posteriors to directly compare to the prior.

3. Results

The results discussed in this section can be reproduced, and the MCMC data examined, using code provided for the reader's convenience.45

3.1. Detection

Using the setup for a detection run as described in Section 2.4.1, we find no evidence for a GW signal from 3C 66B. We calculate a Savage–Dickey Bayes factor of ${{ \mathcal B }}_{10}\,=0.74\pm 0.02$. Therefore, there is no evidence for the detection of a GW signal in the data. The posterior for this run is plotted in Figure 2.

Figure 2.

Figure 2. Posterior for the detection analysis described in Section 3.1 (blue). The vertical orange region describes the area of parameter space where a signal with the parameters found by I10 would lie. While the upper end of the parameter space is ruled out, there is clearly no value that is preferred by the sampler.

Standard image High-resolution image

3.2. Upper Limits

As no GW signal is detected from 3C 66B, we set upper limits on the chirp mass using the procedure described in Section 2.4.2. Using the constant-value frequency prior at 60.4 nHz (corresponding to the 1.05-year orbital period), we set a 95% upper limit of (1.65 ± 0.02) × 109 M for the ${ \mathcal M }$ of the SMBHB in 3C 66B. This value corresponds to a strain of (2.47 ± 0.05) × 10−14. To compare, the expected strain of the model in I10 is $({7.2}_{-5.8}^{+6.8})\times {10}^{-15}$. As can be seen in Figure 3, while we achieve a factor of 4.3 improvement over the limit set by J04, we cannot rule out the I10 model. The posterior distribution of samples does include a peak at about 1 ×109 M, which is within the error region for the chirp mass calculated from I10. However, this peak is not statistically significant and is able to be traced to a single pulsar, J1909−3744. By examining the posterior distributions constructed from samples corresponding to this peak, we find structure in the GW phase posterior at J1909−3744 that does not occur for any other pulsar. This likely occurs due to covariances between the model and sinusoidal behavior caused by noise processes in the data, as a real GW signal would be recovered by more than one pulsar. Therefore, this peak in the posterior is not indicative of a signal, and our upper limit can be considered robust. We will note that the upper limit listed can be calculated for the reader's preferred distance using the transformation described in Section 2.1.

Figure 3.

Figure 3. The chirp mass posterior histogram is plotted in blue, with a vertical line depicting the 95% upper limit. Shown in orange is the chirp mass upper limit of I10, with the shaded region representing the error on the value. With these methods, the I10 mass estimate is impossible to rule out. We also note that the peak in the posterior at $1\times {10}^{9}\,{M}_{\odot }$ is not statistically significant.

Standard image High-resolution image

3.3. Frequency Prior Testing

As described in Section 2.4.3, we also performed tests to quantify how much our upper limits might improve if we constrain (through electromagnetic observation) the orbital frequency of the target. While for 3C 66B, the orbital frequency is assumed to be known to within small errors, for other targets, a frequency may not be known or be only poorly constrained. This test provides a sense of how well the period must be constrained to provide effective sensitivity gains for a GW search.

Using the three scenarios described above, we are able to characterize the change in the reweighted upper limit between the setups. The result of the log-uniform prior search over the entire frequency band is summarized in Figure 4. The white area represents the area of ${ \mathcal M }$fGW parameter space ruled out in this analysis. From the uniformity of the samples over the parameter space, it is clear that there are no sampling issues. This is due to the improved sampling methods described in Section 2.3. The weighted 95% upper limit is plotted for each frequency bin, allowing us to quantify for which frequencies we are the most sensitive to 3C 66B. We note that for the very lowest frequencies, the upper limit is dependent upon the choice of prior, as the search cannot rule out any of the prior range.

Figure 4.

Figure 4. The 2D histogram of samples in the log-uniform prior setup. Also plotted is the weighted 95% upper limit for each frequency bin (blue) from the scenario 3 setup. The white area indicates the section of parameter space ruled out by our search. It is clear from the uniform distribution of samples across all frequency and mass channels that all sampling issues have been resolved. This uniform distribution also makes clear that there is no indication of a signal at the distance and sky location of 3C 66B. We only plot the upper half of the parameter space in ${ \mathcal M }$ to resolve more detail. Below ${\mathrm{log}}_{10}{ \mathcal M }=8.5$, all sampling is uniformly distributed, identical to the upper half of the figure. For comparison, the scenario 1 weighted upper limit (orange triangle) and I10 chirp mass estimate (red star) are also shown.

Standard image High-resolution image

In addition to the three runs described above, it was also possible to infer the upper limit that would be derived from a run with a frequency prior width between those of the three separate runs. To do this, we bin the samples in the scenario 3 (widest ${f}_{\mathrm{GW}}$ prior) run to keep only a certain range of frequencies and recalculate the weighted upper limit for this subset. These bins increase symmetrically in log space about the value of fGW reported by S03, from a log space full width of 0 dex (a constant) until the upper bound reaches fGW = 100 nHz. After this, only the lower bound expands to reach a full log space full width of 2 dex (essentially, 2 orders of magnitude in linear space). The weighted upper limits calculated from these binned samples are plotted in Figure 5.

Figure 5.

Figure 5. Chirp mass upper limits plotted with respect to frequency prior width (blue). Also shown as horizontal lines are previous upper limits set by S03 (red), J04 (green), and I10 (orange). Shaded regions describe error bars on the quoted limit. It is clear that none of these upper limits rule out that of I10. However, this figure accentuates the fact that when a period is known to less than 1 order of magnitude of precision, the limits on the target's mass improve by nearly 1 order of magnitude; that is, while the tightest prior produces the lowest upper limit, moderately wide priors also produce similar results, indicating that perfect orbital models would not be necessary to perform such a search on other systems. It is not until the prior spans approximately an order of magnitude that sensitivity is lost. Also plotted for comparison are the weighted upper limits for each of the three separate runs.

Standard image High-resolution image

Also plotted in Figure 5 are the upper limits from the three individual runs. From the consistency of these points with the calculated curve, it is clear that this technique is robust. Additionally, this shows the feasibility of searching over fGW, as the results are consistent with those calculated for both an individual frequency and a small range.

As can be seen in Figure 5 and Table 2, there is nearly an order of magnitude difference in the upper limits derived from frequency-varied runs of different prior widths. Because the upper limits at the very lowest frequencies are dependent upon the prior choice, the difference seen here is a lower limit. However, from the curve calculated from binned samples, we see that this increase does not begin until about 1 order of magnitude in frequency space about the I10 value is included. This implies that a targeted search such as this is worthwhile even without exact orbital information, as long as the frequency is known to within an order of magnitude.

3.4. Test of a Specific Binary Model

To directly test our sensitivity to a GW from the model of 3C 66B proposed in I10, we directly test priors as described in Section 2.4.4. In Figure 6, we can compare the prior and posterior for both fGW and ${ \mathcal M }$. These distributions are quantified in Table 3, where the errors on the posterior values are calculated with the percentiles of the posterior distribution corresponding to 1σ error bars. The values of fGW are consistent with those of the prior, but for ${ \mathcal M }$, we are able to significantly lower the upper bound on the value, effectively ruling out part of the high-mass region of the model.

Figure 6.

Figure 6. Posteriors (blue) and priors (orange) for the direct test of the model presented in I10. Vertical bars mark the 16.86th, 50th, and 84.13th percentiles of each, to represent the 1σ error bars.

Standard image High-resolution image

Additionally, we report the information gained between the posterior and the prior as described in Section 2.4.4. The differences in the distributions for fGW produce a K–L divergence of 0.0096, while those of the ${ \mathcal M }$ distributions produce a K–L divergence of 0.2597. While neither of these values is large, it is clear that much more information is gained about the chirp mass of 3C 66B from this model test than the GW frequency.

4. Discussion

To provide context for the upper limit on 3C 66B set in this work, we can compare to the limits set in A19, which do not have the benefit of electromagnetic constraints (i.e., a "blind" search). This comparison will allow us to estimate the improvement in sensitivity gained by including electromagnetic data over a typical blind search. By comparing our strain upper limit of (2.47 ± 0.05) × 10−14 to the sensitivity curve in Figure 3 of A19, where the strain upper limit at the nearest searched frequency is 5.3 × 10−14 nHz, we observe that we have gained a factor of 2.1 in sensitivity by holding the source position fixed in our search. Note that a much greater improvement comes from knowing the binary candidate's period, as demonstrated in Figures 4 and 5.

With the framework developed in Hazboun et al. (2019b), we can construct detection sensitivity curves to estimate the PTA that will be required to detect or rule out the mass model presented in I10. The hasasia (Hazboun et al. 2019a) package46 allows us to construct these detection sensitivity curves using a straightforward matched filter statistic and simulate PTA data with control over the number of pulsars, observing cadence, timing precision, and data length. Using this software to estimate an idealized signal-to-noise ratio (S/N; see Equation (79) in Hazboun et al. 2019b), assuming the parameters in I10 and using the pulsar noise parameters in Arzoumanian et al. (2018b), we obtain S/N = 0.87. We used this software to extend the baseline of the existing 11-year NANOGrav data set by adding new data to the existing pulsars with a timing precision and cadence that match recent data. We also augmented the PTA, adding new pulsars with timing precisions and cadences similar to those already in the array; we added pulsars for each projected year at a rate comparable to the current growth rate of NANOGrav, which has been approximately seven pulsars per year for the past three data sets.

We find that NANOGrav should be able to detect or rule out the existence of an SMBHB in 3C 66B with the I10 mass within 5–8 years from the end of the data set considered here. However, while hasasia allows us to calculate the PTA's sensitivity to a CW at a specific sky location, it is unable to set other parameters (such as luminosity distance) as known due to electromagnetic information about the GW source, as is done in this work. As is discussed above, including source parameters that are electromagnetically derived to reduce the parameter space of the GW search allows for increased sensitivity. Because of this, using electromagnetic information will likely allow us to accelerate this estimated timeline. To more reliably estimate this timescale, detailed simulation work will be necessary to quantify the improvement made by including electromagnetic information over typical searches.

Because the sensitivity of the array depends heavily on the observing baseline of each pulsar, the inclusion of additional data can help tremendously. Data of this sort are accessible through the IPTA (Perera et al. 2019), and follow-up analyses of 3C 66B by the international community could prove fruitful. This timeline for the PTA sensitivity required to confirm or deny 3C 66B as an SMBHB will be reduced with the more rapid addition of pulsars to the array, e.g., by adding more than seven per year. This improvement will be accelerated if the newly included pulsars are near the sky location of 3C 66B, as, currently, there are few pulsars in the array that are very sensitive to 3C 66B. To accomplish this, pulsar searches should be undertaken near the sky locations of potential PTA targets to begin improving our sensitivity more rapidly. Some pulsars in this area of the sky can be included through use of data provided by the IPTA (Perera et al. 2019), showing once again that an international effort to detect 3C 66B could be worthwhile.

In addition to the results for GWs from 3C 66B, our work has many implications for the detection prospects of other binary candidates. As discussed in Section 3.3 and shown in Figure 5, for 3C 66B, it was not until we widened our prior to span an order of magnitude in frequency space on either side of the target frequency that sensitivity was lost. For similar candidates, particularly those at similarly high orbital frequencies, we presume that this result will hold. Therefore, as long as the sky location and luminosity distance of a potential target are known, a search of this type is worth attempting if at least an estimate of an orbital period can be obtained. We will caution that this improvement will differ depending on the sky location of the source, and that the amount of frequency space that can be effectively searched with this method will be larger for higher-frequency sources. As can be seen in Figure 4, it is the inclusion of samples at low frequencies that raises the upper limit. However, typical errors on binary periods are quite a bit smaller than the limit suggested here, meaning that this method should prove useful for most binary candidates. This method will also account for any frequency error due to unaccounted-for frequency evolution of the SMBHB, which, in the case of a detection, would provide important constraints for evolutionary models.

5. Conclusions

In this work, we present a new method for performing multimessenger searches for individual SMBHBs, using 3C 66B as a test case. The galaxy 3C 66B was first identified as a binary candidate by S03 and first visited by PTAs in J04, which ruled out the proposed binary model. In the intervening 15 years, a revised model was published by I10, and PTA data and analysis methods have greatly improved. We used the NANOGrav 11-year data set, as well as the collaboration's flagship GW detection package, enterprise, to search for GWs from 3C 66B. Here we are able to limit 3C 66B's chirp mass, at 95% confidence, to (1.65 ± 0.02) × 109 M, a factor of 4.3 smaller than the limit set in J04. However, we are unable to rule out the existence of a binary corresponding to the revised model proposed in I10.

In addition to directly placing a limit on the chirp mass of 3C 66B for the published orbital period, we are able to quantify how much this multimessenger approach increases our sensitivity over a typical "blind" PTA search. We have conducted a search on real data that includes GW frequency as a free parameter, and from this analysis, we learn that by including frequency constraints from electromagnetic binary source measurements to restrict the prior, we can gain approximately an order of magnitude in sensitivity when compared to a frequency-blind search spanning the whole PTA band. However, this drop in sensitivity does not occur until a relatively wide range of frequencies is searched, meaning that this approach will be useful even for candidates with relatively poor constraints on their orbital periods.

Author contributions. We list specific contributions to this paper as follows. C.A.W. led the work on this paper, ran the GW searches, and led the development of the manuscript. J.S., S.R.T., S.J.V., and S.B.S. provided guidance throughout the project and key development of the project motivation and scientific interpretation. J.A.E., S.R.T., P.T.B., S.J.V., and C.A.W. designed and implemented the Bayesian search algorithms in enterprise. J.S.H. performed and interpreted the S/N simulations with hasasia. R.D.E. performed initial literature reviews on 3C 66B. N.J.C., J.S.H., D.L.K., M.T.L., T.J.W.L., M.A.M., C.M.F.M., and D.J.N. contributed valuable scientific comments. The NANOGrav data are the result of the work of dozens of people over the course of more than 13 years. Z.A., K.C., P.B.D., M.E.D., T.D., J.A.E., E.C.F., R.D.F., E.F., P.A.G., M.L.J., M.T.L., R.S.L., M.A.M., C.N., D.J.N., T.T.P., S.M.R., P.S.R., R.S., I.H.S., K.S., J.K.S., and W.Z. developed the 11-year data set. All authors are key contributing members of the NANOGrav collaboration.

Acknowledgments. S.B.S. and C.A.W. were supported in this work by NSF award Nos. 1458952 and 1815664. The NANOGrav collaboration is supported by NSF Physics Frontier Center award No. 1430284. C.A.W. acknowledges support from West Virginia University through the Outstanding Merit Fellowship for Continuing Doctoral Students. S.B.S. is a CIFAR Azrieli Global Scholar in the Gravity and the Extreme Universe program. This research made use of the Super Computing System (Spruce Knob) at WVU, which is funded in part by the National Science Foundation EPSCoR Research Infrastructure Improvement Cooperative Agreement No. 1003907, the state of West Virginia (WVEPSCoR via the Higher Education Policy Commission), and WVU. We acknowledge the use of Thorny Flat at WVU, which is funded in part by National Science Foundation Major Research Instrumentation Program (MRI) award No. 1726534 and WVU. NANOGrav research at UBC is supported by an NSERC Discovery Grant and Discovery Accelerator Supplement and the Canadian Institute for Advanced Research. M.V. and J.S. acknowledge support from the JPL RTD program. S.R.T. was partially supported by an appointment to the NASA Postdoctoral Program at JPL, administered by Oak Ridge Associated Universities through a contract with NASA. J.A.E. was partially supported by NASA through Einstein Fellowship grants PF4-150120. Portions of this work performed at NRL are supported by the Chief of Naval Research. The Flatiron Institute is supported by the Simons Foundation. Portions of this research were carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Data for this project were collected using the facilities of the Green Bank Observatory and the Arecibo Observatory. The Green Bank Observatory is a facility 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 by the University of Central Florida in alliance with Yang Enterprises, Inc., and Universidad Metropolitana. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This work made use of the online cosmology calculator tool (Wright 2006). We also acknowledge use of numpy (Oliphant 2006), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), and astropy (Price-Whelan et al. 2018). This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.

Footnotes

  • 40 

    The NED is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration and can be accessed at https://ned.ipac.caltech.edu/.

  • 41 

    This definition is occasionally written as the negative of the right side of the equation here, e.g., ${s}_{+,\times }(t)-{s}_{+,\times }\left({t}_{p}\right)$, as in A19. This is resolved with a change of convention in the definition of the GW antenna pattern, as we have done here; thus, all results are consistent between these works.

  • 42 
  • 43 
  • 44 

    Note that our restricted priors might not always be Gaussian; in some cases, electromagnetic observations of a source may produce a model that contains greater complexity than Gaussian error bars. In such cases, non-Gaussian priors must be used. The functionality exists in enterprise for studies that would require such a setup. As an example, if cyclic flux variability is observed, the period of variability might represent the fundamental orbital frequency, a harmonic, or even a resonance, requiring a multivalued prior. In our analysis, the reported errors on binary masses from I10 were asymmetric; thus, for some analyses, our chirp mass prior required an asymmetric distribution.

  • 45 
  • 46 
Please wait… references are loading.
10.3847/1538-4357/ababa1