GMRT DISCOVERY OF PSR J1544+4937: AN ECLIPSING BLACK-WIDOW PULSAR IDENTIFIED WITH A FERMI-LAT SOURCE

The Large Area Telescope (LAT; Atwood et al. 2009) on board the Fermi Gamma-ray Space Telescope has discovered a large number of γ-ray point sources, of which many are unidentified or even unassociated with any known potential counterpart (Ackermann et al. 2012). The LAT can localize most of these sources well enough that they can be covered in a single pointing with the large primary beams of radio telescopes at low frequencies, allowing them to be searched efficiently. Targeted radio searches of unassociated LAT point sources by the Fermi Pulsar Search Consortium (PSC) have resulted in the discovery of 43 radio millisecond pulsars (MSPs; Ray et al. 2012). MSPs are thought to evolve from normal pulsars in binary systems via transfer of angular momentum from companions. Thus, the majority of MSPs are naturally expected to be in binaries (∼83% being the binary fraction for MSPs in the Galactic field¹⁵). Binary systems where the pulsar wind evaporates the companion are one way to form isolated MSPs. Such systems where the interaction is ongoing are called black-widow (BW) pulsars. Many exhibit long eclipses (∼10% of the orbital period, apparently larger than the companion's Roche lobe) that are believed to be caused by the material blown from the very low mass companion (M_c ≪ 0.1M_☉) by the pulsar wind. There were two such eclipsing BW systems in the Galactic field known before the launch of Fermi—PSR B1957+20 (Fruchter et al. 1988) and PSR J2051−0827 (Stappers et al. 1996). The BW pulsars are found to have higher values of spin-down energy-loss rate ($\dot{E} \sim 10^{34}$ erg s⁻¹) compared to other MSPs, making these systems good candidates for pulsed γ-ray emission (Roberts 2011). Among 43 new MSPs found in Fermi-directed searches, there are at least 10 BWs (Ray et al. 2012). This Letter describes the discovery and follow-up study of an eclipsing BW MSP, J1544+4937, with the Giant Metrewave Radio Telescope (GMRT).

As a part of the PSC search effort, we observed mid- and high-Galactic-latitude unassociated Fermi point sources with the GMRT at 607 MHz. The GMRT Software Back-end (Roy et al. 2010) produces simultaneous incoherent and coherent filter-bank outputs of 512 × 0.0651 MHz sampled every 61.44 μs. The wider incoherent beam of the GMRT (40' at 607 MHz) can easily cover error circles associated with the Fermi sources. In addition, a coherent beam that is three times more sensitive and narrower (15 at 607 MHz using the central core of the GMRT) can be useful if the pulsar happens to be near the pointing center.

One of the targets was Fermi J1544.2+4941, a γ-ray source from an unpublished internal source list created by the LAT Collaboration using 18 months of data in preparation for the 2FGL catalog (Nolan et al. 2012). The source location (J2000) from that analysis (and used for our telescope pointing) was R.A. = 236074, decl. = 49695, with a 95% confidence error circle of radius 95. This source is very weak, with a likelihood test statistic (TS; Mattox et al. 1996) of 26.2 in the 18 month analysis, and did not make the significance cut to be included in the 2FGL catalog itself.

We processed the data on an IUCAA HPC cluster with Fourier-based acceleration search methods using PRESTO (Ransom et al. 2002). We investigated trial dispersion measures (DMs) ranging from 0 pc cm⁻³ up to 350 pc cm⁻³. A linear drift of up to 200 Fourier-frequency bins for the highest summed harmonic was allowed. The powerline, 50 Hz, and its subsequent harmonics were excised. Using parameters of 32 MHz bandwidth, 10% duty-cycle, incoherent array gain of 2.3 K Jy⁻¹, for 30 minutes of observing, we estimate the search sensitivity as (92 K + T_sky)/(335 K) mJy for a 5σ detection at 607 MHz. Considering |b| > 5°, where T_sky ∼ 10–45 K, our search sensitivity is 0.3–0.4 mJy.

In a 30 minute pointing on 2011 February 1, toward Fermi J1544.2+4941 we discovered a binary MSP of period 2.16 ms with significant acceleration of 2.25 m s⁻² at a DM of 23.2 pc cm⁻³.

We localized J1544+4937 with an accuracy of 5'' (positions listed in Table 1) using continuum imaging for the full GMRT array followed by multi-pixel beamforming (Roy et al. 2012), which allowed us to have sensitive follow-up studies using the coherent array. We estimate a flux of 5.4 mJy at 322 MHz, and a spectral index of −2.3. We started the regular timing campaign for J1544+4937 in 2011 April at 322 MHz with the same coherent filter bank. With the derived position from the multi-pixel search and an a priori binary model predicted by Bhattacharyya & Nityananda (2008), we obtained phase-connected time-of-arrivals (TOAs) from TEMPO,¹⁶ using the JPL DE405 solar system ephemeris (Standish 2004). The binary timing model used is ELL1 (Lange et al. 2001), since J1544+4937 is in a very low eccentricity system. This MSP is in a very compact binary with an orbital period of 2.9 hr. We derive a minimum companion mass (for 90° orbital inclination) of 0.017 M_☉ using the Keplerian mass function, assuming a pulsar mass of 1.4 M_☉. J1544+4937 is eclipsed for about 13% of the orbit at 322 MHz (Section 6). The best-fit timing model (MJD 55680.927–56332.90) is obtained excluding the TOAs around the eclipse phase (0.05−0.35). We achieved a post-fit rms timing residual of 6.9 μs from 652 days of timing (Figure 1). There are still unmodeled residuals, which can be partially absorbed by proper motion fit. However, the inclusion of proper motion reduces the LAT detection significance, indicating that more timing data are required to improve the model. We estimate a precise DM equal to 23.2258(11) pc cm⁻³ by a timing fit using 322 and 607 MHz TOAs from non-eclipsing binary phases. Ephemeris, position, and derived parameters are listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1. Parameters of J1544+4937

Parameter	Valuea
Interferometric positionb
Right ascension (J2000) ...	15h44m04166 ± 03
Declination (J2000) ...	+49°37'5745 ± 47
Offset from survey beam center ...	43
Parameters from radio timing
Right ascension (J2000) ...	15h44m0448722 (2)
Declination (J2000) ...	+49°37'552545 (2)
Position epoch (MJD) ...	51544.0
Pulsar period, P (ms) ...	2.15928839043289 (5)
Pulsar frequency, f (Hz) ...	463.11553585462 (1)
Frequency derivative, $\dot{f}$ (Hz s−1) ...	−6.29 (1) × 10−16
Period epoch (MJD) ...	56007.0
Dispersion measure, DM (cm−3 pc) ...	23.2258 (11)
Binary model ...	ELL1
Orbital period, Pb (days) ...	0.1207729895 (1)
Projected semi-major axis, x (lt-s) ...	0.0328680 (4)
Epoch of ascending node passage, TASC (MJD) ...	56124.7701121 (2)
Span of timing data (MJD) ...	652
Number of TOAs...	280
Post-fit residual rms (μs) ...	6.9
Reduced chi-square ...	2.7
Derived parameters
Mass function, f (M☉) ...	0.0000026132
Min companion mass, mc (M☉) ...	0.017
DM distancec (kpc) ...	1.2
Flux density at 322 MHz (mJy) ...	5.4
Flux density at 607 MHz (mJy) ...	1.2
Spectral index ...	−2.3
Surface magnetic field, Bs (108 G) ...	0.805 (1)
Spin-down luminosity, $\dot{E}$ (1034 erg s−1) ...	1.150 (8)
Characteristic age, τ (Gyr) ...	11.65 (3)
γ-ray parametersd
Photon flux (>0.1 GeV, cm−2 s−1) ...	1.6(8) × 10−9
Energy flux (>0.1 GeV, erg cm−2 s−1) ...	2.1(6) × 10−12
Luminosity, Lγ/fΩ (>0.1 GeV, erg cm−2 s−1) ...	3.6 × 1032
Efficiency, ηγ/fΩ (>0.1 GeV) ...	0.03

Notes. ^aErrors in the last digit are in parentheses. ^bRoy et al. (2012). ^cCordes & Lazio (2002). ^dPhase averaged.

Download table as:  ASCIITypeset image

The radio timing position of J1544+4937 is 43 from the LAT localization of Fermi J1544.2+4941. This is well within the radius of the 95% confidence error circle, suggesting an association between the pulsar and the LAT source, which we sought to confirm with detection of pulsations. In addition, the γ-ray detectability metric $\dot{E}^{1/2}$/d² for J1544+4937 is 7.4 × 10¹⁶ erg^1/2 kpc⁻² s^−1/2, which is comparable to other γ-ray detected MSPs (Figure 12 of Abdo et al. 2010).

Because of the low significance of the LAT source, a phase-averaged spectral analysis of the region is required to optimize the sensitivity of the H-test for pulsed significance using photon probability weighting (Kerr 2011). We performed a binned likelihood spectral analysis of the region using the Fermi Science Tools.¹⁷ We selected LAT data with reconstructed energies in the range 100 MeV to 20 GeV collected between 2008 August 4 and 2013 February 6. We used the P7SOURCE_V6 instrument response functions and excluded events at a zenith angle >100° and times when the LAT was in the SAA or the rocking angle was greater than 52°. We modeled¹⁸ a region of radius 15° using a source list from the 2FGL catalog with an additional source added to the location of the pulsar. We modeled the pulsar with an exponentially cutoff power-law spectrum. Because it is so faint, we fitted only for the power-law index and normalization, keeping the cutoff energy fixed at 2.5 GeV, a value typical of other MSPs (Abdo et al. 2010). The spectral model included both the isotropic and Galactic diffuse contributions with normalizations free in the fit. We detected the source with a TS of 38.2, corresponding to ∼6σ. The spectral index and cutoff energy are very poorly constrained for this weak source, so to estimate the uncertainty on the integrated flux we repeated the spectral fit with the cutoff frozen at 0.7, 2.5, and 5.0 GeV, and with the index frozen at 1.3 (the average value for LAT-detected MSPs) and the cutoff free. We find that the uncertainty due to the spectral shape is about 3 × 10⁻¹⁰ cm⁻² s⁻¹ for the photon flux and 2 × 10⁻¹³ erg cm⁻² s⁻¹ for the energy flux. These are smaller than the statistical errors for this faint source and are added in quadrature to compute the errors on the fluxes reported in Table 1. Using the fitted energy flux (G) and the DM distance (d), we compute the γ-ray luminosity L_γ = 4πf_ΩGd² and efficiency $\eta = L_\gamma /\dot{E}$, where the beaming factor f_Ω is assumed to be 1 (Waters et al. 2009).

Using this model for the continuum emission from the region, we search for pulsed γ-ray emission by selecting events within 2° of the source and computed the probability that they originated from the source with the Fermi Science Tool gtsrcprob. Using these probabilities as weights and pulse phases computed using the radio timing model, we computed a weighted H-test of 37.1, corresponding to a detection significance of 5.1σ and confirming that J1544+4937 is indeed a γ-ray pulsar and identifying the LAT source as this MSP. The LAT weighted light curves in two energy bands, phase aligned to the radio profiles, are shown in Figure 2. For this faint source, it is difficult to determine the peak multiplicity, but the γ-ray emission seems to be mainly between phases 0.4 and 0.6, relative to the single radio peak. Fitting a single Gaussian peak to the light curve gives a peak at phase 0.54 ± 0.02 with a width of 0.25 ± 0.05.

Zoom In Zoom Out Reset image size

We analyzed five Swift/X-ray Telescope (XRT) observations (totaling 9.3 ks) of this field from 2013 March 19–24. No X-ray counterpart is detected, providing a flux upper limit of ∼1.7 × 10⁻¹⁴ erg cm⁻² s⁻¹ for 0.3–10 keV.

We have checked for an optical counterpart for J1544+4937 using archival SDSS data. This field was imaged (Adelman-McCarthy et al. 2008) on MJD 52046.294, at binary phase 0.57. No optical counterpart was detected. The strongest upper limits were g' > 22.5, r' > 22.0, and i' > 21.5 (90% CL). We also obtained a 300 s Hα image with the MiniMo camera on the 3.5 m WIYN telescope on MJD 55975.521 (S. Brownsberger & R. W. Romani 2013, in preparation). No stellar counterpart or extended emission was seen with an Hα flux upper limit of 2.5 × 10⁻¹⁶ erg cm⁻² s⁻¹ (90% CL) for a compact <5'' nebula.

Breton et al. (2013) have obtained optical photometry of the heated companions of a number of Fermi MSPs. They find a typical heating efficiency η ∼ 0.15. For systems with a very strong wind η approaches unity (e.g., J1810+1744, Breton et al. 2013; J1311−3430, Romani et al. 2012). Breton et al. (2013) give an estimate for the irradiation temperature $T_{{\rm Irr}}=(\eta {\dot{E}}_{{\rm SD}}/4\pi \sigma a^2)^{1/4}$, which is 4370 K considering η = 0.15 and 6360 K considering η = 0.6 for J1544+4937. If we assume that the deep, variable radio eclipses imply a near Roche-lobe filling companion, and consider a typical inclination angle = 60°, an unheated (backside) companion temperature of 2500 K and η = 0.15, we can use the binary light-curve synthesis program "Icarus" (Breton et al. 2012) to predict magnitudes at phase ϕ_B = 0.57 of g' ≈ 25.2, r' ≈ 23.9, and i' ≈ 23.1 (magnitudes at maximum are 23.0, 21.9, and 21.5, respectively). Thus, the observed SDSS limits are not constraining. However, if the irradiation efficiency is higher, the fluxes can be detectable; for example, with η = 0.6 we expect g' ≈ 22.4, r' ≈ 21.8, and i' ≈ 21.6 at ϕ_B = 0.57. These are comparable to our observed magnitude limits, so higher efficiencies are ruled out unless the Roche-lobe filling factor is small or the source distance is larger. These estimates also depend weakly on the observer inclination and the secondary unheated temperature and composition.

Our Hα limit for J1544+4937 corresponds to <0.15 of the highest surface brightness ∼5'' patch of the B1957+20 bow shock. The average flux ratio is even more constraining, with an upper limit of ∼0.03 of the total bow-shock flux. However, given the small fraction of MSPs that show optical bow shocks (likely due to the small filling factor of the neutral interstellar medium), the non-detection at such a large distance from the Galactic plane (b ∼ 50°) is not surprising.

Figure 3 presents the timing residuals (and electron column densities) around the eclipse phase at 322 and 607 MHz, with 90 s time resolution. The effect of the eclipses is generally seen from 0.18 to 0.31 orbital phase (eclipse zone hereafter) at 322 MHz. The eclipses are centered at binary phase 0.24 with a duration of around 22 minutes. We estimate the radius of the companion's Roche lobe, R_L (Eggleton 1983),

$$\begin{equation} {R_{L}}=\frac{0.49 a q^{2/3}}{0.6q^{2/3}+ \ln (1+q^{1/3})} \sim 0.13\, R_{\odot }, \end{equation} \tag{ 1 }$$

where q = m_c/m_p is the mass ratio of the companion and the pulsar, and a is the separation of the companion from the pulsar (a ∼ 1.2 R_☉ for J1544+4937, indicating an extremely compact binary). The opaque portion of the companion's orbit is 0.98 R_☉, much larger than R_L of 0.13 R_☉, so the volume occupied by the eclipsing material is well outside the companion's Roche lobe, and thus is not gravitationally bound to the companion. This confirms that this binary is a BW, where the pulsar is ablating its companion, creating a significant amount of intrabinary material that obscures the pulsar's emission.

Zoom In Zoom Out Reset image size

Our sample consisted of six eclipses at 322 MHz where the pulsar emission was fully obscured by companion and its wind and two 607 MHz observing sessions covering the full orbit (Figure 3). At 607 MHz we detect the MSP throughout the 322 MHz eclipse zone. However, we observe a flux fading at 607 MHz near the superior conjunction (orbital phase 0.24). Significant delays in pulse arrival times are observed at 607 MHz during the eclipse zone at 322 MHz. Maximum delay in pulse arrival time at 607 MHz near the eclipse superior conjunction is around 300 μs, which corresponds to an increase in DM of 0.027 pc cm⁻³ and an added electron density, N_e of 8 × 10¹⁶ cm⁻².

Thompson et al. (1994, T94 hereafter) elucidate a collection of eclipse mechanisms. According to them, eclipsing due to refraction of the radio beam demands an order of magnitude higher group delay (∼few tens of milliseconds) than we observe for J1544+4937 (250 μs at 322 MHz egress). Using N_e ∼ 8 × 10¹⁶ cm⁻² observed during superior conjunction at 607 MHz and absorption length about twice the size of the eclipse zone, according to Equation (11) of T94, we find that free–free absorption is possible (absorption optical depth τ_ff > 1) if the plasma temperature T ⩽ 4 × f_cl^2/3 K, where f_cl^2/3 is the clumping factor. This demands either a very low temperature or a very high value of the clumping factor, both of which are not physically achievable. Eclipsing by pulse smearing (due to the increase of N_e along the line of sight) can be ruled out, as the excess electron column density inferred from 607 MHz predicts 373 μs smearing of pulses at 322 MHz near superior conjunction (considering incoherent dedispersion), which is less than one-fifth of pulse period. Since J1544+4937 has a narrow main pulse, and no significant profile evolution is apparent at the eclipse boundary, pulse broadening due to scattering can reduce the detectability but cannot explain the eclipse. In addition, since J1544+4937 is relatively weak, nearby, and has a shallower spectrum than B1957+20, the expected induced Compton scattering optical depth is much less than one (Equation (26) of T94). Another eclipse mechanism considered by T94 is cyclotron–synchrotron absorption of the radio waves by non-relativistic/relativistic electrons, which requires a magnetic field in the vicinity of the companion. We calculate a magnetic field B ∼ 11 G and corresponding cyclotron absorption frequency ∼31 MHz (Equations (35) and (37) of T94, assuming a moment of inertia of 10⁴⁵ g cm²). Thus, 322 and 607 MHz will correspond to 10th and 20th harmonics of the cyclotron resonance. For a fixed temperature, the optical depth for cyclotron absorption drops with harmonics, which may explain the lack of absorption seen at 607 MHz. A modified model is proposed by Khechinashvili et al. (2000) based on kinematic treatment of cyclotron damping, assuming white-dwarf companions with reasonably strong surface magnetic fields. Further observations over a wider radio spectrum may help to investigate the frequency-dependent degree of damping predicted by this model.

We observe a temporal variation of the ingress phase at 322 MHz. For two eclipses, the ingress phase is shifted to 0.16 (from 0.18 for the other four eclipses), whereas there are no apparent shifts in the egress phase. This corresponds to an increase of the opaque portion by 0.15 R_☉ and N_e > 1.2 × 10¹⁶ cm⁻² at the ingress boundary. Such asymmetric increase of eclipse duration may indicate that our line of sight is probing a wind zone where there is systematic outflow of eclipse material.

We also observe strong phase modulations and additional short-duration absorptions at ingress and egress, in time-series data at higher resolution. The durations of these features are in general around 10–20 s, and hence they are not seen in Figure 3. However, in one of the observing epochs these modulations lasted longer—phase modulation of duration 100 s, followed by a short-duration absorption of ∼180 s, then regular emission resumes for 500 s, after which the eclipse starts (Figure 4). Fragmented blobs of plasma randomly oriented around an eclipsing zone, obscuring radiation from the pulsar, can explain these short-duration absorptions.

Zoom In Zoom Out Reset image size

We report the GMRT discovery of an eclipsing BW pulsar, at the position of an unassociated LAT source, Fermi J1544.2+4941. This is the first Galactic field MSP discovered at the GMRT. The detection of pulsed γ-rays from this pulsar demonstrates it as the source powering Fermi J1544.2+4941. Due to the limited significance of the source in γ-rays, additional data are needed before conclusions on peak multiplicity and system geometry can be drawn. The implied efficiency (∼3%) of converting spin-down energy into γ-rays is typical of LAT-detected MSPs. This is the first discovery of a radio MSP in a LAT source fainter than the 2FGL catalog limit (Ray et al. 2012). Since the radio pulsar is relatively bright, this provides strong justification to continue these searches as new LAT sources are revealed in analyses of longer data sets. The radio flux is uncorrelated with the γ-ray flux (Ackermann et al. 2012), so even faint new LAT sources can harbor bright radio MSPs.

Eclipsing BW pulsars have the potential of providing information on the evolutionary connection between the low-mass X-ray binaries and isolated MSPs. With long-term monitoring of this pulsar, we aim to estimate $\dot{P}_{b}$ and its higher derivatives, which can provide an estimate of the lifespan of the system. Bates et al. (2011) noted that for BW systems the measured value of $\dot{E}$/a² is an order of magnitude higher than for other MSP binaries, indicative of greater energy flux needed to ablate the companion. For J1544+4937 we calculate $\dot{E}/a^2 \sim 1.5 \times 10^{33}$ erg lt-s⁻² s⁻¹, which is similar to other BW systems in the Galactic field. Dual frequency observations presented in this paper suggest that cyclotron absorption by the plasma formed via interaction of the pulsar wind with ablated material can obscure the pulsed emission. However, exploring the radio spectrum on either side to probe the reduced/increased opaqueness of the stellar wind during the eclipse phase may provide better insight into the plausible eclipse mechanism.

The Fermi-LAT Collaboration acknowledges support from a number of agencies and institutes for both development and the operation of the LAT as well as scientific data analysis. These include NASA and DOE in the United States, CEA/Irfu and IN2P3/CNRS in France, ASI and INFN in Italy, MEXT, KEK, and JAXA in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the National Space Board in Sweden. Additional support from INAF in Italy and CNES in France for science analysis during the operations phase is also gratefully acknowledged. We acknowledge the support of the telescope operators of the GMRT, which is run by the National Centre for Radio Astrophysics of the Tata Institute of Fundamental Research. We thank the Swift team at the Pennsylvania State University, especially Abe Falcone. We acknowledge help of C. Cheung in interpreting the XRT data. We thank D. Thompson and T. Johnson for their comments and R. Breton for a discussion of heating fluxes.