Multi-messenger Approaches to Supermassive Black Hole Binary Detection and Parameter Estimation. II. Optimal Strategies for a Pulsar Timing Array

A PTA regularly monitors a number of millisecond pulsars (MSPs), which are characterized by rapid and extremely stable rotation periods. It is therefore sensitive to fluctuations in the arrival times of the radio pulses induced by GWs emitted by astrophysical sources including SMBHBs. We first simulate ten realizations of the MSP population using the pulsar population synthesis package PsrPopPy ⁹ (Bates et al. 2014) in the snapshot mode, where simulated MSPs are generated by randomly drawing properties from distributions constrained by earlier studies. The overall population size is constrained by halting generation when the population contains the same number of MSPs (28) that would have been discovered by the Parkes Multibeam Survey (Manchester et al. 2001). We assume that the MSPs are distributed radially from the Galactic Center according to a Gaussian distribution with a standard deviation of 7.5 kpc and exponentially above and below the galactic plane with a mean scale height of 500 pc (Lorimer 2012). We assume that the MSP spin periods are distributed according to Lorimer (2012) and that their luminosities (at 1.4 GHz) are distributed according to a log-normal distribution with $\langle {\mathrm{log}}_{10}L\rangle =-1.1$ and ${\sigma }_{{\mathrm{log}}_{10}L}=0.9$ (Faucher-Giguere & Kaspi 2006), where L is in units of mJy kpc². Aggarwal & Lorimer (2022) found that MSP spectral indices are best described by a broken power law with a break at 300 MHz, but since all of our simulated observations occur above 300 MHz, we simply use the high frequency part of their model: Gaussian with a mean spectral index of −1.47 and a standard deviation of 0.7. We assume that the MSP pulse widths are distributed log₁₀ normal with a mean of 0.05 P^0.9 and a standard deviation of 0.3. Finally, dispersion measures (DM) and scattering timescales are computed by the NE2001 electron density model of the Milky Way (Cordes & Lazio 2002).

In order to construct a PTA from the simulated MSP population, we simulate a timing program with the proposed Deep Synoptic Array-2000 (DSA-2000; Hallinan et al. 2019), where an anticipated 25% of the time will be used for pulsar timing. We simulate timing observations of each MSP using the full DSA-2000 collecting area operating in the frequency range 0.7–2 GHz and assume an integration time of 60 minutes per pulsar per observation epoch. We assume that DSA-2000 will have a system temperature of 25 K, a gain of 10 K Jy⁻¹, and will be able to observe declinations from −30° to 90°. If a model pulsar is within the decl. limits, its pulse duty cycle is <90%, and it is sufficiently bright to be detected, we then include it in our mock PTA and compute the uncertainty on its pulse time of arrival σ_TOA, which is dependent on both the pulsar and telescope parameters. To do so, we use the FrequencyOptimizer ¹⁰ software package (Lam et al. 2018), which defines the TOA uncertainty as ¹¹

$\begin{eqnarray*}&&{\sigma }_{\mathrm{TOA}}^{2}={\sigma }_{\widehat{\mathrm{DM}}}^{2}+{\sigma }_{\delta {t}_{C}}^{2}+{\sigma }_{\mathrm{DM}(\nu )}^{2}+{\sigma }_{\mathrm{tel}}^{2}.\end{eqnarray*}$

Here, ${\sigma }_{\widehat{\mathrm{DM}}}$ is the standard error on the infinite-frequency TOA when fitting the timing residuals for the DM; it contains sources of white noise including template-fitting errors, scintillation errors, and jitter errors, where the jitter errors are computed using the scaling relationship shown in Figure 7 of Lam et al. (2019)

$\begin{eqnarray*}&&{\sigma }_{J}={10}^{-0.218}{\delta }^{1.216}\displaystyle \frac{{P}^{3/2}}{\sqrt{{t}_{\mathrm{int}}}},\end{eqnarray*}$

which is a function of the pulsar spin period P, pulse duty cycle δ, and integration time t_int. ${\sigma }_{\delta {t}_{C}}$ is the systematic error due to chromatic variations in the pulse width caused by interstellar scattering and dispersion. σ_DM(ν) is the uncertainty in the DM estimation between two frequency channels. σ_tel contains uncertainties due to radio frequency interference, incorrect gain calibration, and instrumental self-polarization. NANOGrav considers σ_TOA < 1 μ s to be the minimum criterion for the inclusion of an MSP in the PTA, so we adopt that cutoff here. The results are ten mock PTAs with an average of ∼180 MSPs with σ_TOA < 1 μ s. The sky distribution and distribution of σ_TOA for one realization are plotted in the first panel of Figure 1 and the left panel of Figure 2 (all_sky), respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we consider two variants of our mock timing program: (1) an all-sky PTA consisting of only pulsars with σ_TOA < 100 ns (all_sky_best) and (2) a targeted campaign in an area of the sky consisting of ∼10 pulsars within a ∼10° radius, where the GW source is chosen to be located within that sky area (targeted). Since the antenna pattern functions increase with a closer (but not perfect) alignment between the source and the pulsar (see the definitions in Ellis et al. 2012), a PTA campaign which targets pulsars near an EM-selected source may be more advantageous for detecting its GW signal. We plot the pulsars in one realization of all_sky_best as large filled circles in the second panel of Figure 1, and those in targeted as stars in the third panel. Their respective σ_TOA distributions are shown in the last two panels in Figure 2. The properties of a sample of simulated MSPs are shown in Table 1, and ten realizations of the mock population are available in machine-readable form online.

Table 1. Sample of Simulated Pulsars

P	$\dot{P}$	DM	tscatter	W50	l	b	R.A.	Decl.	S1400	L1400	α	dtrue	X	Y	Z	σTOA	σtel	σδDM	σW	σJ	Ared	γred	Aeff/Aeff,tot
(ms)	(s s−1)	(pc cm−3)	(ms)	(ms)	(deg)	(deg)	(deg)	(deg)	(mJy)	(mJy kpc2)		(kpc)	(kpc)	(kpc)	(kpc)	(μs)	(μs)	(μs)	(μs)	(μs)	(μs yr1/2)
6.64	2.19E-21	615.09	7.42	0.80	−0.39	1.85	264.39	−28.29	9.11E-04	0.08	−0.91	9.18	−0.06	−0.67	0.30	−1.00	−1.00	−1.00	−1.00	1.09	6.50E-04	−5.60	1.00
3.19	1.86E-19	177.42	2.01	0.06	6.82	−9.80	279.87	−27.63	0.004	0.15	−2.08	5.86	0.69	2.77	−1.00	1757.16	0.01	1757.16	5.68	0.05	0.09	−5.60	1.00
3.27	9.73E-21	336.38	3.86	0.11	−15.20	−7.19	264.06	−45.61	0.002	0.33	−2.45	14.57	−3.79	−5.45	−1.82	−2.00	−2.00	−2.00	−2.00	0.11	0.02	−5.60	1.00
4.42	1.37E-19	112.95	0.00	0.48	33.95	15.35	269.57	7.82	8.00E-04	0.03	−2.70	6.26	3.37	3.49	1.66	−1.00	−1.00	−1.00	−1.00	0.53	0.02	−5.60	1.00
14.04	4.74E-21	638.35	10.42	0.50	−3.47	0.52	263.76	−31.60	0.012	0.71	−3.04	7.63	−0.46	0.89	0.07	−2.00	−2.00	−2.00	−2.00	1.00	8.17E-04	−5.60	1.00
25.52	7.16E-20	1415.46	1143.50	0.79	12.99	−0.46	273.89	−17.89	2.26E-04	0.08	−0.18	19.35	4.35	−10.36	−0.16	−1.00	−1.00	−1.00	−1.00	2.12	1.95E-04	−5.60	1.00
7.00	7.94E-21	552.88	1514.56	0.21	54.88	0.25	293.04	19.55	2.54E-06	7.14E-04	−2.13	16.76	13.71	−1.14	0.07	−1.00	−1.00	−1.00	−1.00	0.29	0.00	−5.60	1.00
14.44	2.79E-19	943.11	500.37	0.29	−43.43	−0.23	221.07	−60.08	0.002	0.54	−1.69	15.14	−10.41	−2.50	−0.06	−2.00	−2.00	−2.00	−2.00	0.58	0.02	−5.60	1.00
28.56	5.91E-21	215.86	0.04	1.29	4.94	10.38	259.74	−19.19	0.006	0.52	−0.80	9.23	0.78	−0.54	1.66	101.11	0.04	101.11	27.34	3.62	3.04E-04	−5.60	1.00
1.98	6.71E-21	372.40	1.68	0.34	−2.35	4.21	260.91	−28.64	0.002	0.16	−0.64	8.46	−0.35	0.07	0.62	−1.00	−1.00	−1.00	−1.00	0.25	0.01	−5.60	1.00

Notes. P—pulse period; $\dot{P}$—period derivative; DM—NE2001 DM; t_scatter–scattering timescale ; W₅₀—intrinsic pulse FWHM; and l, b—galactic longitude and latitude. Note that in PsrPopPy, −180° < l < + 180°. R.A., decl.—R.A. and decl.; S₁₄₀₀—flux at 1.4 GHz; L₁₄₀₀—luminosity at 1.4 GHz; α—spectral index; d_true—true distance from Earth; and X, Y, Z—galactic Cartesian positions. The Sun is at X = 0, Y = 8.5 kpc, Z = 0 in this coordinate system. σ_TOA—total TOA uncertainty; σ_tel—telescope noise; σ_δDM—DM estimation uncertainty: ${\sigma }_{\delta \mathrm{DM}}=\sqrt{{\sigma }_{\widehat{\mathrm{DM}}}^{2}+{\sigma }_{\delta {t}_{C}}^{2}+{\sigma }_{\mathrm{DM}(\nu )}^{2}}$; σ_W—white noise errors; σ_J—jitter noise; A_red, γ_red—red noise power spectrum amplitude and spectral index, respectively; and A_eff/A_{eff, tot}—fraction of the total collecting area used in observation.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Then, for each pulsar sample in each realization, we construct mock PTA observations which have a baseline of 15 yr, which is approximately the length of the current PTAs. The cadence of observations is chosen to be either monthly or weekly, to imitate the current "monthly" and "high-cadence" campaigns of NANOGrav (Alam et al. 2021).

We consider a binary system with component masses m₁ and m₂. The so-called chirp mass ${ \mathcal M }$ is given by ${ \mathcal M }\equiv {({m}_{1}{m}_{2})}^{3/5}/{({m}_{1}+{m}_{2})}^{1/5}$. We assume the system is in a circular orbit whose GW frequency is related to its orbital frequency f_GW = 2f_orb and is often expressed as ω = π f_GW. The source is located at a luminosity distance d_L; since in this work we will only consider nearby sources, we ignore the effect of redshift: (1 + z) ≈ 1. For simplicity, we only consider the "Earth term" in the timing residual and do not consider the "pulsar term," ¹² and hence only the combination ${{ \mathcal M }}^{5/3}/{d}_{{\rm{L}}}$ can be constrained. It is often convenient to redefine the amplitude as $\alpha ={{ \mathcal M }}^{5/3}/{d}_{{\rm{L}}}/{\omega }^{1/3}$.

The sky position of the source is written as polar and azimuthal angles ϕ and θ, respectively, (corresponding to R.A. and decl. in units of radians: ϕ = R. A., θ = π/2 − decl. ). We assume the frequency evolution of the binary over the course of ∼15 yr is negligible (i.e., a monochromatic signal). We further assume the BHs are non-spinning, since PTAs are largely insensitive to the effects of spin in an SMBHB system (see Sesana & Vecchio 2010).

Such a binary system is then described by the following parameters: p = {α, ω, ϕ, θ, i, ψ, Φ₀}, where i, ψ, and Φ₀ are the inclination, GW polarization angle, and initial orbital phase, respectively. For simplicity, we adopt fixed values of ψ = π/2, Φ₀ = π/2, and an intermediate inclination of i = π/4 for each of the sources, since these are the less astrophysically interesting parameters for our study. We then follow the usual definitions in Ellis et al. (2012) and Aggarwal et al. (2019) to compute the Earth term timing residual Res corresponding to each pulsar which has coordinates ϕ_psr and θ_psr.

The first source (S1) which will be "observed" by our mock PTA is the circular and non-spinning analog of the well-known SMBHB candidate OJ 287 (e.g., Sillanpaa et al. 1988; Lehto & Valtonen 1996; see, e.g., Dey et al. 2019; Valtonen et al. 2021 for recent reviews). Adopting the binary parameters in Dey et al. (2018), we estimate a GW amplitude of α ≈ 3 ns. We note that while it is possible to compute a more accurate waveform for this complex binary candidate (see, e.g., Valtonen et al. 2021), in this work we only seek to obtain estimates of its detection and parameter measurement prospects.

As can be seen in Figure 1, S1 is located in an area of the sky where few mock pulsars are nearby (e.g., within a ∼10 deg radius), thereby making a future targeted pulsar timing campaign with DSA-2000 unlikely. Prompted by this, we consider a mock source (S2) which (1) has similar binary parameters but is located at the part of the sky where it may be observed by a targeted campaign and (2) represents a generic SMBHB candidate, whereas S1 is based on an actual candidate. We additionally consider S3, which shares the same binary parameters as S2 but is at a less advantageous location for pulsar searching and timing (and hence the possibility of a targeted campaign). On the other hand, it is located outside the Galactic plane and suffers from less reddening and extinction effects. Thus, S3 represents sources which could be easily detected or observed by most EM observatories. The parameters of the three sources are listed in Table 2, where the S1 parameters are based on Dey et al. (2018).

Table 2. SMBHB Parameters

Source	R.A.	decl.	${ \mathcal M }$	fGW	dL	α
	deg	deg	$\mathrm{log}{M}_{\odot }$	log Hz	log Mpc	ns
S1	133.7036	20.1085	9	−8.28	3.2	3.4
S2	300	30	9	−8	2.5	13.9
S3	250	−10	9	−8	2.5	13.9

Note. Ψ = π/2, Φ₀ = π/2, and i = π/4 for each source. S1 represents an OJ 287-like SMBHB, whereas S2 and S3 are mock SMBHBs.

Download table as:  ASCIITypeset image

To summarize, we observe the three sources with two all-sky PTAs with either weekly or monthly cadence over 15 yr. However, an all_sky campaign at a weekly cadence is not under consideration here since the total observing time would exceed the 25% available for pulsar timing. For S2, we additionally consider a targeted campaign at both cadences. Those mock programs are also summarized in Table 3.

The Fisher matrix is often used to assess the parameter estimation performance of an experiment in fields such as cosmology (e.g., Albrecht et al. 2006) and GW astrophysics (e.g., Sesana & Vecchio 2010; McGrath & Creighton 2021). While it can lead to a less accurate prediction of the performance of the experiment (see, e.g., Vallisneri 2008) than the mock data analysis approach, it is significantly less computationally expensive and thus can be applied to a large number of experimental setups and realizations. In this section, we describe how we apply the Fisher matrix formalism to predict the prospects of mock PTAs as described in Section 2.1 of measuring the parameters of the three SMBHBs as depicted in Section 2.2.

We begin by first computing the Fisher matrix for each pulsar (indexed by a)

$\begin{eqnarray*}&&{{ \mathcal F }}_{{ij}}^{a}=\displaystyle \sum _{b}\displaystyle \frac{1}{{\sigma }^{2}}\displaystyle \frac{\partial \mathrm{Res}b}{\partial {p}_{i}}\displaystyle \frac{\partial {\mathrm{Res}}_{b}}{\partial {p}_{j}},\end{eqnarray*}$

where i and j are the indices of the binary parameters, b denotes each observation, σ represents the σ_TOA of the a-th pulsar as discussed in Section 2.1, ¹³ and Res is the Earth term timing residual corresponding to the a-th pulsar as described in Section 2.2. Here we fix the values of ϕ and θ to mimic the search for a CW signal where the source has been pinpointed (which is a reasonable assumption if the candidate is identified electromagnetically), which means we do not compute the partial derivatives with respect to ϕ or θ, and hence i, j = 1...5 which correspond to the parameters {α, ω, i, ψ, Φ₀}.

We then compute the Fisher matrix for multiple pulsars at various locations {ϕ_psr, θ_psr} (i.e., a PTA) by summing the matrices

$\begin{eqnarray*}&&{ \mathcal F }=\displaystyle \sum _{a}{{ \mathcal F }}^{a}.\end{eqnarray*}$

Instead of directly inverting the matrix ${ \mathcal F }$, we first ensure it is well-conditioned, by normalizing the matrix by the factor $\sqrt{{{ \mathcal F }}_{{ii}}{{ \mathcal F }}_{{jj}}}$ so that the diagonal elements are 1 and the off-diagonal elements are of order unity. We then use singular value decomposition to invert the normalized matrix: ${{ \mathcal C }}_{N}={{ \mathcal F }}_{N}^{-1}$. We then divide ${{ \mathcal C }}_{N}$ by the normalization factor to obtain the final covariance matrix ${ \mathcal C }$, where ${ \mathcal C }={{ \mathcal F }}^{-1}$. Finally, we check the accuracy of the inversion by multiplying the covariance matrix by the original Fisher matrix and confirm that its difference from the identity matrix is lower than an error threshold: max($| {\mathbb{I}}-{ \mathcal F }{ \mathcal C }| )\lt {10}^{-3}$. We can then obtain the measurement uncertainty σ_i of parameter p_i from the i-th diagonal element ${\sigma }_{i}^{2}={{ \mathcal C }}_{{ii}}$. In this work, we focus on the uncertainties of the two more astrophysically interesting parameters, α and ω.

Finally, we estimate the total S/N by summing the contribution from each pulsar: ${{\rm{S}}/{\rm{N}}}^{2}={\sum }_{a}\,{{\rm{S}}/{\rm{N}}}_{a}^{2}$, where SNR_a = $\sqrt{{\sum }_{b}{\left(\tfrac{{\mathrm{Res}}_{b}}{\sigma }\right)}^{2}}$. In this work, we focus our analysis in the strong signal regime (defined in this work as S/N > 5) so that σ_i , instead of being a lower limit, approaches the actual measurement uncertainty. Therefore we only record the value of σ_i if S/N is greater than 5 in a realization.

We first investigate whether S1 can be observed by a future all-sky PTA at a either weekly or monthly cadence, i.e., our simulated all_sky and all_sky_best, by performing a Fisher matrix analysis for each simulated PTA. At a monthly cadence, only 2/10 of the realizations of an all_sky PTA exceed our S/N threshold of 5. Their S/Ns slightly decrease in all_sky_best, which is expected because the pulsars in all_sky_best are a subset of the ones in all_sky (solid line in Figure 3). If only the best pulsars (all_sky_best) are monitored at a weekly cadence, 9/10 of the realizations have S/N > 5, and the mean value increases to ∼9 (dashed line).

Zoom In Zoom Out Reset image size

We proceed to compare the measurement uncertainties between the GW- and EM-based methods, where the GW-based uncertainties are computed using the Fisher matrix method in Section 2.3, and the EM-based uncertainties are based on the reported uncertainties of m₁, m₂, and P_orb from Dey et al. (2018). The measurement uncertainty of the GW amplitude (quantified by Δα/α) is at the ∼100% level (middle panel), far exceeding the EM-based measurement uncertainty of ∼0.2% for this binary candidate. By contrast, ω can be constrained at the few percent level despite the modest S/N (right panel); it is however still greater than the EM-based uncertainty of ∼0.1%.

Thus, we tentatively conclude that a future PTA with DSA-2000 which resembles our all_sky or all_sky_best could detect the GW signal from OJ 287 with ∼15 yr of data at >monthly cadence. However, evidence for its binarity based only on GW observations will be modest, unless GW data (especially the more precise measurement of f_GW) are interpreted alongside EM observations. Further, the stochastic GWB has a fiducial power-law spectral shape which increases in amplitude at low GW frequencies; therefore the GWB would have a significant effect on the detectability of sources emitting at (particularly) lower frequencies. See Section 3.4 for a discussion of this effect on the three mock sources considered in this paper.

It is worth noting that the all_sky campaign only increases the S/N by <1 compared to all_sky_best, despite monitoring ∼6 times the number of pulsars. This suggests that if the total number of observations is a limited resource, it may be best spent on the highest quality pulsars, either at the same cadence or a higher cadence. We explore this further in the following sections.

In this section, we consider the detectability and parameter measurement uncertainties of the mock source S2, which has a larger GW amplitude (α ≈ 14 ns) and is located in a part of the sky where more pulsars may be found and monitored. We first compute the expected results of a monthly, 15 year-long campaign. As we show in the left panel of Figure 4, the highest S/N is reached in the all_sky campaign; it decreases in all_sky_best and decreases further in targeted. This can be understood by the fact that targeted has comparable pulsar timing noise as all_sky, but ∼10 times fewer pulsars in the array. Similarly, all_sky_best has much fewer pulsars than all_sky and hence a lower S/N, despite a higher median pulsar timing quality.

Zoom In Zoom Out Reset image size

As we also show in the last two panels of Figure 4, targeted results in poorer constraints on α and ω, consistent with its lowest S/N. Nonetheless, these GW-based parameter measurements are better than the constraints obtained from typical EM observations for the following reasons. (Recall that S2 represents a generic EM binary candidate and therefore its parameters are assumed to be measured from standard methods.) First, an observed variability period caused by modulated accretion in the binary may not directly correspond to its intrinsic orbital period, since this relationship may be mass ratio-dependent and the observed period may be a few times the orbital period (see, e.g., D'Orazio et al. 2013; Farris et al. 2014; Noble et al. 2021); therefore, an EM-based frequency measurement would have an (underestimated) uncertainty of 100%, which is ∼2 orders of magnitude larger than the GW-basement measurement (right panel in Figure 4). Second, standard BH mass measurements suffer from a systematic uncertainty of ∼0.3 dex (e.g., Shen 2013); combined with the aforementioned uncertainty on ω, this translates to a Δα/α of approximately 120%, which is a few times higher than the GW-based measurement (middle panel in Figure 4). These precise GW parameter measurements can allow meaningful comparisons with EM-based measurements as a robust test of the binary model, breaking the degeneracies arising from EM-only observations, and testing the predictions of SMBHB theory.

It is however important to note that all_sky boasts the highest S/N simply due to the fact that it is timing ∼6 times more pulsars (than all_sky_best) and that the combined contribution of the bottom 85% of the pulsars to the S/N only increases the S/N by ∼1. Therefore it would be interesting to investigate whether all_sky_best or targeted is a more efficient approach to achieve a similar scientific outcome.

To do so, we compare the results of the high-cadence, weekly all_sky_best and targeted campaigns with the (monthly) all_sky campaign. As shown in Figure 4, we predict that all_sky_best results in the highest S/N and best parameter measurements; this is despite the fact that all_sky_best costs ∼40% less total telescope time than all_sky. Likewise, a high-cadence targeted campaign requires only of order a third of the time as all_sky, but achieves a similar S/N. Furthermore, both high-cadence campaigns can measure α and ω at precise levels, which are better than typical EM measurements as discussed previously. Therefore, we conclude that all_sky_best and targeted are both possible, if not preferable, alternatives to all_sky.

Finally, we examine the detection and parameter estimation prospects of S3. As we show in Figure 5, the expected outcomes of all_sky and all_sky_best follow the same trend; namely, all_sky results in a higher S/N for the same cadence, but the high-cadence all_sky_best campaign outperforms all_sky at a lower cost of total telescope time. Further, despite the slightly lower S/N than S2 (consistent with the difference in sky location), the parameters of S3 can still be measured to moderate to high precision.

Zoom In Zoom Out Reset image size

As discussed earlier, the sky location of S3 may be more representative of that of a likely EM candidate: it is located where it is more likely to be found by EM observatories, including all-sky surveys; on the other hand, it is not at such a fortuitous location that it happens to be in close proximity to a large number of (high quality) MSPs. Fortunately, its detection prospects are nonetheless good in both all-sky campaigns and are only mildly impacted by its "inferior" location. This suggests that, contrary to conventional wisdom, dedicated timing campaigns targeting pulsars near likely sources may not be necessary, efficient, or possible for CW observations.

Recently, a red noise process which has a common spectrum across pulsars has been observed in regional PTAs and the combined International Pulsar Timing Array data set (Arzoumanian et al. 2020; Chen et al. 2021; Goncharov et al.2021; Antoniadis et al. 2022). While the PTAs have not confirmed this common red noise process as the GWB, it already needs to be accounted for in more recent CW searches (e.g., Arzoumanian et al. 2023) and will have an effect on PTA's sensitivity to CWs, especially at lower frequencies. Further, some level of red noise is present in many, if not most, MSPs, which would further compromise the observability of CW sources.

While our Fisher matrix methodology does not handle red noise, in this section, we seek to understand the potential effect of the GWB and pulsar red noise on the detectability of CW sources using simulated PTA sensitivity curves and speculate their implications for binary parameter estimation.

The GWB has a power-law power spectrum of the form

$\begin{eqnarray*}&&{h}_{c}(f)={A}_{\mathrm{GWB}}{\left(\displaystyle \frac{f}{{\mathrm{yr}}^{-1}}\right)}^{\alpha },\end{eqnarray*}$

where we adopt the median value of the amplitude from the NANOGrav 12.5 yr analysis: A_GWB = 1.92 × 10⁻¹⁵ (Arzoumanian et al. 2020) and α = − 2/3 (Phinney 2001).

The power spectrum of the pulsar red noise is given by (e.g., Alam et al. 2021)

$\begin{eqnarray*}&&P(f)={A}_{\mathrm{RN}}^{2}{\left(\displaystyle \frac{f}{{\mathrm{yr}}^{-1}}\right)}^{{\gamma }_{\mathrm{RN}}},\end{eqnarray*}$

where A_RN and γ_RN are the amplitude and the index, respectively. We simulate A_RN and γ_RN for each mock pulsar using the model developed by Shannon & Cordes (2010). They determined that the natural log of the measured, red timing noise after a second-order polynomial fit is normally distributed as $\mathrm{ln}({\sigma }_{\mathrm{TN},2})\sim { \mathcal N }\left(\mathrm{ln}({\hat{\sigma }}_{\mathrm{TN},2}),\delta \right)$. ${\hat{\sigma }}_{\mathrm{TN},2}$ is the post-fit red noise scaled to the spin period P and period derivative $\dot{P}$ given by

$\begin{eqnarray}&&{\hat{\sigma }}_{\mathrm{TN},2}={10}^{15}{C}_{2}{P}^{-\alpha }{\left(\displaystyle \frac{\dot{P}}{{P}^{2}}\right)}^{\beta }{T}_{\mathrm{yr}}^{\gamma }\,\mathrm{ns},\end{eqnarray} \tag{ 1 }$

where C₂, α, and β are estimated over the pulsar population and T_yr is the timescale of the residuals in years. We assume that the spin period derivatives are uncorrelated with spin period and are distributed log-normally with $\langle {\mathrm{log}}_{10}\dot{P}\rangle \,=-19.9$ and ${\sigma }_{{\mathrm{log}}_{10}\dot{P}}=0.5$. We adopt the updated values for C₂, α, β, δ, and γ, estimated for a larger sample of MSP red noise measurements, in Table 4 of Lam et al. (2017), row "MSP_10,PPTA + NANO." The red noise index, γ_RN =− (2γ + 1), is assumed to be −5.6 for all MSPs. By combining Equation (1) and Equation (15) of Lam et al. (2017) for the measured red noise in terms of A_RN and γ_RN, we can directly draw A_RN for each model pulsar where

$\begin{eqnarray*}&&f({A}_{\mathrm{RN}})=\displaystyle \frac{1}{\sqrt{2\pi {\delta }^{2}}}\ \exp \left[-\left(\displaystyle \frac{{\left(\mathrm{ln}({\hat{A}}_{\mathrm{RN}})-\mathrm{ln}({A}_{\mathrm{RN}})\right)}^{2}}{2{\delta }^{2}}\right)\right],\end{eqnarray*}$

and

$\begin{eqnarray*}&&{\hat{A}}_{\mathrm{RN}}={\left(-(1+{\gamma }_{\mathrm{RN}})\right)}^{1/2}{10}^{15}{C}_{2}{P}^{-\alpha }{\left(\displaystyle \frac{\dot{P}}{{P}^{2}}\right)}^{\beta }\mu {\rm{s}}\ {\mathrm{yr}}^{1/2}.\end{eqnarray*}$

To simulate the PTA sensitivity to CWs in the presence of these noise terms, we use the software package hasasia ¹⁴ (Hazboun et al. 2019). We first compute the sensitivity curve by using the sky positions of simulated all_sky pulsars and only include white timing noise (see Section 2.1). Following the previous sections, we assume a 15 yr PTA with a monthly cadence. The resulting sensitivity is shown as the solid curve in Figure 6. For a monochromatic source (which we assume throughout the paper) emitting at f₀ with a strain amplitude h₀, the expected S/N is ${h}_{0}\sqrt{{T}_{\mathrm{obs}}/{S}_{\mathrm{eff}}({f}_{0})}$, where S_eff is the effective strain–noise power spectral density for the PTA and is related to the sensitivity curve ${h}_{c}(f)=\sqrt{{{fS}}_{\mathrm{eff}}(f)}$ (Hazboun et al. 2019). Therefore, the resulting S/Ns for S1 and S2/S3 are 2.2 and 9.3, respectively, consistent with our results in Section 3 that the sources are detectable at the intermediate to high S/N level at monthly cadence. ¹⁵

Zoom In Zoom Out Reset image size

Next, we recompute the sensitivity in the presence of the GWB. The result is shown as the dashed curve. As expected, the GWB severely decreases the PTA sensitivity to CWs, especially at frequencies lower than ∼2 × 10⁻⁸ Hz where it decreases by up to an order of magnitude. Consequently, S1 would not be detectable above the GWB (S/N = 0.2), while S2 and S3 may still be detectable, albeit with a lower S/N (2.7).

Finally, we include the pulsar red noise and show the resulting sensitivity curve as the dashed–dotted curve. The effect is less pronounced, with the sensitivity decreasing by less than a factor of a few at all frequencies, resulting in an of S/N = 0.2 (2.1) for S1 (S2/S3).

We further anticipate that PTA's ability to estimate binary parameters would also decrease in the presence of red noise. This is based on the fact that at a given frequency, source detectability and parameter measurability closely track each other (see Figures 3–5), i.e., a lower S/N corresponds to a larger measurement uncertainty.

While our Fisher matrix method cannot handle red noise and hence our conclusions in this work regarding source detectability and observing strategies are drawn assuming only white noise (i.e., the solid curve in Figure 6), future work would be able to make more accurate predictions of the detectability of CW sources by taking into account additional noise, especially the GWB. In practice, the (relatively small) effect of pulsar red noise could be mitigated by only including pulsars with low red noise levels in the PTA, while the effect of the astrophysical GWB cannot removed.