PSR J0437-4715: The Argentine Institute of Radioastronomy 2019-2020 Observational Campaign
V. Sosa Fiscella, S. del Palacio, L. Combi, C.O. Lousto, J.A. Combi, G. Gancio, F. García, E. Gutiérrez, F. Hauscarriaga, P. Kornecki, F.G. López Armengo, G.C. Mancuso, A.L. Müller, A. Simaz Bunzel
DDraft version November 30, 2020
Typeset using L A TEX twocolumn style in AASTeX63
PSR J0437 − V. Sosa Fiscella,
1, 2
S. del Palacio, L. Combi,
2, 3
C. O. Lousto, J. A. Combi,
2, 1
G. Gancio, F. Garc´ıa,
2, 4
E. Guti´errez, F. Hauscarriaga, P. Kornecki, F. G. L´opez Armengol, G. C. Mancuso, A. L. M¨uller, andA. Simaz Bunzel Facultad de Ciencias Astron´omicas y Geof´ısicas, Universidad Nacional de La Plata, Paseo del Bosque, B1900FWA La Plata, Argentina Instituto Argentino de Radioastronom´ıa (CCT La Plata, CONICET; CICPBA; UNLP),C.C.5, (1894) Villa Elisa, Buenos Aires, Argentina. Center for Computational Relativity and Gravitation, School of Mathematical Sciences, Rochester Institute of Technology, 85 LombMemorial Drive, Rochester,New York 14623, US Kapteyn Astronomical Institute, University of Groningen, P.O. BOX 800, 9700 AV Groningen, the Netherlands (Received October 5, 2020; Revised November 25, 2020; Accepted November 30, 2020)
Submitted to ApJABSTRACTWe present the first year data set of high-cadence, long-duration observations of the bright millisecondpulsar J0437 − µ s. We characterize the white and red timing noise in IAR’s observations, we quantify theeffects of scintillation, and we perform single pulsar searches of continuous gravitational waves, settingconstraints in the nHz– µ Hz frequency range. We demonstrate IAR’s potential for performing pulsarmonitoring in the 1.4 GHz radio band for long periods of time with a daily cadence. In particular, weconclude that the ongoing observational campaign of the millisecond pulsar J0437 − Keywords: editorials, notices — miscellaneous — catalogs — surveys INTRODUCTIONThe Argentine Institute of Radio astronomy (IAR)is equipped with two single-dish 30-meter antennas –dubbed A1 and A2 – capable of performing daily ob-servations of pulsars in the southern hemisphere at 1.4GHz. These antennas were recently refurbished to ob-tain high-quality timing observations as described inGancio et al. (2020).Pulsar Monitoring in Argentina (PuMA), is a scien-tific collaboration dedicated to pulsar observations fromthe southern hemisphere. As part of IAR’s observatorydeveloping stage, accurate timing observations of the Corresponding author: C. O. [email protected] In 2019 the antennas A1 and A2 were renamed “Varsavsky” and“Bajaja”, respectively. http://puma.iar.unlp.edu.ar millisecond pulsar (MSP) J0437 − − S = 150 . d = 156 . ± .
25 pc) pulsars. It has a short period( P = 5 .
758 ms) and it is one of the most massive pul-sars known to date ( m = 1 . ± .
07 M (cid:12) ; Reardon 2018).This pulsar is in a binary system and in an almost cir-cular orbit of period 5.74 days. The secondary star is alow mass ( ∼ . (cid:12) ), helium white dwarf, with strongvisible emission (Danziger et al. 1993). In the interstel-lar region, an optical bow shock was also reported byBell et al. (1993). In addition, J0437 − a r X i v : . [ a s t r o - ph . GA ] N ov Sosa Fiscella et al. although in this wavelength its spectrum is consistentwith that of a black body (Lorimer & Kramer 2012) andpulsed emission was not seen (Kargaltsev et al. 2004).Because of its proximity to Earth, J0437 − − − − − − (cid:46) µ s.The paper is organized as follows: Sect. 2 introducesthe observations and the reduction methods. In Sect. 3we describe the observations in terms of their signal-to-noise ratio (S/N) and its relation to interstellar scintil-lation. In Sect. 4 we present the timing results and westudy the influence of the S/N and bandwidth ( BW )on the timing analysis; further details on this analy- sis are provided in the Appendix A. In Sect. 5 we usethe ENTERPRISE software to perform a noise analysisof the observations and estimate the contribution of agravitational-wave background. Finally, in Sect. 6 wepresent the main conclusions of our study. OBSERVATIONSAs described in more detail in Gancio et al. (2020),the design of the antennas allows to observe a sourcecontinuously during 220 min. Their receivers are notcurrently refrigerated and have a system temperature of T sys ∼
100 K. The back-end is based on two SDRs whichacquire raw samples with a maximum rate of 56 MHz perboard. A1 uses these two digital plates in consecutiveradio frequencies with a total bandwidth of 112 MHz ina single polarization mode; while A2 uses those digitalplates in one per polarization, thus covering a band-width of 56 MHz. Those characteristics are summarizedin Table 1. In November 2019, A1’s receiver front-endwent into commissioning. The electronics and systemswere verified and improved, resulting in a slightly highersensitivity and the recovery of the second polarization.Nonetheless, the observations were retaken with the pre-vious configuration to have a homogeneous data set.
A1 A2Number of observations 170 (145*) 197 (171*)MJD start – MJD finish 58596.7 – 58999.6Total observation time [h] 391 (372*) 393 (381*)Central frequency [MHz] 1400, 1415, 1428 1428Bandwidth ( BW ) 112 MHz 56 MHzPolarization modes 1 2Frequency channels ( n chan ) 64/128 64Time resolution [ µ s] 73.14Phase bins ( n bin ) 512/1024 Table 1.
Parameters of the observations analyzed in thiswork. Values marked with (*) correspond to the restricteddata set used in Sect. 5 (observations lasting more than 40min that achieve a S / N >
40 and σ TOA < µ s). In this work we present the analysis of a data set of 170observations with A1 and 197 with A2 over an intervalof 13 months, from April 23rd 2019 to May 30th 2020.This includes days with multiple observations (89 dayswith 2 observations, 24 days with 3, and 1 day with 4).The observations add up to over 390 hs of observationwith each antenna (Table 1), leading to an observationefficacy of 0.26 for both antennas. This efficacy is aimedto be improved in a future considering that (i) A1 under-went maintenance between October 8th and November29th 2019, (ii) an unusually loud source of local radio
AR timing analysis of PSR J0437 − rfiClean (Maan et al. 2020, in prep.)gave better results than the rfifind task in PRESTO toclean RFIs. We therefore ran both softwares in all A1observations carried out from November 2019 onwards.The observations, stored in filterbank format , werefolded and de-dispersed with PRESTO (Ransom et al.2003; Ransom 2011) using n bins = 512 or 1024 phasebins and n chan = 64 frequency channels for A2 ob-servations and n chan = 64 or 128 for A1 observations.The data were folded using the timing flag of the task prepfold and the parameter ( .par ) file provided byIPTA , “Combination B” with edits adapted to the IARsite. We then calculated the time of arrival (TOA) ofthe pulses using the pat package in PSRCHIVE (Hotanet al. 2004) with a Fourier Phase Gradient (PGS) match-ing template fitting (Taylor 1992). The template wasobtained applying a smoothing wavelet algorithm to abest profile; a more detailed discussion of the templateselection is provided in Appendix A.2. The TOAs inthis data set were fixed of clock systematics on April22nd 2019 (MJD 58595), when we reached an accuracyof < µ s (we refer to Gancio et al. 2020, for details onclock settings). ANALYSIS OF THE OBSERVATIONS3.1.
Signal-to-noise ratio of the observations
In order to characterize the S/N of the observations weuse the functions getDuration and getSN of the Pythonpackage
PyPulse (Lam 2017). In Fig. 1 we show theS/N of each observation as a function of their duration.The mean S/N of observations with A1 is 151 and withA2 is 105, with mean observing times of 147 min and116 min, respectively. However, we note that these num-bers are affected by many short and low-quality obser-vations.When we restrict our analysis to observations withS / N >
50, the mean S/N for observations with A1 in- https://github.com/ymaan4/rfiClean http://sigproc.sourceforge.net/ In Appendix A.3 we show that the number of phase bins doesnot affect the posterior analysis as long as n bins ≥ http://ipta4gw.org//data-release/ https://github.com/mtlam/PyPulse
10 50 100 150 200 t obs [minutes] S / N A1A2
Figure 1.
S/N of the observations of each antenna as afunction of their t obs . We also plot f ( t obs ) = a √ t obs , where a = 13 . − / for A1 and a = 11 . − / for A2. creases to 166 and with A2 to 122, with mean observingtimes of 162 min and 124 min, respectively. We summa-rize these and other values in Table 2.We observe a positive correlation between S/N and t obs , fitting to a S / N ∝ √ t obs as expected (Lorimer &Kramer 2012): S/N = (cid:112) n P t obs BW (cid:18) T peak T sys (cid:19) (cid:112) W ( P − W ) P , (1)where P is the pulsar period and W its width, T peak its maximum amplitude, T sys is the noise temperatureof the system, t obs is the observing time, and n P thenumber of polarizations observed.We collect the observation per S/N for each antennaand display them as histograms in Fig. 2. We observe adistribution for A1 with a mean higher than the corre-sponding distribution for A2, perhaps due to the broaderband sensitivity of A1.We collect the observations into sets of S / N >
1, 50,80, 110, 140 and 170, corresponding to roughly the po-sition of the larger bins in the A2 histogram. In Table 2we specify the number of observations, mean durationand mean S/N for each of these sets.3.2.
Scintillations
In what follows we assume that the expected S/Nscales ∝ √ t obs and that additional variations in the S/Nare due to scintillation. We note that the observationsdescribed in Sec. 2 lack of absolute flux calibrations andthus possible variations in T sys are not accounted for.Moreover, RFIs are also variable and their mitigationleads to variations in the effective bandwidth of each Sosa Fiscella et al.
S/N > >
50 S/N >
80 S/N >
110 S/N >
140 S/N > N (cid:104) S / N (cid:105) (cid:104) t obs (cid:105) N (cid:104) S / N (cid:105) (cid:104) t obs (cid:105) N (cid:104) S / N (cid:105) (cid:104) t obs (cid:105) N (cid:104) S / N (cid:105) (cid:104) t obs (cid:105) N (cid:104) S / N (cid:105) (cid:104) t obs (cid:105) N (cid:104) S / N (cid:105) (cid:104) t obs (cid:105) A1 170 151 147 159 160 155 150 166 160 120 183 166 96 197 180 59 223 187A2 197 105 116 164 120 121 128 136 146 88 153 178 58 168 192 22 192 194A1+A2 367 127 130 323 140 140 278 152 154 208 170 171 154 186 182 81 214 191
Table 2.
Number of observations N , mean S/N, and mean t obs expressed in minutes per S/N subset per antenna. S/N N u m b e r o f o b s e r v a t i o n s A1A2
Figure 2.
Histogram of the observations for each antenna,A1 and A2, according to its S/N. observation, so additional dispersion in the S/N vs. t obs relation is also expected.To quantify the variations due to scintillation we builda projected pulse S/N as S / N proj = S / N (cid:112) t max /t obs ,with t max = 217 min. Given that short observationshave a large uncertainty in their determined S/N, weonly use observations with t obs >
20 min (which isroughly half of the scintillation timescale). Fig. 3 showsa histogram of the projected pulse S/N for A1 and A2.The line shows the estimated probability density func-tion (PDF) from scintillation (Cordes & Chernoff 1997) f S ( S | n ISS ) = ( Sn ISS /S ) n ISS S Γ( n ISS ) exp (cid:18) − Sn ISS S (cid:19) Θ( S ) , (2)where n ISS is the number of scintles, S is the meanvalue of the signal S (i.e., S = (cid:104) S / N (cid:105) ), and Θ isthe Heaviside step function. We calculate n ISS by fit-ting the normalized data for each antenna. We obtain n ISS = 2 . ± .
31 for A1 and n ISS = 2 . ± .
25 for A2,with S = 127 .
27 for A1 and S = 87 .
16 for A2. Thebin size is determined using Knuth’s rule (Knuth 2006)algorithm provided in astropy (Astropy Collaboration We normalize the number of observations in each S/N bin by thetotal of observations of each antenna. et al. 2013, 2018), though we confirm that the obtainedvalues do not depend on the binning by repeating theanalysis for different bin sizes.In addition, we make use of the long duration of theobservations that is significantly larger than the typ-ical scintillation timescale for J0437 − t min = 2000 s and t min = 5000 s and repeat the previous analysis. In thiscase we obtain larger values of n ISS ∼ n ISS in Eq. 2) under the null hy-pothesis that the sample is drawn from the referencedistribution. The null hypothesis can be rejected at agiven confidence level α if the resulting p -value is lowerthan 1 − α . The p -values obtained are summarized inTable 3. For α = 0 . p -value.However, the fits to the splitted observations of A2 failthis test, suggesting that, for short observations withA2, Eq. 2 may not be entirely valid or that the estimateof the projected S/N becomes unreliable. A1 A2 n ISS error p n
ISS error p No split 2.67 0.31 0.38 2.17 0.25 0.24Split t min = 5000 s 6.33 0.54 0.90 5.53 1.04 0.009Split t min = 2000 s 5.50 0.36 0.70 4.63 0.43 0.004 Table 3.
Adjusted values of n ISS for each set of observationsand the KS test p -value for each fitting. We compare our values of n ISS , with theoretical es-timations following Lam & Hazboun (2020). We scalethe scintillation parameters given at the frequency of1.5 GHz by Keith et al. (2013) to match our observa-tions centered at 1.4 GHz and obtain the scintillationbandwidth ∆ ν d = 740 MHz and scintillation timescale AR timing analysis of PSR J0437 − Projected S/N C o un t s S = 127.27 n ISS = 2.67
Projected S/N C o un t s S = 87.16 n ISS = 2.17
Projected S/N C o un t s S = 199.4 n ISS = 5.5
Projected S/N C o un t s S = 146.7 n ISS = 4.63
Figure 3.
Histograms of projected pulse S/N for PSR J0437 − t min = 2000 s. The line showsthe estimated scintillation distribution from fitting n ISS in Eq. 2. ∆ t d = 2290 s. We calculate n ISS via the usual formula n ISS ≈ (cid:18) η t T ∆ t d (cid:19) (cid:18) η ν BW ∆ ν d (cid:19) (3)where η t and η ν are filling factors ∼ .
2. The estimated n ISS for T = 220 min are 2.22 for A1 ( BW = 112 MHz)and 2.18 for A2 ( BW = 56 MHz). We confirm thatthe value obtained with A2 is consistent with the ex-pectations, although for A1 it is larger than expected,perhaps due to additional factors affecting the variabil-ity observed.Gwinn et al. (2006) found two scintillation scales ob-serving J0437 − t d , = 5727 s and ∆ t d , = 515 s, leadingto n ISS , = 1 .
46 for both antennas and n ISS , = 6 . n ISS , = 6 .
35 for A2. The later values areclose to the ones displayed in Table 3 for the split obser-vations, consistent with the shorter observations being more sensitive to the shorter-scale scintillations. Notealso that those scintillations scales have been observedto vary notably between epochs (Smirnova et al. 2006). TIMING ANALYSISHere we discuss the timing-error dependence on threeparameters: (i) the S/N of the observations, (ii) thenumber of bins used in the reduction of the observations,(iii) the BW of the observations. In addition we studyand quantify other sources of systematic errors.4.1. Timing Residuals
We compute the timing residuals of the TOAs us-ing
Tempo2 (Hobbs et al. 2006) and its Python wrap-per, libstempo (Vallisneri 2020), with the timing modelgiven in the file
J0437 − provided by IPTA Sosa Fiscella et al. and adapted to IAR observatory . Tempo2 returns: i)the MJD, residual, and template-fitting error ( σ TOA )of each observation, and ii) the timing model param-eters, the weighted errors of the residuals (RMS), andthe χ = χ /n free of the timing model fit to the resid-uals. The χ test considers a good fit when χ ∼ χ (cid:29) • Assume a certain systematic error in the computa-tion of the TOAs due, for instance, to instrumen-tal errors such as observation timestamp , reduced BW , hidden RFIs, etc. • Define a criteria to discard the outliers, for in-stance by vetting residuals above a certain value.
100 200 300 400 500
MJD - 58500 − − − R e s i du a l s [ µ s ] A1A2
Figure 4.
Timing residuals for the complete data set for A1and A2.
Fig. 4 shows the timing residuals of the observationstaken with each antenna. The values of the χ fromthe fits are greater than 1, indicating the presence ofoutliers or underestimated errors.To account for possible systematic errors, we adopt asimplified approach in which we add quadratically acommon σ sys to all the σ TOA , producing a total error σ = σ + σ . (4) In this .par we also included four
JUMPs to account for the differ-ent central frequencies of the observations, and the correspondingantenna (A1/A2; see Table 1). In Sec. 5 we compare the results of this simplified model withthose obtained using a standard and more refined white noisemodel as in Arzoumanian et al. (2016).
We calculate the value of σ sys that leads to χ = 1,obtaining σ sys ∼ . µ s for the observations with A1and σ sys ∼ . µ s for the observations with A2. Werecompute the RMS using the corrected errors in theresiduals by adding σ sys as in Eq. 4. We obtain RMS =0 . µ s for A1 and RMS = 1 . µ s for A2.In order to determine the effect of the outliers mea-surements we set a 3- σ criteria, but since the σ tot itselfdepends on the assumed value of σ sys , we apply the fol-lowing iterative process:1. Given an initial σ ( i )sys (as obtained previously), toeach TOA we assign an error σ ( i ) 2tot = σ + σ ( i ) 2sys .2. If the residual of an observation is such that | δt | > σ ( i )tot , then this observation is discarded as an out-lier.3. If the residual is such that | δt | ≤ σ ( i )tot , then wekeep this observation and its TOA error is giventhe new value σ (i+1)tot 2 = σ TOA2 + σ ( i +1)sys 2 , (5)where σ ( i +1)sys is chosen such that when the newresiduals are computed we get χ = 1. In prac-tice, the process converges after 1–2 iterations.In this way, we eliminate all the outliers in our dataset (5 observations for A1 and 24 for A2) and obtainrefined values of the systematic errors σ sys ∼ . µ sfor A1, σ sys ∼ . µ s for A2, and σ sys ∼ . µ s forA1+A2. 4.2. Timing versus S/N
We study the timing residuals for each S/N subset foreach antenna; these are shown in Fig. 13. By filteringout the low S/N observations, those with large residualsare eliminated. Thus we conclude that outliers tend tohave low S/N; we note, however, that some low S/Nobservations 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 χ (cid:29)
1. We interpret this as indicative of unaccountedsystematic errors and we perform the procedure detailedin Sec. 4.1 to find the values of σ sys that lead to χ ≈ / N >
50, we obtain σ sys = 0 . µ s for A1, 0 . µ s for A2, and 0 . µ s forA1+A2. We note that these values change if we do notremove the 3- σ outliers, leading to σ sys = 0 . µ s forA1, 0 . µ s for A2, and 0 . µ s for A1+A2.In Fig. 5 we display the values of σ sys and RMS foreach subset of observations with their 1- σ error bars ( ∼
68% confidence limits). The error bars for σ sys are com-puted as the values σ sys , min that yield χ ( n free , α/ AR timing analysis of PSR J0437 − σ sys , max that yield χ ( n free , − α/ α =0 . (cid:38) . µ s, though values a bit higher ( ≈ . µ s) areobtained for A2 when low-S/N ( < / N >
140 we get RMS ≈ . µ s for A1 and 0 . µ s for A2, which is a slightimprovement over those reported in Gancio et al. (2020)(0 . µ s for A1 and 0 . µ s for A2).We also obtain a consistent value of σ sys ≈ . µ s.There is a systematic trend of lower σ sys towards in-creasing S/N (Fig. 5), though with a small significance(close to or below 1- σ level). We conclude that the sys-tematic errors of both IAR’s antennas are of the order of0 . . µ s when accounting for outliers and S/N effects.Finally, the values of σ sys and RMS are smaller forA1 than for A2 for each subset of S / N min ; this behaviorsubsists at the same (cid:104) S / N (cid:105) (Fig. 5). Given that the maindifferences between the two antennas are BW and n P ,we explore those dependencies in detail to understandthe reason(s) behind the improved timing precision ofA1. 4.3. Timing versus bandwidth
In Sec. 4.2 we found that the A1 observations havea lower timing RMS than the corresponding from A2.Since both antennas differ in their BW (112 MHz for A1and 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 thosehardware differences on errors. For this analysis we usethe 6 subsets of observations defined by their S/N inSec. 3.1. We split the A1 observations into two subin-tervals of BW = 56 MHz using pat for the scrunching with the options -j "T { n } " , where n = 2 is the numberof subintervals. For the observations with A2 we wouldlike to split the two polarizations separately; however,this is not possible as these observations only store thesum of both polarization modes. From the radiometerequation (Lorimer & Kramer 2012) σ sys ∝ T sys √ n P BW , (6)we see that the errors scale with n − / ; hence, we multi-ply the errors of A2 residuals by a factor √ n P = 1 (assuming we are not strongly af-fected by the polarization of the source).In this way, for each subset of S / N min we have fivegroups of observations: three from A1 (one with BW =112 MHz and two reduced to 56 MHz), and two from A2(with errors modeled to 2 and 1 polarization modes). We
120 136 153 168 183 192 214 223 h S / N i . . . . . . σ s y s S/N > > > >
140 S/N >
120 136 153 168 183 192 214 223 h S / N i . . . . . . R M S [ µ s ] S/N > > > >
140 S/N > Figure 5.
Top: σ sys for the (cid:104) S / N (cid:105) of each subset of obser-vations for each antenna and their corresponding error bars.Bottom: RMS from recomputed TOAs errors augmented by σ sys and their corresponding error bars. modeled S / N( σ TOA ) and then computed (cid:104) S / N (cid:105) for eachof these subset from the σ TOA of their observations.The resulting RMS are plotted in Fig. 6. For a givenvalue of S / N min , the higher frequency sub-band of theA1 observations have lower RMS than the lower fre-quency sub-band, which in turn are similar to the A2observations in one polarization. The inclusion of allthe BW for A1 or both polarizations of A2 show consis-tently lower RMS. These results can be interpreted asdue to1. RFIs affecting more the lower frequency sub-band.2. Effects of differential scintillation.3. Dispersion effects being better modeled at higherfrequencies.In conclusion, we found that the main difference in thetiming errors between the antennas can be attributed to Sosa Fiscella et al. the difference in BW , followed by n P and an increase ofS/N in the selection of the observations. An increase in BW seems paramount to improve the timing errors.
80 100 120 140 160 180 200 h S / N i . . . . . . R M S [ µ s ] S/N > > > > > >
170 A1 (112MHz)A1 - inf. sub-band (56MHz)A1 - sup. sub-band (56MHz)A2 - 2 polarizationsA2 - 1 polarization
Figure 6.
RMS of the timing residuals with the A1 observa-tions scrunched to BW = 56 MHz and those observed withA2 reduced to one polarization. We also reproduce the fullA1 and A2 original residuals. Timing versus observation length (with split ofobservations)
The RMS values improve both with longer observationtimes and with a higher number of TOAs (Lorimer &Kramer 2012; Wang 2015). Here we investigate whetherits possible to improve the overall timing by splitting thelong ( >
200 min) observations into multiple subintegra-tions, producing various TOAs from each observation.In this way, we obtain additional data points at the ex-pense of lower timing precision in each of them.We start with a set of 268 unsplit observations with t obs >
75 min and σ TOA < . µ s. First we calculatetheir RMS (dashed line in Fig. 7). Then we systemati-cally split these observations for different values of theminimum duration of the subintervals considered, from10 min to 75 min (the shorter the subinterval, the largerthe number of TOAs obtained). We plot the RMS asa function of t min in Fig. 7, and specify the total num-ber of points obtained from splitting in each case. Wesee that the RMS diminishes monotonously as t min in-creases, showing that this method is not suitable forimproving the timing of our observations. This is mostlikely a sign of the S/N being a major factor affecting our current timing precision; for t min <
70 min it is alsopossible that jitter affects the TOAs.
10 15 20 25 30 35 40 45 50 55 60 65 70 75 t min [minutes] . . . . . . . . . . . R M S N =268 N =2776 N =2140 N =1692 N =1380 N =1166 N =1002 N =862 N =768 N =663 N =600 N =557 N =516 N =435 N =412 Figure 7.
Timing obtained when splitting observations.The total number of points after splitting is detailed as N .The horizontal dashed line corresponds to no splitting. Theprojected crossing of curves occurs at about 90 minutes.5. NOISE ANALYSISIn the following sections, we analyze: i) the whitenoise in our data set, which is needed to estimate thesystematic timing errors; ii) the red noise, which is cor-related in time and has a larger amplitude at low fre-quencies; iii) the GWB at µ Hz frequencies, which is pro-duced from a variety of sources that we cannot iden-tify individually. For this purpose, we use the software
ENTERPRISE (Enhanced Numerical Toolbox Enabling aRobust PulsaR Inference SuitE), a pulsar-timing anal-ysis code which performs noise analysis, gravitational-wave (GW) searches, and timing-model analysis (Elliset al. 2019).
ENTERPRISE uses the timing model, previ-ously fit with
Tempo2 , as the basis to construct a designmatrix centered around the timing parameters. Thisis then used to find the maximum-likelihood fit for thewhite- and red-noise parameters.5.1.
White-noise analysis
As described by Alam et al. (2020), the white noise ismodeled using three parameters: • EQUAD accounts for sources of uncorrelated andsystematic (Gaussian) white noise in addition tothe template-fitting error in the TOA calculations. • EFAC is a dimensionless constant multiplier tothe TOA uncertainty from template-fitting errors.It accounts for possible systematics that lead tounderestimated uncertainties in the TOAs.
AR timing analysis of PSR J0437 − • ECORR describes short-timescale noise processesthat have no correlation between observing epoch,but are completely correlated between TOAs thatwere obtained simultaneously at different observ-ing frequencies. This parameter accounts for wide-band noise processes such as pulse jitter (Os(cid:32)lowskiet al. 2011; Shannon et al. 2014).Considering σ TOA the template-fitting error of a givenobservation, the resulting white-noise model is modelledby the noise covariance matrix (Lentati et al. 2014) σ νν (cid:48) ,tt (cid:48) = δ tt (cid:48) [ δ νν (cid:48) (cid:0) EFAC σ + EQUAD (cid:1) + ECORR ] , (7)where t and ν are the time and frequency of the obser-vation, respectively.Given that we have multiple TOAs per day (seeSec. 2), we need to consider an ECORR contribution.We then incorporate all these noise components andtiming-model parameters (as specified in Sec. A.1) into ajoint likelihood using ENTERPRISE . We sample the poste-rior distribution using the sampler
PTMCMCSampler (Ellis& van Haasteren 2017), setting uniform prior distribu-tions.Firstly, we investigate the consistency between theanalysis with
ENTERPRISE and the independent analy-sis we presented in Sect. 4.2. With this end, we usethe same set of observations as in the aforementionedanalysis while we fix the value EFAC = 1 and we ex-clude the ECORR parameter from our analysis, so thatthe Gaussian white noise EQUAD becomes equivalentto the parameter σ sys in Eq. 4. We obtain a noto-rious agreement between the values of EQUAD and σ sys : when removing 3- σ outliers (Sect. 4.2) we obtainEQUAD ≈ . µ s, fully consistent with the systematicerror of σ sys ≈ . µ s that we found in Sect. 4.2 forA1+A2 and S/N >
50; without removing the outliers,the results are EQUAD ≈ . µ s and σ sys ≈ . µ s,which again are fully consistent.Secondly, we repeat the previous analysis now takingboth EFAC and EQUAD as free parameters. By doingso, we obtain EFAC = 2 . +0 . − . and log EQUAD = − . +0 . − . (EQUAD ≈ . µ s) as the best-fit param-eters. The quoted error bars correspond to the 1- σ ( ≈ corner.quantile from the corner.py Python module (Foreman-Mackey 2016) andtaking the 16th and 84th percentiles. We present a cor-ner plot for these parameters in Fig. 8. In this plot wealso show confidence intervals considering that the rele-vant 1- σ contour level for a 2D histogram of samples is1 − e − . ∼ .
393 (39 . ∼ ∼ . . . . . EFAC − . − . − . − . l og E Q UA D − . − . − . − . log EQUAD
Figure 8.
White noise
ENTERPRISE timing analysis forJ0437 − Red-noise analysis
The red noise is assumed to be a stationary Gaus-sian process, which is parameterized with a power-lawmodel in frequency, such that the spectral power densityis given by (see Hazboun et al. 2020) P ( f ) = A π (cid:18) ff ref (cid:19) Γ rn yr , (8)where f is a given Fourier frequency in the powerspectrum, f ref is a reference frequency (in this case,1 yr − ), A rn is the amplitude of the red noise at thefrequency f ref , and Γ rn is the spectral index. We takea prior on the red noise amplitude that is uniform onlog ( A rn [yr / ]) ∈ [ − . , − rn ∈ [0 , . f ∈ [1 /T span , /T span ], where T span is the span of thepulsar’s data set (in this case, 1.1 yr).In the following analysis we use a total of 319 obser-vations obtained with A1 and A2 between April 2019and June 2020 that meet the criteria t obs >
40 min,S/N >
40, and σ TOA < µ s. Details of this data set aresummarized in Table 2.0 Sosa Fiscella et al.
We analyze the data sets of each antenna both in-dependently and altogether. As described in Sec. 5.1,ECORR accounts for noise that is correlated betweenobservations that were obtained simultaneously at dif-ferent frequencies. Since such observations are not avail-able for a single antenna, we exclude the ECORR pa-rameter from the analysis of the individual data sets.However, we do include this parameter when analyzingthe A1+A2 data set in order to profit from the simul-taneous observations at different frequencies. A cornerplot for these parameters and their errors is shown inFig. 9.The fitted values to the white- and red-noise parame-ters for the different data sets are presented in Table 4.Complementary, we explore the possibility of splittinglong-duration observations into two subintegrations of t min = 75 min in order to sample shorter timing frequen-cies. The adjusted values for this case, also presented inTable 4, are consistent within 1 σ to the ones obtainedwithout the splitting. We therefore conclude that split-ting long observations does not improve the timing anal-ysis, in line with the conclusion from Sect. 4.4.We obtain EQUAD ≈ . µ s in all cases. The valueof Γ rn is less constrained and consistent within 0.5–1.5,while the amplitude is A rn ≈ × − . 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 is due to the relativelyshort baseline of our current data set. EFAC log EQUAD Γ rn log A rn A1 2 . +0 . − . − . +0 . − . . +0 . − . − . +0 . − . A2 2 . +0 . − . − . +0 . − . . +0 . − . − . +0 . − . A1+A2 2 . +0 . − . − . +0 . − . . +0 . − . − . +0 . − . A1+A2* 2 . +0 . − . − . +0 . − . . +0 . − . − . +0 . − . Table 4.
Adjusted values for the white and red noise pa-rameters. Values marked with (*) were obtained by splittingthe observations as described in Sec. 4.4.
Gravitational-wave analysis
We now embark on setting the first bounds to theGW amplitude from massive binary black holes usingobservations from IAR. In doing so, we aim to exploitthe high cadence of these observations.5.3.1.
Gravitational wave analysis: stochastic background
The contribution of the GWB coming from an en-semble of supermassive black-hole binaries or primordial fluctuations during the
Big Bang is modeled similarly tothat of the red noise (Eq. 8). Any GWB component ismodeled as a single stationary Gaussian process with apower-law timing-residual spectral density P ( f ) = A π (cid:18) ff ref (cid:19) Γ gwb yr . (9)The analysis is nearly identical to the red noise anal-ysis described in Sec. 5.2. The prior on the GWBamplitude is taken uniform on log ( A gwb [yr / ]) ∈ [ − . , − gwb ∈ [0 , . t obs ≥
75 min (seeSec. 4.4). The best-fitted values to each GWB param-eter and for each set of observations are presented inTable 5.Using the split observations, we get a higher ca-dence at the cost of worsening the S/N (and thereforethe TOAs precision) of each data point. Our resultsshow consistent values of A gw ≈ (3 ± × − andΓ gw ≈ . ± . gw and slightly higher valuesof A gw , though these differences are not significant asthey are within 1 σ of the values obtained without thesplitting.While the amplitude we find is consistent with ex-pected bounds for the stochastic background, Γ gw fallsshort from the expected 13/3 for a stochastic GW back-ground of SMBH binaries (Siemens et al. 2013), possiblydue to our current relatively short observational baselineof ≈ . rn ∈ [0 ,
7] andan extra red-noise process with Γ gwb set to 4 .
33. Wealso fix all of the white-noise parameters to the valuesobtained in Sec. 5.1. In Fig. 10 we show a corner plotof the fit to the joint A1+A2 data sets. We find valuesof A rn ≈ (4 ± × − , consistent with our previousresults, and Γ rn ≈ . ± .
1. Such uncertainties maybe attributed to the short time span.5.3.2.
Gravitational-wave analysis: continuous source
A single supermassive binary black-hole system pro-duces “continuous” GWs because the system does notevolve notably over the few years of a pulsar-timing dataset. We used the Python package
Hasasia (Hazboun
AR timing analysis of PSR J0437 − . . . E F A C − . − . − . − . − . l og E Q UA D − . . . . Γ r n − . − . − . − . − . log ECORR − . − . − . − . − . l og A r n . . . EFAC − . − . − . − . − . log EQUAD − . . . . Γ rn − . − . − . − . − . log A rn Figure 9.
ENTERPRISE timing analysis of red noise for 1.1 yr of observations of J0437 − gwb . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . log A gwb − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . − . +0 . − . Table 5.
Best-fit values to the GWB parameters. Sosa Fiscella et al. − − − − − l og A r n . . . . Γ rn − . − . − . − . l og A g w b − − − − − log A rn − . − . − . − . log A gwb Figure 10.
ENTERPRISE gravitational wave analysis forJ0437 − gwb = 4 . et al. 2019) to calculate the single-pulsar sensitivitycurve of our data set of J0437 − h c ( f ) = (cid:112) f S ( f ) (10)where S is the strain-noise power spectral density forthe pulsar. This is related to the power spectrum of theinduced timing residuals of Eq. 9 by (see Jenet et al.2006) P ( f ) = 112 π f h c ( f ) (11)The white- and red-noise parameters that were ad-justed in Secs. 5.1 and 5.2 using ENTERPRISE are loadedinto the package in order to account for these effects inthe calculations. Since our observations have a timebaseline of T obs = 1 . / (10 T obs ) ∼ . × − Hz and 1 / (1 day) ∼ . × − Hz.The resulting sensitivity curve is shown in Fig. 11.It is readily seen that there is a loss of sensitivity ata frequency of (1 yr) − , caused by fitting the pulsar’sposition, and at a frequency of (PB) − ∼ µ Hz (withPB the orbital period), caused by fitting the orbital pa-rameters of the binary system. The additional spikesseen at frequencies higher than (PB) − correspond toharmonics of the binary orbital frequency. In addition, the sensitivity at lower frequencies is re-duced by: i) the fit of a quadratic polynomial to theTOAs required to model the pulsar spin-down, and ii)the fitting of ‘jumps’ to connect the timing residualsobtained with different backends (Yardley et al. 2010).The frequency dependence ( ∼ f − / ) at low frequenciesis evidence of a fit to a quadratic spin-down model forthe pulsar spin frequency. As a result, the minimum ofthe sensitivity curve should be attained at a frequencyof 1 /T obs . However, given that the T obs of our dataset is close to one year, this feature coincides with theloss of sensitivity at (1 yr) − . We expect to obtain awell-defined minimum at ≈ /T obs in a future by accu-mulating more observations and achieving a significantlylonger time baseline.For completeness, we tested the significance of thered-noise contribution by calculating a sensitivity curvewithout this component. The curve was essentially in-sensitive to those changes in the priors. This is expected,since the injection of red noise should lead to a flat sensi-tivity curve around the minimum (Hazboun et al. 2019),though in our case it is coincident with the spike at(1 yr) − .For comparison, we used ENTERPRISE to perform afixed-frequency MCMC at four different frequencies. Weobtained a posterior distribution for log h gw at eachof these frequencies with a mean value in great agree-ment with the curve obtained with Hasasia , as shownin Fig. 11. − − − − − Frequency [Hz] − − − − − C h a r a c t e r i s t i c S t r a i n , h c J0437-4715 f − / /T obs / (10 T obs )1 / (1 day)1 / (1 year)1 / PB ENTERPRISE
Figure 11.
Sensitivity curve for J0437 − /T obs ,the dotted red line to 1 /T yr and the dotted purple line to1 / PB. The black crosses correspond to the mean values ofthe log h gw distributions obtained using ENTERPRISE . AR timing analysis of PSR J0437 − z ≈ h ∼ − at the time of ar-rival to our Galaxy, roughly a million years ago (Loustoet al. 2017). CONCLUSIONSWe presented the first detailed analysis of the obser-vational campaign towards the bright MSP J0437 − • The number of phase bins used in the reductiondoes not have an impact on the timing precisionas long as n bins ≥ • The S/N of the individual observations plays a cru-cial role in determining the timing precision. Inparticular, to achieve a timing precision < µ s,observations with S / N >
140 are required, a con-dition that is currently fulfilled by ∼ / ∼ / • Splitting long observations into shorter intervalsdoes not improve the timing precision, most likelydue to current limitations in the S/N for short ob-servations. • A1 slightly outperforms A2, probably due to itslarger bandwidth configuration. • The systematic errors of the observations are σ sys ≈ . µ s, although this value is likely to beS/N-limited. The RMS of the data set is ≈ . µ s • The white-noise analysis performed with
ENTERPRISE indicates that the error bars aretypically underestimated by a factor ∼ • We placed upper limits to the GWB in thetens nHz to sub- µ Hz frequency range. Although the current sensitivity is not sufficient for plac-ing physically-interesting constraints, the ongo-ing campaign –together with incoming hardwareupgrades– is likely to significantly improve in thenext 5–10 years (see also Lam & Hazboun 2020).In particular, observations lasting over 3 h arepromising for exploring GW signals with frequen-cies above 0.1 µ Hz by splitting them into hour-scale subintegrations.Ongoing and future hardware upgrade of IAR’s anten-nas, such as installing larger bandwidth boards, promiseto expand IAR’s observational capabilities and improveits achievable timing precision. Such upgrades wouldallow to reduce the systematical errors of the antennasand to include (sub)daily high-precision timing of otherMSPs of interest, such as PSR J2241 − Facility:
IAR
Software:
ENTERPRISE (Ellis et al. 2019),PRESTO (Ransom et al. 2003; Ransom 2011),PSRCHIVE (Hotan et al. 2004), PyPulse (Lam 2017),Hasasia (Hazboun et al. 2019)APPENDIX4
Sosa Fiscella et al. A. DETAILS OF THE ANALYSISA.1.
Reduction of observations
To de-disperse and fold the observations we use thesoftware
PRESTO (Ransom et al. 2003; Ransom 2011). Ithas a variety of tools for the reduction of observations.The processed data is stored in a .pfd file that containsthe 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 model-ing the variation of the pulsar period. These coefficientsallow to determine the period of pulsation in a topocen-tric reference system and are necessary to compute thetiming residuals.If the observation is in the file obs.fil and the maskin the file m.mask , then the command-line used has thefollowing syntax, prepfold -nsub 64 -n 1024 -timingJ0437 − (A1)where the option -timing indicates prepfold to gener-ate a file .polycos based on the pulsar parameters Thisprocess is currently automatized through local Pythonscripts. A.2. Templates
Considering that A1 and A2 have different configura-tions (number of polarizations and bandwidth; see Ta-ble 1), it is possible that slight differences arise in theintegrated profile seen by each antenna. We thereforestudy whether the template used has a significant im-pact in the timing residuals.To create each template we chose observations with n bins = 1024 phase bins and n chan = 64 frequencychannels. We select for each antenna data the highestS/N observation and extract the noise by using the task psrsmooth in the package psrchive . This choice of tem-plates seems adequate since J0437 − jit-ter of the pulsar (Liu et al. 2012). The selected tem-plates for each antenna correspond to the profiles with n bins = 1024 phase bins in Fig. 12. The relative errorbetween them is below 5% near the peak, with largerrelative differences towards the wings, but those do nothave major influence in the determination of the TOAs.In the preliminary timing analysis of PSR J0437 − .
25 0 .
50 0 . . . . . . . N o r m a li ze d fl u x A1 n bin = 64n bin = 128n bin = 256n bin = 512n bin = 1024 .
25 0 .
50 0 . n bin = 64n bin = 128n bin = 256n bin = 512n bin = 1024
64 128 256 512 1024 n bins . . . . . . R M S [ µ s ] A1 64 128 256 512 1024 n bins A2 Figure 12.
Top:
Templates for each antenna for differentvalues of n bins . Bottom: RMS found for each subset per n bins , and its corresponding 1- σ error bars. assumption was valid to the current level accuracy, pro-ducing an RMS = 0 . µ s residual for A2 observationswith the use of either template . Notwithstanding this aposteriori verification, we consistently use different tem-plates for A1 and A2 throughout this work.A.3. Timing versus number of phase bins
In order to study the effect of the number of phasebins ( n bins ) used in the folding of observations on thetiming residuals, we have taken data folded originallywith n bins = 1024, and processed with the routine bscrunch of the psrchive package for Python, to gen-erated copies of the observations and their correspond-ing templates for each antenna, but with values of n bins = 512 , , ,
64 and 32. Through this process of scrunching we obtained 6 sets of observations for eachantenna only differing by their n bins . Fig. 12 shows theeffect of the n bins on the templates for each antenna.While for n bins = 32 we lose temporal resolution, thedifferences beyond n bins ≥
256 are almost negligible toour precision.Next we compute the timing residuals for each n bins subset. Interestingly, only for n bins ≤
64 the timing er-rors are too large; for n bins ≥
128 the derived TOAs arevery consistent, being the size of the error bar the maindifference (with smaller error bars obtained for larger n bins ). The RMS of the residuals for each subset afteradjusting σ sys as a function of n bins is shown in Fig. 12. AR timing analysis of PSR J0437 − n bins significantlyfor 64 to 256 bins , showing that we cannot attain goodtiming for n bins ≤
64 and need at least 256 bins to obtaina precision higher than 1 µ s for a pulsar like J0437 − µ s at 1400 MHz. A.4. Timing vs S/N
Here we present additional figures that support thehypothesis that our timing studies are limited due tothe S/N of the observations. This effect has a largerimpact for A2, as can be seen in Fig. 13.REFERENCES
Alam, M. F., Arzoumanian, Z., Baker, P. T., et al. 2020,arXiv e-prints, arXiv:2005.06495.https://arxiv.org/abs/2005.06495Arzoumanian, Z., Brazier, A., Burke-Spolaor, S., et al.2016, ApJ, 821, 13, doi: 10.3847/0004-637x/821/1/13Arzoumanian, Z., Baker, P. T., Blumer, H., et al. 2020,arXiv e-prints, arXiv:2009.04496.https://arxiv.org/abs/2009.04496Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33,doi: 10.1051/0004-6361/201322068Astropy Collaboration, Price-Whelan, A. M., Sip˝ocz, B. M.,et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4fBecker, W., & Tr¨umper, J. 1993, Nature, 365, 528,doi: 10.1038/365528a0Bell, J. F., Bailes, M., & Bessell, M. S. 1993, Nature, 364,603, doi: 10.1038/364603a0Cordes, J. M., & Chernoff, D. F. 1997, ApJ, 482, 971,doi: 10.1086/304179Danziger, I. J., Baade, D., & della Valle, M. 1993, A&A,276, 382Ellis, J., & van Haasteren, R. 2017,jellis18/PTMCMCSampler: Official Release, 1.0.0,Zenodo, doi: 10.5281/zenodo.1037579Ellis, J. A., Vallisneri, M., Taylor, S. R., & Baker, P. T.2019, ENTERPRISE: Enhanced Numerical ToolboxEnabling a Robust PulsaR Inference SuitE.http://ascl.net/1912.015Ferdman, R. D., van Haasteren, R., Bassa, C. G., et al.2010, Classical and Quantum Gravity, 27, 084014,doi: 10.1088/0264-9381/27/8/084014Foreman-Mackey, D. 2016, The Journal of Open SourceSoftware, 1, 24, doi: 10.21105/joss.00024Gancio, G., Lousto, C. O., Combi, L., et al. 2020, A&A,633, A84, doi: 10.1051/0004-6361/201936525Gwinn, C. R., Hirano, C., & Boldyrev, S. 2006, A&A, 453,595, doi: 10.1051/0004-6361:20054280Hartnett, J. G., & Luiten, A. N. 2011, Reviews of ModernPhysics, 83, 1, doi: 10.1103/RevModPhys.83.1 Hazboun, J., Romano, J., & Smith, T. 2019, Journal ofOpen Source Software, 4, 1775, doi: 10.21105/joss.01775Hazboun, J. S., Romano, J. D., & Smith, T. L. 2019,PhRvD, 100, 104028, doi: 10.1103/PhysRevD.100.104028Hazboun, J. S., Simon, J., Taylor, S. R., et al. 2020, ApJ,890, 108, doi: 10.3847/1538-4357/ab68dbHobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006,MNRAS, 369, 655, doi: 10.1111/j.1365-2966.2006.10302.xHotan, A. W., van Straten, W., & Manchester, R. N. 2004,PASA, 21, 302, doi: 10.1071/AS04022Jenet, F., Hobbs, G., van Straten, W., et al. 2006,Astrophys. J., 653, 1571, doi: 10.1086/508702Johnston, S., Lorimer, D. R., Harrison, P. A., et al. 1993,Nature, 361, 613, doi: 10.1038/361613a0Kargaltsev, O., Pavlov, G. G., & Romani, R. W. 2004,ApJ, 602, 327, doi: 10.1086/380993Kaspi, V. M., Taylor, J. H., & Ryba, M. F. 1994, ApJ, 428,713, doi: 10.1086/174280Keith, M. J., Coles, W., Shannon, R. M., et al. 2013,MNRAS, 429, 2161, doi: 10.1093/mnras/sts486Knuth, K. H. 2006, arXiv e-prints, physics/0605197.https://arxiv.org/abs/physics/0605197Lam, M. T. 2017, PyPulse: PSRFITS handler.http://ascl.net/1706.011Lam, M. T., & Hazboun, J. S. 2020, arXiv e-prints,arXiv:2007.00260. https://arxiv.org/abs/2007.00260Lentati, L., Alexander, P., Hobson, M. P., et al. 2014,MNRAS, 437, 3004, doi: 10.1093/mnras/stt2122Liu, K., Keane, E., Lee, K., et al. 2012, Mon. Not. Roy.Astron. Soc., 420, 361,doi: 10.1111/j.1365-2966.2011.20041.xLorimer, D. R., & Kramer, M. 2012, Handbook of PulsarAstronomy (Cambridge University Press)Lousto, C. O., Zlochower, Y., & Campanelli, M. 2017,ApJL, 841, L28, doi: 10.3847/2041-8213/aa733cMatsakis, D. N., Taylor, J. H., & Eubanks, T. M. 1997,A&A, 326, 924Os(cid:32)lowski, S., van Straten, W., Hobbs, G. B., Bailes, M., &Demorest, P. 2011, MNRAS, 418, 1258,doi: 10.1111/j.1365-2966.2011.19578.x Sosa Fiscella et al. − − R e s i du a l s [ µ s ] S/N > > − − R e s i du a l s [ µ s ] S/N >
80 S/N >
100 200 300 400 500
MJD - 58500 − − R e s i du a l s [ µ s ] S/N >
100 200 300 400 500
MJD - 58500
S/N > − − R e s i du a l s [ µ s ] S/N > > − − R e s i du a l s [ µ s ] S/N >
80 S/N >
100 200 300 400 500
MJD - 58500 − − R e s i du a l s [ µ s ] S/N >
100 200 300 400 500
MJD - 58500