Serendipitous Discovery of Three Millisecond Pulsars with the GMRT in Fermi-directed Survey and Follow-up Radio Timing

We report the discovery of three millisecond pulsars (MSPs): PSRs J1120−3618, J1646−2142, and J1828+0625 with the Giant Metrewave Radio Telescope (GMRT) at a frequency of 322 MHz using a 32 MHz observing bandwidth. These sources were discovered serendipitously while conducting the deep observations to search for millisecond radio pulsations in the directions of unidentified Fermi Large Area Telescope (LAT) γ-ray sources. We also present phase coherent timing models for these MSPs using ∼5 yr of observations with the GMRT. PSR J1120−3618 has a 5.5 ms spin period and is in a binary system with an orbital period of 5.6 days and minimum companion mass of 0.18 M ⊙, PSR J1646−2142 is an isolated object with a spin period of 5.8 ms, and PSR J1828+0625 has a spin period of 3.6 ms and is in a binary system with an orbital period of 77.9 days and minimum companion mass of 0.27 M ⊙. The two binaries have very low orbital eccentricities, in agreement with expectations for MSP-helium white dwarf systems. Using the GMRT 607 MHz receivers having a 32 MHz bandwidth, we have also detected PSR J1646−2142 and PSR J1828+0625, but not PSR J1120−3618. PSR J1646−2142 has a wide profile, with significant evolution between 322 and 607 MHz, whereas PSR J1120−3618 exhibits a single peaked profile at 322 MHz and PSR J1828+0625 exhibits a single peaked profile at both the observing frequencies. These MSPs do not have γ-ray counterparts, indicating that these are not associated with the target Fermi LAT pointing emphasizing the significance of deep blind searches for MSPs.


Introduction
Millisecond pulsars (MSPs) are rapidly rotating neutron stars (rotational periods of a few tens of milliseconds) with very small spin-down rates (P 10 18 < - ). They are thought to have acquired their high rotational rate by accretion of matter, and thereby transfer of angular momentum, from a binary companion (Alpar & Cheng 1982;Bhattacharya & van den Heuvel 1991), which is now supported by observational evidence (Archibald & Stairs 2009;Papitto et al. 2013;Roy et al. 2015). MSPs are still a small population compared to the normal pulsars and much diversity in their intrinsic characteristics as well as evolutionary history are yet to be explored. About 15% of the ∼3200 known pulsars are MSPs, either in the Galactic disk or in globular clusters. 11 Even though estimates for the Galactic population of MSPs range from 40,000 to 90,000 objects (Faucher-Giguere & Kaspi 2006), MSP discoveries are hindered by their intrinsic faintness, multipath propagation of the radio signals (known as "scattering"), and dispersive smearing of the radio pulses, particularly at low radio frequencies. In addition, the majority of these objects are in binary systems; for the most compact of them, the orbital motion will cause a variation of their spin periods, which makes the discovery process for binary MSPs more difficult and computer-intensive. Pointed surveys allow the concentration of telescope time and computing resources in particular areas of the sky that are more likely to yield new discoveries.
A reasonable fraction of MSPs have large enough spin-down luminosities and small enough distances to be detectable in high-energy γ-ray data from the Fermi Large Area Telescope (LAT; 12 Atwood et al. 2009) in areas of the sky where there were previously no sources known. Some of these MSPs were found in searches for pulsations in the arrival times of the gamma-ray photons (e.g., Clark et al. 2018). For a few others, the identification and precise location of companions at optical wavelengths, together with precise orbital period information obtained from their brightness variations, has allowed the determination of P, P  , and the orbital parameters from the γ-ray photons (e.g., Pletsch et al. 2012;Nieder et al. 2020). However, the majority of MSPs are in binary systems without optically identifiable companions. In this case, the search for pulsations in γ-rays becomes prohibitively expensive, given the vast parameter space (for positional, spin, and orbital parameters) that must be searched. In these cases, the best hope for determining these parameters is the discovery of an associated radio pulsar. Pointed surveys of these unassociated Fermi sources with radio telescopes, which are coordinated by the Fermi Pulsar Search Consortium (PSC; Ray et al. 2012) have resulted in the discovery of 110 radio MSPs so far. The experience to date is that, because of the steep radio spectra of most MSPs, observing capabilities at lower radio frequencies (300-600 MHz) are the optimal choice to look for nearby undiscovered MSPs, especially at high Galactic latitudes, where the radio background and dispersive smearing and scattering are not so large. Because these pulsars were missed in earlier surveys, they will necessarily be faint, hence the need for a sensitive telescope.
The Giant Metrewave Radio Telescope (GMRT 13 ) is a multielement aperture synthesis telescope consisting of 30 antennas each of 45 m diameter, having a maximum baseline length of 25 km (Swarup et al. 1997). Among the PSC, the GMRT is one of the few with low frequency receivers, being one of the most sensitive instruments below 1 GHz. Its capabilities for pulsar searches have been demonstrated by the discovery of 36 pulsars in targeted and blind searches at an encouraging pulsar-per-square degree discovery rate (e.g., Ray et al. 2012;Bhattacharyya et al. 2013Bhattacharyya et al. , 2016Bhattacharyya et al. , 2019. Being a large array of small telescope pulsar searches with the GMRT has multiple advantages like wide field of view with incoherent beam (FWHM of 80′ at 322 MHz, and 40′ at 607 MHz); high sensitivity coherent beam (four to five times incoherent array beam; good for follow-up timing observations), and rapid precise localization (∼10″) using the imaging capability, even on search observations .
This sensitivity warranted a search for new pulsars in unassociated Fermi γ-ray point sources with the GMRT. This survey has, to date, discovered four MSPs, which show γ-ray pulsations . In this work we present the discovery of three new MSPs in this GMRT survey. Even though all three MSPs were discovered in error circles of LAT sources from the 1FGL catalog (Abdo et al. 2010), none have associated sources in the 4FGL−DR3 catalog (an incremental update to the 4FGL catalog using 12 yr of LAT data; Abdollahi et al. 2020;Abollahi et al. 2022). This suggests that the MSPs are likely unassociated with the original Fermi γ-ray point sources targeted in the survey, i.e., they are serendipitious discoveries. As such, they give an indication of the number of MSPs we might discover if a survey of the depth described by  were extended to the whole sky.
In Section 2, we provide details on the discoveries. Results from follow-up timing studies of the discovered pulsars are reported in Section 3. Section 4 presents a discussion on serendipity of the discoveries. Finally Section 5 presents a summary of the paper.

Discovery of Three MSPs
As mentioned above, the discoveries described in this paper were made by a GMRT survey of Fermi γ-ray point sources that had no association with any previously known astronomical objects. These sources were observed with the GMRT at 322 and 607 MHz with a 32 MHz observing bandwidth with GMRT software backend (Roy & Gupta 2010). Details on the survey strategy, target list, observations, and data processing are described in Bhattacharyya et al. (2013. The properties of the discoveries presented in this work are listed in Table 1. We estimated the pulsed flux density of the pulsars at the discovery epoch using the radiometer equation for an incoherent array with 1.6 K Jy −1 gain and by extrapolating the sky temperature at 322 MHz from the all-sky 408 MHz image of Haslam et al. (1982). Thus flux density estimation is done in the same way as reported in . Pulse profiles for these MSPs at 322 and 607 MHz are presented in Figure 1.
PSR J1120−3618 was discovered in a 30 minute observing run with the GMRT at 322 MHz. This pulsar was not detected with the GMRT 607 MHz observations at multiple epochs (4 epochs each of 1 hr duration), which could be due to a relatively steep spectra; we estimate a limiting spectral index of <−3.0. We also note that scintillation might have played a major role in the nondetection as was also observed by Camilo et al. (2015) for the MSPs discovered with the Parkes in the Fermi-directed surveys. With ongoing wide bandwidth observations using the upgraded GMRT (Gupta et al. 2017) we will be able to probe this in detail and will report that in future publications. The profile of this MSP is single peaked at 322 MHz. There is a hint of two different components but due to dispersion smearing the profile components are not well resolved at 322 MHz. The nearest γ-ray source, 4FGL J1117.7−3650, is 0°. 73 away, with 95% confidence-level semimajor and semiminor error-ellipse axes of 0°. 0895 and Notes. The numbers in the parentheses are uncertainties in preceding digit. a Ray et al. (2012). b Flux density at 322 MHz without primary beam correction. c Flux density at 607 MHz without primary beam correction. d Error on spectral index is calculated considering a typical 10% uncertainty in the flux measurement. e 5σ nondetection limit at 607 MHz for 30 minutes of GMRT observations using 17 antennas in phased array.
0°. 0633, respectively. The 4FGL−DR3 source is associated with NVSS J111758−364918. Incidentally another MSP, J1124−3653 was discovered while targeting the same Fermi sources with the Green Bank Telescope (GBT) (P. Bangale et al. 2021, in preparation) and was also detected by us. PSR J1124−3653 is the true counterpart of the gamma-ray source (1FGL J1124.4−3654, now 4FGL J1124.0−3653) that was originally targeted by our survey, while PSR J1120−3618, the subject of this work, is unrelated to either of the 4FGL sources. PSR J1646−2142 was discovered in a 60 minute observing run with the GMRT at 322 MHz. We observe significant evolution of its average profile between 322 MHz to 607 MHz: The pulse profile at 607 MHz seems to be wider than the 322 MHz profile, which is contrary to the radius-to-frequency mapping generally observed for pulsars (Lorimer & Kramer 2004). We aim to closely investigate the observed change in the profile shape with the wide band observations. The nearest γray source, 4FGL J1646.7−2154, is 0°. 23 away, with 95% confidence-level semimajor and semiminor error-ellipse axes of 0°. 0792 and 0°. 0690, respectively. The 4FGL source does not have an association in the 4FGL−DR3 catalog. PSR J1828+0625 was discovered in a 30 minute observing with the GMRT at 322 MHz. The pulse profile is single peaked for this MSP at both frequencies. The nearest γ-ray source, 4FGL J1830.1 + 0617, is 0°. 43 away, with 95% confidencelevel semimajor and semiminor error-ellipse axes of 0°.0437 and 0°.0350, respectively. The 4FGL−DR3 source is associated with TXS 1827 + 062.
We localized the newly discovered MSPs in the image plane with the GMRT interferometric array with an accuracy of approximately ±10″ (half of the typical synthesized beam used in the image made at 322 MHz) using a coherently dedispersed gated imaging beam former ( Table 2; . Once the MSPs are localized in the image plane, we use the sensitive coherent array (four to five times more incoherent array for the GMRT) for follow-up observations with a smaller field of view but with enhanced sensitivity. We determine a 10σ detection significance of 0.3 mJy using the coherent array with the central core of the GMRT having 17 antennas (i.e., gain of ∼7 K Jy −1 ).

Timing Study
The precise astrometric localization of these three MSPs derived with the GMRT interferometric array allowed us to conduct dense follow-up observations of these sources with the coherent array of the GMRT. For J1646−2142 and J1828 +0625 these observations were made with the 322 and 607 MHz receivers, for PSR J1120−3618 we used only the 322 MHz receiver. The highest signal-to-noise ratio profiles, obtained from dedispersing and folding the data (see again details in , were used as templates for extracting pulse times of arrival (TOAs). We used the standard pulsar timing software TEMPO 14 to derive timing models that can describe the spin and (in two cases) orbital motion with enough precision to fold the data. After an initial dense campaign, which was important for establishing an initial timing model, we started relatively sparse observing of the MSPs, in order to extend these timing models. These lowcadence observations eventually extended to about 5 yr.
Using the manual method described in Section 3 of Freire & Ridolfi (2018), we were able to obtain a phase-connected timing solution for the isolated pulsar in our sample, PSR J1646−2142. For the others, the low flux densities and relatively blunt profiles resulted in a rather low timing precision, which makes the determination of a timing solution difficult, especially given the fact that we must also fit for orbital parameters. For this reason, we used an earlier version of the phase-connection algorithm proposed in Section 4 of Figure 1. Radio pulse profiles of PSR J1120−3618, J1646−2142, and J1828+0625 from the GMRT observations. The blue curves show the 322 MHz radio profiles, and the red dashed curves show the 607 MHz radio profiles, when available with a bandwidth of 32 MHz using the GMRT legacy system. PSR J1120−3618 is not detected with 607 MHz observations with the GMRT. Profiles are broadened due to significant dispersion smearing, which was left unaccounted with the incoherent dedispersion. Freire & Ridolfi (2018). 15 Similar methods were used, for the same reason, for the three MSPs found by .
The timing solutions thus derived are presented in Table 2. The parameters are presented in Dynamic Barycentric Time; we used the DE440 solar system ephemeris (Park et al. 2021). The orbital motion of our two binaries is described using the ELL1 orbital model (Lange et al. 2001); this has the advantage of eliminating the correlation between two Keplerian parameters, the longitude of periastron (ω) and the time of passage through periastron (T 0 ); this correlation is especially large for low-eccentricity orbits such as those of these two pulsars. The orbits are described instead by the time of passage through the ascending node (T asc ), which can be measured precisely even for circular orbits, and the two Lagrange-Laplace parameters, ò 1 = e sin ω and ò 2 = e cos ω. The disadvantage of this description is that it is not exact: The geometric orbital delay is expressed in successive terms proportional to xe n , where n is an integer. Recently, the TEMPO implementation of this model has been updated to include terms of order xe 2 (Zhu et al. 2019). These models can describe adequately all binaries where the largest ignored term (now of amplitude xe 3 ) is much smaller than the uncertainty in x, which includes the vast majority of binary pulsars in the Galaxy and certainly the two new binary MSPs presented in this work. Using the timing solutions in Table 2, we obtain the timing residuals (TOA minus the prediction of the timing solution prediction for that TOA) in Figure 2. The postfit rms of the timing residuals are around 23, 14, and 21 μs, respectively. We aim to reduce these residual rms and improve the timing models with an ongoing study of these MSPs with the more sensitive wide band observations using the upgraded GMRT (Gupta et al. 2017), which has a coherent dedispersion capability.
We now describe some of the timing results in more detail. PSR J1120−3618 is an MSP in a 5.6 days binary system in a low-eccentricity (e = 3.0(4) × 10 −5 ) orbit. Assuming a pulsar mass of 1.4 M e and orbital inclinations of 90°, 60°, and 25°its mass function yields companion masses of 0.18, 0.21, and 0.47 M e , respectively. We have not observed any evidence of eclipsing for this binary system. According to the Tauris & Savonije (1999) model, for the orbital period of this system the mass of the helium white dwarf (WD) companion will be ∼0.2 M e ; this suggests that the orbital inclination is close to the median of 60°.
PSR J1828+0625 is in a wide (P b = 78 days) orbit, which also has a low eccentricity (e = 9.789(12) × 10 −5 ). Using the 322 and 607 MHz observations allowed us to determine a DM of 22.4165 pc cm −3 ; we also detect a slow decrease of the DM with time. For this MSP we also detect significant proper motion of 15 mas per year, translating to a transverse velocity (V T ) of ∼77 km s −1 . With the same pulsar mass and inclination assumptions as above, its mass function yields companion masses of 0.27, 0.31, and 0.73 M e , respectively. No evidence of eclipsing is observed from this wide binary system yet. For its orbital period, the expectation of the Tauris & Savonije (1999) model is a companion mass of ∼0.3 M e ; this suggests that the orbital inclination is also close to the median of 60°.
Finally, our timing results confirm that PSR J1646−2142 is an isolated MSP. One of the aspects of these three MSPs is that their spin period derivatives (P  ) are relatively small when compared to the known MSP population. For this reason, it is important to take into account kinematic contributions to this quantity, in order to obtain an estimate of the intrinsic spin period derivatives, P int  . To do this, we follow the analytic equations and quantities used in Section 7 of Guo et al. (2021). In the case of PSR J1828+0625, we can estimate the intrinsic spindown since we know the proper motion. This spin-down yields a characteristic age (τ c ) of about 25 Gyr. For the other pulsars, we do not yet know the proper motion, so we cannot estimate one of the kinematic effects: the Shklovskii effect (Shklovskii 1970). However, since the latter effect is always positive, we can calculate an upper limit for the P int  , which then yields upper limits for the spin-down energy (E  ) and surface magnetic field (B 0 ), and lower limits for τ c . These estimates are also presented in Table 2. From this, it is clear that all pulsars have characteristic ages in excess of 10 Gyr. If we assume that P int  must be positive, we can derive, from the equation for the Shklovskii effect and the assumed distance d, upper limits on the proper motion of each system: where P int,max  is the upper limit for P int  in Table 2. These are 11 mas yr −1 for PSR J1120−3618 and 22 mas yr −1 for PSR J1646−2142. Such proper motions are not detectable given current timing: The 2σ uncertainties of the total proper motion of both systems are ∼20 and ∼80 mas yr −1 , respectively.

Discussion on Serendipity of Discoveries
The three MSPs reported in this paper were discovered while following up the rank-ordered Fermi γ-ray point sources with the GMRT. In such targeted surveys, specific parts of the sky with prior information of the presence of pulsars are searched deeper with higher sensitivities than what is possible with conventional blind surveys. The lack of a γ-ray counterpart in later Fermi LAT catalogs indicates that these three MSPs are chance discoveries, not associated with the original Fermi targets. Smith et al. (2019) have demonstrated that it is possible to detect γ-ray pulsations from pulsars that are not found as significant point sources. Therefore, we cannot rule out the possibility that gamma-ray pulsations might be detected from these MSPs in a careful analysis of the LAT data, which has not yet been performed, to our knowledge.
As part of Fermi-targeted searches with the GMRT, we followed up 375 unassociated sources listed in the Appendix of . Out of the seven discoveries, four MSPs were associated with the target Fermi pointing reported in Bhattacharyya et al. (2013 whereas the remaining three MSPs reported in this paper are unassociated to the Fermi target pointing. Considering size of the beam as around 1.4 deg 2 , this will translate to 0.005 MSP per square degree discovery rate. This indicates a population of weak MSPs waiting to be discovered with deep enough blind searches. We note the two other MSPs that were serendipitously discovered in Fermi-directed radio searches are J1103−5403 (with the Parkes telescope; Keith et al. 2012) and J1551−0658 (with the GBT; P. Bangale et al. 2021, in preparation). Blind surveys with other telescopes are also discovering MSPs at a very high rate. This is exemplified by the discovery of 25 new MSPs by the GBNCC survey (McEwen et al. 2020). Considering 0.063 pulsar detections per square degree and the discovery of 670 pulsars reported by McEwen et al. (2020), we arrive at a 0.002 MSP per square degree discovery rate for the GBNCC survey. A factor of two higher MSP per square degree discovery rate reported in this paper is most likely due to small number statistics.
As can be inferred from the flux density listed in Table 1, all three MSPs are weak, the strongest one having S 322 ∼ 3 mJy. A lower flux density is observed at ∼400 MHz only for 6% of the known MSPs. Thus it is evident that relatively longer integration time with the Fermi-directed targeted searches reaching up to 0.3 −0.9 mJy at 322 MHz and 0.3−0.4 mJy at 607 MHz with a 32 MHz observing bandwidth of the GMRT in an observing duration of 30 minutes (as detailed in  made the discoveries possible. MSPs are about one order of magnitude less luminous and less efficient radio-emitters compared to ordinary pulsars (Kramer et al. 1998), who commented that relatively faint and far away MSPs are likely being missed by the search efforts. While comparing the luminosity distribution of ordinary pulsars and MSPs within a distance of 1.5 kpc, Kramer et al. (1998) found that the difference is less prominent. Table 2 of  lists ongoing blind and targeted surveys with different radio telescopes. We note that the TRAPUM-UHF 16 survey with the MeerKAT, SUPERB 17 with the Parkes, GPPS 18 with FAST, and GHRSS 19 with the GMRT will be reaching up to 0.2 mJy while surveying the sky blindly for pulsars. Given these surveys, our results imply that there should be an increased rate of discovery of such low luminosity pulsars in the coming decade. However, we note that the corresponding targeted surveys with these telescopes will reach sensitivities up to a factor of 10 deeper; thus probing an even fainter population of pulsars and MSPs.
The uGMRT will also contribute to this effort, especially at the lower frequencies where its sensitivity is highest. It is now equipped with higher instantaneous bandwidth, with effective bandwidth at least four to six times higher than the 32 MHz bandwidth of the legacy GMRT observing system with which the Fermi-directed searches were performed. Thus one will be able to reach similar (if not improved) search sensitivity with the uGMRT with around 10 minutes of integration. Thus along with the sensitive surveys conducted by the other telescopes, the ongoing blind survey (GHRSS survey; Bhattacharyya et al. 2019) reaching sensitivity up to 0.2 mJy and the recently restarted Fermi-directed targeted survey with sensitivity up to 0.1 mJy or better using the uGMRT, promises to discover a good fraction of the weaker MSPs probing an unexplored part of MSP luminosity distributions. This will in turn provide inputs for investigation of the birth history of MSPs (Faucher-Giguere & Kaspi 2006).

Summary
In this paper we report the discovery of three MSPs: PSRs J1120−3618, J1646−2142, and J1828+0625 with the GMRT. Follow-up phase coherent timing for ∼5 yr with the GMRT was conducted for these MSPs. While PSR J1646−2142 is an isolated MSP, the other two MSPs, PSRs J1120−3618 and J1828+0625, are in binary systems, with orbital periods of 5.6 and 77.9 days, respectively. The two binaries are likely to be MSP-helium WD systems. The median companion masses estimated from their mass functions are consistent with the predictions of the Tauris & Savonije (1999) model for their orbital periods. Their orbital eccentricities are also similar to those observed for other MSP-helium WD systems, and are consistent with the prediction of Phinney (1992) for their orbital periods. Finally, their spin periods are also typical of what one finds among the MSPs in this population.
One of the characteristics of the MSPs presented in this work are their small values for P int  , which yield relatively small surface magnetic fields, and very large characteristic ages: PSR J1828 + 0625 has a characteristic age of ∼25 Gyr, for PSR J1646−2142 and J1120−3610, the lower limits are 11.6 and 54 Gyr, respectively. These small P int  are consistent, but at a lower end of the distribution observed among MSP-helium WDs and isolated MSPs; This suggests that these MSPs went through a very extensive and prolonged recycling, and that their spin periods have been little changed since accretion stopped, this is especially true for PSR J1120−3618. These pulsars were discovered while following up the rankordered Fermi γ-ray point sources. However, we conclude that these are chance discoveries and the MSPs are not associated to the originally targeted Fermi γ-ray points sources. The lack of γ-ray point-source detection of these MSPs might be due to their relatively small values E  . About 40% of the pulsars with E d 2  similar to that of PSR J1828 + 0625 are not detected, the number of nondetections increases for lower spin-down powers (Guillemot et al. 2016). Additional information from such pulsars, including detailed searches for gamma-ray pulsations without requiring a point-source detection (Smith et al. 2019), is important for understanding the conditions for the occurrence of γ-ray emission.