Constraints on individual supermassive binary black holes using observations of PSR J1909-3744

TL;DR

This study uses pulsar timing observations to search for gravitational waves from supermassive binary black hole candidates, setting upper limits on their chirp masses and GW strain amplitudes, and exploring methods to improve detection sensitivity.

Contribution

It provides new upper limits on GW emissions from specific SMBBH candidates and analyzes how orbital parameter knowledge can enhance search sensitivity.

Findings

No significant GW signals detected.

Upper limits on chirp masses are established.

Improved upper limits on GW strain amplitude by a factor of 2 to 7.

Abstract

We perform a search for gravitational waves (GWs) from several supermassive binary black hole (SMBBH) candidates (NGC 5548, Mrk 231, OJ 287, PG 1302-102, NGC 4151, Ark 120 and 3C 66B) in long-term timing observations of the pulsar PSR J1909$-$3744 obtained using the Parkes radio telescope. No statistically significant signals were found. We constrain the chirp masses of those SMBBH candidates and find the chirp mass of NGC 5548 and 3C 66B to be less than $2.4 \times 10^9\,\rm M_{\odot}$ and $2.5 \times 10^9\,\rm M_{\odot}$ (with 95% confidence), respectively. Our upper limits remain a factor of 3 to 370 above the likely chirp masses for these candidates as estimated from other approaches. The observations processed here provide upper limits on the GW strain amplitude that improve upon the results from the first Parkes Pulsar Timing Array data release by a factor of 2 to 7. We investigate how information about the orbital parameters can help improve the search sensitivity for individual SMBBH systems. Finally, we show that these limits are insensitive to uncertainties in the Solar System ephemeris model.

Tables

Table 1. Table 1 : Source information and upper limits. Column 1: name of the source; column 2: red shift; column 3: the observed orbital period; column 4: eccentricity; column 5: total mass; column 6: mass ratio; column 7: chirp mass; column 8: upper limit on chirp mass based on PPTA DR1 data set; column 9: upper limit on chirp mass based on PPTA P18 data set; column 10: ratio between l ​ i ​ m ​ i ​ t P ​ 18 𝑙 𝑖 𝑚 𝑖 subscript 𝑡 𝑃 18 limit_{P18} and M c subscript 𝑀 𝑐 M_{c} ; column 11: strain amplitude; column 12: evolution timescale; column 13: references.

Name	$z$	$P$	$e$	$M_{total}$	$q$	$M_c$	$limi{t}_{DR1}$	$limi{t}_{P18}$	ratio	$h$	$\tau$	Ref.
		[yr]		[${10}^8{M}_{\odot }$]		[${10}^8{M}_{\odot }$]	[${10}^9{M}_{\odot }$]	[${10}^9{M}_{\odot }$]		[${10}^{-17}$]	[Myr]
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
NGC 5548	0.0172	14.1	0.13	2.7	1.0	1.2	7.6	2.4	20	6.4	33	1
Mrk 231	0.0422	1.2	-	1.5	0.026	0.16	15	6.0	370	0.45	1.4	2
OJ 287	0.3056	12.1	0.700	184	0.008	10.2	230	88	86	11.8	0.83	3
PG 1302-102	0.2784	5.2	-	3.0	-	1.3	20	11	85	0.88	1.3	4
NGC 4151	0.0033	15.9	0.42	0.56	0.27	0.19	14	2.0	110	1.4	980	5
Ark 120	0.0327	$\sim 20$	-	2.6	-	1.1	43	15	140	2.2	100	6
3C 66B	0.0213	1.05	-	19	0.58	7.9	3.4	2.5	3.2	730	$1.4\times {10}^{-3}$	7

Equations

$$\begin{array}{rl}s(t,\hat{\Omega })=& (F^{+}(\hat{\Omega })\cos 2\psi +F^{\times }(\hat{\Omega })\sin 2\psi )△s_{+}(t)\\ & +(F^{+}(\hat{\Omega })\sin 2\psi -F^{\times }(\hat{\Omega })\cos 2\psi )△s_{\times }(t),\end{array}$$

$$\begin{array}{rl}{F}^{+}=& \frac{1}{4(1-\cos \theta )}\{(1+{\sin }^2\delta ){\cos }^2{\delta }_p\cos [2(\alpha -{\alpha }_p)]\\ & -\sin 2\delta \sin 2{\delta }_p\cos (\alpha -{\alpha }_p)+{\cos }^2\delta (2-3{\cos }^2{\delta }_p)\},\\ F^{\times }=& \frac{1}{2(1-\cos \theta )}\{\cos \delta \sin 2{\delta }_p\sin (\alpha -{\alpha }_p)\\ & -\sin \delta {\cos }^2{\delta }_p\sin [2(\alpha -{\alpha }_p)]\},\end{array}$$

$$△s_{+,\times }(t)=s_{+,\times }(t)-s_{+,\times }(t_p),$$

$$t_p=t-d_p(1-\cos \theta )/c,$$

$$\begin{array}{rl}{s}_{+}(t)=& \sum _n-(1+{\cos }^2\iota )[a_n\cos (2\gamma )-b_n\sin (2\gamma )]\\ & +(1-{\cos }^2\iota )c_n,\\ s_{\times }(t)=& \sum _{n}2\cos \iota [b_n\cos (2\gamma )+a_n\sin (2\gamma )],\end{array}$$

$$\begin{array}{rl}{a}_n=& -\zeta {\omega }^{-1/3}[J_{n-2}(ne)-2e{J}_{n-1}(ne)+(2/n)J_n(ne)\\ & +2e{J}_{n+1}(ne)-J_{n+2}(ne)\sin [nl(t)],\\ b_n=& \zeta {\omega }^{-1/3}\sqrt{1-e^2}[J_{n-2}(ne)-2J_n(ne)\\ & +J_{n+2}(ne)]\cos [nl(t)],\\ c_n=& (2/n)\zeta {\omega }^{-1/3}{J}_n(ne)\sin [nl(t)].\end{array}$$

$$\frac{\mathrm{d}f}{\mathrm{d}t}$$

$$\frac{\mathrm{d}e}{\mathrm{d}t}$$

$$\frac{h_{DR1}}{h_{P18}}\sim \sqrt{\frac{N_{P18}}{N_{DR1}}}=\sqrt{\frac{393}{125}}\sim 2$$

$$\sqrt{(s|s)}\ge 5.$$

$$\sqrt{\frac{1}{2}N\frac{A^2}{{\sigma }^2}}\ge 5,$$

$$A_{th}=\sqrt{\frac{50}{N}}\sigma ,$$

$$f_M(M_{95})=f_h(h_{95})=A_{0.05}=A_{th}$$

$$\begin{array}{rl}s=& ϵ(1+{\cos }^2\iota )(F^{+}\cos 2\psi +F^{\times }\sin 2\psi )\sin (nl(t))\\ & +ϵ(2\cos \iota )(F^{+}\sin 2\psi -F^{\times }\cos 2\psi )\cos (nl(t)),\end{array}$$

$$ϵ=\frac{(GM{)}^{5/3}}{c^4{D}_L(2\pi f{)}^{1/3}}.$$

$$h=\frac{2(GM{)}^{5/3}(2\pi f{)}^{2/3}}{c^4{D}_L},$$

$$h_{95}\propto {ϵ}_{95}\propto M_{95}^{5/3}.$$

$$\begin{array}{rl}{A}^2=\Vert s{\Vert }^2=& {ϵ}^2[(1+{\cos }^2\iota {)}^2(F^{+}\cos 2\psi +F^{\times }\sin 2\psi {)}^2\\ & +(2\cos \iota {)}^2(F^{+}\sin 2\psi -F^{\times }\cos 2\psi {)}^2].\end{array}$$

$$A=2ϵ\sqrt{F^{+2}+F^{\times 2}}=(1+\cos \theta )ϵ.$$

$$h_{95}\propto {ϵ}_{95}\propto A_{th}\propto \frac{1}{\sqrt{N}}$$

Full text

1 Introduction

Pulsar timing arrays (PTAs) can be used to detect and study gravitational waves (GWs) in the nanohertz frequency band (Detweiler 1979; Foster & Backer 1990). The primary source of GWs in this band are thought to be inspiralling supermassive binary black holes (SMBBHs), formed in the aftermath of galaxy mergers (Begelman et al. 1980). The separations of these binary systems are generally at sub-parsec scales, which are usually too small to be resolved with current electromagnetic instruments. This makes identification of SMBBHs and measurements of their physical properties (such as the black hole mass and orbital parameters) particularly challenging. Thus far there have been a number of SMBBH candidates identified through several lines of approaches, although the nature of these candidates are still under investigation (e.g., Sundelius et al. 1997, Bon et al. 2012,Graham et al. 2015a, Yan et al. 2015, Li et al. 2016).

There has recently been significant progress in searches for GWs from individual SMBBHs in PTA data sets (Yardley et al. 2010, Zhu et al. 2014, Babak et al. 2016, Aggarwal et al. 2018). Although such projects have not detected GWs yet, the first detection is likely in the next one or two decades. For instance, Mingarelli et al. (2017) claimed that GWs from at least one nearby SMBBH will be detected within 10 years if the GW background (GWB) can be successfully isolated. Once detected, the SMBBH can also be followed up by electromagnetic observations, thus enabling a multi-messenger view of the black hole system (Burke-Spolaor 2013). Furthermore, the detection of SMBBH systems can yield direct information about the masses and spins of the black holes (Mingarelli et al. 2012), which would shed light on the formation and evolution of supermassive black holes. The precisely measured distance of PSR J0437$-$4715 by Reardon et al. (2016) would allow the polarization angle of GWs generated by OJ 287 to be measured to better than $8^{\circ}$ (Chen & Zhang 2018). Chen and Zhang (2018) also showed that measurement uncertainties of the spin and quadrupole scalar from future PTAs are comparable or even smaller than those in the optical measurements.

Lommen & Backer (2001) using pulsar observations to place upper limits on the mass ratio of hypothetical SMBBHs in the Magorrian et al. (1998) sample of galaxies that were as small as 0.06 for orbital periods of $\sim$ 2000 days. This work was continued and updated by Schutz & Ma (2016), which constrained the mass ratios of hypothetical SMBBHs in nearby galaxies using PTA limits (Zhu et al. 2014; Babak et al. 2016). Sudou et al. (2003) detected apparent periodic variations in the radio core position of the nearby elliptical galaxy 3C 66B with a period of 1.05 yr and suggested the presence of a SMBBH with a total mass of $5.4\times 10^{10}\,\rm M_{\odot}$ and a mass ratio of 0.1. Jenet et al. (2004) subsequently placed limits on the mass and eccentricity of the proposed SMBBH and ruled out the adopted system with 95% confidence.

In this work, we use data from long-term timing observations of the pulsar, PSR J1909$-$3744, from the Parkes Pulsar Timing Array project to place limits on the chirp mass of a few well-known SMBBH candidates using a method similar to that in Jenet et al. (2004). Our chosen SMBBH candidates, listed below, have all been monitored in detail and show evidence for a SMBBH system:

• •

The well-known BL Lac object OJ 287 exhibits a $\sim$12-year periodicity in the 100-year long optical light curve (Sundelius et al. 1997).

• •

NCG 4151 (Bon et al. 2012) exhibits periodic variations in the H$\alpha$ light and radial velocity curves.

• •

NGC 5548 (Li et al. 2016) shows evidence for periodicites in both the long-term continuum light curves and emission line profile variations.

• •

PG 1302$-$102 exhibits periodicity in optical light curve, which is well modeled by the relativistic Doppler boosting of the secondary mini-disk (Graham et al. 2015a)111There are several modern time-domain surveys that have released candidate samples for SMBBHs through periodicity searches (Graham et al. 2015b; Charisi et al. 2016; Liu et al. 2016). Most of these candidates are fairly distant so that their GW emissions are clearly too weak to detect with the current PTA datasets. In addiation, due to the limited temporal baselines (generally $\sim$1.5 cycles of the periods), the reported periodicities may suffer high false alarm probabilities (Vaughan et al. 2016). We thereby only select PG 1302$-$102 reported by Graham et al. (2015a), the most significant candidate found in the Catalina Real-Time Transient Survey..

• •

Ark 120 (Li et al. 2019) exhibits a wave-like pattern in the optical emissions with a total temporal baseline over four decades and the H$\beta$ integrated flux series varies with a similar behaviour.

• •

Mrk 231 was found to have a flux deficit at UV wavelength, which could signify a SMBBH system (Yan et al. 2015).

• •

3C66B: since Jenet et al. (2004) ruled out the original model of Sudou et al. (2003), Iguchi et al. (2010) performed follow-up observations and still found evidence of a SMBBH system, but with significantly lower SMBBH mass estimates than were previously assumed. Such low masses imply that the expected GW emission would be below the detection thresholds of previous all-sky searches using PTA data (see, e.g., Zhu et al. 2014).

Table 1 also lists the basic properties of the above candidates. Some of the SMBBH candidates mentioned above have published estimates of the full orbital parameters for the postulated SMBBH systems (e.g., NGC 5548, 3C 66B and OJ 287). Many (but not all) of the previous PTA studies usually did not use such information when placing upper limits on the chirp mass of such systems. Here, we investigate in detail how the extra information, e.g., the eccentricity of the orbit and the inclination angle of the binary orbit with respect to the line of sight, can help improve the search sensitivity. Recently, it has been shown that the GWB constraints are sensitive to the choice of the Solar System ephemeris (SSE) adopted in the data analysis (Arzoumanian et al. 2018). We explore how the choice of the SSE influences the analysis results presented here.

The organization of this paper is as follows. Section 2 describes the expected signature of GW emission from a general binary system. The observations of PSR J1909$-$3744 are described in Section 3. Section 4 presents the limit on the chirp mass for our SMBBH candidates. We also discuss here how further information about the orbital parameters (if they were available) would influence the limit and how the limit changes with those choice of SSE model used. We conclude in Section 5.

2 THE SIGNATURE OF A SMBBH

We adopt a quadrupole formalism for estimating the signature of GW signals in pulse arrival times (ToAs) from pulsar observations. The pulsar timing residuals induced by a single GW source can be written as

[TABLE]

where $\hat{\Omega}$ is a unit vector defining the direction of GW propagation, $\psi$ is the GW polarization angle. The two functions $F^{+}$ and $F^{\times}$ are the geometric factors which only depend on the GW source position relative to a given pulsar (Lee et al. 2011)

[TABLE]

where $\cos\theta=\cos\delta\cos\delta_{p}\cos(\alpha-\alpha_{p})+\sin\delta\sin\delta_{p}$ with $\theta$ being the angle between the GW source and pulsar direction with respect to the observer, $\delta(\delta_{p})$ and $\alpha(\alpha_{p})$ are the declination and right ascension of the GW source (pulsar), respectively.

The GW-induced pulsar timing residuals can be expressed as the combination of two terms - the Earth term $s_{+,\times}(t)$ and the pulsar term $s_{+,\times}(t_{p})$:

[TABLE]

[TABLE]

where $d_{p}$ is the pulsar distance. $s_{+}(t)$ and $s_{\times}(t)$ are source-dependent functions, and at Newtonian order take the same forms for GWs emitted by SMBBHs in eccentric orbits as (Taylor et al. 2016):

[TABLE]

where

[TABLE]

Here, $\zeta=(GM)^{5/3}/c^{4}D_{L}$ is the amplitude parameter, where $M$ is the binary chirp mass defined as $M^{5/3}=m1m2(m1+m2)^{-1/3}$ with $m1$ and $m2$ being the binary component masses, $D_{L}$ is the luminosity distance of the binary, $e$ is the eccentricity, and $\omega=2\pi f$. $f$ is the orbital frequency; $l(t)$ is the mean anomaly; $\gamma$ is the initial angle of periastron; and $\iota$ is the inclination angle of the binary orbit with respect to the line of sight. In the following, $f$ and $M$ refer to the observed redshifted values, such that $f_{r}=f(1+z)$ and $M_{r}=M/(1+z)$, where $f_{r}$ and $M_{r}$ are rest frame values, and $z$ is the cosmological redshift of the binary. We assume no binary evolution over typical pulsar timing baselines. The Earth term $s_{+,\times}(t)$ and pulsar term $s_{+,\times}(t_{p})$ are computed by assuming that the binary’s mean orbital frequency $f$ and eccentricity $e$ remain constant over the total time span of our observations of a given pulsar. Meanwhile, we assume binary evolution over the lag time between the Earth and pulsar term signals. $e$ and $f$ for the pulsar term are obtained by solving the coupled differential equations (Peters 1964)

[TABLE]

3 TIMING OBSERVATIONS

We use data from long-term timing observations of the pulsar PSR J1909$-$3744, which is one of the best high-time-precision pulsars. We have processed two, related data sets. The first (known as “PPTA DR1”) is the initial PPTA data release (Manchester et al. 2013) and we make use of the data collection for this data set provided by Madison et al. (2016). We also process the more up-to-date data set that was used to constrain ultralight scalar-field dark matter by Porayko et al. (2018) (hereafter we call this the “P18” data set).

The DR1 data set was acquired between March 2005 and March 2011 for 20 pulsars. The P18 data set consists of observations for 26 pulsars collected between February 5, 2004 and January 31, 2016. It includes the DR1 data along with some earlier data for some pulsars. The observing cadence is normally once every two to three weeks. At most observing epochs, each pulsar was observed in three radio bands (10, 20 and 50 cm) with a typical integration time of one hour.

Pulsar times of arrival were fitted with a timing model using the standard TEMPO2 (Hobbs et al. 2006) software package. Typical parameters included in this fit are the five astrometric parameters (sky position, proper motion, and parallax), spin frequency, spin-down rate, dispersion measure, and (when applicable) binary orbital parameters. Additionally, constant offsets or jumps were fitted among ToAs collected with different receiver/backend systems. For the DR1 data set, we used TT(BIPM2011) as the reference time standard and adopted the JPL DE421 solar system ephemeris (SSE) model. For the more up-to-date data set, we used TT(BIPM2015) and adopted the JPL DE418 SSE model.

The pulsars in the data sets are not equally sensitive to GW signatures. For the work that we describe here we simply make use of a single pulsar. This significantly simplifies the analysis procedures and provides a baseline result for more detailed studies. We note that selecting a single pulsar dataset can lead to biased bounds (which we discuss later), but this work is primarily a proof of concept and we are currently very unlikely to detect the GW signature from a SMBBH even with the full sensitivity obtained using all pulsars in an array. Such a detailed analysis will be presented by the PPTA team when the second data release has been completed. We therefore use observations of PSR J1909$-$3744. Similar to Zhu et al. (2014) and Madison et al. (2016), ToAs of the best band (i.e. where the lowest rms timing residuals are seen and 10 cm observing band for PSR J1909$-$3744) have been selected after correcting for dispersion measure variations.

4 RESULTS AND DISCUSSION

The timing residuals from PSR J1909$-$3744 were searched for the signature of GWs, using a Lomb-Scargle periodogram (LP; Lomb 1976; Scargle 1982). There have recently been lots of PTA data analysis methods for single GW sources developed (Yardley et al. 2010, Ellis et al. 2012, Ellis 2013, Zhu et al. 2014, Wang et al. 2014, Zhu et al. 2015, Wang et al. 2015, Zhu et al. 2016, Taylor et al. 2016), ranging from frequentist to Bayesian techniques (e.g., Ellis 2013, Taylor et al. 2016), and from Earth-term-only approaches to a coherent inclusion of both Earth terms and pulsar terms (e.g., Wang et al. 2014, Wang et al. 2015, Zhu et al. 2016). Although Lomb-Scargle periodogram is sub-optimal, this simplifies the analysis procedures and provides a baseline result for more detailed studies. And we defer analysis using optimal methods including ellipticity and the evolution of binary orbit during the period of observation to a future study. The timing residuals and Lomb-Scargle power are shown in the top and bottom panel of Figure 1 for the two data sets respectively. The LP for DR1 data set is shown in the bottom left panel of Figure 1. The periodogram power was determined for 250 frequencies ranging from 1/22.0 $\rm yr^{-1}$ to 11.4 $\rm yr^{-1}$, with a resolution of 1/22.0 $\rm yr^{-1}$. This corresponds to a frequency-oversampling factor of four as the time span of the DR1 data is 5.5 yrs. In order to determine the statistical significance of a given peak, we generate synthetic data sets of Gaussian deviates with the same sampling as in the real data and find the largest periodogram power for each such data set. Using the distribution of largest periodogram power of $10^{5}$ synthetic data sets, the periodogram power of a significant peak is 0.18 assuming a False-Alarm probability (FAP) of 0.001 . Hence there is no significant peak in this LP.

The LP for P18 data set is shown in the bottom right panel of Figure 1. The periodogram power was calculated for 786 frequencies ranging from 1/48.0 $\rm yr^{-1}$ to 16.4 $\rm yr^{-1}$, with a resolution of 1/48.0 $\rm yr^{-1}$. This corresponds to a frequency-oversampling factor of 4 as the time span of P18 data is 12 yrs. Similar to PPTA DR1 data set. The periodogram power of a significant peak is 0.067 according to Monte Carlo simulations. Hence there is no significant peak in this LP.

4.1 CONSTRAINTS FROM PULSAR TIMING

Since our analysis find no significant signal, we choose to set upper limits on the chirp masses of a few well-known SMBBH candidates. A Monte Carlo analysis similar to that in Jenet et al. (2004) was then used to place 95% upper limits on the chirp mass of a few well-known SMBBH candidates. Before laying out our computational procedure, we wish to state that the 95% upper limit on the chirp mass implies that the true value of the chirp mass is less than the upper limit with 95% probability.

The distance to PSR J1909$-$3744 is $1.140\pm 0.012$ kpc (Reardon et al. 2016). So the distribution of $d_{p}$ is taken as $(d_{p}/\rm{kpc})\in\mathcal{N}(1.140,0.012)$. The orbital frequency, $f$, and the SMBBH distance are calculated using values in Table 1. The eccentricity, $e$, uses the value in Table 1 or taken as 0 when it is not available in electromagnetic observations (denoted as ‘-’ in Table 1). Also the mass ratio, $q$, is taken as 1 when it is not available from the electromagnetic observations.

The distribution for the other parameters are as follows: $l_{0}\in U[0,2\pi],\gamma\in U[0,2\pi],\iota\in U[0,\pi],\psi\in U[0,\pi]$. For 3C 66B, we use the reported inclination angle in Iguchi et al. (2010), i.e. $\iota\in U[0.08,0.11]$ . For NGC 5548, the full orbital parameters were reported with small error bars by Li et al. (2016) and a schematic for the geometry of the SMBBH system can be found therein. The distribution for the parameters are as follows: $(f/\rm{Hz})\in U[2.15\times 10^{-9},2.43\times 10^{-9}],e\in U[0,0.31],l_{0}\in U[4.52,5.22],\gamma\in U[3.23,3.84],\iota\in U[0.36,0.48],\psi\in U[0,\pi]$. In top panel of Figure 1, we show 5 realizations of pre-fit timing residuals induced by NGC 5548 with the pulsar term included for each data set. For these 5 realizations, we use the upper limit chirp mass based on PPTA P18 data set, i.e. $2.4\times 10^{9}\,\rm M_{\odot}$, to calculate timing residuals. The other parameters of each line are as follows:

red line: $(f/\rm{Hz})=2.40\times 10^{-9},e=0.108,l_{0}=5.18,\gamma=3.33,\iota=0.421,\psi=0.58,d_{p}/\rm{kpc}=1.1524$ 2. 2.

orange line: $(f/\rm{Hz})=2.36\times 10^{-9},e=0.016,l_{0}=4.79,\gamma=3.34,\iota=0.364,\psi=3.08,d_{p}/\rm{kpc}=1.1384$ 3. 3.

green line: $(f/\rm{Hz})=2.33\times 10^{-9},e=0.231,l_{0}=5.21,\gamma=3.64,\iota=0.422,\psi=2.48,d_{p}/\rm{kpc}=1.1453$ 4. 4.

blue line: $(f/\rm{Hz})=2.38\times 10^{-9},e=0.0243,l_{0}=4.83,\gamma=3.52,\iota=0.455,\psi=0.89,d_{p}/\rm{kpc}=1.1588$ 5. 5.

magenta line: $(f/\rm{Hz})=2.36\times 10^{-9},e=0.218,l_{0}=4.96,\gamma=3.76,\iota=0.476,\psi=1.46,d_{p}/\rm{kpc}=1.1405$

Timing residuals induced by other SMBBH candidates can be found in Appendix B.

The corresponding waveform was generated using equations in Section 2. The waveform was then added to the ToA data. Next we computed the residuals using the TEMPO2 (Hobbs et al. 2006) package. The residuals were then analyzed using the LP method described above. If a significant peak was found (see above), then the signal was considered to be detected. A total of 1000 waveforms were tested for each $M$. We adjust $M$ until $95\%$ of the injected signal can be detected.

In Table 1, we tabulate the upper limits on the chirp mass. For example, 95$\%$ confidence level upper limit on the chirp mass of NGC 5548 is $7.6\times 10^{9}\,\rm M_{\odot}$ and $2.4\times 10^{9}\,\rm M_{\odot}$ for PPTA DR1 and P18 data set respectively. And 95$\%$ confidence level upper limit on the chirp mass of 3C 66B is $2.5\times 10^{9}\,\rm M_{\odot}$ for PPTA P18 data set. Our upper limits remain a factor of 3 to 370 above the reported values for these candidates.

4.2 COMPARISON BETWEEN DR1 AND P18

The DR1 data set for PSR J1909$-$3744 has a time span of 5.5 yrs with 125 data points and P18 data set has a time span of 12.0 yrs with 393 data points. Note that the time span of DR1 data set is even shorter than the period of the gravitational wave expected of NGC 5548, NGC 4151, OJ 287 and Ark 120, which reduce their detectabilities in such a data set significantly.

The corresponding upper limit on the chirp mass for the P18 data set is typically a factor of 2 to 3 more stringent than that of DR1 data set. The corresponding upper limit on the strain amplitude for the P18 data set is typically a factor of 2 to 7 more stringent than that of DR1 data set. This is caused by two effects: (1) the P18 data set has more data points than DR1 data set, and (2) the fitting process is known to significantly reduce a PTA’s sensitivity to GWs at the lowest frequencies (because of the fit for pulsar spin period and its first time derivative).

We estimate the two effects respectively. According to the scaling relation 19 (see Appendix A for derivation),

[TABLE]

So the number of data points contributes about a factor of 2 and the fit for pulsar spin period and its first derivative contributes a factor of about 1 to 3. Specifically, the fitting process contributes a factor of about 3 for low frequency GW source NGC 5548 , OJ 287 and Ark 120. For NGC 4151, the DR1 upper limit is much worse, this is because low frequency combined with high eccentricity can reduce the sensitivity significantly. The fitting process can have a huge effect on GW detection of low frequency sources. A factor of 3 worse caused by the fitting process would require 9 times longer observation to compensate. Longer time span of PTA data is critical for detection of low frequency and eccentric SMBBHs.

4.3 INFLUENCE OF EXTRA INFORMATION

In this section, we discuss how extra information about the orbital parameters of the SMBBH would influence the upper limits on the strain amplitude. As we study the relative influence of the orbital parameters, so we just consider the Earth term of the signal. We take NGC 5548 as an example and run the following cases for both DR1 and P18 data sets:

fiducial case: We know everything about the GW source. $(f/\rm{Hz})=2.15\times 10^{-9},e=0,l_{0}=0,\gamma=0,\iota=0,\psi=0$ 2. 2.

phase effect: How varying phase influence the limit. $(f/\rm{Hz})=2.15\times 10^{-9},e=0,l_{0}\in U[0,2\pi],\gamma=0,\iota=0,\psi=0$ 3. 3.

polarization angle effect: How varying polarization angle influence the limit. $(f/\rm{Hz})=2.15\times 10^{-9},e=0,l_{0}=0,\gamma=0,\iota=0,\psi\in U[0,\pi]$ 4. 4.

inclination effect: especially, how would edge on inclination influence the limit. $(f/\rm{Hz})=2.15\times 10^{-9},e=0,l_{0}=0,\gamma=0,\iota=0.5\pi,\psi\in U[0,\pi]$ 5. 5.

eccentricity effect: how would large eccentricity influence the limit. $(f/\rm{Hz})=2.15\times 10^{-9},e=0.7,l_{0}=0,\gamma=0,\iota=0,\psi\in U[0,\pi]$

The results are shown in Figure 2. Comparing case 1 and case 2, different phases do not influence upper limits. For case 3, polarization angle effect is equivalent to phase effect when $\iota=0$, so the result of case 3 is the same as case 2. Comparing case 3 and 4, different inclinations can cause upper limits on the strain amplitude to vary by a factor about 10. We get a worse upper limit when $\cos\iota\sim 0$. Comparing case 3 and 5, increasing the eccentricity leads to a worse upper limit. This is because the LS power is less when the eccentricity gets larger.

4.4 INFLUENCE OF THE SOLAR SYSTEM EPHEMERIS

SSE errors can mimic a GWB signal (Tiburzi et al. 2016). Recently, it was shown that the GWB constraints are sensitive to the SSE model used (Arzoumanian et al. 2018). The GWB constraints are very different if they use different SSE models. If they take the SSE errors into account, the constraints would converge. However, this removes most evidence for the presence of a GWB.

To investigate how SSE model used would influence the upper limits on the strain amplitude, we repeated the fiducial case in Section 4.3 using each of the JPL DE418, JPL DE421, JPL DE430 and JPL DE436 SSEs for both DR1 and P18 data sets. The results are shown in Figure 3. It indicates that the choice of SSE does not influence the upper limits on the chirp mass of NGC 5548. This is consistent with our expectation due to the intrinsic difference between the two detection experiment. The GWB ‘signal’ is a kind of noise in the PTA data sets, thus its detection is more prone to inaccurate or incomplete characterization of various sources of noises, including that of SSE.

5 CONCLUSION

We have performed a search for GWs from several SMBBH candidates in long-term timing observations of the pulsar PSR J1909$-$3744 as part of the PPTA. No statistically significant signals were found. We constrained the chirp masses for these SMBBH candidates. For example, the 95$\%$ confidence level upper limit on the chirp mass of NGC 5548 is $2.4\times 10^{9}\,\rm M_{\odot}$. The corresponding limit for 3C 66B is $2.5\times 10^{9}\,\rm M_{\odot}$. Our upper limits remain a factor of 3 $\sim 370$ above the reported values for these candidates. The corresponding upper limit on the strain amplitude for PPTA P18 data set is about a factor of 2 to 7 more stringent than that of DR1 data set. The number of data points contributes about a factor of 2 and the fit for pulsar spin period and its first derivative contributes about a factor of 1 to 3. Longer time span of PTA data is critical for detection of low frequency and eccentric SMBBHs. We investigated how information about the orbital parameters would influence the upper limit. This may suggest strategies for future PTA single GW source searches. If a PTA’s goal is to detect a known single GW source, we should consider the full orbital parameters to decide which source as our detection target.

Although it was shown that the GWB constraints are sensitive to the SSE model used (Arzoumanian et al. 2018) very recently, we find that our upper limits on the chirp mass of NGC 5548 is not sensitive to the SSE model used. But whether this result only holds for LP method or single pulsar analysis is not known. We do not attempt to answer this question in this paper but leave it to a future study.

We note that there are ongoing large time-domain surveys that have yielded hundreds of quasar and AGN candidates with possible periodic variability. The list of best candidates is growing. Future work will place constraints on other known nearby and new candidate SMBBH systems.

In this work, we just used one pulsar, PSR J1909$-$3744, which has the second best timing precision in the PPTA data set. PSR J0437$-$4715 has the best timing precision, but the timing residuals of PSR J0437$-$4715 has a strong red noise component and are not suitable for the LP method. Shannon et al. (2015) found no evidence for red noise for PSR J1909$-$3744, so PSR J1909$-$3744 is suitable for the LP method. Current PTA’s sensitivity to GWs is often dominated by several best-timing-precision pulsars. For example, Babak et al. (2016) pointed out that the array is heavily dominated by PSR J1909-3744, contributing more than 60 percent of $(\rm S/N)^{2}$, followed by PSR J1713+0747 at about 20 percent. Moreover, Shannon et al. (2015) found that their limit placed using PSR J1909$-$3744 is the same as the joint limit, indicating that PSR J1909$-$3744 is the dominant pulsar in the sample. Although the bounding techniques for single GW and GWB are different, PSR J1909$-$3744 contributes about 30 % of the total S/N of the full PPTA data set and the upper chirp mass limits given by the one pulsar analysis are meaningful even compared with a full PPTA analysis. However, if the opening angle between the pulsar and the SMBBH is large, it is not sensitive to the SMBBH. In this case, the constraints can be bad, and the full PPTA data set with a better sky coverage can be used to place better constraints. Nevertheless, the one pulsar analysis provides valuable insights into problems such as the influence of the SSE, etc. We defer analysis using the full PPTA data set to future work. We will also discuss how adding pulsars to present PTAs would influence the constraints.

Acknowledgements.

We thank the anonymous referee for very useful comments on the manuscript. This work is supported by National Key R$\&$D Program of China No. 2017YFA0402600, by the CAS International Partnership Program No.114- A11KYSB20160008, by the CAS Strategic Priority Research Program No. XDB23000000, and by NSFC grant No. 11725313 and 11690024. Y.R.L. acknowledges financial support from the NSFC (11573026) and the Youth Innovation Promotion Association CAS. J.M.W acknowledges financial support from the NSFC (11833008), the National Key R&D Program of China (2016YFA0400701), and the Key Research Program of Frontier Sciences of the Chinese Academy of Sciences (QYZDJ-SSWSLH007). The Parkes radio telescope is part of the Australia Telescope, which is funded by the Commonwealth of Australia for operation as a National Facility managed by the Commonwealth Scientific and Industrial Research Organisation (CSIRO). This work is supported by Aliyun Cloud. The authors wish to acknowledge crucial guidance from George Hobbs and many useful discussions with Xingjiang Zhu, Zhoujian Cao and Youjun Lu.

Appendix A scaling relation

$\rm SNR\geq 5$ is usually required for a definitive detection, i.e.

[TABLE]

This formular can be approximated as

[TABLE]

where $N$ is the total number of data points, $A$ is the amplitude of the timing residuals induced by the SMBBH, $\sigma$ is the noise rms. The thresold of the amplitude to be detected is

[TABLE]

i.e. if $A$ is larger than $A_{th}$, then we can detect the signal. We can detect the signal at 95 confidence level when $A_{0.05}=A_{th}$, where $A_{0.05}$ is 0.05 fractile of $A$ . 95% upper limit of M($M_{95}$) and h($h_{95}$) is determined by $A_{0.05}$ through the following formular

[TABLE]

So the functions $f_{M}$ and $f_{h}$ are essential to determine $M_{95}$ and $h_{95}$. Assuming the fiducial parameters: $e=0,l_{0}=0,\gamma_{0}=0$, the Earth term pulsar timing residuals can be written as

[TABLE]

where

[TABLE]

The strain amplitude is given by

[TABLE]

so

[TABLE]

The amplitude $A$ is calculated by

[TABLE]

Further assuming $\iota=0$, the amplitude can be simplified as:

[TABLE]

Combining the relations above, we have

[TABLE]

Appendix B Gw-induced timing residuals

The parameters of each line are as follows:

red line: $(f/\rm{Hz})=2.42\times 10^{-9},e=0.061,l_{0}=4.66,\gamma=3.49,\iota=0.361,\psi=1.64,d_{p}/\rm{kpc}=1.1322$ 2. 2.

orange line: $(f/\rm{Hz})=2.25\times 10^{-9},e=0.035,l_{0}=4.97,\gamma=3.48,\iota=0.389,\psi=1.06,d_{p}/\rm{kpc}=1.1330$ 3. 3.

green line: $(f/\rm{Hz})=2.28\times 10^{-9},e=0.023,l_{0}=5.10,\gamma=3.78,\iota=0.450,\psi=0.76,d_{p}/\rm{kpc}=1.1382$ 4. 4.

blue line: $(f/\rm{Hz})=2.30\times 10^{-9},e=0.102,l_{0}=4.80,\gamma=3.45,\iota=0.447,\psi=3.89,d_{p}/\rm{kpc}=1.1562$ 5. 5.

magenta line: $(f/\rm{Hz})=2.41\times 10^{-9},e=0.216,l_{0}=5.14,\gamma=3.55,\iota=0.431,\psi=2.72,d_{p}/\rm{kpc}=1.1481$

The parameters of each line are as follows:

red line: $(f/\rm{Hz})=2.53\times 10^{-8},e=0,l_{0}=3.39,\gamma=5.98,\iota=1.22,\psi=0.60,d_{p}/\rm{kpc}=1.1348$ 2. 2.

orange line: $(f/\rm{Hz})=2.53\times 10^{-8},e=0,l_{0}=4.33,\gamma=0.67,\iota=1.41,\psi=0.93,d_{p}/\rm{kpc}=1.1469$ 3. 3.

green line: $(f/\rm{Hz})=2.53\times 10^{-8},e=0,l_{0}=0.34,\gamma=3.87,\iota=0.78,\psi=1.80,d_{p}/\rm{kpc}=1.1460$ 4. 4.

blue line: $(f/\rm{Hz})=2.53\times 10^{-8},e=0,l_{0}=5.02,\gamma=3.94,\iota=2.73,\psi=0.06,d_{p}/\rm{kpc}=1.1332$ 5. 5.

magenta line: $(f/\rm{Hz})=2.53\times 10^{-8},e=0,l_{0}=1.52,\gamma=3.48,\iota=0.11,\psi=2.00,d_{p}/\rm{kpc}=1.1208$

The parameters of each line are as follows:

red line: $(f/\rm{Hz})=2.62\times 10^{-9},e=0.7,l_{0}=2.85,\gamma=0.11,\iota=0.29,\psi=2.33,d_{p}/\rm{kpc}=1.1378$ 2. 2.

orange line: $(f/\rm{Hz})=2.62\times 10^{-9},e=0.7,l_{0}=0.27,\gamma=4.60,\iota=1.54,\psi=1.68,d_{p}/\rm{kpc}=1.1526$ 3. 3.

green line: $(f/\rm{Hz})=2.62\times 10^{-9},e=0.7,l_{0}=1.67,\gamma=4.57,\iota=1.79,\psi=2.96,d_{p}/\rm{kpc}=1.1408$ 4. 4.

blue line: $(f/\rm{Hz})=2.62\times 10^{-9},e=0.7,l_{0}=5.62,\gamma=0.38,\iota=0.70,\psi=0.73,d_{p}/\rm{kpc}=1.1522$ 5. 5.

magenta line: $(f/\rm{Hz})=2.62\times 10^{-9},e=0.7,l_{0}=3.10,\gamma=3.41,\iota=2.79,\psi=1.91,d_{p}/\rm{kpc}=1.1493$

The parameters of each line are as follows:

red line: $(f/\rm{Hz})=6.15\times 10^{-9},e=0,l_{0}=2.93,\gamma=3.07,\iota=0.14,\psi=0.81,d_{p}/\rm{kpc}=1.1332$ 2. 2.

orange line: $(f/\rm{Hz})=6.15\times 10^{-9},e=0,l_{0}=4.92,\gamma=2.62,\iota=2.83,\psi=2.35,d_{p}/\rm{kpc}=1.1477$ 3. 3.

green line: $(f/\rm{Hz})=6.15\times 10^{-9},e=0,l_{0}=1.11,\gamma=5.12,\iota=2.00,\psi=2.91,d_{p}/\rm{kpc}=1.1448$ 4. 4.

blue line: $(f/\rm{Hz})=6.15\times 10^{-9},e=0,l_{0}=6.14,\gamma=2.19,\iota=1.89,\psi=2.59,d_{p}/\rm{kpc}=1.13926$ 5. 5.

magenta line: $(f/\rm{Hz})=6.15\times 10^{-9},e=0,l_{0}=5.83,\gamma=6.09,\iota=0.21,\psi=1.13,d_{p}/\rm{kpc}=1.1401$

The parameters of each line are as follows:

red line: $(f/\rm{Hz})=2.01\times 10^{-9},e=0.42,l_{0}=2.39,\gamma=0.55,\iota=0.63,\psi=2.68,d_{p}/\rm{kpc}=1.1243$ 2. 2.

orange line: $(f/\rm{Hz})=2.01\times 10^{-9},e=0.42,l_{0}=2.54,\gamma=3.57,\iota=0.66,\psi=0.09,d_{p}/\rm{kpc}=1.1252$ 3. 3.

green line: $(f/\rm{Hz})=2.01\times 10^{-9},e=0.42,l_{0}=3.87,\gamma=5.82,\iota=2.85,\psi=2.72,d_{p}/\rm{kpc}=1.1406$ 4. 4.

blue line: $(f/\rm{Hz})=2.01\times 10^{-9},e=0.42,l_{0}=1.96,\gamma=2.36,\iota=2.58,\psi=2.19,d_{p}/\rm{kpc}=1.1263$ 5. 5.

magenta line: $(f/\rm{Hz})=2.01\times 10^{-9},e=0.42,l_{0}=2.62,\gamma=5.63,\iota=0.54,\psi=1.29,d_{p}/\rm{kpc}=1.1497$

The parameters of each line are as follows:

red line: $(f/\rm{Hz})=1.53\times 10^{-9},e=0,l_{0}=5.75,\gamma=1.78,\iota=3.08,\psi=0.26,d_{p}/\rm{kpc}=1.1425$ 2. 2.

orange line: $(f/\rm{Hz})=1.53\times 10^{-9},e=0,l_{0}=5.06,\gamma=3.44,\iota=0.10,\psi=1.34,d_{p}/\rm{kpc}=1.1259$ 3. 3.

green line: $(f/\rm{Hz})=1.53\times 10^{-9},e=0,l_{0}=3.71,\gamma=5.47,\iota=0.60,\psi=1.10,d_{p}/\rm{kpc}=1.1415$ 4. 4.

blue line: $(f/\rm{Hz})=1.53\times 10^{-9},e=0,l_{0}=1.58,\gamma=1.44,\iota=0.29,\psi=0.04,d_{p}/\rm{kpc}=1.1476$ 5. 5.

magenta line: $(f/\rm{Hz})=1.53\times 10^{-9},e=0,l_{0}=2.86,\gamma=5.37,\iota=1.34,\psi=1.80,d_{p}/\rm{kpc}=1.1432$

The parameters of each line are as follows:

red line: $(f/\rm{Hz})=3.02\times 10^{-8},e=0,l_{0}=5.43,\gamma=2.03,\iota=0.11,\psi=0.01,d_{p}/\rm{kpc}=1.1276$ 2. 2.

orange line: $(f/\rm{Hz})=3.02\times 10^{-8},e=0,l_{0}=1.06,\gamma=3.95,\iota=0.09,\psi=1.50,d_{p}/\rm{kpc}=1.1428$ 3. 3.

green line: $(f/\rm{Hz})=3.02\times 10^{-8},e=0,l_{0}=3.41,\gamma=6.00,\iota=0.11,\psi=0.67,d_{p}/\rm{kpc}=1.1432$ 4. 4.

blue line: $(f/\rm{Hz})=3.02\times 10^{-8},e=0,l_{0}=2.05,\gamma=1.08,\iota=0.10,\psi=2.13,d_{p}/\rm{kpc}=1.1660$ 5. 5.

magenta line: $(f/\rm{Hz})=3.02\times 10^{-8},e=0,l_{0}=0.21,\gamma=3.87,\iota=0.08,\psi=1.05,d_{p}/\rm{kpc}=1.1316$