Decade-long Timing of Four GMRT Discovered Millisecond Pulsars

Millisecond pulsars (MSPs) are fast rotating neutron stars that exhibit exceptional clock-like behavior over time. Due to their rotational stability, we can probe the effects due to the interstellar medium along the line of sight as well as the effects that are intrinsic to them (e.g., Foster & Cordes 1990). Pulsar timing (e.g., Lorimer & Kramer 2004) is a technique that utilizes the time-of-arrival (ToA) information of the pulses coming from MSPs, thus revealing precise astrometric, rotational, and binary models as well as propagation effects.

The long-term timing of MSPs helps to reduce covariance between system model parameters, resulting in more precise and accurate estimation of these parameters. It can reveal higher-order effects in the system, such as the higher-order derivatives of the spin frequency, temporal variations in orbital motion in the case of binaries, and a change of the electron–proton plasma along the line of sight (e.g., Lorimer & Kramer 2004). Better modeling of pulsars by taking these higher-order effects into account will result in more accurate predictions for their ToAs, making them better interstellar clocks.

Detweiler (1979) introduced the utilization of long-term pulsar timing of MSPs to identify an isotropic stochastic gravitational wave (GW) background of cosmic origin. Its primary contributor is expected to be an ensemble of merging galaxies having supermassive black hole binaries at their centers (Burke-Spolaor et al. 2019). A pulsar timing array (PTA) is an array of well-timed MSPs with varying angular separations in the sky designed to detect the effect of nanohertz (nHz) GWs on the cross-correlation of timing residuals from pairs of MSPs (Hellings & Downs 1983). The current major PTAs are the North American Nanohertz Observatory for Gravitational Waves (NANOGrav; Jenet et al. 2009), the Parkes Pulsar Timing Array (PPTA; Manchester 2006), the European Pulsar Timing Array (EPTA; Stappers et al. 2006), the Chinese Pulsar Timing Array (CPTA; Lee 2016), the MeerKAT Pulsar Timing Array (MPTA; Bailes et al. 2016), and the Indian Pulsar Timing Array (InPTA; Joshi et al. 2018).

Recently, Agazie et al. (2023; NANOGrav), Reardon et al. (2023; PPTA), Antoniadis et al. (2023; EPTA+InPTA), and Xu et al. (2023; CPTA) have shown Hellings and Downs curve signatures in cross-correlation statistics of the residuals of a set of MSPs. For the four PTAs, the signal-to-noise ratio (S/N) of the common correlated signal among the pairs of MSPs ranges from 3 to 5. The S/N is defined here as ρ/σ (Jenet et al. 2005). ρ measures the similarity between the cross-correlation distribution of residuals (as a function of angular separation between the MSPs) and the Hellings and Downs function (Hellings & Downs 1983). σ is the standard deviation of the cross-correlation distribution. According to Agazie et al. (2023), the next step will be to merge data sets from NANOGrav, PPTA, EPTA, and InPTA, which will consist of around 80 MSPs and a temporal baseline of up to 24 yr. This is expected to boost the detection significance of the common correlated signal. A high S/N cross-correlation curve will clarify the source of its origin, and a larger number of pulsars may reveal anisotropy in the GW background (e.g., Hotinli et al. 2019), polarization structures (e.g., Sato-Polito & Kamionkowski 2022), and so on. Along with incorporating data sets from multiple PTAs, the individual PTAs intend to improve their sensitivity to this signal by increasing the number of MSPs and the observational timing baseline.

Siemens et al. (2013) describe an equation that shows the dependency of the S/N of the cross-correlation statistics on the number of MSPs included in the GW detection experiment, the observational cadence, the timing baseline, the timing precision achieved for individual MSPs, and the amplitude of GWs. It shows that in the intermediate signal regime (Siemens et al. 2013) the number of MSPs included in the experiment has the greatest influence on the detection significance of the common correlated GW signal, implying that identifying new MSPs which are suitable for PTAs should be of greater importance.

The Giant Meterwave Radio Telescope (GMRT) is an interferometer with 30 parabolic antennas and a longest baseline of 25 km. Each antenna dish has a 45 m diameter and two orthogonal polarizations ³ (Swarup et al. 1997). The legacy GMRT system (Roy et al. 2010) is equipped with the GMRT Software Backend (GSB) with a maximum instantaneous observational bandwidth of 33 MHz within a frequency range of 120 to 1460 MHz. The upgraded GMRT (uGMRT; Gupta et al. 2017; Reddy et al. 2017) with the GMRT Wide-band Backend (GWB) became operational in 2017, which can have a maximum instantaneous observational bandwidth of 400 MHz in the same frequency range as legacy GMRT. The GMRT High Resolution Southern Sky (GHRSS) survey (Bhattacharyya et al. 2016, 2019; Singh et al. 2022; Sharma et al. 2023; Singh et al. 2023, 2023) and Fermi-directed survey (Bhattacharyya et al. 2013; Roy et al. 2015; Bhattacharyya et al. 2021, 2022) have so far discovered 12 MSPs using legacy GMRT and uGMRT systems. For the timing campaign, we have been monitoring four of these regular MSPs that are bright in the lower-frequency bands of the GMRT. These four MSPs were observed using the legacy GMRT system from 2011 to 2017, and after 2017 they have been observed using the uGMRT system. Phase-coherent timing has been obtained for the GSB observations for these four MSPs and was reported in Bhattacharyya et al. (2019) and Bhattacharyya et al. (2022). Sharma et al. (2022) reported the phase-coherent timing with GWB observations for the same MSPs, which accounts for 2−4 yr of observational timing baseline.

In this paper, we present the timing of GWB observations with an increased timing baseline and combine it with legacy GMRT timing data sets, resulting in around a decade-long baseline for the four MSPs. Sections 2 and 3 provide the specifics of the observations and the timing technique used in this study, respectively. Section 4 reports the results from our timing analysis, a comparison between GSB and GWB timing results, and a comparison of the timing residuals for the four GMRT MSPs with other PTAs. Section 5 summarize the findings from the decade-long timing of these four MSPs, including the suitability of GMRT-discovered MSPs and the decade-long timing data set for PTAs.

For phased array (PA) GWB observations in uGMRT band 3 (300–500 MHz) and band 4 (550–750 MHz), 70% and 80% of the array is phased, providing gains of 7 and 8 K Jy⁻¹, respectively. An array with similar gains as GWB is phased for PA beam observations in legacy GMRT band 3 (306–339 MHz) and band 4 (591–624 MHz; Bhattacharyya et al. 2019). Table 1 lists the configurations for the observations in various receiver backends, dedispersion modes, frequency bands, time resolution, and antenna count. Intrachannel smearing has been corrected for using coherent dedispersion mode data sets with known DM values for the observed MSPs. The residual smearing for incoherent dedispersion mode data sets is determined by 2048 (GSB) or 4096 (GWB) channels over a bandwidth of 33 and 200 MHz, respectively.

We observed four GMRT-discovered MSPs, which are presented in Table 2. The table shows the mean observation time, the GMRT backend utilized, the S/N in various frequency bands, the number of epochs in various observing modes and bands, and the timing baseline for the observed MSPs. The flux density values and profiles of these MSPs at various frequencies are reported in Sharma et al. (2022). The spin periods and DMs are listed in the timing ephemerides table (Table 4) in Section 4.

The four MSPs were selected from the set of GMRT-discovered MSPs based on S/N (provided in the same table) and the timing precision reported in Bhattacharyya et al. (2022), Sharma et al. (2022), and Bhattacharyya et al. (2019). The four MSPs are bright in GMRT band 3 and band 4, with the exception of J1120−3618, which has a significantly lower detection significance in band 4. These MSPs have been observed roughly with a monthly cadence for the last 7.7–11.0 yr.

Following observations made with the GWB, the GMRT raw PA beam data sets are converted to filterbank format, which is then incoherently dedispersed and folded with the PREPFOLD command of PRESTO (Ransom 2011). For folding purposes, we make use of the ephimeris presented in Bhattacharyya et al. (2022) and Bhattacharyya et al. (2019). The folded data cubes (PFDs) produced by the PREPFOLD command are then converted to FITS format using the PAM command of PSRCHIVE (van Straten et al. 2012). Following the same profile bin configuration as Sharma et al. (2022), we created 128 and 64 bins for the coherent and incoherent dedispersed profiles, respectively, in FITS files for all the MSPs. The FITS files are averaged to a single time integration, 64/128 bins, and 16 frequency subbands. These FITS files are then used for further timing analysis, as described in Section 3.

For all GSB observations, phase-coherent timing has been established in Bhattacharyya et al. (2022), and Bhattacharyya et al. (2019). The GSB ToAs from these works are being used directly for this study for aiding long-term timing.

This section outlines the procedure that we used to create templates for the MSPs, derive the DM from each epoch of observation, and produce ToAs for the timing analysis.

Procedure for creating the templates: we collected all high-S/N FITS files in a specified frequency band/mode for a particular MSP. We use ppalign module of the PulsePortraiture ⁴ (Pennucci et al. 2016; Pennucci 2019, and Pennucci et al. 2014) software to create a two-dimensional template using the FITS files. The ppalign module aligns the selected set of FITS files iteratively. This module is equipped with a robust alignment algorithm that takes into consideration not only a constant temporal phase offset between different epoch observations but also rotates each frequency subband profile of each epoch by an offset proportional to ν⁻². Here, ν is the frequency of the corresponding frequency subband profile. This results in an aligned set of FITS files which are then averaged while keeping the frequency resolution to obtain a two-dimensional template (function of frequency and phase bins).

For each MSP the resultant two-dimensional template has one time integration, 64/128 profile bins, and 16 frequency subbands. The two-dimensional template is then frequency averaged to provide a single profile with 64/128 bins. This profile is then used as a template for the particular dedispersion mode and frequency band for that MSP. Note that we use separate templates for the two dedispersion modes and different uGMRT receiver bands for each MSP.

Deriving the DM: we considered all the FITS files of a single MSP in a single band and a single dedispersion mode at a time. The template corresponding to this band, mode, and MSP is used to generate frequency subband ToAs from the FITS files (with one time integration, 64/128 bins, and 16 frequency subbands), resulting in 16 subband ToAs per epoch. We use the PAT command of PSRCHIVE to generate ToAs. TEMPO2 (Hobbs & Edwards 2012) is used to fit for the DM using these subband ToAs for each epoch, keeping the other parameters fixed while fitting.

ToAs for timing: we averaged all FITS files corresponding to a specific MSP, band, and mode in frequency so that they have one time integration, 64/128 bins, and one frequency subband. The corresponding template for this set is used to generate a single ToA for each epoch observation. Again, we create ToAs using the PAT command of PSRCHIVE. These band-averaged ToAs are used for timing purposes.

We have measured the temporal variation of the DM for each of the four MSPs from the band-3 GWB data sets, and show them in each pulsar subsection below (Figures 1–4). We did not include the DM values from GWB band 4 and the GSB data sets in these plots since their error bars (GWB band 4: ∼5 × 10⁻³ pc cm⁻³; GSB band 3/4: ∼10⁻¹ pc cm⁻³ or higher) are significantly larger than those from the GWB band-3 data sets (∼5 × 10⁻⁴ pc cm⁻³). Due to the very large DM error bars in the GSB data sets, we did not consider the DM variations with time for any of the MSPs to treat the GSB and GWB data sets equally when doing timing. Instead, we used a global DM for the combined timing of the band 3+band 4 and GSB+GWB data sets (named GSB+GWB from now on) of each MSP.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3 shows the number of ToAs corresponding to individual frequency bands and receiver backends for each MSP, as well as the median ToA uncertainty and the rms of the postfit timing residuals obtained from the individual GSB, GWB, and GSB+GWB timings for these MSPs. We show the postfit residuals for each MSP, corresponding to GSB+GWB timing, in each pulsar subsection below (Figures 5–8). Table 4 shows the timing ephimerides and derived parameters obtained from the phase-coherent timing of GSB+GWB.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. J1120−3618

J1120−3618 is a binary MSP with a spin period of ∼5.56 ms, an orbital period of ∼5.66 days, and a DM of ∼45.12 pc cm⁻³. The timing precision obtained in the individual GWB and GSB timings is 6.1 and 27.0 μs respectively, each spanning more than 4 yr of baseline. We achieved an overall timing precision of 10.5 μs after combining the GSB and GWB timing data sets (spanning 10.3 yr), which is clearly limited by the larger error bars of ToAs in the GSB timing data.

This MSP's line of sight has been crossing an electron-rich medium for the last 4–5 yr, as indicated by its DM trend seen in GWB band 3. Its DM value has increased by 8.6 × 10⁻³ pc cm⁻³ from 2018 February to 2023 January while the median DM error bar for this MSP in GWB band 3 is 1.0 × 10⁻³ pc cm⁻³.

We measure a total proper motion for this MSP, for the first time, of 6.0(6) mas yr⁻¹ from GSB+GWB timing, and individually 5.2(6) mas yr⁻¹ in R.A. and 2.8(7) mas yr⁻¹ in decl. The proper motion for J1120−3618 falls within the typical proper motion range for known MSPs [few tens of mas yr⁻¹; ATNF pulsar catalog ⁵ (Manchester et al. 2005)].

We estimate the distance ⁶ to this MSP using its DM, position, and the galactic electron density models NE2001 (Cordes & Lazio 2002) and YMW16 (Yao et al. 2017), which result in a distance of 1.75 and 0.95 kpc, respectively. Using these estimates and the proper motion value for this MSP, we calculate transverse velocities of 49 km s⁻¹ (for 1.74 kpc; NE2001) and 27 km s⁻¹ (for 0.95 kpc; YMW16). These transverse velocity values are significantly lower than the maximum expected space velocity of the MSPs (∼270 km s⁻¹) in the scenario of spherically symmetric supernova explosions, suggesting that this MSP did not experience an asymmetric kick during its recoiling phase (Tauris & Bailes 1996).

The timing fit for J1120−3618 yields a spin-period derivative (P1) of 9.4(3) × 10⁻²² s s⁻¹, an order of magnitude less than the other three MSPs in this work. It has an intrinsic spin-period derivative (Toscano et al. 1999) of ∼5.5 × 10⁻²² s s⁻¹ (∼0.6 P1; determined using Equation 2 of Toscano et al. 1999; its proper motion, and the distance measurements from the YMW16 model). The P1 value of J1120−3618 is the fifth smallest of all MSPs currently known.

Figure 9 shows a comparison of fitted model parameters from the GWB, GSB, and GSB+GWB timing for this MSP. Commonly fitted model parameters (listed in Figure 9) obtained from GSB+GWB timing are on average 1.1 and 3.6 times more precise than individual GWB and GSB timings, respectively. The model parameters from the GSB+GWB timing are consistent with individual GWB or GSB timings within ±1σ, where σ is the error bar of the fitted model parameter from either GWB or GSB.

Zoom In Zoom Out Reset image size

4.2. J1646−2142

J1646−2142 is an isolated MSP with a spin period of ∼5.85 ms and a DM of ∼29.74 pc cm⁻³. Timing precision of the GWB and GSB timings are 7.3 and 15.7 μs, respectively. With the combined GSB and GWB timing data, we achieved a timing precision of 11.1 μs for a span of 11.0 yr.

The median DM uncertainty for this MSP in GWB band 3 and band 4 is 6.5 × 10⁻⁴ and 5.1 × 10⁻³ pc cm⁻³, respectively. We see an overall DM change of less than ±3 σ_DM from 2017 July to 2022 October, where σ_DM is the median DM uncertainty in GWB band 3.

The GSB+GWB timing fit resulted in an insignificant detection of proper motion; therefore, we excluded its fitting in the timing. Figure 10 shows the differences in the model parameters and their precision in the GSB+GWB timing against the individual GWB and GSB timings. We see an average improvement in model parameter precision of 5.7 and 2.4 times in the GSB+GWB timing when compared individually to the GWB and GSB timings, respectively. The model parameters from the GSB+GWB timing are consistent with the individual GSB or GWB timings within ±1σ, where σ is the error bar of the fitted model parameter from either GSB or GWB.

Zoom In Zoom Out Reset image size

4.3. J1828+0625

J1828+0625 is a binary MSP with a spin period of ∼3.63 ms, an orbital period of ∼77.92 days, and a DM of ∼22.42 pc cm⁻³. We achieved timing precision of 6.2 and 20.9 μs in individual timing of the GWB and GSB data sets, respectively. The combined GSB+GWB timing resulted in a timing precision of around 11.8 μs for a span of 10.9 yr.

The median DM uncertainty of this MSP in GWB band 3 and band 4 is 8.1 × 10⁻⁴ and 5.2 × 10⁻³ pc cm⁻³, respectively. We find that the DM variation for this MSP is confined within ±3 σ_DM from 2019 June to 2022 September, where σ_DM is the median DM uncertainty in GWB band 3.

For this MSP, we estimated a total proper motion of 14.4(2) mas yr⁻¹, which is well within ±2σ of the values reported in Bhattacharyya et al. (2022). The increased timing span with precise GWB ToAs allows us to derive the proper motion with at least five times higher precision than the values reported in Bhattacharyya et al. (2022) using only GSB data.

The distance estimates for this MSP using the NE2001 and YMW16 models are 1.12 and 1.00 kpc, respectively. Its proper motion and these distances give a transverse velocity of 77 km s⁻¹ (for 1.12 kpc; NE2001) and 68 km s⁻¹ (for 1.00 kpc; YMW16), which is comparable to the transverse velocities of J1120-3618 and suggests the absence of an asymmetric kick during recoiling. The timing fit for this MSP resulted in a precise P1 value of 4.6882(9) × 10⁻²¹ s s⁻¹. Using its transverse velocity and distance (for YMW16), the intrinsic spin period derivative is estimated to be ∼2.7 × 10⁻²¹ s s⁻¹ (∼0.6 P1).

Figure 11 compares the model parameters obtained from the GSB+GWB timing to the individual GSB or GWB timings. The model parameters obtained from the GSB+GWB timing are on average 10.5 times more precise than the individual GWB timings while being 4.5 times more precise than the individual GSB timings. The model parameters resulting from the GSB+GWB timing are in agreement with the parameters obtained from the individual GWB or GSB timings within ±1σ, where σ is the error in model parameter obtained from the individual timings of the GWB or GSB data sets.

Zoom In Zoom Out Reset image size

4.4. J2144−5237

J2144−5237 is a binary MSP with a spin period of ∼5.04 ms, an orbital period of ∼10.58 days, and a DM of ∼19.55 pc cm⁻³. The timing precision achieved in individual timing of GWB and GSB data is 9.0 and 17.1 μs, respectively, whereas the combined GSB+GWB timing resulted in a timing precision of 10.6 μs over a 7.7 yr span.

The median DM uncertainty for this MSP in GWB band 3 and band 4 is 4.3 × 10⁻⁴ and 5.8 × 10⁻³ pc cm⁻³, respectively. We see small variations in the DM from 2017 July to 2023 January, which are contained within ±3 σ_DM, where σ_DM is the DM uncertainty in GWB band 3.

Individual GSB and GWB timings resulted in unreliable proper motion values with large error bars. We find the proper motion for this MSP (9.8(5) mas yr⁻¹) for the first time using the combined GSB+GWB timing.

The NE2001 and YMW16 models (together with DM and position) gave distance estimates of 0.80 and 1.66 kpc for this MSP, respectively. The resulting transverse velocities are 37 km s⁻¹ (for 0.80 kpc; NE2001) and 77 km s⁻¹ (for 1.66 kpc; YMW16), respectively. These numbers are comparable to the transverse velocities of MSPs J1120-3618 and J1828-0625, and suggest a symmetric supernova explosion for this MSP as well. The timing fit for this MSP resulted in a P1 value of 9.055(3) × 10⁻²¹ s s⁻¹, from which the intrinsic spin-period derivative is calculated to be ∼5.9 × 10⁻²¹ s s⁻¹ (∼0.6 P1).

Figure 12 compares the model parameters obtained from the GSB+GWB timing with those obtained from the individual GWB and GSB timings. The model parameters obtained from the GSB+GWB timing are on average 3.4 and 5.7 times more precise than the parameters from the individual GWB and GSB timings, respectively. Also, the model parameter values obtained from the GSB+GWB timing are within ±1σ of the model parameters obtained from the individual GWB and GSB timings, where σ is the error in model parameter obtained from individual timings of the GWB or GSB data sets.

Zoom In Zoom Out Reset image size

4.5. Comparison with MSPs Currently Included in PTAs

The International Pulsar Timing Array (IPTA; Verbiest et al. 2016; Perera et al. 2019) combines timing data from three major PTAs, NANOGrav, PPTA, and EPTA, with the goal of improving the timing precision for individual MSPs. In the IPTA second data release, Perera et al. (2019) reported timing analyses of 65 MSPs, providing the model parameters and timing residuals obtained for each MSP. In Figure 13, we plotted the rms of the timing residuals obtained for the 65 IPTA MSPs and four GMRT MSPs from GSB+GWB timing.

Zoom In Zoom Out Reset image size

The four GMRT-discovered MSPs fall inside the timing precision range of the IPTA MSPs. We note that the rms timing residuals obtained from the individual GWB timings are (on average) 1.5 times less than the rms timing residuals obtained from the GSB+GWB timing. Therefore, considering solely the GWB timing precision makes them even more promising candidates for PTA experiments. Moreover, Tarafdar et al. (2022; InPTA) reported the rms timing residuals of 14 well-timed PTA MSPs (observed with uGMRT) from GWB timing observations for a span of ∼4 yr. The timing precision for these bright PTA MSPs at the uGMRT ranges between 0.76 and 23.4 μs (Figure 13; bottom panel). This clearly illustrates that the timing precision of the four GMRT-discovered MSPs (based on GWB observations) is well within the range of PTA MSPs over a similarly sized timing span.

In this work, we compared the timing results of four GMRT-discovered MSPs using GWB (spanning 3.1–5.5 yr), GSB (spanning 1.8–6.1 yr), and GSB+GWB (covering 7.7–11.0 yr) combined observations. The rms of timing residuals obtained from the GSB+GWB timing is on average ∼1.6 times higher than the individual GWB timings but ∼1.8 times (on average) lower than the residual's rms obtained from the GSB timing. The fitted model parameters in the GSB+GWB timing are on average five and four times more precise than those of the individual GWB and GSB timings, respectively. The model parameters from the GSB+GWB timing are consistent with the individual GWB and GSB timings within ±1σ, where σ is the error in model parameter obtained from the individual timing of the GWB or GSB observations.

We presented DM variations for the four MSPs derived from GWB observations. We find that the line of sight for J1120−3618 crosses an electron-rich medium. We see an ∼9σ_DM increase in the DM value for this MSP during the GWB observational period of 4.2 yr, where σ_DM is the median DM uncertainty. For the remaining three MSPs, the change in DM values is contained within ±3σ_DM.

We detected proper motion for the first time for J1120−3618 (6.0(6) mas yr⁻¹) and J2144−5237 (9.8(5) mas yr⁻¹). For J1828−0625, we detect a proper motion of 14.4(2) mas yr⁻¹, which is at least five times more precise than that reported in Bhattacharyya et al. (2022). The transverse velocities for the three MSPs, calculated from these proper motion values and the distance estimates derived from NE2001 and YMW16, fall within the range expected for MSPs (e.g., Hobbs et al. 2005).

Finally, we compared the timing precision achieved for the four GMRT MSPs to those of 65 IPTA MSPs reported in the second IPTA data release (Perera et al. 2019) and with 14 MSPs reported in first data release of InPTA (Tarafdar et al. 2022). The GMRT MSPs fall within the timing precision range of the IPTA MSPs and InPTA-observed MSPs, making them potential candidates for PTA experiments. The decade-long timing data of the four GMRT MSPs utilized in this study may be useful in the ongoing global effort to aid in the detection significance of the common correlated signal recently detected by multiple PTAs.

We gratefully acknowledge the Department of Atomic Energy, Government of India, for its assistance under project No. 12-R&D-TFR-5.02-0700. The GMRT is run by the National Centre for Radio Astrophysics—Tata Institute of Fundamental Research.