A Multimessenger Search for the Supermassive Black Hole Binary in 3C 66B with the Parkes Pulsar Timing Array

TL;DR

This study searches for gravitational waves from a suspected supermassive black hole binary in galaxy 3C 66B using pulsar timing data, setting upper limits and proposing a combined electromagnetic and gravitational-wave analysis approach.

Contribution

It introduces a new methodology for joint likelihood analysis of electromagnetic and gravitational-wave data for SMBHBs, and assesses the candidate in 3C 66B with updated limits.

Findings

No definitive detection of the SMBHB signal.

Upper limits on chirp mass and strain amplitude are established.

Partially rules out some electromagnetic-based parameter space.

Abstract

A subparsec supermassive black hole binary (SMBHB) at the center of the galaxy 3C 66B is a promising candidate for continuous gravitational-wave searches with pulsar timing arrays (PTAs). In this work, we search for such a signal in the third data release of the Parkes Pulsar Timing Array. Matching our priors to estimates of binary parameters from electromagnetic observations, we find a log Bayes factor $\ln B = - 0.0027(7)$, highlighting that the source can be neither confirmed nor ruled out. We place upper limits at $95\%$ credibility on the chirp mass $M < 6.90 \times 10^{8}\ M_{\odot}$, and on the characteristic strain amplitude $\textrm{log}_{10}(h_0)< -14.44$. This partially rules out the parameter space suggested by electromagnetic (EM) observations of 3C 66B. We also independently reproduce the calculation of the chirp mass with the 3 mm flux monitor data from the unresolved core of 3C 66B. Based on this, we outline a new methodology for constructing a joint likelihood of EM and gravitational-wave data from SMBHBs. Finally, we suggest that targeted searches may allow firmly established SMBHB candidates to be treated as standard sirens, for complementary constraints on the Universe expansion rate.

Tables

Table 1. Table 1: Summary of the Properties of the Supermassive Black Hole Binary Candidate in 3C 66B.

	Value	Reference
$ℳ$ ($M_{\odot }$)	${7.9}_{-4.5}^{+3.8}\times {10}^8$	Iguchi et al. (2010)
$f_{\text{GW}}$ (nHz)	$60.4\pm 1.73$	Sudou et al. (2003)
$D_{\text{L}}$ (Mpc)	$94.17\pm 0.20$	van den Bosch et al. (2015)
R.A. (HH:MM: SS.SSS)	02:23:11.4110	Hunt et al. (2021)
Decl. (DD:MM: SS.SSS)	+42:59:31.385	Hunt et al. (2021)
$h_0$	${7.3}_{-5.8}^{+6.8}\times {10}^{-15}$	Iguchi et al. (2010)

Table 2. Table 2: Table Describing the Pulsar Distances, the Methods Used to Obtain Them, and the Source from Which They Were Obtained.

PSR	Distance [pc]	Measurement	Reference
	(pc)
J0030+0451	${329}_{-5}^{+6}$	VLBI astrometry	Ding et al. (2023)
J0125-2327	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J0437-4715	$156.96\pm 0.11$	Shklovskii effect	Shklovskii (1970); Reardon et al. (2024)
J0613-0200	$125\pm 54.69$	Timing	Verbiest et al. (2009)
J0614-3329	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J0711-6830	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J0900-3144	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J1017-7156	$256.41\pm 78.90$	Timing	Ng et al. (2014)
J1022+1001	${720.98}_{-14.55}^{+21.31}$	VLBI astrometry	Deller et al. (2019)
J1024-0719	${1072}_{-49}^{+67}$	Optical astrometry	Moran et al. (2023)
J1045-4509	$303.03\pm 174.47$	Timing	Verbiest et al. (2009)
J1125-6014	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J1446-4701	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J1545-4550	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J1600-3053	$5000\pm 3750$	Timing	Verbiest et al. (2009)
J1603-7202	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J1643-1224	${4763.36}_{-104.89}^{+99.06}$	VLBI astrometry	Ding et al. (2023)
J1713+0747	${1052.63}_{-55.40}^{+66.48}$	VLBI astrometry	Chatterjee et al. (2009)
J1730-2304	${510}_{-30}^{+30}$	VLBI astrometry	Ding et al. (2023)
J1744-1134	$416.67\pm 17.36$	Timing	Verbiest et al. (2009)
J1824-2452A	$10000\pm 50000$	VLBI astrometry	Ding et al. (2023)
J1832-0836	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J1857+0943	$909.09\pm 165.29$	VLBI astrometry	Verbiest et al. (2009)
J1902-5105	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J1909-3744	$1265.82\pm 32.05$	Timing	Verbiest et al. (2009)
J1933-6211	$1000\pm 200$	…	Ellis and Van Haasteren (2017)
J1939+2134	${2604.17}_{-311.96}^{+325.52}$	VLBI astrometry	Ding et al. (2023)
J2124-3358	$322.58\pm 57.23$	Timing	Verbiest et al. (2009)
J2129-5721	$2358.49\pm 489.50$	Optical astrometry	Jennings et al. (2018)
J2145-0750	${623.83}_{-3.50}^{+24.51}$	VLBI astrometry	Deller et al. (2019)
J2241-5236	$1000\pm 200$	…	Ellis and Van Haasteren (2017)

Equations

$$\{\theta ,\varphi ,f_{\text{GW}},{\Phi }_0,\psi ,\iota ,\mathcal{M},(D_{\text{L}}\wedge h_0)\}$$

$$h_0=\frac{2(G\mathcal{M}{)}^{5/3}(\pi f_{\text{GW}}{)}^{2/3}}{D_{\text{L}}{c}^4},$$

$$f_{\text{GW}}^{-8/3}(t)=\frac{(8\pi {)}^{8/3}}{5}(\frac{G\mathcal{M}}{c^3}{)}^{5/3}(t-t_{\text{c}}),$$

$$\mathcal{L}(\mathbf{\delta t}|\mathbf{\theta })=\frac{\exp (-\frac{1}{2}(\mathbf{\delta t}-\mathbf{\mu }{)}^{\text{T}}{\mathbf{C}}^{-1}(\mathbf{\delta t}-\mathbf{\mu }))}{\sqrt{\det \left(2\pi \mathbf{C}\right)}},$$

$$\begin{array}{r}\mathbf{S}(t)=A_{\text{EM}}\sin (2\pi f_{\text{EM}}t+{\varphi }_{\text{EM}})+\\ +a_{\text{EM}}(t-t_0{)}^3+b_{\text{EM}}(t-t_0{)}^2+\\ +c_{\text{EM}}(t-t_0)+d_{\text{EM}}\end{array}$$

$$H_0=\frac{\nu }{D}.$$

$$\hat{\Omega }=(\varphi ,\theta )=-\text{sin}\theta \text{cos}\varphi \hat{x}-\text{sin}\theta \text{sin}\varphi \hat{y}-\text{cos}\theta \hat{z},$$

$$h_{ab}(t,\hat{\Omega })=e_{ab}^{+}(\hat{\Omega })h_{+}(t,\hat{\Omega })+e_{ab}^{\times }(\hat{\Omega })h_{\times }(t,\hat{\Omega }).$$

$$\hat{n}=$$

$$\hat{p}=$$

$$\cos \psi \cos \theta \sin \varphi +\sin \psi \cos \varphi ,-\cos \psi \sin \theta ),$$

$$\hat{q}=$$

$$\sin \psi \cos \theta \sin \varphi -\cos \psi \cos \varphi ,-\sin \psi \sin \theta ).$$

$$e_{ab}^{+}={\hat{p}}_a{\hat{p}}_b-{\hat{q}}_a{\hat{q}}_b,$$

$$e_{ab}^{\times }={\hat{p}}_a{\hat{q}}_b+{\hat{q}}_a{\hat{p}}_b.$$

$$F^A(\hat{\Omega })\equiv \frac{1}{2}\frac{{\hat{u}}^a{\hat{u}}^b}{1+\hat{\Omega }\cdot \hat{u}}{e}_{ab}^A(\hat{\Omega }),$$

$$\delta t(t,\hat{\Omega })=F^{+}(\hat{\Omega })\delta t_{+}^{\prime }(t)+F^{\times }(\hat{\Omega })\delta t_{\times }^{\prime }(t),$$

$$\delta t_{+,\times }^{\prime }(t)=\delta t_{+,\times }(t_{\text{p}})-\delta t_{+,\times }(t),$$

$$t_p=t-\frac{L}{c}(1+\hat{\Omega }\cdot \hat{u}),$$

$$\delta t_{+}^{\prime }(t)=$$

$$\delta t_{\times }^{\prime }(t)=$$

$$\mathcal{M}\equiv \frac{(m_1{m}_2{)}^{3/5}}{(m_1+m_2{)}^{1/5}}$$

$$\begin{array}[]{l}\Phi(t)=\Phi_{0}+\frac{1}{32}(\frac{G\mathcal{M}}{c^{3}})^{-5/3}[\omega_{0}^{-5/3}-\omega(t)^{-5/3}],\\ \omega(t)=\omega_{0}(1-\frac{256}{5}(\frac{G\mathcal{M}}{c^{3}})^{5/3}\omega_{0}^{8/3}t)^{-3/8},\end{array}$$

$${\mathcal{M}}_{\text{o}}$$

$$f_o$$

$$\{\begin{array}{l}\beta =\sqrt{(1-r{)}^2+({\beta }_{\text{app}}r{)}^2},\\ \iota =\text{atan2}({\beta }_{\text{app}}r,1-r),\end{array}$$

$$r=\frac{P_{\text{obs}}}{(1+z)P_{\text{k}}}.$$

$${\beta }_{\text{bl,0}}=\frac{1-\beta \cos \iota }{\sin \iota }\frac{S_{\text{r}}^{\frac{1}{3+\alpha }}-1}{S_{\text{r}}^{\frac{1}{3+\alpha }}+1}.$$

$$\{\begin{array}{l}a=\sqrt[3]{\frac{G{M}_{\text{tot}}{P}_{\text{k}}^2}{(2\pi {)}^2}},\\ R=\frac{v_{\text{bl,0}}{P}_{\text{k}}}{2\pi }.\end{array}$$

$$\{\begin{array}{l}\mathcal{M}=M_{\text{tot}}[x(1-x){]}^{3/5},\\ x\equiv \frac{R}{a}=v_{\text{bl,0}}\sqrt[3]{\frac{P_{\text{k}}}{2\pi G{M}_{\text{tot}}}}.\end{array}$$

Full text

gravitational waves — pulsars: general — methods: data analysis

††software: ptmcmcsampler (Ellis and van Haasteren, 2019), enterprise (Ellis et al., 2019), enterprise_warp at github.com/bvgoncharov/enterprise_warp, enterprise_extensions (Taylor et al., 2021).

I Introduction

Pulsar timing arrays (PTAs) are experiments for long-term monitoring of millisecond pulsar pulse arrival times with the primary goal of detecting nanohertz-frequency gravitational waves (GWs). Recent PTA data show evidence for a nanohertz-frequency gravitational wave background (GWB; Goncharov et al., 2021b; Arzoumanian et al., 2020a; Chen et al., 2021; Antoniadis et al., 2022; Agazie et al., 2023a; EPTA Collaboration et al., 2023; Reardon et al., 2023). The background origin is consistent with the superposition of adiabatically inspiraling supermassive black hole binaries  (SMBHBs; Goncharov et al., 2025). Provided confident detection of the background, the detection of a continuous gravitational wave (CGW) from an individual SMBHB is the next milestone for PTAs. Such an observation would decisively resolve the final parsec problem, which represents difficulties in modeling the evolution of SMBHBs to subparsec separations (Milosavljević and Merritt, 2003). When the host galaxy of CGW source is known, the detection unlocks multimessenger astronomy and cosmology with PTAs (Schutz, 1986; Holz and Hughes, 2005; Charisi et al., 2022).

One of the most promising sources of CGWs is the SMBHB candidate in the center of galaxy 3C 66B. By using long baseline interferometry, Sudou et al. (2003) found evidence for elliptical binary orbital motion in this radio galaxy’s core and inferred a GW frequency of $f_{\textrm{GW}}=60.4\pm 1.73$ nHz. Further observations by Iguchi et al. (2010) found that the binary has a chirp mass $\mathcal{M}=7.9^{+3.8}_{-4.5}\times 10^{8}\ M_{\odot}$, primarily through measuring the orbital radius of the electromagnetic (EM)-emitting black hole and the semi-major axis of an SMBHB. Galaxy 3C 66B is located at a redshift $z=0.02126$ (van den Bosch et al., 2015), which corresponds to a luminosity distance $D_{\text{L}}=94.17^{+0.20}_{-0.20}$ Mpc (Planck Collaboration et al., 2020).

Based on the above, the expected strain amplitude $h_{0}$ of a CGW from 3C 66B is $7.3^{+6.8}_{-5.8}\times 10^{-15}$. This is the highest known strain amplitude of CGW sources in the PTA frequency band, and is within reach for PTA experiments. A summary of the expected properties of the SMBHB in 3C 66B is shown in Table 1. Because the properties in the table are inferred from EM observations, we will also refer to this model of the SMBHB in 3C 66B as the “EM model”.

It is worth noting that Iguchi et al. (2010) revised the chirp mass of the SMBHB in 3C 66B inferred in the earlier work of Sudou et al. (2003). The initial model of the SMBHB at the center of 3C 66B proposed by Sudou et al. (2003) suggested a very high chirp mass of $1.3\times 10^{10}\ M_{\odot}$. Jenet et al. (2004) have ruled out this model based on 7 yr of timing from PSR B1855+09 (Kaspi et al., 1994). Assuming a circular SMBHB in 3C 66B, they placed an upper limit of $7\times 10^{9}\ M_{\odot}$ on its chirp mass. This upper limit does not exclude the model by Iguchi et al. (2010), which motivates further studies.

Following the revision of the chirp mass of the SMBHB candidate by Iguchi et al. (2010), a number of searches for a CGW from 3C 66B with PTAs have been performed. The authors of Arzoumanian et al. (2020b) performed a search for an SMBHB in a circular orbit with the 11 yr dataset of the North American Nanohertz Observatory for Gravitational Waves (NANOGrav). An upper limit on the chirp mass of $1.65\times 10^{9}\ M_{\odot}$ was placed. A subsequent search was then performed for both circular and eccentric SMBHBs using the NANOGrav 12.5 yr dataset. The search for a circular binary (Arzoumanian et al., 2023) constrained the mass to $1.41\times 10^{9}\ M_{\odot}$. The search, which included eccentricity in the SMBHB system (Agazie et al., 2024), placed their tightest constraints on the upper limit of the chirp mass with the value of $1.81\times 10^{9}\ M_{\odot}$. The results from the 12.5 yr NANOGrav dataset are all still compatible with the Iguchi et al. (2010) model, yet claim no detection, which motivates further searches for the source in PTA datasets. Similar to the results shown later in our work, the NANOGrav 15 yr dataset places a limit on the chirp mass at $95\%$ credibility as $9.54\times 10^{8}~M_{\odot}$ (Agarwal et al., 2025). On this same dataset, Tian et al. (2025) found that for an eccentric model, the upper limit on the chirp mass is of $1.82\times 10^{9}~M_{\odot}$.

In this work, we perform a search for a CGW from an SMBHB with a circular orbit in 3C 66B using the third data release (DR3) of the Parkes Pulsar Timing Array  (PPTA; Zic et al., 2023). It is the first targeted search for 3C 66B with PPTA data. We expect the search based on these data to provide better constraints on a CGW from 3C 66B, thanks to an increased time span exceeding 18 yr in 13 pulsars (Mingarelli et al., 2017). As motivated by Liu and Vigeland (2021), the targeting is performed by setting the sky position as a fixed parameter, as well as setting priors for other parameters based on the EM observations summarized in Table 1. We employ both Bayesian and frequentist methods to calculate the significance of the hypothesis that an SMBHB with subparsec separation resides in 3C 66B. Finally, we discuss the current feasibility of cosmology with PTAs and produce constraints on the Hubble–Lemaître constant, assuming that the EM model is valid.

The rest of the paper is organized as follows. In Section II, we introduce the dataset from DR3 of the PPTA. In Section III, we introduce the signal model, the Bayesian analysis methodology, and our noise models. In Section IV, we present our results and discuss the model selection, upper limits on the chirp mass, upper limits on the characteristic strain amplitude, upper limits as a function of CGW frequency, and the frequentist analysis results. In Section V, we discuss new methods for a joint-likelihood analysis of 3C 66B and whether targeted searches in PTAs can be used for cosmology experiments. In Section VI, we present our conclusions.

II Data

The PPTA DR3 dataset (Zic et al., 2023) is based on monitoring of $N_{\text{psr}}=32$ millisecond radio pulsars. The total observing span of PPTA DR3 is $T_{\text{obs}}=18$ yr. The first observations took place on 2004 February 6, and the last observations took place on 2022 March 8. The observations were conducted using the 64 m Parkes ‘Murriyang’ radio telescope in Australia, and observations were taken for each pulsar approximately every 3 weeks. The most recent 3 yr of observations in DR3 are carried out using the ultra-wide bandwidth, low-frequency receiver (UWL; Hobbs et al., 2020). In this work, only the data from 31 pulsars are used. PSR J1741+1351 is excluded from the analysis due to the fact that it was only added to the PPTA observations after the commissioning of the UWL and was only observed with low priority. With only 16 observations, this pulsar is not sensitive to any CGW signal and will therefore not contribute to the results.

As with any PTA data, PPTA DR3 is a set of pulse arrival time observations, including observation details such as radio frequency and backend-receiver information, as well as the pulsar ephemeris, known as the timing model. A timing model predicts pulse arrival times based on the pulsar’s specific spin, astrometry, properties of the interstellar plasma along the line of sight, and the binary parameters of the pulsar when applicable. We use the DE440 ephemeris (Park et al., 2021) to obtain pulse arrival times in an inertial frame centered at the solar system barycenter (SSB).

Additionally, we follow Iguchi et al. (2010) by employing the same millimeter-wavelength flux variation data for 3C 66B. There are two components of this data. The first component of the data is collected at $93.7$ GHz using the Nobeyama Millimeter Array of the National Astronomical Observatory of Japan. The second component of the data is collected at $86.2$ GHz using the Plateau de Bure Interferometer of the Institut de Radioastronomie Millimétrique.

III Methodology

III.1 Signal Model

We use the CGW signal model for nonspinning SMBHBs in circular111Iguchi et al. (2010) point out that 3C 66B may be in a circular orbit. orbits at zeroth post-Newtonian order described in Aggarwal et al. (2019) and Arzoumanian et al. (2020b, 2023). More details of the model are provided in Appendix A. The CGW model for metric perturbations at SSB (the Earth term) is based on eight parameters:

[TABLE]

where $(\theta,\phi)$ is the sky position, $f_{\text{GW}}$ is the GW frequency, $\Phi_{0}$ is the orbital phase, $\psi$ is the GW polarization angle, $\iota$ is the binary inclination angle, $\mathcal{M}$ is the chirp mass, $D_{\text{L}}$ is the luminosity distance and $h_{0}$ is the gravitational strain amplitude. At zeroth post-Newtonian order, we compute the amplitude as

[TABLE]

unless it is provided as a free parameter (hence, $D_{\text{L}}\wedge h_{0}$). The full CGW model is based on metric perturbations at both the SSB and each pulsar (pulsar terms). This requires additional $2N_{\text{psr}}$ parameters $\{L_{i},\Phi_{i}\}$, $i=(1,N_{\text{psr}})$. They are the pulsar distance and the orbital phase at the pulsar, respectively. We provide more information on the pulsar distances in Appendix C. Searching only for the Earth term is sufficient for detection and upper limits, whereas including the pulsar term enhances the search sensitivity and remains important for parameter estimation when the signal is detected.

We model the evolution of CGW frequency on the timescales of the light travel time between the pulsars and Earth. However, for simplicity, we neglect frequency evolution on the timescale of PPTA DR3. This approach is referred to as the stationary phase approximation. To demonstrate the validity of this approach, we employ an equation for the evolution of CGW frequency as a function of time at zeroth post-Newtonian order from (Abbott et al., 2017a):

[TABLE]

where $f_{\text{GW}}$ is the GW frequency, $G$ is Newton’s constant, $c$ is the speed of light, and $t_{\text{c}}$ is the time of coalescence. We present the expected frequency evolution of 3C 66B at SSB in Figure 1. Based on the EM model for 3C 66B, a circular SMBHB is only expected to evolve by about $0.04$ $\text{nHz}\ \text{yr}^{-1}$. In the 20 yr since the initial frequency was determined by Sudou et al. (2003), the frequency may have increased by $0.86$ nHz. This width does not exceed the frequency resolution of PPTA DR3 of $T^{-1}_{\text{obs}}=1.75$ nHz. Therefore, according to the model provided by Sudou et al. (2003), the stationary phase approximation is sufficient for our analysis.

III.2 Bayesian Inference

We perform parameter estimation and model selection by computing posteriors, $\mathcal{P}(\bm{\theta}|\bm{\delta t})=\mathcal{Z}^{-1}\mathcal{L}(\bm{\delta t}|\bm{\theta})\pi(\bm{\theta})$, where $\mathcal{L}(\bm{\delta t}|\bm{\theta})$ is the likelihood of data time series vector $\bm{\delta t}$ given a vector of model parameters $\theta$ for signals and noise, and $\pi(\bm{\theta})$ is the prior. The constant $\mathcal{Z}$ is the integral of the numerator over $\bm{\theta}$. It is also referred to as the evidence or the fully marginalized likelihood. We neglect $\mathcal{Z}$ in our parameter estimation and obtain posterior samples using parallel-tempered Markov-Chain Monte-Carlo (Ellis and van Haasteren, 2019). The model selection follows from computing the Bayes factor, $\mathcal{B}=\mathcal{Z}/\mathcal{Z}_{\varnothing}$, where $\mathcal{Z}_{\varnothing}$ is the evidence of the noise model and $\mathcal{Z}$ is the evidence of the signal model, which includes both the noise and the CGW. Instead of computing $\mathcal{Z}$ separately for the two models, we compute the Bayes factor directly using the product-space sampling method (Hee et al., 2016).

We use the multivariate Gaussian likelihood (van Haasteren et al., 2009)

[TABLE]

where vector $\bm{\mu}$, a function of $\bm{\theta}$, represents a prediction of deterministic contributions to pulse arrival times $\bm{\delta t}$, which includes a CGW model, as well as deterministic noise transients such as exponential dips in the pulse arrival times (Goncharov et al., 2021a). The form of our covariance matrix $\bm{C}$ and the methods for computing it are described in detail in Arzoumanian et al. (2016). In short, $\bm{C}=\bm{N}+\bm{T}\bm{B}\bm{T}^{\text{T}}$, where the block diagonal matrix $\bm{N}$ represents temporally uncorrelated “white” noise, $\bm{T}$ is the design matrix for reduced-rank modeling (Lentati et al., 2013), and $\bm{B}$ is the prior covariance matrix of reduced-rank coefficients $\bm{b}$. The likelihood is marginalized over $\bm{b}$, which determines the time series realization $\bm{T}\bm{b}$. Thus, the design matrix $\bm{T}$ is a mapping from the parameter space $\bm{b}$ to the time series. It is made up of three blocks. Block $\bm{M}$ maps timing model parameters to the time series. Block $\bm{F}$ maps Fourier sine and cosine amplitudes of time-correlated “red” noise to the time series. Block $\bm{U}$ maps delays in the arrival times at different observing epochs due to pulse jitter to all arrival times. Coefficients $\bm{b}$ are assumed to be samples from the zero-mean Gaussian distribution described by the covariance matrix $\bm{B}$. For the coefficients of the pulsar timing model, we assume an uninformative prior with a standard deviation of $10^{40}$ in respective units, which is much greater than the typical coefficient value.

III.3 Noise Model

Noise modeling in PTAs is crucial for correct astrophysical interpretability of the GW searches (Goncharov and Sardana, 2025; Hazboun et al., 2020; Larsen et al., 2024; Reardon et al., 2023). Broadly, modeled PTA noise is categorized into temporally uncorrelated (white) noise, temporally correlated (red) noise with the power spectral density modeled as a power law, as well as noise transients. We use the noise model from Reardon et al. (2023). In addition to this, we model the contribution of the common-spectrum noise process associated with the GWB (Goncharov et al., 2022). It is shown in Zhao et al. (2025) that not modeling this process leads to false positives. Agazie et al. (2023b) suggest that further modeling of the Hellings and Downs (1983) correlations may be necessary for the data where they are evident. We do not model these correlations in this work; however, it may be required for future CGW searches with the PPTA data. Furthermore, we impose a minimum allowed value of the red noise power spectral density of $T_{\text{obs}}\times 10^{-20}~\text{s}^{3}$, which is well below the noise in PPTA DR3. This is done to avoid numerical errors in the likelihood.

III.4 Analysis of the 3 mm Data

We perform a Bayesian reanalysis of the 3C 66B 3 mm flux measurements rather than using the original frequentist results from Iguchi et al. (2010). This allows a direct comparison between the GW posteriors and the EM posterior. We employ the univariate Gaussian likelihood, $\mathcal{L}(\bm{S}|\bm{\theta}_{\text{EM}})$, where $\bm{S}$ is a vector of flux measurements in units of mJy. The general form is the same as Equation 4, except that the covariance matrix is diagonal with elements $\bm{\sigma}_{\text{EM}}^{2}=\bm{\sigma}_{S}^{2}e_{\text{f}}^{2}+e_{\text{q}}^{2}$. Here, $\bm{\sigma}_{S}$ is the flux measurement uncertainty. White noise parameters $(e_{\text{f}},e_{\text{q}})$ are the error factor and the error added in quadrature, respectively. The analog of $\bm{\mu}$ from Equation 4 for the EM data is the sine cubic function from Iguchi et al. (2010):

[TABLE]

In total, $\bm{\theta}_{\text{EM}}$ is the set of parameters ($e_{\text{f}}$, $e_{\text{q}}$, $A_{\text{EM}}$ [mJy], $f_{\text{EM}}$ [Hz], $\phi_{\text{EM}}$, $t_{0}$ [s], $a_{\text{EM}}$ [mJy$\text{s}^{-3}$], $b_{\text{EM}}$ [mJy$\text{s}^{-2}$], $c_{\text{EM}}$ [mJy$\text{s}^{-1}$], $d_{\text{EM}}$ [mJy]).

The chirp mass is calculated from $\bm{\theta}_{\text{EM}}$ and Equation 5 as follows. Iguchi et al. (2010) suggest that the sine component is due to the Doppler boosting of the flux from the core of 3C 66B due to binary motion around the center of mass, whereas the cubic component may be responsible for the nonthermal radiation in the knots of the jet in the core. From the frequency of the sine component $f_{\text{EM}}$, the SMBHB orbital frequency $f_{\text{k}}$ from Sudou et al. (2003), and the apparent boost in the relativistic jet $\beta_{\text{app}}$ (Iguchi et al., 2010), we compute the jet outflow boost $\beta$ and the inclination angle $\iota$ between the jet axis and the line of sight. The flux ratio from Doppler boosting (Rieger and Mannheim, 2000; De Paolis et al., 2002, 2004) is $\bm{S}_{\text{r}}=\max(\bm{S})/\min(\bm{S})$, where $\max(\bm{S})=\bar{S}+A_{\text{EM}}$ and $\min(\bm{S})=\bar{S}-A_{\text{EM}}$ are computed from the mean flux $\bar{S}$. From the above, assuming a power-law relation between the Doppler factor and $\bm{S}_{\text{r}}$ with spectral index $\alpha$, we calculate the SMBHB orbital radius. Considering that the jet originates from one of the black holes in a binary, but not two, and from the total mass of the binary $M_{\text{tot}}$, which can be computed based on the galaxy velocity dispersion relation, the chirp mass $\mathcal{M}$ is obtained. To sum up, $\mathcal{M}=\mathcal{G}(\theta_{\text{EM}},\bm{\varphi})$, where $\bm{\varphi}=(f_{\text{k}}$, $\beta_{\text{app}}$, $M_{\text{tot}}$, $\alpha)$ and $\mathcal{G}$ is a function with the aforementioned operations. For more details, please refer to Appendix B and Iguchi et al. (2010).

IV Results

IV.1 Model Selection

We first perform model selection by calculating the Bayesian evidence in favor of an SMBHB at 3C 66B against the null hypothesis. For this calculation, we match our prior on the CGW parameters to the results of EM searches. This approach allows for precise confirmation or exclusion of the EM model. In particular, for the $\log_{10}$ chirp mass, we use a skewed Normal prior with hyperparameters for the skewness, mean, and the standard deviation which have values of $-3.08$, $9.12~[\log_{10}M_{\odot}]$ and $0.42~[\log_{10}M_{\odot}]$, respectively. We found empirically that this prior has the same unequal uncertainty levels as constraints from the EM data in Table 1. Similarly, we fix the luminosity distance based on the redshift of 3C 66B, assuming the $\Lambda$CDM cosmological model based on Planck Collaboration et al. (2020). We use an exponential-Normal prior for $\log_{10}f_{\text{GW}}$, with mean and variance from Table 1. We perform this calculation with only the Earth term of the CGW. Finding $\ln\mathcal{B}=-0.0027(7)$, we claim no evidence for an SMBHB at 3C 66B; nevertheless, we cannot rule out the EM model. Finding no evidence, we then calculate upper limits on the parameters of a potential SMBHB in 3C 66B.

Exponentiating our result, the Bayes factor is $\mathcal{B}=0.9973(6)$, whereas Agarwal et al. (2025) find a value of $\mathcal{B}=0.293(3)$. This disagreement is due to noise fluctuations, which are capable of yielding Bayes factors of up to $\sim 1$ according to the null distribution (see Figure 1, left, in Agarwal et al., 2025). The use of different pulsars, different instruments, as well as different noise models can influence the fluctuations. Additionally, the Bayes factor calculated in this study is based on the Earth term, whereas the Bayes factor in Agarwal et al. (2025) is calculated based on both the Earth term and the pulsar term.

IV.2 Limits on the Supermassive Black Hole Binary Chirp Mass in 3C 66B

We perform parameter estimation of the log-10 chirp mass. As in the previous calculation, we assume the luminosity distance as a fixed parameter based on the known redshift of 3C 66B and Planck 2018 results (Planck Collaboration et al., 2020). The difference is that now we use a broad uniform prior $\pi(\log_{10}\mathcal{M})=\mathcal{U}(7.5,9.5)~[\log_{10}M_{\odot}]$. This provides a more EM-independent perspective consistent with the decreasing probability of finding an SMBHB toward high $\mathcal{M}>10^{9}~M_{\odot}$ in nature at $z=0$ (Hopkins et al., 2007), defined by the SMBHB mass function. We also use a more convenient prior on $\pi(\log_{10}(f_{\text{GW}}))=\mathcal{N}(-7.22,0.01)~[\log_{10}\text{Hz}]$, which is still consistent with the EM model.

The results are shown in Figure 2. When only the Earth term is modeled in our analysis, we find an upper limit of $\mathcal{M}<8.25\times 10^{8}~M_{\odot}$ at 95% credibility. When including the pulsar term, we find a more constraining limit, $\mathcal{M}<6.90\times 10^{8}~M_{\odot}$. For comparison, we show the value from the EM model from Iguchi et al. (2010), $(7.9^{+3.8}_{-4.5})\times 10^{8}$ $M_{\odot}$. Our limits are not below the prediction of the EM model. Therefore, we cannot rule it out. However, the tightest upper limit calculated rules out 58% of the uncertainty region of the EM-model. For consistency with Arzoumanian et al. (2023) and Agarwal et al. (2025), we also obtain a posterior by reweighting our posteriors to a uniform prior in linear space $\pi(\mathcal{M})=\mathcal{U}(10^{7.5},10^{9.5})~[M_{\odot}]$. While not matching astrophysical expectations for the SMBHB mass function, this presents the most conservative limits. With this prior, we report an upper limit at $95\%$ credibility to be $\mathcal{M}<11.06\times 10^{8}~M_{\odot}$ for the Earth term and $\mathcal{M}<10.63\ \times 10^{8}~M_{\odot}$ when the pulsar term is included. We find a 16% improvement in the chirp mass upper limit when including the pulsar term. This difference is on the larger end of the 5-15% improvement found in previous studies (Charisi et al., 2024; Zhu et al., 2016), which compare cases where the signal is detected with a signal-to-noise ratio between 5 and 15. Upper limits use information from the tail of the distribution, whereas detections use the bulk of the distribution, which may explain the larger improvement when including the pulsar term in this work. Our limits are a factor of 1.55 lower compared to Arzoumanian et al. (2020b), a factor of 1.33 lower compared to Arzoumanian et al. (2023), and similar to Agarwal et al. (2025).

We also report the result of our reanalysis of the 3 mm data from Iguchi et al. (2010). Following the methodology in Section III.4, we obtain a posterior on the log-10 chirp mass of the SMBHB in 3C 66B. It is shown as the black shaded area in Figure 2. The posterior is skewed toward higher masses, consistently with the uncertainties reported in Iguchi et al. (2010) and presented in Table 1.

IV.3 Limits on the Gravitational-wave Strain from 3C 66B

Instead of calculating the strain amplitude $h_{0}$ of the CGW from 3C 66B based on the luminosity distance using Equation 2, we now treat it as a free parameter. This way, the distance to 3C 66B is not involved in the calculation. Similar to the chirp mass prior, we use the prior $\pi(\log_{10}h_{0})=\mathcal{U}(-17.0,-12.5)$, which is more consistent with the decreasing probability of an SMBHB with a high $h_{0}$. The other priors are as in Section IV.2. The posterior on $h_{0}$ is shown in Figure 3. With only the Earth term, we find an upper limit at 95% credibility to be $\log_{10}h_{0}<-14.12$. When including the pulsar term, we find a more constraining limit, $\log_{10}h_{0}<-14.44$. Using the more conservative prior, $\pi(h_{0})=\mathcal{U}(10^{-17.0},10^{-12.5})$, we find an upper limit for the Earth term of $\log_{10}h_{0}<-13.86$ and an upper limit of $\log_{10}h_{0}<-14.20$ when including the pulsar term. As in Section IV.2, the range of $h_{0}$ of the EM model from Iguchi et al. (2010) is partially ruled out in all cases.

It is also of interest to compare the strain amplitude of a CGW signal from 3C 66B with the expectation for the GWB at the EM-predicted frequency of the binary. Goncharov et al. (2025) propose a fiducial prior on the characteristic strain amplitude of the background, $h_{\text{c}}(f)=A(f~\text{yr})^{-2/3}$, characteristic of SMBHBs. Where $A$ is the strain amplitude at the gravitational wave frequency $f$. The prior is based on black hole masses inferred from kinematic observations and galaxy velocity dispersions. Based on this prior, we calculate $h_{0,\text{GWB}}=2h_{\text{c}}/(2^{2/3}\sqrt{32/5})$, to revert the rms averaging over the SMBHB inclination and polarization assumed for the background (see Equation 8 in Burke-Spolaor et al. (2019)). The best-fit $h_{\text{c}}(f=\text{yr}^{-1})$ of the GWB from the data of the European Pulsar Timing Array by Goncharov et al. (2025) is within $3\sigma$ interval of the prior, above the prior mean. Thus, the prior is consistent with observations.

The prior from Goncharov et al. (2025) marginalized over $\pi(f_{\text{CGW}})$ is shown as the solid black line in Figure 3, elucidating a connection between the EM model of the CGW candidate 3C 66B and the expected GWB. The comparison suggests the following:

• •

the EM model of 3C 66B corresponds to a loud source in the expected GWB;

• •

a fraction of the parameter space of the EM model for 3C 66B, which we have ruled out, is also not supported by the predicted GWB amplitude; and

• •

our limits on $f_{\text{CGW}}$ of 3C 66B are, on average, higher than the expected GWB because we are limited by white noise at the frequency of 3C 66B.

The value of $h_{0,\text{GWB}}$ can also be interpreted as the equivalent expected CGW amplitude of the background if it is dominated by one source (3C 66B) at the 3C 66B CGW frequency. This scenario is possible, given the scaling of the number of SMBHB sources contributing to the background with frequency. This may also lead to a deviation from the power-law model for $h_{\text{c}}(f)$ at the frequency of 3C 66B, depending on the number density of SMBHBs (Phinney, 2001; Agazie et al., 2025). It is also worth noting that Zhu et al. (2019) find that the EM model for 3C 66B is in tension with the observed GWB amplitude based on the expected SMBHB merger rate.

IV.4 Upper Limits as a Function of the Continuous Gravitational-Wave Frequency

Following Arzoumanian et al. (2020b), for a broader perspective on how PPTA DR3 limits on $\mathcal{M}$ depend on $f_{\text{GW}}$ for any SMBHB at the position of 3C 66B, we show Figure 4. For this calculation, we use the same prior as in Section IV.2, except the prior on the CGW frequency is now $\pi(\log_{10}f_{\text{GW}})=\mathcal{U}(-9.0,-7.0)~[\log_{10}\text{Hz}]$. We also only model the Earth term of the CGW signal. The figure highlights that the PPTA is well suited to target 3C 66B, as the value suggested by the EM model is at the frequency to which our targeted search is most sensitive. Consistent with Figure 2, when modeling the pulsar term and using astrophysically motivated priors, we are close to ruling out the best-fit value of the EM model. We observe the highest density of posterior samples at around $f_{\text{GW}}=\text{yr}^{-1}$, which corresponds to uncertainties associated with measuring pulsar positions. There is also a minor overdensity of posterior samples at $f_{\text{GW}}=(2~\text{yr})^{-1}$, which corresponds to uncertainties associated with measuring pulsar parallaxes. Given that the posterior is not uniform at these frequencies, some of the uncertainties may be systematic. We find that the chirp mass upper limits are constrained to lower values in the targeted search as compared to the results from Figure 7 in the PPTA DR3 all-sky search (Zhao et al., 2025).

IV.5 Frequentist Search

To confirm the nondetection of a CGW, we perform a frequentist targeted search based on the $\mathcal{F}$ statistic (Ellis et al., 2012), a frequentist measure of evidence for a CGW in our data. In particular, we chose to calculate the $\mathcal{F}_{\text{e}}$ statistic, which is based on only modeling the Earth term of the CGW signal. This statistic is a function of CGW frequency and sky position, which we set to the 3C 66B EM model values. We find a value of $\mathcal{F}_{\text{e}}=4.90$, corresponding to the false alarm probability $\text{FAP}=0.50$. Ellis et al. (2012) suggest a threshold for FAP to be $10^{-4}$ to claim a detection. The trials factor for our search is one. The background distribution of the $\mathcal{F}_{\text{e}}$ statistic is the chi-squared distribution with four degrees of freedom and is displayed in Figure 5 along with the value of the $2\mathcal{F}_{\text{e}}$ statistic. Based on this null distribution, we find a $p\text{-value}=0.04$ which corresponds to about $2.1\sigma$. Thus, the $\mathcal{F}_{\text{e}}$ statistic value we find is consistent with noise.

V Discussion

V.1 Multimessenger Joint Likelihood

Here, we lay out the groundwork for the joint analysis of CGW and EM data, which may be explored in future work. Because EM and CGW data are independent, the total likelihood of these data is a product of the likelihoods: $\mathcal{L}(\bm{\delta t},\bm{S}|\bm{\theta},\bm{\theta}_{\text{EM}})=\mathcal{L}(\bm{\delta t}|\bm{\theta})\mathcal{L}(\bm{S}|\bm{\theta}_{\text{EM}})$, where $\bm{\theta}=(\bm{\theta}_{\text{GW}},\mathcal{M})$. Effectively, the common parameter in both the CGW model and the EM model from Iguchi et al. (2010) is the chirp mass $\mathcal{M}$. Therefore, the posterior is $\mathcal{P}(\bm{\theta}_{\text{GW}},\bm{\theta}_{\text{EM}},\varphi|\bm{\delta t},\bm{S})=\mathcal{L}(\bm{\delta t},\bm{S}|\bm{\theta}_{\text{GW}},\bm{\theta}_{\text{EM}},\varphi)\pi(\bm{\theta}_{\text{GW}})\pi(\bm{\theta}_{\text{EM}})\pi(\bm{\varphi})$. Here, we also express $\mathcal{L}(\bm{\delta t}|\bm{\theta}_{\text{GW}},\theta_{\text{EM}},\varphi)=\mathcal{L}(\bm{\delta t}|\bm{\theta}_{\text{GW}},\mathcal{M})=\mathcal{G}(\theta_{\text{EM}},\bm{\varphi})$.

V.2 Feasibility of Standard Siren Cosmology using Targeted Searches with Pulsar Timing Arrays

Standard sirens are compact binaries that emit GW and EM radiation simultaneously, allowing for measurement of the Universe’s expansion rate (Schutz, 1986). The redshift of the host galaxy of a standard siren from EM information enables the determination of its recessional velocity $v$, whereas the proper distance $D$ is measured based on the GW signal. The Hubble–Lemaître constant, which quantifies the Universe’s expansion at the present time, is

[TABLE]

The primary challenge in the standard siren approach is localizing the EM counterpart of a GW. To this date, the first observed binary neutron star inspiral, GW170817 (Abbott et al., 2017a), is the only GW signal for which an EM counterpart is established, for galaxy NGC 4993. The sky position uncertainty for GW170817 is $28\ \text{deg}^{2}$ at 90% credibility. For PTAs, the sky position uncertainty of all-sky CGW searches is expected to be less precise. Moreover, the signals are not transient, which makes identification of the host galaxy even more difficult. For example, IPTA DR3 simulations (116 pulsars, 20 yr baseline) from Petrov et al. (2024) found that for an ideal case where the $\text{SNR}=15$, the 90% credible area would range from $29\ \text{deg}^{2}$ to $241\ \text{deg}^{2}$. However, for a realistic case, with $\text{SNR}=8$, the 90% credible area increases from $287\ \text{deg}^{2}$ to $530\ \text{deg}^{2}$. Petrov et al. (2024) conclude that without better CGW localization, host galaxy searches are impractical. Truant et al. (2025) have additionally simulated a 30 yr Square Kilometer Array PTA observation and found that out of the $\sim 35$ CGWs detected, only 35% have a 90% credible area $<$ 100 $\text{deg}^{2}$. They conclude that observations with an SNR $>$ 25 will have fewer than 100 host galaxy candidates. Even in simulations which look at a very high SNR, such as those performed by Zhu et al. (2016), with a 15 yr baseline and 30 pulsars, a strong signal with $\text{SNR}=100$ yields a localization of $\sim 40$ $\text{deg}^{2}$. Therefore, standard siren cosmology based on all-sky CGW searches is not feasible with currently operating PTAs. Nevertheless, the use of standard sirens has been explored with simulated future PTA datasets, mostly using dark sirens, and results suggest that precision cosmology may be possible (Wang et al., 2025; Yan et al., 2020; Xiao et al., 2025).

We propose that by taking advantage of targeted searches for SMBHB candidates, we may be able to revisit the prospect of standard siren cosmology with PTAs.

Current standard siren methods mostly focus on follow-ups from all-sky searches (Jin et al., 2025), which do not offer a priori knowledge of the source position. In contrast, targeted searches based on EM candidates assume the known sky location of the host galaxy. It is important to note, however, that targeted searches do not localize a source within a large uncertainty region, rather they test if the CGW model is consistent with the location of the candidate. This means that source misspecification is possible if there exist multiple hosts that overlap in parameter space. Further studies to determine the extent of the host identification problem and possible solutions are ongoing.

As a demonstration, let us consider a toy example of obtaining constraints on the Hubble–Lemaître constant, provided precise knowledge of the sky position and redshift (Hunt et al., 2021; van den Bosch et al., 2015) and the (lack of) GW emission from 3C 66B. Because $z=0.02\ll 0.1$ and because we do not expect stringent constraints on $D$, we assume $D=D_{\text{L}}$ and use the linear approximation $v=cz$ for simplicity. To develop high-precision cosmology with PTAs as part of future work, it is necessary to use full expressions for $D$ and $v$, while also correcting the recession velocities based on peculiar velocities of galaxies (Abbott et al., 2017b). In any case, the key requirement of a PTA in cosmology, based on targeted CGW searches, is to constrain the luminosity distance of an SMBHB candidate.

Therefore, to obtain a posterior on $H_{0}$, we place constraints on the luminosity distance $D_{\text{L}}$ to 3C 66B with the PPTA DR3 data. We impose a prior $\pi(D_{\text{L}})=\mathcal{U}(10^{-2},10^{3})$ Mpc. We chose a prior on $\log_{10}f_{\text{GW}}$ as in Section IV.2, and we chose a prior on $\log_{10}\mathcal{M}$ as in Section IV.1. Based on the posterior on $\pi(D_{\text{L}})$ and Equation 6, we obtain a posterior on $H_{0}$ in Figure 6. The values $H_{0}$ are shown along the bottom horizontal axis. The values of $\pi(D_{\text{L}})$ corresponding to $H_{0}$ based on Equation 6 and the known redshift of 3C 66B are shown along the top horizontal axis. Assuming that the EM model is correct, since we do not detect a GW signal, we rule out the nearest luminosity distances. This, in turn, limits the highest values of $H_{0}$. We find that the PTA constraints on $H_{0}$ based on the targeted search generally recover the prior, other than the largest values of $H_{0}$. Our choice of $\pi(D_{\text{L}})$ is presented as a toy model, which is why we do not currently report the $95\%$ credible upper-limit value for $H_{0}$. Because we do not claim any detection of a CGW signal from 3C 66B, we emphasize that these constraints are subject to the validity of the CGW model, where 3C 66B is a true SMBHB. In other words, the EM model for 3C 66B must be correct in order for the constraints to be meaningful. Nevertheless, a strong detection of a CGW from an SMBHB in a targeted search, especially with the resolved pulsar term, may prove to be useful in constraining $H_{0}$, and investigations into this matter are underway.

VI Conclusion

Using PPTA DR3, we conduct a targeted multimessenger search for an SMBHB at the center of 3C 66B. We search for a CGW signal by performing model selection, where we compare a model with a CGW signal to the null hypothesis. With priors carefully matched to the EM model, we find a log Bayes factor of $\ln\mathcal{B}=-0.0027(7)$. This suggests we can neither confirm nor rule out the EM model. We also find a frequentist statistic, $\mathcal{F}_{e}=5.04$, with an associated $\text{FAP}=0.50$, which confirms nondetection of the signal. Based on the background estimation, we find that the value of $\mathcal{F}_{e}=5.04$ corresponds to $p=0.04$, about $2.1\sigma$.

Following nondetection, we place upper limits on the mass as $\mathcal{M}<6.90\times 10^{8}~M_{\odot}$ and strain amplitude $\log_{10}(h_{0})<-14.44$, at 95% credibility.

For both the mass and the strain, we rule out over 50% of the parameter space of the EM model. Using conservative priors, we still rule out a portion of the EM model, finding a $95\%$ upper limit on the chirp mass of $\mathcal{M}<9.55\times 10^{8}~M_{\odot}$ and for the strain amplitude we find a value of $h_{0}<-14.20$. Finally, we demonstrate a new methodology for a joint-likelihood approach for the analysis of 3C 66B as well as a proof of principle, showing that it may be possible to use firmly established SMBHB candidates in targeted searches with PTAs as standard sirens, complementary to high-precision cosmology.

VII Acknowledgments

Murriyang, CSIRO’s Parkes radio telescope, is part of the Australia Telescope National Facility (https://ror.org/05qajvd42) which is funded by the Australian Government for operation as a National Facility managed by CSIRO. We acknowledge the Wiradjuri people as the Traditional Owners of the Observatory site. This research has made use of the NASA/IPAC Extragalactic Database, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. Parts of this work were supported by the Australian Research Council Centre of Excellence for Gravitational Wave Astronomy (OzGrav) through grant CE230100016

Appendix A A Continuous Gravitational Wave from a Supermassive Black Hole Binary

Let us start by considering a plane wave traveling from an SMBHB to the SSB to which the pulse arrival times are referenced. A unit vector that points from the CGW source to the SSB is

[TABLE]

where $(\hat{x},\hat{y},\hat{z})$ is the Cartesian coordinate basis with the north celestial pole in the $\hat{z}$-direction, and the vernal equinox is in the $\hat{x}$-direction. The polar angle $\theta$ and azimuthal angle $\phi$ of the source are related to the equatorial coordinates as $\text{R.A.}=\phi$ and $\text{decl.}=\pi/2-\theta$.

The wave is characterized by the metric strain rank-$2$ tensor in the transverse traceless gauge, indexed with $(a,b)$:

[TABLE]

Components $h_{+,\times}$ are the polarization amplitudes, and $e^{+,\times}_{ab}$ are the polarization tensors. Rotating the Cartesian coordinate system such that it is oriented toward a CGW source, we obtain new basis vectors in $(\hat{x},\hat{y},\hat{z})$ as in Arzoumanian et al. (2023):

[TABLE]

The polarization angle $\psi$ is the angle between the natural polarization basis (at source) and the reference polarization basis (at SSB). With this, polarization tensors are expressed as

[TABLE]

The response of a PTA to the plane wave is described by the antenna pattern function (Romano and Cornish, 2017)

[TABLE]

where $\hat{u}$ is the unit vector pointing from the SSB to the pulsar. This function describes the pulsar’s geometrical sensitivity to the CGW source (Taylor et al., 2016). Pulsar timing delays $\delta t$ are then expressed as

[TABLE]

where we implied that $\delta t_{+,\times}$ is the time integral of $h_{+,\times}$, and $\delta t^{\prime}_{+,\times}$ is the difference between the delay experienced at the SSB (Earth) and the signal induced at the pulsar. This difference is

[TABLE]

where $t$ is the time when the CGW passes the SSB and $t_{\text{p}}$ is the time when the CGW passes the pulsar. One can then geometrically relate $t$ and $t_{\text{p}}$, as

[TABLE]

where $L$ is the distance to the pulsar.

For an SMBHB, at zeroth post-Newtonian order, and choosing a coordinate system where $\psi=0$, we arrive at our model of the CGW projected onto the PTA detector:

[TABLE]

where $\iota$ is the SMBHB inclination angle, $D_{\text{L}}$ is the luminosity distance, $\omega(t)=\pi f_{\text{GW}}$ is the angular orbital frequency222The SMBHB angular orbital frequency is half the CGW angular frequency $2\pi f_{\text{GW}}$., $\Phi(t)$ is the phase, and

[TABLE]

is the chirp mass, where $m_{1,2}$ are the SMBHB component masses. The last observation in the PPTA dataset (2022 March 8, or MJD 59646) is used as the reference time. The orbital phase and frequency of the SMBHB are then defined as

[TABLE]

with $\Phi_{0}$ being the initial orbital phase and $\omega_{0}$ the initial orbital frequency; one may notice that this equation is equivalent to Equation 3.

The chirp mass and frequency are redshifted for CGW sources at cosmological distances, such that the observed (“$o$”) values are related to intrinsic source (“$s$”) values as

[TABLE]

For 3C 66B, $z=0.02\ll 1$, so given our measurement uncertainty we assume $\mathcal{M}=\mathcal{M}_{o}=\mathcal{M}_{s}$ and $f=f_{o}=f_{s}$.

Appendix B Analysis of the 3 mm Data

Although key details of the analysis of the 3 mm data are provided in Iguchi et al. (2010), here we provide explicit equations that we used in our reanalysis. Throughout, we switch between the velocity $v$ and the boost parameter $\beta=v/c$, which may also have subscripts “bl,0” and “app” for the absolute binary component velocity around the center of gravity and the apparent relativistic jet velocity, respectively. First, we point out that Equations 1 and 2 from Iguchi et al. (2010) can be expressed directly in terms of $(\beta,\iota)$, as

[TABLE]

where $r$ is the ratio

[TABLE]

Next, $(\beta,\iota)$ are required in the calculation of the absolute value of the rotational velocity of one of the binary components, which has a jet, around the center of mass. This velocity causes Doppler boosting of the flux. We solve Equation 4 from Iguchi et al. (2010) for this velocity, obtaining

[TABLE]

Considering Keplerian motion in a circular orbit, we use the following expressions of the orbital radius $R$ and the orbital separation $a$,

[TABLE]

Introducing $x\equiv R/a$, one finds that $x=m/M_{\text{tot}}$ with $m$ the mass of the lighter black hole for the case where a heavier black hole emits a jet (smaller radius). Similarly, one finds $x=M/M_{\text{tot}}$ with $M$ the mass of the heavier black hole for the case where the lighter black hole emits a jet (larger radius). The total mass, $M_{\text{tot}}$, is obtained by using the relationship between the central stellar velocity dispersion and black hole masses (Ferrarese and Merritt, 2000; Merritt and Ferrarese, 2001). Using this relationship Noel-Storr et al. (2007) obtain a total mass for 3C 66B of $1.9^{+1.0}_{-0.8}\times 10^{9}\ M_{\odot}$. Apparently, for the calculation of the chirp mass, these two cases reduce to the same answer due to a symmetry in $x(1-x)$ in the final expression

[TABLE]

The value of $x(x-1)$ is also known as the symmetric mass ratio.

Appendix C Pulsar Distances

Pulsar distances are required in the calculation of the expected signal when modeling the pulsar term of the CGW signal. Pulsar distances are measured using very long baseline interferometry (VLBI) astrometry, optical astrometry, timing, and the Shklovskii (1970) effect. A list of the priors used for each pulsar, as well as the technique used to find this distance, are listed in Table 2. We impose a Gaussian prior on pulsar distance in accordance with the estimates found in the literature. When no such estimates are available for a given pulsar, we impose a sufficiently conservative prior $\mathcal{N}(1000,200)~\text{[pc]}$. Distance estimates based on pulsar dispersion measures are not used due to the large uncertainties in the electron density model (Yao et al., 2017).