Multimessenger pulsar timing array constraints on supermassive black hole binaries traced by periodic light curves
DDraft version September 25, 2020
Typeset using L A TEX twocolumn style in AASTeX62
Multimessenger pulsar timing array constraints on supermassive black hole binaries traced by periodic light curves
Chengcheng Xin, Chiara M. F. Mingarelli,
2, 3 and Jeffrey S. Hazboun Columbia University, Department of Astronomy, 550 West 120th Street, New York, NY, 10027, USA Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY, 10010, USA Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT 06269-3046, USA University of Washington Bothell, 18115 Campus Way NE, Bothell, WA 98011, USA (Revised September 25, 2020)
ABSTRACTSupermassive black hole binary systems (SMBHBs) emitting gravitational waves may be tracedby periodic light curves. We assembled a catalog of 149 such periodic light curves, and using theirmasses, distances, and periods, predicted the gravitational-wave strain and detectability of each binarycandidate using all-sky detection maps. We found that the International Pulsar Timing Array (IPTA)provides almost uniform sky coverage – a unique ability of the IPTA – and by 2025 will improveNANOGrav’s current minimum detectable strain by a factor of 6, and its volume by a factor of 216.Moreover, IPTA will reach detection sensitivities for three candidates by 2025, and 13 by the end of thedecade, enabling us to constrain the underlying empirical relations used to estimate SMBH masses. Wefind that we can in fact already constrain the mass of a binary in Mrk 504 to
M < . × M (cid:12) . Wealso identify 24 high-mass high-redshift galaxies which, according to our models, should not be able tohost SMBHBs. Importantly the GW detection of even one of these candidates would be an essentiallyeternal multimessenger system, and identifying common false positive signals from non-detections willbe useful to filter the data from future large-scale surveys such as LSST. Keywords: supermassive black holes – gravitational waves – quasars INTRODUCTIONGravitational waves (GWs) emitted by supermassiveblack hole binary (SMBHB) systems with periods ofyears to decades are expected to be the most powerfulsources of GWs in the Universe. Searches for individ-ual SMBHBs can be aided by looking for periodic lightcurves in Active Galactic Nuclei (AGN), which tracethe SMBHB activity and serve as electromagnetic (EM)counterparts to these GW sources (e.g. Kauffmann &Haehnelt 2000; Goulding et al. 2018).With decades of precision pulsar timing data, Pul-sar Timing Arrays (PTAs) are ready and able to makethe first detection of low-frequency GWs (Burke-Spolaoret al. 2019). Ultra-stable radio millisecond pulsarsare the best clocks in nature – deviations in their ar-rival times manifesting over years to decades can there-
Corresponding author: Chengcheng [email protected] fore signal the presence of passing low-frequency GWs(Hellings & Downs 1983).Identifying the host galaxy of a SMBHB system wouldyield an essentially eternal multimessenger system —the black holes are millions of years from coalescence,and would provide a treasure trove of science: directevidence that the infamous final parsec problem is solved(Begelman et al. 1980), clues as to how it was solved bythe quantity of gas, stars, and residual eccentricity inthe system (Milosavljevi´c & Merritt 2001; Sesana 2013),as well as probe General Relativity, e.g. O’Beirne et al.(2019).Here we assembled a catalog of 149 periodic lightcurves which may trace SMBHB activity (Graham et al.2015b; Charisi et al. 2016; Liu et al. 2019; Lehto & Valto-nen 1996; Arzoumanian et al. 2020), though in the endwe restrict our analysis to the Graham et al. (2015b)sample since we found that the other candidate binariesare too high-frequency to be detected.Using published binary parameters, we predict theGW strain of each binary candidate, and for the first a r X i v : . [ a s t r o - ph . GA ] S e p Xin, Mingarelli, Hazboun time predict their time to detection by constructingall-sky detection maps to simulate the IPTA in 2025(IPTA2025), Phase 1 of the Square Kilometre Array(SKA1), and its second phase (SKA2) using hasasia (Hazboun et al. 2019b), from hereon H19.These public data, simulations, and open-source codeswill help to inform where to search for new pulsarsto accelerate GW detection, and lead to targeted GWsearches – making these searches more sensitive by atleast factor of 2 (Arzoumanian et al. 2020), more if theGW frequency is known. GRAVITATIONAL-WAVE SIGNALSPreviously, the periodic light curves from Grahamet al. (2015b), from hereon G15, were used by Sesanaet al. (2018) to test the binary hypothesis: they showedthat if all the AGN light curves were really SMBHBs,that the inferred cosmic population of SMBHBs wouldcreate a GW background currently in tension with PTAupper limits (Arzoumanian et al. 2018b). Kelley et al.(2019) also carried out simulations which predict that,statistically, five of the
CRTS sources may contain gen-uine binaries.We now have new tools which we can apply directlyto the
CRTS sample for the first time: hasasia – apublic code which can create continuous GW (CGW)detection sky maps, and new results from Arzouma-nian et al. (2020) which show that targeted searches canboost CGW sensitivity by at least a factor of two. In-deed, it is only now possible to rigorously investigate thedetectability of each individual binary candidate’s GWsignal.The GW strain for a SMBHB is determined via h = 2 M / c ( πf ) / D L , (1)where M / c = q/ (1 + q ) M / is the chirp mass, M = M + M is the total binary mass, q = M /M ≤ D L is the luminosity distance to the sourceand f is the observed GW frequency. Many of the meanBH mass estimates are reported in previous literaturethat identified periodic sources, e.g. G15, Charisi et al.(2016); Liu et al. (2019). While G15 computed the strainin their Figure 10, they did not include errors on thetotal mass M (their error bars in Figure 10 are fromvarying 0 . ≤ q ≤ . h , since h ∝ M / ,and is therefore even more important than the error onthe mass ratio q .We compute the SMBH masses as in G15, using thewidths of H β and Mg-II spectral emissions (Shen et al.2008). We redo this calculation directly from SDSS spec- tra since G15 do not report mass errors which are cru-cial for determining the uncertainties in h . In order toget a handle on these errors, we generate 160 ,
000 MonteCarlo (MC) samples to create a distribution of the mass, M , and also sample in the ratio q from uniform distri-butions over the range 0.25 to a maximum of 1. Sincethe BH mass error dominates the uncertainty in h , thechoice of the lower limit on q is of less consequence assmaller values of q simply extend the lower error baron h . We interpret the peak of the distribution as themean, h and the error bars are taken as one standarddeviation (1 σ ) in the mass distribution. While only 37 CRTS sources have the available SDSS spectra neededfor this calculation, this method is generic and can beapplied to a larger sample when spectra are available.We report all our mass calculations and their error esti-mates in Table 3.The GW frequencies are obtained via reported binaryorbital periods P , i.e. f GW = 2 /P . The luminositydistance D L is computed via the redshift of the source,assuming cosmological parameters, H = 0 . − s − and Ω m = 0 . A ∼ × − (Arzoumanian et al. 2020). Potential in-dividual GW amplitudes in the datasets we investigatehere are much lower than this, and at best, at the levelof some of the most conservative GW background mod-els (Ryu et al. 2018; Bonetti et al. 2018; Sesana et al.2018; Zhu et al. 2019). We therefore study the impactof this stochastic background signal on the detection ofthe individual SMBHB signals. In Figure 1 we illus-trate the effect of a stochastic GW background on thesensitivity to single sources as a separate noise source,though the background is obviously tied in with individ-ual sources, especially at the detection threshold. Thenumber of resolvable binaries in a single frequency bin –the so-called confusion limit for PTAs – varies in the lit-erature from 4 (Babak & Sesana 2012) to 2 N − N is the number of pulsars in thePTA. Ever evolving data analysis techniques (B´ecsy &Cornish 2020) will allow individual sources to be pulledfrom the GWB as they become significant in more sen-sitive datasets. MA PTA Constraints on SMBHBs − − − Gravitational wave frequency [Hz] − − − − − G r a v i t a t i o n a l w a v e s t r a i n ( h ) SKA1 PSRIPTA2025 PSRSKA1SKA1 w/ GWB − − Gravitational wave frequency [Hz]10 − − − − − − − G r a v it a ti on a l w a v e s t r a i n ( h ) NANOGrav 11 yrIPTA 2025SKA1SKA2OJ2873C66B
Figure 1. Current and Future GW Detection Curves . Left: Individual SKA1 pulsar sensitivities and full PTA sensitivitycurves. Pulsars with short observing time spans considerably increase the sensitivity at higher frequencies. The atypical flatbucket of the curve comes from the addition of low RMS - short time span pulsars to IPTA pulsars. An unresolved GWbackground with A GWB = 1 . × − deteriorates our sensitivity to individual SMBHBs by at most a factor ∼ CRTS sources. Right: Truncated GW frequency regime relevant to
CRTS sources (orange circles with 1 − σ errorbars), OJ287, and 3C66B. The continuous curves are the S/N = 3 detection curves for NANOGrav, IPTA and SKA1 and 2experiments. We only plot candidates with strain h > − , GW strain less than this is unlikely to be detectable by PTAs. Since our catalog of candidate SMBHBs are mid-to-high GW frequency for PTAs, we also compute theirGW frequency evolution, ˙ f ,˙ f = 965 π / M / c f / . (2)The expert reader will notice that small values of ˙ f meanthat the Earth term and the pulsar term lie in the sameGW frequency bin, and will thus double a pulsar’s resid-ual. This has been taken into account in hasasia , H19. ELECTROMAGNETIC COUNTERPARTSHydrodynamical simulations have demonstrated thatbright quasar variability can be modulated periodicallyby the orbital period of the binary SMBH, due to per-turbations by the surrounding gas in the circumbinaryaccretion disk, see e.g. Roedig et al. 2014; Farris et al.2014b; D’Orazio et al. 2016 and references therein.A promising observational approach to identify SMB-HBs in the optical bands is to search for quasars withperiodic variability. One of the first candidate SMBHBs,OJ287, was identified with variable observed luminosityfluctuations in the form of repeating outbursts occurringevery ∼
12 years (Lehto & Valtonen 1996; Laine et al.2020). Another prominent object is an unequal-massSMBHB candidate in quasar PG1302-102. It emergedfrom a systematic search in
CRTS for its quasi-sinusoidaloptical variability (Graham et al. 2015a) possibly due tothe fact that the emissions from its mini-disc are Dopplerboosted (D’Orazio et al. 2015; Xin et al. 2019).Infrared (IR) variability can also be caused by aSMBHB heating its surrounding dust torus, see e.g. An-tonucci 1993; Krolik & Begelman 1988. As such, can- didate binaries may further benefit from follow-up mea-surements by the James Webb Space Telescope,
JWST .In addition to periodicity, spectral signatures acrossoptical, UV and X-ray bands are widely used to dif-ferentiate binary SMBHs from normal AGNs poweredby single SMBHs. Merging binary SMBHs likely havedifferent X-ray spectral profiles from isolated SMBHs,including harder X-ray spectra (Roedig et al. 2014; Far-ris et al. 2014a; Ryan & MacFadyen 2017). Of particu-lar relevance is Saade et al. 2020, who measured X-rayspectra of 7
CRTS candidates within the
Chandra
X-rayenergy range and computed their optical/UV-to-X-rayspectral indices. While their spectra showed no differ-ence from the broader AGN population with a purportedsingle SMBH, they are careful to note this is not entirelyunexpected: thermal X-ray profiles can only distinguishbinaries separated by 100 r g , where r g = GM/c is thegravitational radius (Roedig et al. 2014).Interestingly, Saade et al. 2020 identified HS 1630+2355(also known as FBQS J163302.6+234928) as the onlyAGN in their sample that could host a SMBHB sinceits semi-major axis is a ∼ r g , making it a tantalizingSMBHB candidate, which we highlight moving forward. FORECASTING PTA CAPABILITIESIn addition to cataloging current candidate SMBHBsystems, we estimate when these binaries may bedetected. Using H19’s open source Python package hasasia , we assess the current sensitivity of NANOGravand forecast the sensitivity of future PTA experiments.For resolvable individual signals from SMBHBs this in-volves using a matched filter statistic and building aneffective strain-noise power spectral density (interpreted
Xin, Mingarelli, Hazboun as the sensitivity) using the sky location, detector re-sponse functions, and noise parameters of the pulsars.The detection thresholds are then calculated using theexpectation value of the S/N for a circular binary giventhe sensitivity. There are many subtleties involved withcalculating a PTA’s sensitivity to various GW sources,and the interested reader in encouraged to refer to H19for details necessarily left out here.There are a number of statistics in the literature(Babak & Sesana 2012; Ellis et al. 2012; Taylor et al.2014) developed for single-source (CGW) searches inPTA data sets. Here we use the S/N from Hazbounet al. (2019b), as it was developed with calculations ofgeneric sky sensitivities in mind. In our forecasts wehave used the full RMS errors quoted by current PTAdata release papers, but have avoided explicitly inject-ing any time-correlated (red) noise into these pulsars,since it is unclear how much of the red noise currentlyobserved in PTA pulsars is due to the stochastic back-ground of GWs and how much is due to intrinsic spinnoise. It is also likely that in the SKA era, if red noisemodels are not sufficient for PTA use, then it will befeasible to avoid the small number of millisecond pul-sars (Arzoumanian et al. 2018a) with large red noiseamplitudes.The NANOGrav 11-year detection curve is based onthe noise parameters in Arzoumanian et al. (2018a).The IPTA2025 curve uses the pulsars and noise param-eters from Perera et al. (2019) as a starting place, andbuilds the array by adding 4 pulsars per year with RMSvalues of 300 ±
100 ns. In order to add 4 IPTA pulsarsper year to those in Perera et al. (2019) we use the skypositions of the current pulsars to build empirical dis-tributions for drawing new pulsar positions. The RMSvalues for the pulsars are drawn from a truncated nor-mal distribution with mean and standard deviation of300 ns ±
100 ns, truncated at 10 ns and 600 ns. The ca-dences are pulled from an empirical distribution basedon those in Arzoumanian et al. (2018a); Perera et al.(2019), where any new pulsar with an RMS less than300 ns is observed weekly.The SKA sensitivity estimates are based on projectedIPTA data with a large addition of pulsars distributedaccording to the planned SKA MID and LOW surveys inthe first few years of the SKA1 (Keane et al. 2015). Bothconservative and more optimistic estimates for SKA1and 2 are included here. The SKA LOW survey will con-centrate on regions further than ± ◦ from the galacticplane, while the SKA MID survey will concentrate on re-gions within 10 ◦ of the galactic plane (Keane et al. 2015).We assume that 15% of the millisecond pulsars discov-ered will be suitable for PTAs. All pulsars in IPTA2025 are used in the SKA with extended baselines from IPTADR2 (Perera et al. 2019), as it is assumed that any SKAPTA will be based in large part on the extensive datasets from other PTAs.We offer both optimistic and conservative outlooks forconstraining SMBHB candidate masses in Table 1 (de-tection prospects in Table 2). Our conservative sensi-tivity projections use the same empirical distributionof RMS errors as are used to build IPTA2025, whilethe optimistic sensitivities are built using a distributionmore inline with the pulsars’ RMS being “jitter limited”(Shannon et al. 2014; Lam et al. 2016, 2019), with meanand standard deviation of 100 ns ±
30 ns, truncated at9 ns and 201 ns. In general we leave the RMS values ofthe adopted IPTA pulsars unchanged, except for PSRJ1640+2224 which is extremely close in sky position toHS 1630+2355, and the two best timers in this region ofthe sky, PSR J1909-3744 and PSR J1713+0747. Thesethree pulsars are assumed to be observed by the SKAand their RMS values for the time spans considered hereare set to the jitter values from Lam et al. (2019).All PTA detection curves are presented as S/N = 3thresholds on the GW strain . This is the threshold setby the PTA community for detection of single sourcesin the nanohertz band. These forecasts represent onlya limited number of the many permutations one couldpostulate, hence the Jupyter Notebooks used to makethese projections are available on GitHub. SPECIAL CONSIDERATIONS FORMULTIMESSENGER SIGNALS
When a periodic light curve or another EM tracer fora SMBHB system is identified, there are some concretesteps to take to increase a PTA’s sensitivity to the can-didate GW source.It is our good fortune that the top two potentialSMBHB host galaxies in Table 1 lie either in, orvery close to, both the European PTA (EPTA) andNANOGrav’s most sensitive region of the sky (Babaket al. 2016; Aggarwal et al. 2019). This already improvestheir detection prospects. Moreover, if a pulsar is closeto the sky location of a potential SMBHB system, onealso gains a factor of a few in sensitivity from the pulsardetector response function. We can take concrete stepsto increase the S/N ( ρ ) of a potential GW candidate: ρ ∝ (cid:10) N T c/σ (cid:11) / , where for simplicity the N pulsarshave the same intrinsic properties, T is the length ofthe dataset, c is the cadence of the observation, and σ There is a factor of 2 difference in the definition of h betweenH19 and this manuscript which has been taken into account forall of our calculations. MA PTA Constraints on SMBHBs . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − . × − Figure 2. Minimum strain h for an S/N = 3 detection for current and future PTAs . CRTS
GW candidates aremarked in orange: HS 1630+2355 (orange star), SDSS J164452.71+4307 (orange pentagon) and SDSS J114857.33+1600 (orangetriangle) are shown over the detection threshold of S/N = 3 for the NANOGrav 11-year data set, and the hasasia -modeledIPTA2025, SKA1 and SKA2 PTAs at 11 . is the white noise RMS, e.g. Arzoumanian et al. 2018a.We can therefore increase the S/N by increasing thenumber of pulsars, the pulsar observing cadence c , andspending more time on observing these in an effort todecrease the RMS white noise value σ (Burt et al. 2011;Mingarelli et al. 2017; Lam 2018). Moreover combiningpulsar datasets under the auspices of the IPTA is an ex-cellent way to readily increase a pulsar’s total observingtime T and cadence c .Furthermore, Arzoumanian et al. 2020 show that tar-geted GW searches increase NANOGrav’s GW sensitiv-ity by at least a factor of two – more if there is GWfrequency information.Some things are left to chance: both NANOGrav andEPTA show uneven CGW sensitivity on the sky due tothe largely anisotropic distribution of pulsars, and theiruneven timing properties. If we are fortunate, the sourcewill lie in an area of high sensitivity. Moreover, theantenna beam pattern (or detector response function) of a pulsar in direction ˆ p to a GW is ∝ / (1 + ˆΩ · ˆ p ), whereˆΩ is the direction of GW propagation. This responsefunction has been well-studied, and is clearly maximalfor CGW sources in direction − ˆΩ, since the denominatorbecomes small as − ˆΩ → ˆ p . At exactly ˆΩ · ˆ p = − UNLIKELY BINARIES IN
CRTS
Mingarelli et al. (2017) computed the probability ofa galaxy hosting a SMBHB system emitting nanohertzGWs: P merger = t c T z (cid:90) . dµ ∗ dNdz ( M ∗ , z, µ ∗ ) T z (3)where t c = 5 / πf ) − / M − / c , dN/dz is the cumula-tive galaxy-galaxy merger rate (Rodriguez-Gomez et al. Xin, Mingarelli, Hazboun
Table 1.
Top 13 periodic SMBHB candidates ranked by strain upper limit (UL; h σ ), where the mass can be constrained bySKA2 (11 in CRTS plus OJ287 and 3C66B). For detection claims we use h , Table 2. A further 15 candidates have marginalS/N ∼
3, but are not shown here. The last four columns report the S/N ( ρ ) on the 1 − σ strain value for current and future PTAexperiments. In the last two columns the S/N values are listed as optimistic(conservative) . Importantly, Arzoumanian et al.(2020) find that strain ULs improve by a factor of at least two in targeted GW searches, which we do not take into accounthere. *For 3C66B h σ is the 95% UL from Arzoumanian et al. (2020).Object Name RA Dec f GW log( M ) strain NG ρ , IPTA ρ , SKA1 ρ , SKA2 ρ ,[Hz] [ M (cid:12) ] h σ µ ∗ is the stellar mass ratio of the parent galaxies,and T z is the estimated binary lifetime.Briefly, this probability is the product of the probabil-ity of two factors: the existence of a pair of SMBHs inthe galaxy resulting from a galaxy merger (the integrandin Equation 3), and the probability that this binary isemitting nanohertz GWs, t c /T z . Here the binaries over-come the final parsec problem (Begelman et al. 1980) bystellar hardening (Quinlan 1996).We use the public software Nanohertz GWs (Min-garelli 2017) to compute this probability of all the AGNin our sample. Since the AGN light curve catalog weassembled is not complete, we can only provide a rela-tive ranking of the probability of a given AGN hostinga SMBHB.Systems at higher redshift may contain more gas,hence the SMBHBs may overcome the final parsec prob-lem by gas interactions in addition to stellar hardening,e.g. Sesana (2013); Tiede et al. (2020). The gas andstellar 3-body interactions may also lead to large bi-nary eccentricities, which is a subject of future studyand could also help to solve the final parsec problem. Afollow-up paper will include more such solutions to thefinal parsec problem. RESULTSOf the 111
CRTS sources, mass estimates are availablefor 98. We therefore compute the strain for 136 SMBHBcandidates identified in
CRTS , PTF , Pan-STARRS1 surveys, alongside individual quasars OJ287 and 3C66B, and compare the strains with current and future PTAdetection curves. None of the
PTF and
Pan-STARRS1 candidates are detectable by SKA2 – we therefore focusour attention on
CRTS .A multimessenger signal could be an important testof GR by comparing the change in the observed EMperiod and the change in the GW frequency, akin to theHulse-Taylor binary (Hulse & Taylor 1975). As such, wecompute the GW frequency derivative, ˙ f gw , Equation 2,for all the SMBHB candidates. We find that most of the CRTS candidates would have a GW frequency shift of ∼ − nHz/yr, and would therefore not be detectable. Wereport the list of computed ˙ f values in Table 3. However,we find that for 3C66B ˙ f gw = 0 .
14 nHz/yr. Since Sudouet al. (2003)’s claim that 3C66B is a binary with anorbital period of 1 . ± .
03 years, the source shouldhave evolved from 60.4 nHz to 62.8 nHz over the last17 years if it were indeed a SMBHB system. Today, itwould have a period of P = 2 / . MA PTA Constraints on SMBHBs h σ which is the 1 − σ upper limit valueof the strain, Table 1, and h , which is the maximum aposteriori value of the strain, Table 2. When we referto constraints on the mass we refer to h σ – this is thepoint where we can start to constrain the mass upperlimits – and for detection claims we refer to h . Whilethe mass ratio q will be uncertain, these upper boundson strain are dominated by total mass uncertainty.In Figure 1 we show the strain for the CRTS candi-dates, and for 37 of these we can compute the mass fromSDSS spectra in order to compute the strain error bars,Table 3 (spectra were not immediately available for theremaining candidates). If the upper error bar touchesthe S/N=3 detection curve with h σ , we claim that wecan constrain the binary’s total mass. We also showhow important the new SKA pulsars will be for PTAs– even with short timing baselines, these new pulsarslead to significant improvements in detection prospectsat mid-to-high GW frequencies.The sky and polarization averaged detection curvedoes not however give the complete picture, since CGWsensitivity is also a function of sky position and over-all alignment with well-timed millisecond pulsars. Wetherefore use hasasia to generate future PTA detectionsky maps, Figure 2. We find that the CGW sensitiv-ity is largely smoothed out over the sky by going from asingle PTA (here NANOGrav) to the IPTA, as expectedbut never concretely shown. We can also see a residualeffect the IPTA data have on SKA1’s sky map, inducinga slight preference for GW observations in the northernhemisphere, while this is virtually washed out by SKA2.Importantly, NANOGrav’s most sensitive area im-proves by a factor of 6 in the strain, or a factor of 216in volume, between the 11-yr data and our projectionsfor IPTA2025. Moving from IPTA2025 to SKA1 yieldsa further factor of 4 improvement for the strain in 5years, and going from SKA2 to from SKA1 yields an-other factor of 2 in strain sensitivity. Improvements be-tween SKA1 and SKA2 are modest since many of thepulsars in the SKA MID and LOW surveys will likelybe found by 2034 and beyond, according to Keane et al.(2015).In Table 1 we rank the selected CRTS sources fromhighest 1 − σ strain upper limit to lowest. Using h σ we find that 3C66B, SDSS J164452.71+4307 and HS1630+2355 should have S/N > . If not detected, we can startto constrain their masses, and eventually rule out theseAGN as SMBHB host galaxies. SKA1 will constrain 12,and SKA2 will constrain 28 SMBHB candidates’ masses.We also find that IPTA2025 can detect 2 sources (3 if weinclude the factor of 2 from Arzoumanian et al. 2020),SKA1 can detect 12, and SKA2 can detect 17 sources.Our findings are summarized in Table 1 and Table 2,and Figure 2 and Figure 3.We also note that SDSS J114857.33+1600’s detectionprospects improve by a factor of 5 from IPTA2025 toSKA1 in five years – this is due to the large increase oflow-RMS pulsars added by preliminary SKA MID andLOW surveys (Keane et al. 2015) and timed for only upto four years. This large increase in sensitivity over ashort time is also illustrated and explained by Figure 1.Furthermore we independently compute the probabil-ity of each galaxy in CRTS hosting a SMBHB systememitting nanohertz GWs, as in Mingarelli et al. (2017).We find that 24 high- z AGN should not be able to hostSMBHB systems at all, Table 4, since there has notbeen enough time for such binary systems to form underthe assumption that they undergo a dynamical frictionphase followed by a stellar hardening phase. The de-tection of GWs from one of these unlikely host galaxieswould imply significant gas interactions and/or binaryeccentricity are accelerating these high z mergers – aplausible scenario and a subject of future study – or thatthese periodic light curves are not truly tracing binaryactivity.HS 1630+2335 and Mrk 504 show periodicity in boththe optical and X-ray bands ( Figure 4) making themvery interesting SMBHB candidates. Quite serendip-itously HS 1630+2335 is closely aligned with PSRJ1640+2224, which is already a target in NANOGrav’shigh-cadence timing campaign. PSR J1640+2224’sresiduals have decreased by a factor of three since thefirst NANOGrav data release, making this pulsar a keytool to detect GWs from HS 1630+2335, or to rule itout as a true SMBHB. According to our calculations,the S/N based on h σ of HS 1630+2335 in the NG 11-yrdata is ∼ .
9. As a simplifying assumption, we assumeall the signal comes from PSR J1640+2224 (a detectionwill require multiple pulsars). Moving from the 11-yr to12.5-yr data will quadruple the cadence, and the RMSof the residuals is half of the 11-yr value, the S/N wouldimprove by a factor of ∼
3. Moreover, a dedicated searchfor this source in the NANOGrav 12.5-yr data wouldgive an additional factor of 2 in sensitivity, with the Note that the strain calculations for SDSS J164452.71+4307do not have mass errors, therefore h ≈ h σ Xin, Mingarelli, Hazboun
Table 2.
Top 13 periodic SMBHB candidates, ranked by the mean strain ( h ), which can be detected or constrained by SKA2(11 in CRTS plus OJ287 and 3C66B). A further 4 candidates have marginal S/N ∼
3, but are not shown here. The last fourcolumns report the S/N ( ρ ) on h – the maximum a posteriori strain value – for current and future PTA experiments. The lasttwo columns list the optimistic(conservative) S/N values. The S/N calculations here do not include the additional factor of twoone achieves from a targeted search (Arzoumanian et al. 2020). *For 3C66B h is calculated in Arzoumanian et al. (2020).Object Name RA Dec f GW log( M ) strain NG ρ , IPTA ρ , SKA1 ρ , SKA2 ρ ,[Hz] [ M (cid:12) ] h . × − . × − . × − . × − . × − Gravitational wave strain, ( h ) Figure 3. Minimum strain for S/N = 3 detection for SKA1.
The candidates HS 1630+2355 (orange star), SDSSJ164452.71+4307 (orange pentagon) and SDSS J114857.33+1600 (orange triangle) are again highlighted over the detectionthreshold for SKA1 at 4 different frequency bins with pulsars as white circles. The orange hexagons show the sources who’s 1 σ masses will be limited in SKA1. The light purple diamonds show the sky position of the remaining CRTS sources from Table 1. albeit optimistic potential of making a detection at h σ with S/N ∼ CRTS . It has a
MA PTA Constraints on SMBHBs Table 3.
List of 36 SMBHB candidates with BH totalmass error estimates obtained with the width of broadlinespectral emissions (Shen et al. 2008), and the rate ofchange in their GW frequencies, in Hz per year.Object Name log( M / M (cid:12) ) ˙ f GW [Hz/yr]HS 1630+2355 9.74 ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± ± Note —A mass uncertainty of ∼ ± +0 . ) or 50% smaller(10 − . ) than the mean BH mass. FBQS J163302.6+234928(HS1630+2335)
Mrk 504MJD
Figure 4.
The optical periodic light curves of theHS1630+2335 (or FBQS J163302.6+234928 in Chandra X-ray) and Mrk 504. Their periodicity (blue) are fitted withminimized sum of squared residuals ≈ ∼ f gw = 16 . h < . × − . Using the source’s distance of160 Mpc, we can now constrain the chirp mass to be M c < . × M (cid:12) . Assuming an equal mass ratio, q = 1, we can limit to the total mass of the binary to be M < . × M (cid:12) .Importantly, we find that the existence of a GW back-ground does not impede the detection of individualSMBHB at higher frequencies probed by CRTS , PTF ,and others. Since the background is predominantly avery low-frequency GW signal, it only decreases our sen-sitivity by as factor of ∼ . ∼
11 nHz, the lowestfrequency considered here. DISCUSSIONUpcoming wide-field time-domain surveys such as theRubin Observatory Legacy Survey of Space and Time0
Xin, Mingarelli, Hazboun
Table 4.
Unlikely SMBHB candidates When computing the probability of each
CRTS candidate hosting aSMBHB system, Equation 3, we found that 24 candidates were not viable binaries, though this may be dueto the limitations of our model. For them to be viable, their parent galaxies would have had to start theirmergers at z > M ) strain NG ρ , IPTA ρ , SKA1 ρ , SKA2 ρ ,[ M (cid:12) ] h (LSST; Ivezi´c et al. 2019), and ambitious spectroscopicsurveys such as the Sloan Digital Sky Survey-V (SDSS-V; Kollmeier et al. 2017) and the Dark Energy Spectro-scopic Instrument (DESI; Flaugher & Bebek 2014) willreveal a sky rich in potential SMBHBs systems (Kelleyet al. 2019). It will therefore be important to understandcommon false-positive signals of binary activity, whichwe may achieve by studying the CRTS sample lookingfor commonalities in the ruled out SMBHB hosts.In general concrete detection strategies can be initi-ated to make a GW detection, or to rule out a givenAGN as a SMBHB host galaxy. To optimize our abil-ity to detect these GWs, we propose the following de-tection strategy for the
CRTS candidates which can bebroadly applied to other such periodic AGN: (i) startor maintain a high-cadence observing program for pul- sars closely aligned with the top candidates in Table 1(or any future candidate in general). If no such pul-sars exist, ones can be added via pulsar searches (Roshiet al. 2019) or targeted radio follow-ups of
Fermi -LATunassociated sources (so far ∼
20% of these have beenadded to NANOGrav ), (ii) an increase in the amount oftime spent observing these pulsars to lower their RMSwhite noise values, (iii) combining pulsar datasets toachieve an immediate increase in cadence, (iv) use tar-geted CGW searches to improve strain limits by a factorof at least two (Arzoumanian et al. 2020).These strategies have not been applied to candidatesin Table 1 and Table 2, and could significantly accelerate https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars MA PTA Constraints on SMBHBs
Fermi is almost immediately valuable for CGWdetection, as we showed in Figure 1 for the SKA.Some particularly interesting candidates to watch are3C66B, which could show GW frequency evolution,HS1630 where the mass can be constrained with up-coming IPTA observations, and Mrk 504 where we wereable to compute a mass upper limit for the first time, M = 3 . × M (cid:12) . In addition to the latter beingrelatively nearby at 160 Mpc, this rather good massconstraint could also be due to Mrk 504’s proximity topulsar J1713+0747 – an exceptional pulsar.Large errors on the SMBH mass estimates from Shenet al. (2008) translate into large strain errors, since h ∝ M / , e.g. Figure 1. This provides PTAs with a uniqueopportunity to constrain the chirp mass M c of thecandidate SMBHBs, which in turn provides new GW-based constraints on the underlying EM-based SMBHmass estimates, e.g. Shen et al. (2008). However, sur-veys like CRTS may select quasars exhibiting red noisewhich only appears to be periodic over a short times-pan (Vaughan et al. 2016). The detection strategies weoutline, however, will allow us to identify true binarySMBHs by confirming and ruling out signals that aresimply red noise, making this toolkit is complementaryto previous methods, i.e. Vaughan et al. 2016 and Kelleyet al. 2019. Zhu & Thrane (2020)’s new Bayesian ap-proach to measuring periodicity is also a promising toolto mitigate this risk and identify truly periodic signals.Lower limits on q can also be determined by the pres-ence of Doppler boosting (D’Orazio et al. 2015), imply-ing a very small mass ratio for the binary. In generalsuch small mass ratios make detecting such potentialbinaries exceedingly difficult. However, in the case ofPG1302-102, with h ∼ × − at ∼
14 nHz (Gra-ham et al. 2015a), detection strategies we describe here,complemented by a targeted search for the additionalfactor of two Arzoumanian et al. (2020), may make itdetectable with SKA2.Some 24 high-mass, high- z AGN in
CRTS are notviable candidates according to our SMBHB evolutionmodels: the host galaxies would have had to merge at z > z = 4 due to limitationsin Illustris (Rodriguez-Gomez et al. 2015). Moreover,the SMBH occupation fraction for the parent galaxieswhich formed the current binaries would have to be 1.0at z >
4, which is also an interesting result. Our modelsdo not, however, include gas and potential binary ec- centricity, both of which would make the binaries evolvemore quickly. A GW detection from one of these 24 hostgalaxies, such as SDSS J164452+4307 (detectable byIPTA in 2025, Table 2), could therefore inform SMBHBevolution models, and may help in understanding theorigins of SMBHs seeds (Volonteri et al. 2008; Tanaka& Haiman 2009) and their occupation fraction.We found that the presence of a stochastic GW back-ground does not impede our detection prospects. Infact, if these AGN really are SMBHB systems, they mayinduce some anisotropy in the GW background (Min-garelli et al. 2013; Mingarelli et al. 2017). Detection ofanisotropy will likely follow that of the isotropic GWbackground in the next few years.So far, no CGWs have been detected at any fre-quency from any compact object. Long-lived CGWsare currently detectable with PTAs, and may eventuallybe complemented by astrometric GW detection (Mooreet al. 2017). Given that SMBHBs in the PTA band willstay in band for tens of megayears, the current periodiclight curves in
CRTS , PTF , and
PAN-STARRS1 offerus a unique opportunity to study SMBHB host galax-ies and their EM emissions. This will be invaluable in-formation for the Laser Interferometer Space Antenna(LISA; Baker et al. 2019), where the angular resolutionis relatively poor, and the signals will only last weeksto months. Understanding the expected EM emissionsfrom SMBHB host galaxies is thus important ground-work to lay to for both PTAs and LISA in this new eraof multimessenger astrophysics. ACKNOWLEDGEMENTSThe authors thank William Garnier for an up-dated SKA1 timeline, and Zolt´an Haiman, MariaCharisi, Michael Lam, Alberto Sesana, and YacineAli-Ha¨ımoud for useful discussions. The Flatiron In-stitute is supported by the Simons Foundation. TheNANOGrav project receives support from National Sci-ence Foundation (NSF) PIRE program award number0968296 and NSF Physics Frontier Center award num-ber 1430284. CX was supported by NASA ADAP grantNNX17AL82G.
Software:
Astropy (Astropy Collaboration et al. 2013,2018), Hasasia (Hazboun et al. 2019a), Healpy (Zonca et al.2019), Matplotlib (Hunter 2007), Numpy (van der Waltet al. 2011), Pandas (McKinney 2010), SciPy (Virtanenet al. 2019), NanohertzGWs (Mingarelli 2017)REFERENCES
Aggarwal, K., Arzoumanian, Z., Baker, P. T., et al. 2019,ApJ, 880, 116 Albareti, F. D., Prieto, C. A., Almeida, A., et al. 2017, TheAstrophysical Journal Supplement Series, 233, 25 Xin, Mingarelli, Hazboun
Antonucci, R. 1993, Annu. Rev. Astron. Astrophys., 31, 473Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al.2018a, ApJS, 235, 37Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2018b,ApJ, 859, 47Arzoumanian, Z., Baker, P. T., Brazier, A., et al. 2020,arXiv:2005.07123Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020,arXiv e-prints, arXiv:2009.04496Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33Astropy Collaboration, Price-Whelan, A. M., Sip˝ocz, B. M.,et al. 2018, AJ, 156, 123Babak, S., & Sesana, A. 2012, PhRvD, 85, 044034Babak, S., Petiteau, A., Sesana, A., et al. 2016, MonthlyNotices of the Royal Astronomical Society, 455, 1665Baker, J., Bellovary, J., Bender, P. L., et al. 2019, arXive-prints, arXiv:1907.06482Baskaran, D., Polnarev, A., Pshirkov, M., & Postnov, K.2008, Phys. Rev. D, 78, 044018B´ecsy, B., & Cornish, N. J. 2020, Classical and QuantumGravity, 37, 135011Begelman, M. C., Blandford, R. D., & Rees, M. J. 1980,Nature, 287, 307Bonetti, M., Sesana, A., Barausse, E., & Haardt, F. 2018,Monthly Notices of the Royal Astronomical Society, 477,2599Boyle, L., & Pen, U.-L. 2012, PhRvD, 86, 124028Burke-Spolaor, S., Taylor, S. R., Charisi, M., et al. 2019,A&A Rv, 27, 5Burt, B. J., Lommen, A. N., & Finn, L. S. 2011, ApJ, 730,17Chamberlin, S. J., & Siemens, X. 2012, PhRvD, 85, 082001Charisi, M., Bartos, I., Haiman, Z., et al. 2016, MonthlyNotices of the Royal Astronomical Society, 463, 2145D’Orazio, D. J., Haiman, Z., Duffell, P., MacFadyen, A., &Farris, B. 2016, MNRAS, 459, 2379D’Orazio, D. J., Haiman, Z., & Schiminovich, D. 2015,Nature, 525, 351Ellis, J. A., Siemens, X., & Creighton, J. D. E. 2012, ApJ,756, 175Farris, B. D., Duffell, P., MacFadyen, A. I., & Haiman, Z.2014a, Monthly Notices of the Royal AstronomicalSociety: Letters, 447, L80Farris, B. D., Duffell, P., Macfadyen, A. I., & Haiman, Z.2014b, ApJ, 783, 12Flaugher, B., & Bebek, C. 2014, in Society ofPhoto-Optical Instrumentation Engineers (SPIE)Conference Series, Vol. 9147, Proc. SPIE, 91470S Goulding, A. D., Greene, J. E., Bezanson, R., et al. 2018,PASJ, 70, S37Graham, M. J., Djorgovski, S. G., Stern, D., et al. 2015a,Nature, 518, 74—. 2015b, Monthly Notices of the Royal AstronomicalSociety, 453, 1562Hazboun, J., Romano, J., & Smith, T. 2019a, The Journalof Open Source Software, 4, 1775Hazboun, J. S., Romano, J. D., & Smith, T. L. 2019b,PhRvD, 100, 104028Hellings, R. W., & Downs, G. S. 1983, ApJL, 265, L39Hulse, R. A., & Taylor, J. H. 1975, ApJL, 195, L51Hunter, J. D. 2007, Computing in Science Engineering, 9,90Ivezi´c, ˇZ., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873,111Kauffmann, G., & Haehnelt, M. 2000, Monthly Notices ofthe Royal Astronomical Society, 311, 576Keane, E., Bhattacharyya, B., Kramer, M., et al. 2015, inAdvancing Astrophysics with the Square KilometreArray (AASKA14), 40Kelley, L., Charisi, M., Burke-Spolaor, S., et al. 2019,BAAS, 51, 490Kollmeier, J. A., Zasowski, G., Rix, H.-W., et al. 2017,arXiv e-prints, arXiv:1711.03234Krolik, J. H., & Begelman, M. C. 1988, ApJ, 329, 702Laine, S., Dey, L., Valtonen, M., et al. 2020, ApJL, 894, L1Lam, M. T. 2018, ApJ, 868, 33Lam, M. T., Cordes, J. M., Chatterjee, S., et al. 2016, ApJ,819, 155Lam, M. T., McLaughlin, M. A., Arzoumanian, Z., et al.2019, ApJ, 872, 193Lehto, H. J., & Valtonen, M. J. 1996, A&A, 460, 207Liu, T., Gezari, S., Gonz´alez, G., & Hynes, R. 2019, ApJ,36, 19McKinney, W. 2010, in Proceedings of the 9th Python inScience Conference, ed. S. van der Walt & J. Millman, 51– 56Milosavljevi´c, M., & Merritt, D. 2001, ApJ, 563, 34Mingarelli, C. 2017, ChiaraMingarelli/nanohertz GWs:First release!, v.v1.0, Zenodo, doi:10.5281/zenodo.838712Mingarelli, C. M., Lazio, T. J. W., Sesana, A., et al. 2017,Nature Astronomy, 1, 886Mingarelli, C. M. F., & Sidery, T. 2014, PhRvD, 90, 062011Mingarelli, C. M. F., Sidery, T., Mandel, I., & Vecchio, A.2013, PhRvD, 88, 062005Moore, C. J., Mihaylov, D. P., Lasenby, A., & Gilmore, G.2017, PhRvL, 119, 261102O’Beirne, L., Cornish, N. J., Vigeland, S. J., & Taylor,S. R. 2019, PhRvD, 99, 124039
MA PTA Constraints on SMBHBs13