The NANOGrav 15 yr data set: Posterior predictive checks for gravitational-wave detection with pulsar timing arrays

We use the NANOGrav 15 yr data set, which contains 67 pulsars that have been timed for more than 3 years, with 16.03 yrs of data between the first and the last time of arrival in the data set [26]. We use the DMX dispersion measure noise model [44] and white noise parameters included in the NANOGrav 15 yr data release [26]. For each pulsar, a best fit timing model is constructed that accounts for deterministic effects like Roemer delay, proper motion, parallax, binary orbits, etc., which is then subtracted from the TOAs to produce a set of timing residuals for each pulsar, $\delta\mathbf{t}$. Stochastic processes like achromatic intrinsic spin-wandering and GW background-induced delays are included in this initial fit as a single “total red noise” contribution, as the first pass analysis is done on a pulsar-by-pulsar basis and so we cannot separate intrinsic pulsar noise from the GW background.

In this subsection, we discuss the full PTA model that is used to search for a GW background. Readers familiar with this already can skip to Sec. III, although later we will make frequent reference to equations introduced in this section. For a more in depth presentation of the PTA analyses, see Refs. [45, 46, 47].

The starting point for the analysis are the timing residuals, $\delta\mathbf{t}$. We characterize stochastic processes like intrinsic pulsar noise and the GW background in the frequency domain using a Fourier matrix $\mathbf{F}$ and associated amplitudes $\mathbf{a}$ [48]. The stochastic processes are covariant with elements of the timing model (specifically the frequency, spin-down, and dispersion measure variations), and so we also introduce deviations from the best-fitting timing model parameters, $\bm{\epsilon}$. We assume these deviations are small, such that changes in $\delta\mathbf{t}$ are linear in changes in $\bm{\epsilon}$ with a design matrix $\mathbf{M}$ made up of derivatives of $\delta\mathbf{t}$ with respect to the timing model parameters. Putting these effects together, we have a model for the residuals

where we have consolidated the frequency domain representation and timing model corrections,

If radio frequency interference is effectively excised and standard pulse profiles are accurate, the resulting noise is dominated by radiometer noise and “pulse profile jitter” which is traditionally assumed to be frequency independent and Gaussian. This leads to a Gaussian likelihood for the timing residuals

where the covariance matrix $\mathbf{N}$ describes the measurement noise of the individual observations, and is block-diagonal. TOAs at different radio frequencies from the same individual observation are correlated with one another due to pulse profile jitter [49], but TOAs from different observations are uncorrelated.

We assume that the GW background and the intrinsic pulsar noise are stationary, and so they can be characterized by the power spectrum of the GW background, correlations between pulsars, and the power spectrum of the intrinsic pulsar noise in each pulsar. The assumption of stationarity for the GW background should hold if the dominant contribution to the background is an ensemble of SMBHBs emitting at roughly constant frequencies. The assumption that intrinsic pulsar noise is stationary is one of expedience that should be tested. Tests on the European Pulsar Timing Array second data release show no signs of non stationarity [50].

Information about the power-law amplitude and spectral index for the intrinsic pulsar noise and the GW background is encoded in the covariance matrix of the sine and cosine amplitudes $\mathbf{a}$ across pulsars. We introduce a set of hyperparameters $\bm{\Lambda}$ to characterize these power laws. We place a Gaussian prior on $\mathbf{b}$,

We use an improper uniform prior on $\bm{\epsilon}$ so that its posterior is determined by the likelihood. This prior is now broadcast across $\mathbf{b}$ parameters for each pulsar. The covariance matrix of the $\mathbf{a}$ coefficients is given by $\bm{\varphi}(\bm{\Lambda})$, which contains blocks corresponding to correlations of the Fourier modes between pulsars. Diagonal blocks encode information about the power spectrum of the total red noise for a given pulsar, including the intrinsic pulsar noise, $\bm{\eta}_{a}(\bm{\Lambda})$ (where the $a$ subscript labels the pulsar) and the GW background spectrum $\bm{\rho}(\bm{\Lambda})$. Off-diagonal blocks between pulsars $a$ and $b$ contain (scaled) contributions from the GW background. Putting all of this together, the covariance matrix for $\mathbf{a}$ is

where $i$ and $j$ label frequencies and $\Gamma_{ab}$ corresponds to the correlations between pulsars. Different angular correlation patterns correspond to different models. In this paper we consider four models. The first states that $\Gamma_{ab}$ follows the Hellings–Downs curve (hd model) that is expected from an isotropic GW background,

where $\theta_{ab}$ is the angle between pulsars $a$ and $b$ on the sky. The second is that $\Gamma_{ab}=\delta_{ab}$, which we call the common uncorrelated red noise, curn, model. We will also consider a mono model that is characterized by monopolar correlations, $\Gamma_{ab}=1$ and a model with dipolar correlations, dip, with $\Gamma_{ab}=\cos\theta_{ab}$. Theoretical models indicate that $\rho_{i}(\bm{\Lambda})$ will roughly take the form of a power law, and past empirical studies suggest that $\eta_{ai}(\bm{\Lambda})$ often follows a power-law as well,

where $A_{\textrm{gw}}$ is the amplitude of the GW background at $f_{\textrm{yr}}=(1\textrm{yr})^{-1}$, $\gamma_{\textrm{gw}}$ is the negative spectral index, $A_{\textrm{rn,a}}$ is the amplitude of intrinsic pulsar noise for pulsar $a$ and $\gamma_{\textrm{rn,a}}$ its associated spectral index. The frequency is given by $f_{i}=i/T$, and $T$ is the time between the first and last TOAs in the data set. For intrinsic pulsar noise we use 30 frequencies, $i\in[1,30]$ and for the GW background we use $i\in[1,14]$. These numbers were chosen based on individual pulsar fitting (for the intrinsic pulsar noise), and a dedicated curn analysis that allows for the common spectrum to “flatten” at high frequencies, where it then becomes indistinguishable from white noise.

The power-law models for the GW background and intrinsic pulsar noise spectra have amplitudes and spectral indices associated with them: $A_{\textrm{gw}}$, $\gamma_{\textrm{gw}}$ for the GW background and, $A_{\textrm{rn,a}}$, and $\gamma_{\textrm{rn,a}}$ for each of the $N_{p}$ pulsars in the array. We collectively denote these parameters as $\bm{\Lambda}$. To reduce the total number of parameters we need to infer, we typically marginalize over the model parameters $\mathbf{b}$, leaving a posterior on the hyperparameters

The covariance matrix is now $\mathbf{C}=\left(\mathbf{N}+\mathbf{T}\mathbf{B}\mathbf{T}^{T}\right)$, and we introduced a prior on the hyperparameters $p(\bm{\Lambda})$. We also note that $p(\mathbf{b}|\delta\mathbf{t},\bm{\Lambda})\propto p(\delta\mathbf{t}|\mathbf{b})p(\mathbf{b}|\bm{\Lambda})=\mathcal{N}(\widehat{\mathbf{b}},\bm{\Sigma})$, which is normal with mean and covariance given by

We estimate the marginalized posterior on $\bm{\Lambda}$ using stochastic sampling methods [51] because of the large dimension of $\bm{\Lambda}$ ($2N_{p}+2$ in the case described above). This yields $N_{s}$ samples $\{\bm{\Lambda}^{s}\}_{s=1}^{N_{s}}$ approximately drawn from the posterior,

Below we use $\delta\mathbf{t}$ to refer to generic timing residuals, $\delta\mathbf{t}^{15\textrm{yr}}$ to refer to residuals from the the 15 yr data set, and $\delta\mathbf{t}^{\textrm{rep}}$ to refer to data replications. We use two models to create sets of data replications to compare to the collected data. Each method proceeds along similar lines:

1. 1.

   Choose $\bm{\Lambda}$ by drawing randomly from $p(\bm{\Lambda}|\delta\mathbf{t}^{15\textrm{yr}})$.

2. 2.

   Draw $\mathbf{b}\sim p(\mathbf{b}|\bm{\Lambda})$. The choice of $p(\mathbf{b}|\bm{\Lambda})$ depends upon the set of replications we are performing. We specify details below when we discuss individual replication sets. This method nominally calls for us to draw from the improper prior on $\bm{\epsilon}$, yielding unusable timing residuals. Therefore, we do not simulate timing model variations, and fix $\bm{\epsilon}\approx 0$.

3. 3.

   Draw $\delta\mathbf{t}^{\textrm{rep}}\sim\mathcal{N}(\mathbf{T}\mathbf{b},\mathbf{N})$ where $\mathbf{b}$ comes from the previous step.

The data replications use different models at each stage. We outline the different data replication sets, their purpose, and what models they use to carry out the procedure described above.

• •

  CURNPosteriorDraws: We create simulated data sets based on the curn model which we index with $s$. We draw $\bm{\Lambda}^{s}\sim p(\bm{\Lambda}|\delta\mathbf{t}^{15\textrm{yr}},\textsc{curn})$, and $\mathbf{b}^{s}\sim\mathcal{N}(0,\mathbf{B}(\bm{\Lambda}^{s})|\textsc{curn})$, Eqs. (5) and (6). The conditioning on curn implies no correlations between pulsars, $\Gamma_{ab}=\delta_{ab}$. We do not simulate timing model variations, i.e., $\bm{\epsilon}=0$. These sets of data replications are compared to the data and recovered model parameters.

• •

  HDPosteriorDraws: We create simulated data sets based on the hd model which we index with $s$. We draw $\bm{\Lambda}^{s}\sim p(\bm{\Lambda}|\delta\mathbf{t}^{15\textrm{yr}},\textsc{hd})$, and $\mathbf{b}^{s}\sim\mathcal{N}(0,\mathbf{B}(\bm{\Lambda}^{s})|\textsc{hd})$. The conditioning on hd implies we include Hellings–Downs correlations between pulsars during simulation. We do not simulate timing model variations, i.e., $\bm{\epsilon}=0$. These sets of data replications are compared to the real data and recovered model parameters to assess how consistent the data are with the hd model.

The curn model is characterized by a lack of spatial correlations between pulsars, $\Gamma_{ab}=\delta_{ab}$. To reject this model, we use the optimal statistic signal-to-noise ratio (SNR) as our test statistic [40, 52, 39, 53]. The SNR, for a given choice of noise parameters, is distributed according to a generalized $\chi^{2}$ distribution [54]; it is large when Hellings–Downs correlations are present and centered around zero when no spatial correlations are present.

The optimal statistic depends upon the total red noise in each pulsar (intrinsic pulsar noise and GW background), which we do not know a priori. Therefore, current analyses average the SNR over the posterior distribution on the noise parameters,

where in the second line we perform a Monte Carlo integral using a finite set of $N_{s}$ posteriors samples [53]. $\overline{\textrm{SNR}}$ is then used as the test statistic for null hypothesis testing. The motivation for using this “noise marginalized optimal statistic” is that we are marginalizing over uncertainty in the noise and signal parameters when we calculate statistical significance. One uses “phase shifts” [55], “sky scrambles” [56], or data replications to estimate the null distribution of $\overline{\textrm{SNR}}$ and calculate a $p$-value for the measured $\overline{\textrm{SNR}}$. In Ref. [1], $\overline{\textrm{SNR}}\approx 5$, which falls at the $3.5-4\sigma$ level in the null distributions, indicating that the curn model does not fully describe the data.

Here, we use a more conservative statistic that gives more weight to low-SNR outliers and lower weight to high SNR outliers than the noise marginalized optimal statistic [28]. The cumulative distribution function is not linear in the SNR, and so we first calculate the $p$-value of the SNR calculated on each $\bm{\Lambda}^{s}$, and then average those $p$-values together. The resulting $p$-value will be less significant than the one calculated on $\overline{\textrm{SNR}}$. Conceptually, this can be thought of as averaging over the risk of rejecting the null hypothesis by placing more weight on the most conservative noise realizations. By contrast, calculating a $p$-value on $\overline{\textrm{SNR}}$ weighs high-SNR (and therefore less conservative noise realizations) equally to low SNR noise realizations.

We compare the value of the optimal statistic on the observed data to its value on data replications from the posterior predictive distribution. We calculate a $p$-value on each $\textrm{SNR}(\delta\mathbf{t}^{15\textrm{yr}};\bm{\Lambda}^{s})$, and average those $p$-values. This final, averaged $p$-value is referred to as a posterior predictive $p$-value or a Bayesian $p$-value because it is marginalized over the posterior predictive distribution for the data [57, 58, 27]. We can do this generically, when we do not know the distribution of the test statistic, by calculating

Here $\Theta$ is the Heaviside function, and $p(\delta\mathbf{t}^{\mathrm{rep}}|\bm{\Lambda},\textsc{curn})$ could be one of the data replications described in Sec. II.3, or it could be a set of “bootstrapped” data replications like sky scrambles or phase shifts.¹¹1We use $\delta\mathbf{t}^{\mathrm{rep}}$ to also denote sky scrambled and phase shifted data sets, in addition to actual simulated data sets. This is somewhat poor notation. When performing sky scrambles and phase shifts the timing residuals themselves are the same as $\delta\mathbf{t}^{15\textrm{yr}}$ and it is the sky position (sky scrambles) or $\mathbf{F}$ (phase shifts) that changes. Nevertheless, we use $\delta\mathbf{t}^{\mathrm{rep}}$ to refer generically to any data sets or schemes used to construct the null distribution, for simplicity. If the analytic distribution for the SNR for a given choice of $\bm{\Lambda}^{s}$ is known, then we do not need to actually perform data replications, and $\int\Theta[\cdot]p(\delta\mathbf{t}^{\mathrm{rep}}|\bm{\Lambda},\textsc{curn})\,\mathrm{d}\delta\mathbf{t}^{\mathrm{rep}}$ is the inverse cumulative distribution function for the SNR. The Bayesian $p$-value reduces to

where in the second line we evaluate the integral numerically using draws from $p(\bm{\Lambda}|\delta\mathbf{t}^{15\textrm{yr}})$. The superscript “rep” indicates that the inverse cumulative distribution function on the measured $\textrm{SNR}(\delta\mathbf{t}^{15\textrm{yr}};\bm{\Lambda})$ is calculated over (theoretical or actual) data replications or sky scrambles.

The probability distribution function for the optimal statistic SNR for fixed $\bm{\Lambda}^{s}$, under the noise model, is a generalized $\chi^{2}$ distribution [54] which we will refer to as GX2 moving forward. In Fig. 1 we show $\textrm{SNR}(\delta\mathbf{t}^{15\textrm{yr}};\bm{\Lambda}^{s})$ for 100 draws along the bottom, and each blue curve is $P[\textrm{SNR}^{\textrm{rep,s}}(\delta\mathbf{t}^{\textrm{rep,s}}\bm{;}\bm{\Lambda}^{s})>\textrm{SNR}(\delta\mathbf{t}^{15\textrm{yr}};\bm{\Lambda}^{s})]$, calculated using the GX2 distribution. The dashed red line gives $p_{B}=7\times 10^{-4}$, which corresponds to $3.2\,\sigma$ significance in favor of rejecting the curn model. By contrast, the significance of the SNR maximum-likelihood draw from $p(\bm{\Lambda}|\delta\mathbf{t})$ calculated using GX2 was $\approx 2\times 10^{-4}$ or $3.5\,\sigma$.

In using the GX2 distribution, we assumed that the noise model is correct. Instead, we can construct replications of our data set that break correlations due to GWs, but preserve any mismodeling in the noise that might cause large SNR values in favor of correlations. The two main methods for doing this are sky scrambles [56] and phase shifts [55]. For each $\bm{\Lambda}^{s}$, we perform 400,000 sky scrambles, where we artificially move the location of the pulsars to different positions on the sky drawn uniformly on the two-sphere and calculate a “new” Hellings–Downs curve using these new positions. Using these sky scrambles, we build a null distribution and calculate significance. We repeat this for 100 draws of $\bm{\Lambda}^{s}$ and average the inverse CDFs as in Eq. (19). Under this procedure, we find $p_{B}=1\times 10^{-4}$, which corresponds to an equivalent Gaussian significance of $3.7\,\sigma$. Using sky scrambles, Ref. [1] found $p=5\times 10^{-5}$ using the traditional procedure of building a null distribution for $\overline{\textrm{SNR}}$. The results are shown in Fig. 2, where the inverse CDFs for sky scrambles (green) in general fall off faster than for the GX2. However, there are some outliers resulting in $p_{B}$ being larger than the $p$-value calculated on $\overline{\mathrm{SNR}}$ using scrambles.

It is unclear why the inverse CDF for sky scrambles generally falls off faster than for the GX2; this is an open area of investigation [59]. In previous work, methods of generating a background distribution from sky scrambling or phase shifting use a “match statistic,” in an attempt to use scrambles or shifts that are quasi-independent of one another and the Hellings–Downs curve [60, 61]. Recently, in Ref. [62], the authors suggested only sky scrambles that produce correlation curves that are independent of one another should be used, where independence is achieved by insisting the match statistic disappear. In Ref. [1], the condition is that the match threshold between any one sky scramble and all others is $\lesssim 0.2$. Here we do not use a match statistic, as our goal is to estimate the probability that the pulsars would be arranged on the sky in such a way that noise fluctuations would produce Hellings–Downs correlations. To test this, we must draw the positions uniformly on the sky. How to produce reliable null distributions for data sets that preserve potentially unmodeled noise is still subject to exploration.

In the previous subsection, we reject the null hypothesis of the curn model at the 3.2 $\sigma$ level. In this section, we use the same scheme to test whether the data are consistent with Hellings–Downs correlations. Given that we only have an analytic null distribution for the optimal statistic, we use a Monte Carlo integral for Eq. (18) to evaluate $p_{B}$ in the presence of a potential signal:

where we average over $N=1,000$ HDPosteriorDraws replications. We show a scatter plot in Fig. 3, where the $x$-axis shows $\mathrm{SNR}(\delta\mathbf{t},\bm{\Lambda}^{s})$ and the $y$-axis shows $\mathrm{SNR}(\delta\mathbf{t}^{\textrm{rep, s}},\bm{\Lambda}^{s})$, and $p_{B}$ corresponds to the fraction of points above the line $y=x$. In this case, we find $p_{B}=0.534$, indicating that the 15 yr NANOGrav data are consistent with data replications that assume Hellings–Downs correlations.

The model for a GW background assumes spatial correlations that follow the HD curve, but other spatial correlations could arise either from statistical fluctuations or due to mismodeling. Monopolar correlations could arise due to an error in the clock at each site, corresponding to a correlated offset that is common to all pulsars. Dipolar correlations could arise due to an error in the effective location and motion of the solar system barycenter [20]. We use the multiple component optimal statistic [52], which estimates the amplitude of monopolar, dipolar, and Hellings–Downs correlations simultaneously, to test whether our estimate of these correlations in the 15 yr data set is consistent with data replications from a pure HD model. The results and methods here follow App. H of Ref. [1]. The analysis is nearly identical, but with more simulations used to calculate $p_{B}$. The results and conclusion are the same as that analysis, but we include it here both for completeness, and because it is strongly related to the rest of the new tests we have performed in this section.

The multiple component optimal statistic simultaneously produces the SNR for all three spatial correlations, SNR${}^{\textrm{MONO}}_{\textrm{MC}}$, SNR${}^{\textrm{DIP}}_{\textrm{MC}}$, and SNR${}^{\textrm{HD}}_{\textrm{MC}}$ where the superscript corresponds to the spatial correlation and the subscript indicates that we are using the multiple component optimal statistic. We again use 1,000 HDPosteriorDraws data replications described in Sec. II and calculate Eq. (20) substituting $\textrm{SNR}^{\textrm{MONO}}_{\textrm{MC}}$ for SNR to produce $p_{B}^{\textrm{MONO}}$. We produce $p_{B}^{\textrm{DIP}}$ and $p_{B}^{\textrm{HD}}$ for dipole and Hellings–Downs correlations defined analogously.

We show similar visualizations to the previous section in Fig. 4. We find $p_{B}^{\textrm{HD}}=0.64$, which again indicates that the HD SNR calculated on the 15 yr data is consistent with what we expect from the pure hd model. Likewise, we find $p_{B}^{\textrm{DIP}}=0.26$, consistent with no dipolar correlations. Finally, we find $p_{B}^{\textrm{MONO}}=0.11$, which is largely consistent with no monopolar correlations.

In this subsection, we summarize the methods outlined in Ref. [29]. Given $p(\bm{\Lambda}|\delta\mathbf{t}^{15\textrm{yr}})$, we calculate $\mathbf{a}$ in two ways. In one method, $\mathbf{a}$ are conditioned on $\delta\mathbf{t}^{15\textrm{yr}}$ and $\bm{\Lambda}$, which we refer to as the “inferred” coefficients (subscript “inf”) because they are drawn from the inferred posterior on $\mathbf{a}$ using information from both the power-law spectrum prior and the real data. In the other method, $\mathbf{a}$ are conditioned only on $\bm{\Lambda}$, which we refer to as “predicted” coefficients (subscript “pre”) because these are the coefficients predicted by the power-law spectrum prior. In both cases, we marginalize over $\bm{\Lambda}$ and $\bm{\epsilon}$. We illustrate the workflow in Fig. 5. The posteriors on the inferred and predicted parameters are formally given by

The first term in the integrand differs between the predicted (no dependence on $\delta\mathbf{t}^{15\textrm{yr}}$) and inferred posteriors (which are conditioned on $\delta\mathbf{t}^{15\textrm{yr}}$). We discuss specifics of how these terms are evaluated below. The second term in the integrand is the posterior on $\bm{\Lambda}$, e.g., power-law amplitudes and spectral indices.

In both Eqs. (21) and (22), we evaluate the posterior with a Monte Carlo integral. We first draw a sample (labeled with “$s$” superscript) $\bm{\Lambda}^{s}\sim p(\bm{\Lambda}|\delta\mathbf{t}^{15\textrm{yr}})$. We then draw from the first term in the integrand. For the inferred coefficients we draw from $p(\mathbf{a},\bm{\epsilon}|\bm{\Lambda}^{s},\delta\mathbf{t}^{15\textrm{yr}})$, which is a Gaussian with mean and covariance that depend on both the prior and the data, and are given by Eqs. (14) and (15). For the predicted coefficients, we draw from just the prior distribution (i.e., no dependence on data), $p(\mathbf{a},\bm{\epsilon}|\bm{\Lambda}^{s})$, which is a zero-mean Gaussian given by Eqs. (5) and (6). More details on this scheme are discussed in Ref. [29].

For each $\bm{\Lambda}^{s}$ we split $\mathbf{a}$ into $\mathbf{a}_{\textrm{rn},a}^{s}$ for intrinsic pulsar noise, and $\mathbf{a}_{\textrm{gw},a}^{s}$ for GWs for each pulsar $a$. We then reconstruct the power spectrum for the intrinsic pulsar noise in each pulsar and the GW power observed by each pulsar. Each pulsar sees a different realisation of the GW background, but $\mathbf{a}_{\textrm{gw},a}^{s}$ are drawn from the same distribution for all pulsars $a$. By contrast, the intrinsic pulsar noise is a different power spectrum for each pulsar, and therefore $\mathbf{a}_{\textrm{rn},a}^{s}$ are drawn from a different distribution for each individual pulsar.

For the inferred coefficients, by conditioning on the data, the power spectrum will deviate from a power law if the true data-generating process differs from a power law. For the predicted spectrum, we obtain different realizations of a power spectrum that are consistent with a power-law. In Sec. IV.2, we discuss results for the inferred and predicted intrinsic pulsar noise and GW spectra for each pulsar, and compare them to an “excess noise” analysis done in Ref. [41]. At each frequency, we use a modified version of the optimal statistic²²2See Sec. IVB2 in Ref. [29] for a discussion. to combine individual-pulsar coefficients to estimate the total the total GW power across the PTA in that frequency bin [63]. We present the results in Sec. IV.3.1.

Finally, for each $\bm{\Lambda}^{s}$, we produce pulsar pair-wise correlations and compare them to the expected Hellings–Downs curve. To do this, we draw coefficients from $p_{\mathrm{inf}}(\mathbf{a}|\delta\mathbf{t}^{15\textrm{yr}},\textsc{hd})$ and $p_{\mathrm{pred}}(\mathbf{a}|\delta\mathbf{t}^{15\textrm{yr}},\textsc{hd})$, and use the optimal statistic to construct pair-wise correlations. In the case of the predicted parameters, this will give us an expected spread on correlation vs. angular separation for a given model. For the inferred parameters, by conditioning on the data, the correlations can deviate from the model. In Sec. IV.3.2, we look at reconstructions of the pair-wise correlations as a function of angular separation, and search for deviations from the Hellings–Downs curve.

Power spectra for each pulsar in the NANOGrav array are explored in Ref. [41], it is therefore worth contrasting the two results. First, Ref. [41] simultaneously estimated the total red noise ($\rho_{i}^{2}+\eta_{ai}^{2}$) in each frequency bin $i$, performing a separate analysis for each individual pulsar $a$. A Savage-Dickey Bayes factor was calculated to estimate the significance of the total red noise at each frequency for each pulsar. Next, both the common red noise and intrinsic pulsar noise were fixed to the maximum likelihood values estimated from a curn analysis that assumes these spectra follow a power-law. Once fixing these parameters in their model, excess noise in each frequency bin for each pulsar was searched for. No evidence for excess noise was found, the power-law model for intrinsic pulsar noise and the common red noise processes are therefore sufficient.

In this work, we instead separate intrinsic pulsar noise and GW background contributions when drawing parameters from the inferred and predicted distributions, and we produce a posterior distribution on both contributions in each frequency bin for each pulsar. This way we are testing both the intrinsic pulsar noise and GW background power-law assumptions at the same time, while constructing full posteriors on the intrinsic pulsar noise and the GW background. This is in contrast to the search for excess noise on top of a power-law common red noise and intrinsic pulsar noise. Another difference is that in this work, individual frequency bin estimates are subject to a prior that follows a power-law, while Ref. [41] used a log-uniform prior on the power in each frequency bin.

In Fig. 6, we show results for two pulsars with strong intrinsic pulsar noise, B1937+21 and J1012+5307 (top two panels), and J1909-3744 which has no measurable intrinsic pulsar noise, but a contribution attributed to the GW background. The green boxes correspond to the estimates of the total red noise power from [41], which was discussed above. The blue and pink boxes correspond to the inferred intrinsic pulsar noise and GW background contributions respectively, orange boxes correspond to the predicted intrinsic pulsar noise, and yellow boxes correspond to the predicted GW background in the bottom panel. Similar plots for each pulsar are included in an electronic supplement.

In the top two panels, the intrinsic pulsar noise (inferred in blue boxes, predicted in orange boxes) typically agrees with the total red noise (green boxes) from Ref. [41]–indicating that the total red noise is dominated by intrinsic pulsar noise. The GW background is significantly below the intrinsic pulsar noise and the total red noise. Additionally, the orange distributions, corresponding to predicted intrinsic pulsar noise, agree with the inferred intrinsic pulsar noise. We quantify this agreement below. In the bottom panel, the total red noise agrees with the GW background, while the intrinsic pulsar noise is significantly lower. This is consistent with the GW background contributing significantly to $\delta\mathbf{t}^{15\textrm{yr}}$ in this pulsar, with no intrinsic pulsar noise contribution. In a few cases the free spectrum total red noise (green boxes) deviates further from the power law than the inferred intrinsic pulsar noise (blue boxes). This is due to the power-law prior used for the inferred intrinsic pulsar noise, which will tend to move those parameters closer to the power law.

In a few situations, e.g., the first and fifth through seventh bins for J1909-3744, the estimated total red noise (green boxes) appears to be lower than the predicted (yellow boxes) and inferred spectra (pink boxes) for the GW background. The low total noise are consistent with the predicted distributions, which follow a $\chi^{2}$ with two degrees of freedom, and therefore have large support at those low values (despite what the box and whiskers show). The inferred distributions broadly agree with the predicted, and do show overlap with the green whiskers; but they do look inflated compared to the total noise. However, simply combining the GW background and intrinsic pulsar noise spectra in each bin may not reproduce the green boxes for a few reasons. First, interference between these two contributions may cancel or amplify the estimated total red noise above or below what we would expect from naively adding their contributions incoherently. Second, a log-uniform prior on the power in each frequency will likely reduce the upper limit on the estimated power in that bin when no red noise detection is made when compared to the upper limit we would set using a prior informed by the power law model

We then quantify excess intrinsic pulsar noise in individual frequency bins for each pulsar. For each pulsar $a$ and frequency bin $i$ we calculate

where $\eta^{2}_{\textrm{inf},ai}$ are drawn from $p_{\rm inf}(\mathbf{a}|\delta\mathbf{t}^{15\textrm{yr}})$, and correspond to the inferred power spectrum, while $\eta_{\textrm{pred},ai}^{2}$ are drawn from $p_{\rm pred}(\mathbf{a}|\delta\mathbf{t}^{15\textrm{yr}})$ and correspond to the predicted power spectrum due to a power law. Each $\eta^{2}$ carries with it an implicit $s$ index, which we have suppressed. We also use the HDPosteriorDraws simulations to calculate

where the superscript “rep” indicates it is the inferred estimate on the power in that frequency calculated on the replicated data. Note that $\eta^{2}_{\textrm{inf},ai}$ and $(\eta_{\textrm{inf},ai}^{2})^{\mathrm{rep}}$ are calculated using the same $\bm{\Lambda}^{s}$. We find that these two methods produce nearly identical results, and so we report results for $p_{B}$ instead of $p_{B}^{\rm sim}$.

For intrinsic pulsar noise across all pulsars and frequencies, we find a minimum of $p_{B}=0.03$, for J1012+5307, $f=51\,\textrm{nHz}$, which is the box that is visibly above the max likelihood curve in the middle panel of Fig. 6, marked with the red arrow. This is not a significant $p$-value, given that we are analyzing 67 pulsars and 30 frequency bins, and so we cannot conclude that this represents a deviation from a power law. This is the same conclusion as Ref. [41]. The minimum and maximum $p_{B}$ for the GW background power spectrum across all pulsars are 0.16 and 0.84.

Deviations from the power-law model may not just take the form of excess noise at individual frequencies. For example, one could have a broken power-law, or excess noise across multiple frequencies that are not individually detectable. We do not develop a statistic to measure this here, as it requires a specific model to compare to the power-law model and there are a broad range of potential models. However, such an analysis should be done in the future.

We show a plot of $p_{B}$ for intrinsic pulsar noise for each pulsar in the top panel of Fig. 7 and for the GW background in each pulsar in the bottom panel. The horizontal axis corresponds to each pulsar, while the vertical axis corresponds to $p_{B}$; each point represents a $p_{B}$ for each pulsar and each frequency bin. When the inferred and predicted power spectrum estimates agree at a given frequency bin, we expect $p_{B}\approx 0.5$. We see that a few noisy pulsars show $p_{B}$ values that stray away from 0.5, including the pulsars in the top two panels of Fig. 6. As stated before, no individual frequency bin shows an extreme value of $p_{B}$, e.g., $p_{B}<0.01$ or $p_{B}>0.99$, meaning we cannot state there are individual bins with excess noise. However, pulsars like B1937+21 and J1012+5307 do show quite a few frequency bins with $p_{B}$ deviating from 0.5, meaning they may benefit from a more flexible model in the future. It is currently prohibitively computationally expensive to estimate intrinsic pulsar noise separately in each frequency bin for each pulsar when we estimate $p(\bm{\Lambda}|\delta\mathbf{t}^{15\textrm{yr}})$. However, with future computational improvements, we may be able to do this for a limited number of pulsars, and this method provides a good starting point for choosing those pulsars.