PSR J0437-4715: The Argentine Institute of Radioastronomy 2019–2020 Observational Campaign

In order to characterize the S/N of the observations we use the functions getDuration and getSN of the Python package PyPulse¹¹ (Lam 2017). In Figure 1 we show the S/N of each observation as a function of its duration. The mean S/N of observations with A1 is 151 and with A2 is 105, with mean observing times of 147 minutes and 116 minutes, respectively. However, we note that these numbers are affected by many short and low-quality observations.

Zoom In Zoom Out Reset image size

When we restrict our analysis to observations with S/N > 50, the mean S/N for observations with A1 increases to 166 and with A2 to 122, with mean observing times of 162 minutes and 124 minutes, respectively. We summarize these and other values in Table 2.

We observe a positive correlation between S/N and t_obs, fitting to a ${\rm{S}}/{\rm{N}}\propto \sqrt{{t}_{\mathrm{obs}}}$ as expected (Lorimer & Kramer 2012):

$$\begin{eqnarray}&&{\rm{S}}/{\rm{N}}=\sqrt{{n}_{{\rm{P}}}\,{t}_{{\rm{o}}{\rm{b}}{\rm{s}}}\,BW}\left(\displaystyle \frac{{T}_{{\rm{p}}{\rm{e}}{\rm{a}}{\rm{k}}}}{{T}_{{\rm{s}}{\rm{y}}{\rm{s}}}}\right)\displaystyle \frac{\sqrt{\left.W(P-W\right)}}{P},\end{eqnarray} \tag{ 1 }$$

where P is the pulsar period and W its width, T_peak is its maximum amplitude, T_sys is the noise temperature of the system, t_obs is the observing time, and n_P is the number of polarizations observed.

We collect the observations per S/N for each antenna and display them as histograms in Figure 2. We observe a distribution for A1 with a mean higher than the corresponding distribution for A2, perhaps due to the broader band sensitivity of A1.

Zoom In Zoom Out Reset image size

We collect the observations into sets of S/N > 1, 50, 80, 110, 140, and 170, corresponding to roughly the position of the larger bins in the A2 histogram. In Table 2 we specify the number of observations, mean duration, and mean S/N for each of these sets.

In what follows we assume that the expected S/N scales $\propto \sqrt{{t}_{\mathrm{obs}}}$ and that additional variations in the S/N are due to scintillation. We note that the observations described in Section 2 lack absolute flux calibrations and thus possible variations in T_sys are not accounted for. Moreover, RFI is also variable and its mitigation leads to variations in the effective bandwidth of each observation, so additional dispersion in the S/N versus t_obs relation is also expected.

To quantify the variations due to scintillation we build a projected pulse S/N as ${\rm{S}}/{{\rm{N}}}_{\mathrm{proj}}={\rm{S}}/{\rm{N}}\sqrt{{t}_{\max }/{t}_{\mathrm{obs}}}$, with ${t}_{\max }=217\,\mathrm{minutes}$. Given that short observations have a large uncertainty in their determined S/N, we only use observations with t_obs > 20 minutes (which is roughly half of the scintillation timescale). Figure 3 shows a histogram of the projected pulse S/N for A1 and A2. The line shows the estimated probability density function from scintillation (Cordes & Chernoff 1997)

$$\begin{eqnarray}&&{f}_{S}(S| {n}_{\mathrm{ISS}})=\displaystyle \frac{{\left({{Sn}}_{\mathrm{ISS}}/{S}_{0}\right)}^{{n}_{\mathrm{ISS}}}}{S{\rm{\Gamma }}({n}_{\mathrm{ISS}})}\exp \left(\displaystyle \frac{-{{Sn}}_{\mathrm{ISS}}}{{S}_{0}}\right)\ {\rm{\Theta }}(S),\end{eqnarray} \tag{ 2 }$$

where n_ISS is the number of scintles, S₀ is the mean value of the signal S (i.e., S₀ = 〈S/N〉), and Θ is the Heaviside step function. We calculate n_ISS by fitting the normalized¹² data for each antenna. We obtain n_ISS = 2.67 ± 0.31 for A1 and n_ISS = 2.17 ± 0.25 for A2, with S₀ = 127.27 for A1 and S₀ = 87.16 for A2. The bin size is determined using the Knuth's rule (Knuth 2006) algorithm provided in astropy (Astropy Collaboration et al. 2013, 2018), though we confirm that the obtained values do not depend on the binning by repeating the analysis for different bin sizes.

Zoom In Zoom Out Reset image size

In addition, we make use of the long duration of the observations, which is significantly larger than the typical scintillation timescale for J0437−4715. We split the observations in segments lasting ${t}_{\min }=2000\,{\rm{s}}$ and ${t}_{\min }=5000\,{\rm{s}}$ and repeat the previous analysis. In this case we obtain larger values of n_ISS ∼ 5.

For each of these fittings we perform a Kolmogorov–Smirnov (KS) test for goodness of fit. This test quantifies the distance between the empirical distribution of the sample (obtained from the projected S/N) and the cumulative distribution function of the reference distribution (obtained from fitting n_ISS in Equation (2)) under the null hypothesis that the sample is drawn from the reference distribution. The null hypothesis can be rejected at a given confidence level α if the resulting p-value is lower than 1−α. The p-values obtained are summarized in Table 3. For α = 0.9 (90% confidence level) we find that the goodness of fit cannot be statistically rejected for complete observations with either A1 or A2, or split observations of A1, all of which have a large p-value. However, the fits to the split observations of A2 fail this test, suggesting that, for short observations with A2, Equation (2) may not be entirely valid or that the estimate of the projected S/N becomes unreliable.

We compare our values of n_ISS with theoretical estimations following Lam & Hazboun (2020). We scale the scintillation parameters given at the frequency of 1.5 GHz by Keith et al. (2013) to match our observations centered at 1.4 GHz and obtain the scintillation bandwidth Δν_d = 740 MHz and scintillation timescale Δt_d = 2290 s. We calculate n_ISS via the usual formula

$$\begin{eqnarray}&&{n}_{\mathrm{ISS}}\approx \left(1+{\eta }_{t}\displaystyle \frac{T}{{\rm{\Delta }}{t}_{{\rm{d}}}}\right)\left(1+{\eta }_{\nu }\displaystyle \frac{\mathrm{BW}}{{\rm{\Delta }}{\nu }_{{\rm{d}}}}\right)\end{eqnarray} \tag{ 3 }$$

where η_t and η_ν are filling factors ∼0.2. The estimated n_ISS for T = 220 minutes are 2.22 for A1 (BW = 112 MHz) and 2.18 for A2 (BW = 56 MHz). We confirm that the value obtained with A2 is consistent with the expectations, although for A1 it is larger than expected, perhaps due to additional factors affecting the variability observed.

Gwinn et al. (2006) found two scintillation scales observing J0437−4715 in 327 MHz. Rescaling those scales to our observing frequency, 1400 MHz, we find timescales of Δt_d,1 = 5727 s and Δt_d,2 = 515 s, leading to n_ISS,1 = 1.46 for both antennas and n_ISS,2 = 6.58 for A1 and n_ISS,2 = 6.35 for A2. The latter values are close to those displayed in Table 3 for the split observations, consistent with the shorter observations being more sensitive to the shorter-scale scintillations. Note also that those scintillations scales have been observed to vary notably between epochs (Smirnova et al. 2006).

We compute the timing residuals of the TOAs using Tempo2 (Hobbs et al. 2006) and its Python wrapper, libstempo (Vallisneri 2020), with the timing model given in the file J0437−4715.par provided by IPTA and adapted to the IAR observatory.¹³ Tempo2 returns: (i) the MJD, residual, and template-fitting error (σ_TOA) of each observation, and (ii) the timing model parameters, the weighted errors of the residuals (rms), and the ${\chi }_{\mathrm{red}}^{2}={\chi }^{2}/{n}_{\mathrm{free}}$ of the timing model fit to the residuals. The χ² test considers a good fit when ${\chi }_{\mathrm{red}}^{2}\sim 1;$ instead, a value of ${\chi }_{\mathrm{red}}^{2}\gg 1$—assuming the timing model is correct—indicates the presence of outliers or an underestimation of the residual errors. In this case we can:

• 1.  
  assume a certain systematic error in the computation of the TOAs due, for instance, to instrumental errors such as observation timestamp, reduced BW, hidden RFI, etc.;
• 2.  
  define a criterion to discard the outliers, for instance by vetting residuals above a certain value.

Figure 4 shows the timing residuals of the observations taken with each antenna. The values of the ${\chi }_{\mathrm{red}}^{2}$ from the fits are greater than 1, indicating the presence of outliers or underestimated errors.

Zoom In Zoom Out Reset image size

To account for possible systematic errors, we adopt a simplified approach¹⁴ in which we add quadratically a common σ_sys to all the σ_TOA, producing a total error

$$\begin{eqnarray}&&{\sigma }_{\mathrm{tot}}^{2}={\sigma }_{\mathrm{TOA}}^{2}+{\sigma }_{\mathrm{sys}}^{2}.\end{eqnarray} \tag{ 4 }$$

We calculate the value of σ_sys that leads to ${\chi }_{\mathrm{red}}^{2}=1$, obtaining σ_sys ∼ 0.67 μs for the observations with A1 and σ_sys ∼ 1.0 μs for the observations with A2. We recompute the rms using the corrected errors in the residuals by adding σ_sys as in Equation (4). We obtain rms = 0.72 μs for A1 and rms = 1.05 μs for A2.

In order to determine the effect of the outliers measurements we set a 3σ criterion, but since the σ_tot itself depends on the assumed value of σ_sys, we apply the following iterative process.

• 1.  
  Given an initial ${\sigma }_{\mathrm{sys}}^{(i)}$ (as obtained previously), to each TOA we assign an error ${\sigma }_{\mathrm{tot}}^{(i)\ 2}={\sigma }_{\mathrm{TOA}}^{2}+{\sigma }_{\mathrm{sys}}^{(i)\ 2}$.
• 2.  
  If the residual of an observation is such that $\left|\delta t\right|\gt 3{\sigma }_{\mathrm{tot}}^{(i)}$, then this observation is discarded as an outlier.
• 3.  
  If the residual is such that $\left|\delta t\right|\leqslant 3{\sigma }_{\mathrm{tot}}^{(i)}$, then we keep this observation and its TOA error is given the new value

  $$\begin{eqnarray}&&{{\sigma }_{\mathrm{tot}}^{({\rm{i}}+1)}}^{2}={{\sigma }_{\mathrm{TOA}}}^{2}+{{\sigma }_{\mathrm{sys}}^{(i+1)}}^{2},\end{eqnarray} \tag{ 5 }$$

  where ${\sigma }_{\mathrm{sys}}^{(i+1)}$ is chosen such that when the new residuals are computed we get ${\chi }_{\mathrm{red}}^{2}=1$. In practice, the process converges after one or two iterations.

In this way, we eliminate all the outliers in our data set (five observations for A1 and 24 for A2) and obtain refined values of the systematic errors σ_sys ∼ 0.50 μs for A1, σ_sys ∼ 0.66 μs for A2, and σ_sys ∼ 0.59 μs for A1+A2.

We study the timing residuals for each S/N subset for each antenna; these are shown in Figure A2. By filtering out the low-S/N observations, those with large residuals are eliminated. Thus we conclude that outliers tend to have low S/N; we note, however, that some low-S/N observations also have small residuals.

We perform a timing analysis for A1, A2, and A1+A2. In all cases—even for large S/N values—we obtain ${\chi }_{\mathrm{red}}^{2}\gg 1$. We interpret this as indicative of unaccounted-for systematic errors and we perform the procedure detailed in Section 4.1 to find the values of σ_sys that lead to ${\chi }_{\mathrm{red}}^{2}\approx 1$. Taking as a reference the case for S/N > 50, we obtain σ_sys = 0.5 μs for A1, 0.66 μs for A2, and 0.59 μs for A1+A2. We note that these values change if we do not remove the 3σ outliers, leading to σ_sys = 0.67 μs for A1, 0.99 μs for A2, and 0.83 μs for A1+A2.

In Figure 5 we display the values of σ_sys and rms for each subset of observations with their 1σ error bars (∼68% confidence limits). The error bars for σ_sys are computed as the values ${\sigma }_{\mathrm{sys},\min }$ that yield ${\chi }_{\mathrm{red}}^{2}({n}_{\mathrm{free}},\alpha /2)$ and ${\sigma }_{\mathrm{sys},\max }$ that yield ${\chi }_{\mathrm{red}}^{2}({n}_{\mathrm{free}},1-\alpha /2)$, with α = 0.32.

Zoom In Zoom Out Reset image size

The timing rms diminishes (i.e., improves) for higher S/N observations. The value of the rms is well-constrained to ≳0.5 μs, though values a little higher (≈0.7 μs) are obtained for A2 when low-S/N ( < 100) observations are included. In particular, for S/N > 140 we get rms ≈ 0.52 μs for A1 and 0.55 μs for A2, which is a slight improvement over those reported in Gancio et al. (2020) (0.55 μs for A1 and 0.81 μs for A2).

We also obtain a consistent value of σ_sys ≈ 0.5 μ s. There is a systematic trend of lower σ_sys toward increasing S/N (Figure 5), though with a small significance (close to or below the 1σ level). We conclude that the systematic errors of both IAR's antennas are of the order of 0.4–0.6 μs when accounting for outliers and S/N effects.

Finally, the values of σ_sys and rms are smaller for A1 than for A2 for each subset of ${\rm{S}}/{{\rm{N}}}_{\min };$ this behavior subsists at the same 〈S/N〉 (Figure 5). Given that the main differences between the two antennas are BW and n_P, we explore those dependences in detail to understand the reason(s) behind the improved timing precision of A1.

In Section 4.2 we found that the A1 observations have a lower timing rms than the corresponding observations from A2. Since both antennas differ in their BW (112 MHz for A1 and 56 MHz for A2), and n_P (1 and 2 for A1 and A2, respectively), we reduce the observations to the same BW and n_P in order to quantify the effect of those hardware differences on errors. For this analysis we use the six subsets of observations defined by their S/N in Section 3.1. We split the A1 observations into two subintervals of BW = 56 MHz using pat for the scrunching with the options -j "F {n}", where n = 2 is the number of subintervals. For the observations with A2 we would like to split the two polarizations separately; however, this is not possible as these observations only store the sum of both polarization modes. From the radiometer equation (Lorimer & Kramer 2012)

$$\begin{eqnarray}&&{\sigma }_{\mathrm{sys}}\propto \displaystyle \frac{{T}_{\mathrm{sys}}}{\sqrt{{n}_{{\rm{P}}}\,\mathrm{BW}}},\end{eqnarray} \tag{ 6 }$$

we see that the errors scale with ${n}_{{\rm{P}}}^{-1/2};$ hence, we multiply the errors in the A2 residuals by a factor $\sqrt{2}$ to simulate a case with n_P = 1 (assuming we are not strongly affected by the polarization of the source).

In this way, for each subset of S/N_min we have five groups of observations: three from A1 (one with BW = 112 MHz and two reduced to 56 MHz), and two from A2 (with errors modeled to two and one polarization modes). We model S/N(σ_TOA) and then compute 〈S/N〉 for each of these subset from the σ_TOA of their observations.

The resulting rms values are plotted in Figure 6. For a given value of S/N_min, the higher-frequency sub-band of the A1 observations has lower rms values than the lower-frequency sub-band, which in turn is similar to the case of the A2 observations in one polarization. The inclusion of all the BWs for A1 or both polarizations of A2 shows consistently lower rms. These results can be interpreted as due to:

• 1.  
  RFI affecting more the lower-frequency sub-band,
• 2.  
  effects of differential scintillation,
• 3.  
  dispersion effects being better modeled at higher frequencies.

Zoom In Zoom Out Reset image size

In conclusion, we found that the main difference in the timing errors between the antennas can be attributed to the difference in BW, followed by n_P and an increase of S/N in the selection of the observations. An increase in BW seems paramount for improving the timing errors.

The rms values improve both with longer observation times and with a higher number of TOAs (Lorimer & Kramer 2012; Wang 2015). Here we investigate whether it is possible to improve the overall timing by splitting the long (>200 minute) observations into multiple subintegrations, producing various TOAs from each observation. In this way, we obtain additional data points at the expense of lower timing precision in each of them.

We start with a set of 268 unsplit observations with t_obs > 75 minutes and σ_TOA < 1.0 μs. First we calculate their rms (dashed line in Figure 7). Then we systematically split these observations for different values of the minimum duration of the subintervals considered, from 10 to 75 minutes (the shorter the subinterval, the larger the number of TOAs obtained). We plot the rms as a function of t_min in Figure 7, and specify the total number of points obtained from splitting in each case. We see that the rms diminishes monotonously as t_min increases, showing that this method is not suitable for improving the timing of our observations. This is most likely a sign of the S/N being a major factor affecting our current timing precision; for ${t}_{\min }\lt 70\,\mathrm{minutes}$ it is also possible that jitter affects the TOAs.

Zoom In Zoom Out Reset image size

As described by Alam et al. (2021), the white noise is modeled using three parameters.

• 1.  
  EQUAD accounts for sources of uncorrelated and systematic (Gaussian) white noise in addition to the template-fitting error in the TOA calculations.
• 2.  
  EFAC is a dimensionless constant multiplier to the TOA uncertainty from template-fitting errors. It accounts for possible systematics that lead to underestimated uncertainties in the TOAs.
• 3.  
  ECORR describes short-timescale noise processes that have no correlation between observing epoch, but are completely correlated between TOAs that were obtained simultaneously at different observing frequencies. This parameter accounts for wide-band noise processes such as pulse jitter (Osłowski et al. 2011; Shannon et al. 2014).

Considering σ_TOA to be the template-fitting error of a given observation, the resulting white-noise model is modeled by the noise covariance matrix (Lentati et al. 2014)

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{\nu \nu ^{\prime} ,{tt}^{\prime} }^{2} & = & {\delta }_{{tt}^{\prime} }[{\delta }_{\nu \nu ^{\prime} }\left({\mathrm{EFAC}}^{2}\,{\sigma }_{\mathrm{TOA}}^{2}+{\mathrm{EQUAD}}^{2}\right)\\ & & +\,{\mathrm{ECORR}}^{2}],\end{array}\end{eqnarray} \tag{ 7 }$$

where t and ν are the time and frequency of the observation, respectively.

Given that we have multiple TOAs per day (see Section 2), we need to consider an ECORR contribution. We then incorporate all these noise components and timing-model parameters (as specified in Appendix A.1) into a joint likelihood using ENTERPRISE. We sample the posterior distribution using the sampler PTMCMCSampler (Ellis & van Haasteren 2017), setting uniform prior distributions.

First, we investigate the consistency between the analysis with ENTERPRISE and the independent analysis we presented in Section 4.2. To this end, we use the same set of observations as in the aforementioned analysis while we fix the value EFAC = 1 and we exclude the ECORR parameter from our analysis, so that the Gaussian white noise EQUAD becomes equivalent to the parameter σ_sys in Equation (4). We obtain a notorious agreement between the values of EQUAD and σ_sys: when removing 3σ outliers (Section 4.2) we obtain EQUAD ≈ 0.57 μs, fully consistent with the systematic error of σ_sys ≈ 0.59 μs that we found in Section 4.2 for A1+A2 and S/N > 50;  without removing the outliers, the results are EQUAD ≈ 0.80 μs and σ_sys ≈ 0.83 μs, which again are fully consistent.

Second, we repeat the previous analysis, now taking both EFAC and EQUAD as free parameters. By doing so, we obtain $\mathrm{EFAC}={2.48}_{-0.30}^{+0.29}$ and ${\mathrm{log}}_{10}\mathrm{EQUAD}=-{6.30}_{-0.07}^{+0.10}$ (EQUAD ≈ 0.5 μ s) as the best-fit parameters. The quoted error bars correspond to the 1σ ( ≈ 68%) confidence limits that were obtained using the lower-level function corner.quantile from the corner.py Python module (Foreman-Mackey 2016) and taking the 16th and 84th percentiles. We present a corner plot for these parameters in Figure 8. In this plot we also show confidence intervals considering that the relevant 1σ contour level for a 2D histogram of samples is 1 − e^−0.5 ∼ 0.393 (39.3%). Values of EFAC ∼ 1 would suggest that observing and timing procedures result in near-true TOA uncertainty estimates; thus, the adjusted values of EFAC ∼ 2.5 indicate that the TOAs error bars are considerably underestimated.

Zoom In Zoom Out Reset image size

The red noise is assumed to be a stationary Gaussian process, which is parameterized with a power-law model in frequency, such that the spectral power density is given by (see Hazboun et al. 2020)

$$\begin{eqnarray}&&P(f)=\displaystyle \frac{{A}_{\mathrm{rn}}^{2}}{12{\pi }^{2}}{\left(\displaystyle \frac{f}{{f}_{\mathrm{ref}}}\right)}^{{{\rm{\Gamma }}}_{\mathrm{rn}}}\,{\mathrm{yr}}^{3},\end{eqnarray} \tag{ 8 }$$

where f is a given Fourier frequency in the power spectrum, f_ref is a reference frequency (in this case, 1 yr⁻¹), A_rn is the amplitude of the red noise at the frequency f_ref, and Γ_rn is the spectral index. We take a prior on the red-noise amplitude that is uniform on ${\mathrm{log}}_{10}({A}_{\mathrm{rn}}\,[{\mathrm{yr}}^{3/2}])\in [-14.5,-12]$, and a prior on the red-noise index that is uniform on Γ_rn ∈ [0, 2.6]. The spectrum is evaluated at 30 linearly spaced frequencies f ∈ [1/T_span, 30/T_span], where T_span is the span of the pulsar's data set (in this case, 1.1 yr).

In the following analysis we use a total of 319 observations obtained with A1 and A2 between 2019 April and 2020 June that meet the criteria t_obs > 40 minutes, S/N > 40, and σ_TOA < 1 μs. Details of this data set are summarized in Table 2.

We analyze the data sets of each antenna both independently and altogether. As described in Section 5.1, ECORR accounts for noise that is correlated between observations that were obtained simultaneously at different frequencies. Since such observations are not available for a single antenna, we exclude the ECORR parameter from the analysis of the individual data sets. However, we do include this parameter when analyzing the A1+A2 data set in order to profit from the simultaneous observations at different frequencies. A corner plot for these parameters and their errors is shown in Figure 9.

Zoom In Zoom Out Reset image size

The fitted values to the white- and red-noise parameters for the different data sets are presented in Table 4. To complement this we explore the possibility of splitting long-duration observations into two subintegrations of ${t}_{\min }=75\,\mathrm{minute}$ in order to sample shorter timing frequencies. The adjusted values for this case, also presented in Table 4, are consistent within 1σ to those obtained without the splitting. We therefore conclude that splitting long observations does not improve the timing analysis, in line with the conclusion from Section 4.4.

Table 4.  Adjusted Values for the White- and Red-noise Parameters

	EFAC	${\mathrm{log}}_{10}\mathrm{EQUAD}$	Γrn	${\mathrm{log}}_{10}{A}_{\mathrm{rn}}$
A1	${2.43}_{-0.23}^{+0.25}$	$-{6.3}_{-0.06}^{+0.06}$	${1.22}_{-0.48}^{+0.13}$	$-{13.88}_{-0.44}^{+0.40}$
A2	${2.82}_{-0.30}^{+0.32}$	$-{6.34}_{-0.09}^{+0.08}$	${1.02}_{-0.42}^{+0.36}$	$-{13.51}_{-0.26}^{+0.20}$
A1+A2	${2.48}_{-0.24}^{+0.26}$	$-{6.32}_{-0.07}^{+0.09}$	${0.97}_{-0.38}^{+0.37}$	$-{13.63}_{-0.27}^{+0.21}$
A1+A2*	${2.76}_{-0.19}^{+0.24}$	$-{6.47}_{-0.14}^{+0.17}$	${0.80}_{-0.37}^{+0.39}$	$-{13.54}_{-0.29}^{+0.43}$

Note. Values marked with (*) were obtained by splitting the observations as described in Section 4.4.

Download table as:  ASCIITypeset image

We obtain EQUAD ≈ 0.5 μs in all cases. The value of Γ_rn is less constrained and consistent within 0.5–1.5, while the amplitude is A_rn ≈ (2–7) × 10⁻¹³. The obtained A_rn lies within the expected order of magnitude, whereas Γ_rn falls below the expected value by at least a factor two (Wang 2015), which we interpret as due to the relatively short baseline of our current data set.

We now embark on setting the first bounds to the GW amplitude from massive binary black holes using observations from IAR. In doing so, we aim to exploit the high cadence of these observations.

5.3.1. GW Analysis: Stochastic Background

The contribution of the GWB coming from an ensemble of supermassive black hole (SMBH) binaries or primordial fluctuations during the big bang is modeled similarly to that of the red noise (Equation (8)). Any GWB component is modeled as a single stationary Gaussian process with a power-law timing-residual spectral density

$$\begin{eqnarray}&&P(f)=\displaystyle \frac{{A}_{\mathrm{gwb}}^{2}}{12{\pi }^{2}}{\left(\displaystyle \frac{f}{{f}_{\mathrm{ref}}}\right)}^{{{\rm{\Gamma }}}_{\mathrm{gwb}}}\,{\mathrm{yr}}^{3}.\end{eqnarray} \tag{ 9 }$$

The analysis is nearly identical to the red-noise analysis described in Section 5.2. The prior on the GWB amplitude is taken uniform on ${\mathrm{log}}_{10}({A}_{\mathrm{gwb}}\,[{\mathrm{yr}}^{3/2}])\in [-14.4,-11]$, whereas the prior on the GWB index is uniform on Γ_gwb ∈ [0, 3.2]. Moreover, we fix EFAC and EQUAD to the values adjusted in Section 5.2 for each data set.

In this analysis we also consider both the original data sets and the data sets obtained by splitting the observations in subintegrations with t_obs ≥ 75 minutes (see Section 4.4). The best-fit values to each GWB parameter and for each set of observations are presented in Table 5.

Table 5.  Best-fit Values to the GWB Parameters

Parameter	A1		A2		A1+A2
	No split	Split	No split	Split	No split	Split
Γgwb	${0.50}_{-0.26}^{+0.25}$	${0.38}_{-0.21}^{+0.20}$	${0.12}_{-0.04}^{+0.04}$	${0.10}_{-0.03}^{+0.03}$	${0.38}_{-0.29}^{+0.28}$	${0.28}_{-0.18}^{+0.17}$
${\mathrm{log}}_{10}{A}_{\mathrm{gwb}}$	$-{13.48}_{-0.23}^{+0.25}$	$-{13.37}_{-0.20}^{+0.20}$	$-{13.33}_{-0.21}^{+0.23}$	$-{13.22}_{-0.17}^{+0.18}$	$-{13.48}_{-0.23}^{+0.24}$	$-{13.41}_{-0.18}^{+0.18}$

Download table as:  ASCIITypeset image

Using the split observations, we get a higher cadence at the cost of worsening the S/N (and therefore the TOA precision) of each data point. Our results show consistent values of A_gw ≈ (3 ± 2) × 10⁻¹⁴ and Γ_gw ≈ 0.3 ± 0.2 for all data sets, both with and without the splitting. In general, splitting the observations leads to slightly lower values of Γ_gw and slightly higher values of A_gw, though these differences are not significant as they are within 1σ of the values obtained without the splitting.

While the amplitude we find is consistent with expected bounds for the stochastic background, Γ_gw falls short from the expected 13/3 for a stochastic GW background of SMBH binaries (Siemens et al. 2013), possibly due to our current relatively short observational baseline of ≈1.1 yr.

In order to account for a background of SMBH binaries, we repeat this analysis including a red-noise model with a uniform prior on the spectral index Γ_rn ∈ [0, 7] and an extra red-noise process with Γ_gwb set to 4.33. We also fix all of the white-noise parameters to the values obtained in Section 5.1. In Figure 10 we show a corner plot of the fit to the joint A1+A2 data sets. We find values of A_rn ≈ (4 ± 3) × 10⁻¹⁴, consistent with our previous results, and Γ_rn ≈ 3.81 ± 2.1. Such uncertainties may be attributed to the short time span.

Zoom In Zoom Out Reset image size

5.3.2. GW Analysis: Continuous Source

A single SMBH binary system produces continuous GWs because the system does not evolve notably over the few years of a pulsar-timing data set. We used the Python package Hasasia (Hazboun et al. 2019a) to calculate the single-pulsar sensitivity curve of our data set of J0437−4715 for detecting a deterministic GW source averaged over its initial phase, inclination, and sky location. The dimensionless characteristic strain is calculated for each sampled frequency as

$$\begin{eqnarray}&&{h}_{c}(f)=\sqrt{f\,S(f)}\end{eqnarray} \tag{ 10 }$$

where S is the strain-noise power spectral density for the pulsar. This is related to the power spectrum of the induced timing residuals of Equation (9) by (see Jenet et al. 2006)

$$\begin{eqnarray}&&P(f)=\displaystyle \frac{1}{12{\pi }^{2}}\displaystyle \frac{1}{{f}^{3}}{h}_{c}{\left(f\right)}^{2}.\end{eqnarray} \tag{ 11 }$$

The white- and red-noise parameters that were adjusted in Sections 5.1 and 5.2 using ENTERPRISE are loaded into the package in order to account for these effects in the calculations. Since our observations have a time baseline of T_obs = 1.1 yr and a nearly daily cadence, we calculate the curve across a frequency range between 1/(10 T_obs) ∼ 2.8 × 10⁻⁹ Hz and 1/(1 day) ∼ 1.2 × 10⁻⁵ Hz.

The resulting sensitivity curve is shown in Figure 11. It is readily seen that there is a loss of sensitivity at a frequency of (1 yr)⁻¹, caused by fitting the pulsar's position, and at a frequency of (PB)⁻¹ ∼ 2μHz (with PB the orbital period), caused by fitting the orbital parameters of the binary system. The additional spikes seen at frequencies higher than (PB)⁻¹ correspond to harmonics of the binary orbital frequency.

Zoom In Zoom Out Reset image size

In addition, the sensitivity at lower frequencies is reduced by (i) the fit of a quadratic polynomial to the TOAs required to model the pulsar spin-down and (ii) the fitting of "jumps" to connect the timing residuals obtained with different backends (Yardley et al. 2010). The frequency dependence ( ∼ f^−3/2) at low frequencies is evidence of a fit to a quadratic spin-down model for the pulsar spin frequency. As a result, the minimum of the sensitivity curve should be attained at a frequency of 1/T_obs. However, given that the T_obs of our data set is close to 1 yr, this feature coincides with the loss of sensitivity at (1 yr)⁻¹. We expect to obtain a well-defined minimum at  ≈ 1/T_obs in future by accumulating more observations and achieving a significantly longer time baseline.

For completeness, we tested the significance of the red-noise contribution by calculating a sensitivity curve without this component. The curve was essentially insensitive to those changes in the priors. This is expected, since the injection of red noise should lead to a flat sensitivity curve around the minimum (Hazboun et al. 2019b), though in our case it is coincident with the spike at (1 yr)⁻¹.

For comparison, we used ENTERPRISE to perform a fixed-frequency Markov chain Monte Carlo procedure at four different frequencies. We obtained a posterior distribution for ${\mathrm{log}}_{10}{h}_{\mathrm{gw}}$ at each of these frequencies with a mean value in great agreement with the curve obtained with Hasasia, as shown in Figure 11.

These first results on GW sensitivity are encouraging, though we still need to achieve a sensitivity of at least a factor 10 higher in order to observe even the most favorable SMBH binary merger events. For instance, the six billion solar mass source of 3C 186 at z ≈ 1 produced a GW of h ∼ 10⁻¹⁴ at the time of arrival to our Galaxy, roughly a million years ago (Lousto et al. 2017).

To de-disperse and fold the observations we use the software PRESTO (Ransom et al. 2003; Ransom 2011). It has a variety of tools for the reduction of observations. The processed data are s stored in a .pfd file that contains the pulse profile for different time and frequency bins. In addition to this profile, PRESTO outputs a .polycos file that contains the coefficients of a polynomial modeling the variation of the pulsar period. These coefficients allow us to determine the period of pulsation in a topocentric reference system and are necessary to compute the timing residuals.

If the observation is in the file obs.fil and the mask in the file m.mask, then the command-line used has the following syntax:

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{c}{\mathtt{prepfold}}\,\text{-}{\mathtt{nsub}}\,{\mathtt{64}}\,\text{-}{\mathtt{n}}\,{\mathtt{1024}}\,\text{-}{\mathtt{timing}}\\ \,\,\,\,{\mathtt{J}}{\mathtt{0437}}\text{-}{\mathtt{4715.}}{\mathtt{par}}\,\text{-}{\mathtt{mask}}\,{\mathtt{m}}.{\mathtt{mask}}\,{\mathtt{obs}}.{\mathtt{fil}}\end{array}\end{array}\end{eqnarray} \tag{ A1 }$$

where the option -timing indicates prepfold to generate a file .polycos based on the pulsar parameters This process is currently automatized through local Python scripts.

Considering that A1 and A2 have different configurations (number of polarizations and BW; see Table 1), it is possible that slight differences arise in the integrated profile seen by each antenna. We therefore study whether the template used has a significant impact on the timing residuals.

To create each template we choose observations with n_bins = 1024 phase bins and n_chan = 64 frequency channels. We select for each antenna data the highest-S/N observation and extract the noise using the task psrsmooth in the package psrchive. This choice of templates seems adequate since J0437−4715 is a very bright pulsar and selected individual observations produce a high enough S/N to create a template. We highlight that the large span of our observations (over 3 hr in many cases) mitigates the impact of the intrinsic jitter of the pulsar (Liu et al. 2012). The selected templates for each antenna correspond to the profiles with n_bins = 1024 phase bins in Figure A1. The relative error between them is below 5% near the peak, with larger relative differences toward the wings, but those do not have a major influence in the determination of the TOAs.

Zoom In Zoom Out Reset image size

In the preliminary timing analysis of J0437−4715 presented in Gancio et al. (2020), we used the same template on both antennas to determine TOAs. We show in our separated analysis of A1 and A2 data that this assumption was valid to the current level accuracy, producing an rms = 0.8 μs residual for A2 observations with the use of either template. Notwithstanding this a posteriori verification, we consistently use different templates for A1 and A2 throughout this work.

In order to study the effect of the number of phase bins (n_bins) used in the folding of observations on the timing residuals, we have taken data folded originally with n_bins = 1024, and processed with the routine bscrunch of the psrchive package for Python, to generated copies of the observations and their corresponding templates for each antenna, but with values of n_bins = 512, 256, 128, 64, and 32. Through this process of scrunching we obtained six sets of observations for each antenna only differing by their n_bins. Figure A1 shows the effect of the n_bins on the templates for each antenna. While for n_bins = 32 we lose temporal resolution, the differences beyond n_bins ≥ 256 are almost negligible to our precision.

Next we compute the timing residuals for each n_bins subset. Interestingly, only for n_bins ≤ 64 are the timing errors too large; for n_bins ≥ 128 the derived TOAs are very consistent, the size of the error bar being the main difference (with smaller error bars obtained for larger n_bins). The rms of the residuals for each subset after adjusting σ_sys as a function of n_bins is shown in Figure A1. The rms decreases with increasing n_bins significantly for 64 to 256 bins, showing that we cannot attain good timing for n_bins ≤ 64 and need at least 256 bins to obtain a precision higher than 1 μs for a pulsar like J0437−4715. This corresponds to a time interval much smaller than the full width at half-maximum of the pulse, that is, 0.3 μs at 1400 MHz.

Here we present additional figures that support the hypothesis that our timing studies are limited due to the S/N of the observations. This effect has a larger impact for A2, as can be seen in Figure A2.

Zoom In Zoom Out Reset image size