A 350 MHz Green Bank Telescope Survey of Unassociated Fermi LAT Sources: Discovery and Timing of 10 Millisecond Pulsars

One of the major results to come out of over a decade of observations with the Fermi Gamma-ray Space Telescope (hereafter Fermi) is that the majority of Galactic high-energy (E ∼ 1 GeV) γ-ray sources so far detected are pulsars. At low Galactic latitudes (b ≲ 10°), it had been suspected since the EGRET era that most discrete sources might be young (τ ≲ 10^5.5 yr), spin-powered pulsars (Roberts 2005). Before the launch of Fermi, the most likely population of Galactic γ-ray sources at higher Galactic latitudes was thought to be nearby (d ≲ 2 kpc) middle-aged pulsars like PSR B1055−52 and Geminga, perhaps associated with the Gould Belt of nearby star-forming regions (Gonthier et al. 2005). Millisecond pulsars (MSPs, which, for the purpose of this paper, are defined as having spin periods P < 10 ms), were considered a possibly important class as well, despite only one marginal detection with EGRET (Kuiper et al. 2000). Radio pulsar surveys targeting EGRET source error boxes (Nice & Sayer 1997; Champion et al. 2005; Crawford et al. 2006) were notoriously unsuccessful, with only three pulsars discovered in these surveys—the MSPs J0751+1807 (Lundgren et al. 1995) and J1614−2230 (Hessels et al. 2005; Demorest et al. 2010), and the young PSR J1028−5819 (Keith et al. 2008)—that were plausibly the counterparts of the γ-ray source targeted (Abdo et al. 2010a). Very deep radio searches of low-latitude extended X-ray sources discovered within EGRET error boxes were somewhat more efficient at discovering potential γ-ray pulsars. Two young and energetic pulsars, PSR J2021+3651 (Roberts et al. 2002a) and PSR J2229+6114 (Halpern et al. 2001) were discovered this way, both of which proved to be powering their respective coincident γ-ray sources (Halpern et al. 2008; Pellizzoni et al. 2009).

With the launch of the Fermi satellite, it was quickly determined that MSPs were efficient γ-ray emitters, with eight MSP detections reported within the first few months of operations (Abdo et al. 2009b). In addition, blind searches of Fermi Large Area Telescope (LAT; Atwood et al. 2009) data soon discovered 16 young pulsars, three of which were at mid-Galactic latitudes (10° < ∣b∣ < 30°; Abdo et al. 2009a). The 46 γ-ray pulsars reported in the First Fermi-LAT Pulsar Catalog (Abdo et al. 2010a) include 15 pulsars at ∣b∣ > 5°. Of these, eight are MSPs, four are young radio pulsars, and three are detected only in γ-rays. It was therefore expected that deep radio searches of unassociated LAT sources at midlatitudes should result in a mix of MSPs and young pulsars.

A sensitive 820 MHz radio search with the 100 m Green Bank Telescope (GBT) of 25 of the brightest LAT sources not associated with previously known sources of potential γ-rays (i.e., energetic pulsars or blazars) discovered no young pulsars. However, out of the eight sources searched at ∣b∣ > 5°, three were found to be binary radio MSPs (Ransom et al. 2011). One of these four, PSR J2214+3000, is in a tight (P_B ∼ 10 hr) orbit around a very-low-mass (minimum companion mass M_c ∼ 0.02 M_⊙) companion and exhibits regular eclipses, becoming at that time only the fourth known "black widow" system discovered in the Galactic field (outside of a globular cluster).

Between the discovery of the first MSP in 1982, the isolated 1.6 ms PSR B1937+21 (Backer et al. 1982), and the launch of Fermi in 2008, psrcat ²¹ (Manchester et al. 2005) lists the discovery of approximately 70 MSPs in the Galactic field. In comparison, since data from Fermi became publicly available in 2009, over 100 MSPs (including the ones reported here) were discovered in the Galactic field by targeting Fermi-LAT sources (Cognard et al. 2011; Keith et al. 2011; Guillemot et al. 2012; Kerr et al. 2012; Ray et al. 2012; Barr et al. 2013; Bhattacharyya et al. 2013; Camilo et al. 2015; Cromartie et al. 2016; Pleunis et al. 2017; Bhattacharyya et al. 2021; Deneva et al. 2021; Wang et al. 2021; Clark et al. 2023a). Surprisingly, only one young γ-ray pulsar has been discovered through targeted radio searches of Fermi sources (Camilo et al. 2012). Surveys such as these are rapidly expanding the known population of MSPs, which may help to understand better the nature of binary and isolated MSPs. Also, continued long-term timing has shown that some of these new MSPs are useful in pulsar timing arrays for gravitational-wave detection (Hobbs et al. 2010; Arzoumanian et al. 2015; Ajello & Atwood et al. 2022). Here we present an overview of the most prolific of the early Fermi targeted searches, reporting the discovery of 10 MSPs and the detection of eight more MSPs that were found in subsequent surveys.

This paper presents our survey results in detail (a brief overview of the new discoveries was published previously in Hessels et al. 2011). Here, we give detailed radio timing solutions for the 10 MSPs discovered in this survey. We also present the detection of pulsed γ-ray emission from nine of these. Six of these were previously reported in the second LAT pulsar catalog (Abdo et al. 2013, hereafter 2PC), and all were included in the third LAT pulsar catalog (Smith et al. 2023, hereafter 3PC). The new discoveries include one isolated MSP and six in close binaries with periods under 1 day. Three of these have very-low-mass companions (M_c ≪ 0.1 M_⊙), and an additional two exhibit extensive eclipses. One of the newly discovered MSPs is a binary which appears to be a chance discovery, well outside the target LAT positional error ellipse. In Section 2 we describe the survey. In Section 3 we describe the subsequent timing program and summarize the survey results as well as describe follow-up radio, X-ray, and γ-ray studies of these sources. In Section 4 we discuss the implications of the new discoveries for different pulsar populations with detailed discussions of the new black widow and redback systems.

Candidate sources for our radio searches were drawn from a preliminary version of the first Fermi-LAT source catalog (1FGL; Abdo et al. 2010b). We considered the fraction of the sky which is visible from the GBT for more than ∼1 hr, or equivalently decl. > −40°. We chose 350 MHz as the observation frequency due to the steep spectra of MSPs (Kramer et al. 1998) and the larger beam size of the GBT at low frequencies. The positional uncertainty of the LAT sources was typically less than the ${35}^{{\prime} }$ FWHM beam of the 350 MHz system. We chose sources well away from the Galactic plane, i.e., ∣b∣ > 5°, where the sky temperature and the effects of interstellar scattering are reduced. We excluded sources that had viable active galactic nucleus (AGN) counterparts (Abdo et al. 2010c) and sources with statistically significant variability as defined in 1FGL, as pulsars had long been recognized as a nonvariable γ-ray population (e.g., McLaughlin et al. 1996; however, see Roberts et al. 2002b; Mayer et al. 2013; Stappers et al. 2014 for discussions of γ-ray variability from pulsar wind and accretion flows; and Allafort et al. 2013 for the discovery of a rare example of γ-ray pulsar variability). We also excluded sources in the Fermi-LAT Bright Source List (Abdo et al. 2009c) since those had largely been surveyed with the GBT at 820 MHz by Ransom et al. (2011). There is no clear relation between γ-ray and radio pulsed flux (3PC), which likely depends to a large degree on beaming geometry, and so targeting a fainter γ-ray population would not be expected to yield a smaller (or larger) fraction of radio pulsar counterparts.

Among the sources that passed these criteria, we generally preferred sources that had an obviously "pulsar-like" spectrum, i.e., a power-law behavior at lower energies but with significant curvature/cutoff at higher energies. However, we did not exclude sources based solely on spectra but balanced spectral desirability with considerations of visibility and observing efficiency within each scheduled observation session. By following this procedure, 49 sources were finally observed throughout 2009. The 1FGL and 3FGL source names, pointing positions, 3FGL γ-ray fluxes, variability indices, observation dates, exposure lengths, and minimum detectable fluxes are listed in Table 1, along with source classifications.

The observations were performed using the 4096-channel 100 MHz bandwidth mode of the (Green Bank Ultimate Pulsar Processing Instrument (GUPPI) backend on the GBT (DuPlain et al. 2008a) with a central frequency of 350 MHz. The integration times were generally 32 minutes, making this about an order of magnitude more sensitive than any prior large-scale survey of the northern sky and 5 to 6 times more sensitive than the Green Bank North Celestial Cap pulsar survey (Stovall et al. 2014), at least for sources near the center of our beam and assuming typical pulsar spectral indices. The sensitivity of the search is obtained by using the modified radiometer equation:

$$\begin{eqnarray}&&{S}_{min}=\displaystyle \frac{{\left(S/N\right)}_{min}\,\left({T}_{{\rm{sys}}}+{T}_{{\rm{sky}}}\right)}{G\sqrt{{n}_{p}{t}_{{\rm{obs}}}\bigtriangleup f}}\sqrt{\displaystyle \frac{W}{P-W}},\end{eqnarray} \tag{ 1 }$$

where (S/N)_min = 8 is the threshold signal-to-noise ratio (S/N) used, G = 2 K Jy⁻¹ is the effective gain of the GBT, n_p = 2 is the number of polarizations summed, △f = 100 MHz is the total observing bandwidth, T_sys = 25 K is the system temperature, T_sky is the sky background temperature, t_obs is the integration time, P is the period of the pulsar, and W is the pulse width. For example, given W = 0.1P, typical for MSPs, a sky background temperature of 50 K, typical for sources off of the Galactic plane at this frequency (Haslam et al. 1982), and an integration time of 32 minutes, the minimum detectable flux is 0.2 mJy. We list the nominal sensitivity toward each pointing position in Table 1, assuming W = 0.1P.

Our actual sensitivity at the beam center was likely somewhat, and perhaps substantially, worse than what would be inferred from the radiometer equation. This equation does not take into account the effect of radio frequency interference (RFI), which results in parts of the bandpass and some intervals in time being unusable. It also assumes full sensitivity across the entire observing bandwidth, whereas, in reality, the sensitivity decreases at the edges of the band, resulting in a smaller effective bandwidth. The assumption of W = 0.1P will not be true for all pulsars. Some may have intrinsically broadened pulse profiles, and pulsars with higher dispersion measures (DMs) may experience appreciable pulse broadening due to interstellar scattering and dispersion. In addition, it is important to note that pulsars with average fluxes greater than the quoted minimum detectable flux could be missed due to interstellar scintillation, which can cause the fluxes of pulsars to vary dramatically from epoch to epoch. Sensitivity is also worse for very-short-orbit binaries (P_b ≲ 5 hr) because we only use an acceleration search over the first-period derivative (see below), which cannot fully correct for the orbital acceleration over a 32 minute observation. And finally, many binary pulsars are in eclipse for a significant fraction of their orbit, and during this survey, we only observed each source one time and so could have missed systems that, on average, are quite bright.

We dedispersed the data over the range from 0 pc cm⁻³ to twice the maximum DM predicted by the NE2001 Galactic electron density model in the direction of each source (Cordes & Lazio 2002); our search was prior to the release of the newer YMW16 electron density model Yao et al. (2017). We performed acceleration searches up to z_max (the maximum Fourier frequency derivative) of 200 to improve sensitivity to short-orbit binaries. Up to eight harmonics were summed in the Fourier power spectra using the standard tools found in PRESTO ²² (Ransom et al. 2002).

In our survey data, we discovered 10 MSPs and detected four more that had recently been discovered by other surveys. A subsequent survey using the GBT at a central frequency of 820 MHz targeted more than 100 LAT sources (Sanpa-Arsa 2016), several of which we had also targeted with our 350 MHz survey. These discoveries are summarized in Table 2. This higher-frequency survey discovered five pulsars within our sample. Subsequent folding of our data using the 820 MHz determined periods and DMs resulted in clear detections of four of these pulsars, although mostly at significance levels below our criteria for a discovery. Of these four, PSRs J0533+6759, J1312+0051, and J1858−2216 have broad pulses, detected with low S/N in our data. The other, PSR J1142+0119, was discovered first in the 820 MHz data. The final 820 MHz discovery, PSR J1921+0137, is not detectable in our survey, likely due to either a flat spectrum, an eclipse, or scintillation.

Table 2. Properties of the Pulsars Observed in This Survey

Name	P	S350	Fγ a	$\dot{E}$ b	DM	Dist. c	PB	Mc d
	(ms)	(mJy)	(10−11 erg cm−2 s−1)	(1034 erg s−1)	(pc cm−3)	(kpc)	(days)	(M⊙)
Pulsars Discovered in Searches of the 350 MHz Data
J0023+0923	3.05	5.4	0.7 ± 0.1	1.3	14.3	1.2	0.14	0.01
J0102+4839	2.96	2.8	1.6 ± 0.1	1.6	53.5	2.3	1.6	0.18
J0307+7443	3.16	1.8	1.6 ± 0.1	2.2	6.3	0.4	37.1	0.20
J0340+4130	3.30	6.5	2.1 ± 0.1	0.7	49.6	1.6	Isolated	N/A
J1124−3653	2.41	2.2	1.3 ± 0.1	1.5	44.9	1.0	0.23	0.02
J1302−3258	3.77	1.3	1.2 ± 0.1	0.4	26.2	1.4	0.78	0.15
J1551−0658 e	7.09	0.7	⋯	0.2	21.6	1.3	5.2	0.20
J1810+1744	1.66	78	2.3 ± 0.1	3.1	39.6	2.4	0.15	0.04
J2129−0429	7.61	2.1	0.8 ± 0.1	2.9	16.9	1.4	0.64	0.37
J2215+5135	2.61	16	1.3 ± 0.1	6.1	69.2	2.8	0.17	0.21
Pulsars Detected in Searches of the 350 MHz Data, Discovered Previously or Concurrently
J1600−3053	3.60	2.7	0.6 ± 0.1	0.7	52.3	2.5	14.3	0.20
J1805+0615	2.13	1.1	0.6 ± 0.1	7.3	64.9	3.8	0.34	0.02
J2043+1711	2.38	0.8	3.0 ± 0.1	1.2	20.7	1.5	1.4	0.17
J2256−1024	2.29	8.3	0.8 ± 0.1	3.4	13.8	1.3	0.21	0.03
Pulsars Discovered in the 820 MHz Survey Detected by Folding the 350 MHz Data
J0533+6759	4.39	2.1	1.0 ± 0.1	0.5*	57.4	2.4	Isolated	N/A
J1142+0119	5.07	0.4	0.6 ± 0.1	8.5	19.2	2.2	1.6	0.15
J1312+0051	4.23	0.4	1.7 ± 0.1	0.4	15.3	1.4	38.5	0.18
J1858−2216	2.38	2.0	0.8 ± 0.1	1.1*	26.6	0.9	46.1	0.21
Pulsars Discovered in the 820 MHz Survey Not Detected in the 350 MHz Data
J0843+67	2.84	⋯	⋯	⋯	20.7	1.5	7.4	0.27
J1921+0137	2.50	⋯	1.6 ± 0.2	4.9*	104.9	5.0	9.9	0.23
Fermi-LAT Pulsars Observed But Not Detected in the 350 MHz Data
J0106+4855	83.16	<0.15	1.9 ± 0.1	2.9*	70.9	3.0	⋯	⋯
J0622+3749	333.21	<0.15	2.0 ± 0.1	2.7*	⋯	⋯	⋯	⋯
J1627+3219	2.18	⋯		⋯	28.1	⋯	0.18	0.02
J2139+4716	282.85	<0.22	2.3 ± 0.2	0.3*	⋯	⋯	⋯	⋯

Notes.

^a γ-ray fluxes in the first group are from this work. Others are taken from the 3FGL catalog. ^b The $\dot{E}$ values have been corrected for the Shklovskii (Shklovskii 1970) effect (except for those marked with an "*") when proper motion information is available, and all have been corrected for the perpendicular and parallel components of acceleration with respect to the Galactic plane (Nice & Taylor 1995). ^c The distance is estimated using the YMW16 electron distribution model (Yao et al. 2017). The distances reported in Hessels et al. (2011) were based on the NE2001 model, as the YMW16 model was not complete then. ^d The minimum companion masses (M_c ) were calculated assuming a pulsar mass of 1.4 M_⊙ and an orbital inclination of 90°. ^e Not associated with the LAT source.

Download table as:  ASCIITypeset image

All of the 10 discoveries are MSPs, with the longest spin period being 7.6 ms. Our survey also included three LAT sources which were later determined to be young, slow pulsars through blind searches of the LAT γ-ray data (Pletsch et al. 2012). Folding our data with the γ-ray period over a plausible range of DMs failed to detect any of these three. These results are summarized in Table 2, and the 350 MHz pulse profiles of the 10 discoveries are shown in Figure 1. The mean flux densities in the table have been calculated using the radiometer equation (Equation (1)) and the appropriate sky temperatures for each source (Haslam et al. 1982). No estimates for statistical uncertainty for variability due to scintillation are available.

Zoom In Zoom Out Reset image size

As described above, our pointing positions were based on an early version of 1FGL. Improved positions and γ-ray spectra were published in the second and third Fermi-LAT catalogs (2FGL and 3FGL; Nolan et al. 2012; Acero et al. 2015), and so we list the associated 3FGL source and the offset from our pointing position (in degrees) in Table 1. In most cases, our pointing positions were consistent with the 3FGL positions, and the entire 3FGL error ellipse was contained within the half-power beamwidth. However, five of our pointings have no 3FGL counterpart within the FWHM beam, presumably due to either improved background modeling or previously confused sources being resolved. Our large beam also allowed for detecting pulsars well outside the LAT error ellipses, as was the case for our chance discovery of PSR J1551−0658. For the other 18 MSPs, the pulsars are apparently associated with the 3FGL sources. Seven 3FGL sources that we searched have now been plausibly associated with AGNs; 3FGL, J0523.5−2529 (Strader et al. 2014) and J1119.9−2205 (Swihart et al. 2022a) are low-mass binary candidates; and 3FGL, J1544.6−1125 is associated with a transitional MSP candidate (Bogdanov et al. 2015). We list these new associations in Table 1.

In the long interval between the survey observations and the finalization of this manuscript, the 4FGL Fermi point-source catalogs (Abdollahi et al. 2020) have been released. Because the γ-ray sources reported in this work are relatively bright, the properties are already well determined in the 3FGL catalog, and we determined that revising all quantities to their 4FGL values is not warranted.

Following each discovery, we began regular timing observations with the GBT as soon as was feasible, which in a few cases was up to 6 months after the initial discovery observation. All GBT observations used the GUPPI spectrometer (DuPlain et al. 2008b), typically at a central frequency of 820 MHz due to receiver availability, and when possible at 350 MHz. There were additional timing observations done at central frequencies of 1500 and 1900 MHz for several of the pulsars. For PSRs J0340+3140 and J1551−0568, the GBT campaigns were supplemented by observations with the Nançay telescope in France at frequencies of 1376, 1408, and 1598 MHz. For PSR J1551−0568, the position outside the Fermi error box was first determined with an imaging/pulsed observation mde with the Giant Metrewave Radio Telescope at a frequency of 322 MHz (Roy & Bhattacharyya 2013). Two of our sources, PSRs J0023+0923 and J0340+4130, have been included as targets in the NANOGrav Pulsar Timing Array and so are being extensively timed (Alam et al. 2021a, 2021b). For more details of the exact observations obtained for each pulsar, see the accompanying files, which contain the pulse time-of-arrival (TOA) data.

We observed each pulsar at least once a month for a span of no less than 1 yr, with at least one densely sampled set of observations over a period of a few weeks to obtain phase-connected solutions. Individual integration times were adjusted according to the pulsar strength, and ranged from 5 minutes to 0.5 hr. In some instances, we performed gridding observations at 820 or 1420 MHz to obtain immediately improved positions that increased the observing efficiency and simplified the determination of timing solutions. In a few cases, X-ray observations from a Neil Gehrels Swift Observatory (Swift) campaign to survey unassociated Fermi-LAT sources systematically provided higher-precision candidate positions (Falcone et al. 2011). Until obtaining enough timing points to yield an adequate orbital solution for each pulsar, we recorded data in search mode, which were analyzed using PRESTO (Ransom et al. 2002). In some cases, once a reasonable initial solution was obtained, we recorded data in full-Stokes (polarimetry) fold mode, with 256 bins per pulse period. A 1 minute observation with a 25 Hz pulsed noise diode was performed before each observation for polarization calibration. These data were analyzed with Psrchive (Hotan et al. 2004). In either case we used 2048 frequency channels typically spanning 100 MHz of band at a central frequency of 350 MHz, 200 MHz of band at 820 MHz, or 800 MHz of band at higher frequencies (1500 and 1900 MHz).

To determine the TOAs, a high-S/N profile at each frequency was used as a reference template. The TOAs were then calculated by cross-correlating each pulse profile with the reference template profile at each frequency (Taylor 1992). A model ephemeris was then fitted to the topocentric TOAs. In the case of a binary, an initial binary solution ephemeris was obtained by fitting a sinusoid to the measured barycentric frequencies. Timing solutions accompanying this paper were obtained with Tempo2 (Hobbs et al. 2006) and generally using the JPL DE421 solar system model (Folkner et al., 2009). Exceptions are noted below.

Six of our discoveries and two of the redetections are in tight binaries with orbital periods P_B < 1 day. Three of the short-period binaries have very-low-mass companions (M_c ≪ 0.1 M_⊙), and so we identify them with the "black widow" class of MSPs. Two of these black widows, PSRs J1124−3653 and J1810+1744, also exhibit eclipses. Two of the other tight binaries, PSRs J2129−0429 and J2215+5135, have minimum companion masses in the M_c ∼ 0.2–0.4 M_⊙ range, exhibit extensive eclipses, and have very strong orbital period variation. We classify them as "redbacks" (Roberts 2011). The final short-period binary, PSR J1302−3258, occasionally seems to eclipse at low frequencies for ∼10% of the orbit but has a companion mass consistent with a white dwarf and shows no other characteristics of being a redback. Three of the other discoveries—PSRs J0102+4839, J0307+7443, and J1551−0658—are in long-period (P_B > 1 day) binaries around companions whose minimum masses are consistent with being white dwarfs. Only two of the detected MSPs are isolated; one of them, PSR J0533+6759, is a weak redetection of a pulsar first discovered at 820 MHz. The other, PSR J0340+4130, is an isolated MSP which we study in this paper.

Because the LAT data span over 10 yr, we found that the timing solutions established solely via these radio pulsar timing campaigns were sometimes inadequate to predict the pulse phase accurately over the full duration. This is particularly true for the compact binaries exhibiting nondeterministic orbital period variations. Consequently, we extended the timing solutions using the Fermi-LAT data. In brief, we used maximum likelihood methods similar to those in Ray et al. (2011) and Kerr et al. (2015). However, these approaches are based on first estimating TOAs from LAT data. Here, we directly optimize the unbinned likelihood using PINT ²³ (Luo et al. 2021) and estimate parameters and their uncertainties using standard methods. Parameters like spin frequency and its derivative and proper motions, whose effect on pulsar timing residuals is cumulative in time, are often best estimated using LAT data, while others benefit from the generally higher precision of the radio TOAs. For each pulsar, we determined the optimal set of parameters to fit using each data set, and we iteratively fit the timing solution to the two data sets independently, yielding a timing solution that adequately describes both data sets and provides optimal parameter estimates. For the two redback pulsars, PSRs J2129−0429 and J2215+5135, the orbital period variations were so substantial that we were unable to obtain a timing solution in this manner. Instead, we used only LAT data and followed the method of Clark et al. (2021) to estimate their orbital period evolution.

The timing solutions ²⁴ for the 10 newly discovered MSPs are presented in Tables 3 and 4. For some pulsars, an "EFAC" term is included to describe additional scatter in the radio data beyond that indicated simply by radiometer noise. Errors reported by Tempo2 and given in these tables are scaled such that the χ² of the residuals is unity. For all of the MSPs, in addition to the measured $\dot{P}$, we give an "intrinsic" $\dot{P}$ which has been corrected for acceleration due to proper motion when available (i.e., the Shklovskii 1970 effect), acceleration toward the Galactic plane, and differential acceleration parallel to the Galactic plane (Nice & Taylor 1995). We also calculate the transverse velocity, given the DM-derived distance. Because of uncertainties in the distance and the Galactic rotation model, there is additional uncertainty in these estimates, which is hard to quantify. We use our best estimate of the intrinsic $\dot{P}$ for all derived quantities ($\dot{E}$, B, and τ) in Tables 3 and 4. The flux densities given are mean values of detections (not including eclipse times) calculated using the modified radiometer equation at 350, 820, and 1500 MHz.

Table 3. Timing and γ-ray Spectral Parameters for the Discovered MSPs

Timing Parameters	PSR J0023+0923	PSR J0102+4839	PSR J0307+7443	PSR J0340+4130	PSR J1124−3653
Fermi-LAT 3 yr Source	3FGL J0023.4+0923	3FGL J0102.8+4840	3FGL J0308.0+7442	3FGL J0340.3+4130	3FGL J1123.9−3653
R.A. (J2000)	00h23m16.s88027(1)	01h02m50.s66874(6)	03h07m55.s886(1)	03h40m23.s288229(6)	11h24m01.s1160(2)
Decl. (J2000)	+09°23'238902(5)	+48°39'427635(5)	+74°43'13426(5)	+41°30'452892(1)	−36°53'19087(4)
Pulsar Period (ms)	3.05020310344768(3)	2.9641124215887(4)	3.1560886884993(7)	3.29933936625009(2)	2.4095726733227(3)
Period Derivative, $\dot{P}$ (10−20 s s−1)	1.14238(3)	1.1406(3)	1.7279(7)	0.70485(8)	0.60151(1)
Reference Epoch (MJD)	55521	55527	55780	56279	55128
DM (pc cm−3)	14.326553(4)	53.5041(6)	6.3430(8)	49.57583(3)	44.8560(2)
Proper Motion in R.A. ${\mu }_{\alpha }\cos \delta$ (mas yr−1)	−11.00(7)	−3.4(4)	−1(1)	−0.77(5)	3.1(4)
Proper Motion in Decl. μδ (mas yr−1)	−8.8(1)	−1.9(5)	3(1)	−3.1(1)	−4.1(6)
Orbital Period (days)	0.13879914308(1)	1.672149564(1)	37.10775488(9)	...	0.226987946(8)
Projected Semimajor Axis (lt-s)	0.03484228(5)	1.8558825(7)	16.108338(8)	...	0.079631(6)
Orbital Eccentricity	<3.4 × 10−5	<2.3 × 10−6	<3.3 × 10−3	...	0.0
Epoch of Ascending Node (MJD)	55186.1134208(1)	55514.5773299(1)	55599.363937(4)	...	55128.586598(3)
Span of Timing Data (MJD)	55130.08−57375.95	54682.78−59050.55	54682.77−59033.74	55186.08−57378.23	55128.59−56854.65
Timing-derived Parameters
Mass Function (10−3 M⊙)	2.35 × 10−3	2.45	3.25	...	0.0105
Minimum Companion Mass (M⊙)	≥0.01	≥0.18	≥0.20	...	≥0.028
Galactic Longitude (deg)	111.3	124.7	131.7	153.7	284.1
Galactic Latitude (deg)	−52.8	−14.2	14.2	−11.1	22.8
DM-derived Distance (YMW16, kpc)	1.2	2.3	0.4	1.6	1.0
Surface Mag. Field, B (108 G)	1.8	1.8	2.3	1.5	1.2
Characteristic Age (Gyr)	4.2	4.1	2.9	7.4	6.3
Spin-down Luminosity, $\dot{E}$ (1034 erg s−1)	1.58922(4)	1.7291(5)	2.1699(9)	0.77477(8)	1.6974(5)
Transverse Velocity (km s−1)	80(16)	43(12)	6(2)	25(5)	25(6)
Acceleration ⊥to Plane (10−14 s s−1)	−5.8(4)	−4.9(3)	2.1(3)	−4.0(2)	4.4(2)
Acceleration ∥to Plane (10−14 s s−1)	−0.7(1)	−0.73(2)	−0.071(8)	2.1(4)	−2.0(4)
Corrected $\dot{P}$ (10−20 s s−1)	0.91(4)	1.096(9)	1.718(2)	0.65(8)	0.55(8)
Corrected $\dot{E}$ (1034 erg s−1)	1.27(6)	1.66(1)	2.158(3)	0.72(1)	1.57(2)
Corrected B (108 G)	1.7	1.8	2.3	1.5	1.2
Observed Parameters
Rotation Measure (rad m−2)	...	−86.3(8)	13(3)	56(1)	...
Flux Density at 350 MHz (mJy)	5.4	2.8	1.8	6.5	2.2
Flux Density at 820 MHz (mJy)	2.4	0.8	0.19	1.6	...
Flux Density at 1400 MHz (mJy)	0.5	...	...	0.5	...
Flux Density at 2000 MHz (mJy)	0.5	...	...	...	...
Spectral-fit Parameters
K (10−12 cm−2 s−1 MeV−1)	3.0(5)	2.4(2)	3.6(6)	2.9(2)	1.8(2)
Photon Index Γ	1.3(2)	1.7(1)	1.0(1)	1.2(1)	1.4(1)
EC (GeV)	1.9(5)	5(1)	1.6(2)	3.6(4)	3.7(7)
F100 (10−8 cm−2 s−1)	0.9(2)	2.0(2)	1.5(1)	1.5(1)	1.2(1)
G100 (10−11 erg cm−2 s−1)	0.7(1)	1.6(1)	1.6(1)	2.1(1)	1.3(1)
TS	470	1349	1920	2448	1149
TScutoff	48	66.5	254	200	85
TSb free	0.9	0.0	0.2	0.8	0.5
Efficiency	0.10	0.62	0.013	0.90	0.10

Note. The numbers in parentheses represent 2σ uncertainties in the last digit as determined by tempo2. Minimum companion masses were calculated assuming a pulsar mass of 1.4 M_⊙. The γ-ray spectral parameters are from fits of a power law with a exponential cutoff shape, given in Equation (3) with b = 1. The first three parameters are as defined in Equation (3): F₁₀₀ and G₁₀₀ give the integrated photon and energy flux above 0.1 GeV, respectively, while the last two parameters are the γ-ray detection significance of the source and significance of an exponential cutoff (as compared to a simple power law). The γ-ray efficiency is estimated as $4\pi {G}_{100}{d}^{2}/\dot{E}$ using YMW16 (Yao et al. 2017) distances assuming spatially uniform emission. The uncertainty on DM produces a negligible <0.0P uncertainty in radio/γ-ray alignment at 820 MHz, with a maximum of <0.02P at 350 MHz for an uncertainty of 0.001 DM units.

Download table as:  ASCIITypeset image

Table 4. Timing and γ-ray Spectral Parameters for the Discovered MSPs

Timing Parameters	PSR J1302−3258	PSR J1551−0658	PSR J1810+1744	PSR J2129−0429	PSR J2215+5135
Fermi-LAT 3 yr Source	3FGL J1302.3−3259	3FGL J1549.7−0658	3FGL J1810.5+1743	3FGL J2129.6−0427	3FGL J2215.6+5134
R.A. (J2000)	13h02m25.s5262(2)	15h51m09.s5279(2)	18h10m37.s28478(5)	21h29m45.s05(8)	22h15m32.s6(2)
Decl. (J2000)	−32°58'36843(7)	−06°58'0783(1)	+17°44'37380(1)	−04°29'0681(8)	+51°35'363(3)
Pulsar Period (ms)	3.770853091420(1)	7.09375570886(1)	1.6627549909837(1)	7.613937470415(5)	2.6096197252726(2)
Period Derivative, $\dot{P}$ (10−20 s s−1)	0.656(1)	2.05(6)	0.4476(1)	32.775(5)	2.8204(3)
Reference Epoch (MJD)	55925	55997	55530	55750	56362
DM (pc cm−3)	26.179(1)	21.578(1)	39.65952(6)	16.8771(1)	69.1947(1)
Proper Motion in R.A. ${\mu }_{\alpha }\cos \delta$ (mas yr−1)	...	...	7.5(2)	12.3(1)	0.3(5)
Proper Motion in Decl. μδ (mas yr−1)	...	...	−3.6(4)	10.1(1)	1.8(6)
Orbital Period (days)	0.78444160(2)	5.20657070(3)	0.148170285(1)	0.63522773(2)	0.172501860(1)
Projected Semimajor Axis (lt-s)	0.92792(4)	4.331018(6)	0.095378(1)	1.85219(2)	0.468131(4)
Orbital Eccentricity	3.2 × 10−5	1.1 × 10−5	4.7 × 10−5	0.0	0.0
Epoch of Ascending Node (MJD)	55924.803600(5)	55999.906196(1)	55529.9597206(2)	55702.111161(7)	55702.111161(7)
Span of Timing Data (MJD)	54683.52−59050.18	55623.28−56371.38	54682.96−59050.59	54682.84−58900.40	54682.84−58900.40
Timing-derived Parameters
Mass Function (10−3 M⊙)	1.39	3.21	0.042	16.9	3.70
Minimum Companion Mass (M⊙)	≥0.15	≥0.20	≥0.04	≥0.37	≥0.21
Galactic Longitude (deg)	305.5	1.5	44.6	48.8	99.8
Galactic Latitude (deg)	29.9	34.8	16.9	−36.8	−4.1
DM-derived Distance (YMW16, kpc)	1.4	1.3	2.4	1.4	2.8
Surface Mag. Field, B (108 G)	1.6	3.9	0.8	15.9	2.7
Characteristic Age (Gyr)	9.0	5.5	5.8	0.36	1.4
Spin-down Luminosity, $\dot{E}$ (1034 erg s−1)	0.483 (1)	0.227 (7)	3.8443 (9)	2.9314 (5)	6.2654 (8)
Transverse Velocity (km s−1)	...	...	96(19)	106(21)	25(9)
Acceleration ⊥to Plane (10−14 s s−1)	5.2(3)	5.3(3)	5.2(3)	−5.5(3)	−3.3(3)
Acceleration ∥to Plane (10−14 s s−1)	−1.2(3)	2.5(6)	−1.5(7)	−0.5(1)	−4.9(6)
Corrected $\dot{P}$ (10−20 s s−1)	0.608(6)	1.92(4)	0.36(1)	32.2(1)	2.770(8)
Corrected $\dot{E}$ (1034 erg s−1)	0.448(5)	0.213(4)	3.1(1)	2.86(1)	6.15(1)
Corrected B (108 G)	1.5	3.7	0.8	15.8	2.7
Observed Parameters
Rotation Measure (rad m−2)	...	...	...	...	...
Flux Density at 350 MHz (mJy)	1.3	0.7	78	2.1	16
Flux Density at 820 MHz (mJy)	2.4	0.3	7	...	0.8
Flux Density at 1400 MHz (mJy)	...	...	1.3	...	0.14
Flux Density at 2000 MHz (mJy)	...	...	0.4	...	0.11
Spectral-fit Parameters
K (10−12 cm−2 s−1 MeV−1)	1.8(2)	...	1.1(6)	1.6(2)	1.7(2)
Photon Index Γ	1.2(2)	...	2.0(1)	1.7(1)	1.3(2)
EC (GeV)	3.1(6)	...	3.8(7)	4(1)	4.2(8)
F100 (10−8 cm−2 s−1)	0.9(1)	< 0.5	4.1(2)	1.0(1)	1.0(1)
G100 (10−11 erg cm−2 s−1)	1.2(1)	...	2.3(1)	0.8(1)	1.3(1)
TS	1132	4	1958	541	887
TScutoff	101	...	72	31	83
TSb free	0.2	...	0.8	0.1	6.2
Efficiency	0.64	...	0.49	0.064	0.19

Note. The numbers in parentheses represent 2σ uncertainties in the last digit as determined by tempo2. Minimum companion masses were calculated assuming a pulsar mass of 1.4 M_⊙. The γ-ray spectral parameters are from fits of a power law with an exponential cutoff shape, given in Equation (3) with b = 1. The first three parameters are as defined in Equation (3): F₁₀₀ and G₁₀₀ give the integrated photon and energy flux above 0.1 GeV, respectively, while the last two parameters are the γ-ray detection significance of the source and significance of an exponential cutoff (as compared to a simple power law). The γ-ray efficiency is estimated as $4\pi {G}_{100}{d}^{2}/\dot{E}$ using YMW16 distances assuming spatially uniform emission. The uncertainty on DM produces a negligible <0.01P uncertainty in radio/γ-ray alignment at 820 MHz, with a maximum of <0.02P at 350 MHz for an uncertainty of 0.001 DM units.

Download table as:  ASCIITypeset image

3.1. Radio Polarization Observations and Results

3.2. X-Ray Observations and Analysis

3.3.  γ-ray Observations and Analysis

To characterize the γ-ray emission and pulsations, we analyzed ∼8 yr (2008 August 4–2016 July 20) of P8R2 Fermi-LAT data. All LAT data analyses were carried out with v10r00p05 of the Fermi Science Tools (STs). ²⁵ We selected SOURCE class events, as defined under the P8R2_V6 instrument response functions (IRFs), which had energies between 100 MeV and 100 GeV, reconstructed directions within 15° of the radio position of each MSP, and zenith angles ≤ 90°. We used the ST gtmktime to exclude times when the LAT was not in nominal science operations mode, when the LAT data were flagged as bad, and during bright LAT-detected solar flares and γ-ray bursts.

We performed a binned maximum likelihood analysis on a 212 × 212 region for each MSP using the P8R2_SOURCE_V6 IRFs, combining events which convert in the front and back sections of the LAT. We calculated binned exposure maps with 30 logarithmic energy bins and spatial pixels that were 01 on a side. The source model for each region included all 3FGL point and extended sources within 25° of each MSP and the diffuse background models. The Galactic diffuse emission was modeled using the gll_iem_v06.fits template, and the isotropic diffuse spectrum, which includes the extragalactic diffuse and residual instrument backgrounds, was modeled using iso_P8R2_SOURCE_V6_v06.txt. The spectral parameters of point sources within 6° of each MSP and with average significance ≥ 10σ, in 4 yr, were left free as well as the normalizations of the Galactic and isotropic diffuse components and point sources within 8° that were flagged as a variable in 3FGL. All other spectral parameters were kept fixed to the 3FGL values.

Since we are using twice as much (8 yr versus 4 yr) and improved (P8 versus P7REP) data compared to the 3FGL catalog, it is possible that our models of the regions may be incomplete. After our initial fits, we looked at the spatial and spectral residuals to see if additional sources needed free parameters or if additional sources needed to be included. For the region around PSR J2215+5135, we did find a few sources outside the 6° radius that were over- or undersubtracted. We refit this region with those additional source normalizations free and found satisfactory fits. In the region around PSR J0340+4130, we saw evidence for a source not included in 3FGL. Using the ST gtfindsrc, we localized this source to (R.A., decl.) = (5988, 5096) with a 95% confidence level error radius of 003, positionally associating this source with the flat-spectrum radio quasar, 4C +50.11, which flared in γ-rays in 2014 January (The Fermi-LAT Collaboration 2014). After including it in the model, we found satisfactory fits.

All of the new MSPs, with the exception of PSR J1551−0658, were found as significant point sources in the 3FGL catalog, with positions consistent with the timing positions within the LAT error ellipses. We fixed them to their timing positions in the region models and modeled their spectra as a power law (Equation (2)) and an exponentially cutoff power law (Equation (3) with b ≡ 1):

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}=K{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-{\rm{\Gamma }}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{dN}}{{dE}}=K{\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-{\rm{\Gamma }}}\exp {\left(-\displaystyle \frac{E}{{E}_{{\rm{C}}}}\right)}^{b}.\end{eqnarray} \tag{ 3 }$$

For both spectral models K is a normalization factor, with units of cm⁻² s⁻¹ MeV⁻¹, E₀ is a scale energy, fixed to the 3FGL pivot energy value for each MSP, and Γ is the photon index. For the cutoff model, E_C is the cutoff energy and b is an exponential index controlling the sharpness of the cutoff. We define a test statistic (TS) as twice the difference in the log likelihood to assess the overall source significance and to test nested models. The power-law model is rejected (large TS_cutoff) in favor of the cutoff model with >5.5σ significance for all detected MSPs. We also performed an analysis with the b parameter free but found no significant improvement (small TS_{b free}, ≲2.5σ) of the fit for any MSP when doing so. The best-fit spectral parameters and TS values for all nine MSPs are given in Tables 3 and 4.

We used the best-fit models and gtsrcprob to calculate spectral weights (probability that each event came from the respective MSP based on the likelihood model) for each event. We then folded the events for each MSP using the ephemerides described in Section 3 using the fermi plugin to tempo2 (Ray et al. 2011) and calculated a weighted H-test (Kerr 2011) for each MSP. We formed γ-ray profiles by generating weighted histograms from events lying within 3° of the source and display them in Figures 2–4 and 7–12. Estimates of the background level (dashed line) and error flags for the intensity in a given phase bin are computed via a bootstrap process described in 2PC. In each panel, we show one or more radio pulse profiles to place the γ-ray emission in context.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This paper presents a 350 MHz GBT survey for pulsars, carried out in 2009, which targeted 49 unidentified Fermi-LAT sources. We detected 18 MSPs in our survey; 10 of these were discoveries unique to our survey. PSR J0340+4130 is an isolated MSP. PSRs J0023+0923, J1124−3653, and J1810+1744 are short-orbital-period binaries with very-low-mass companions (M_c ≪ 0.1 M_⊙), with the latter two showing eclipses at frequencies ≲ 2 GHz, and hence can be classified as black widow systems. PSRs J2129−0429 and J2215+5135 have heavier, nondegenerate (Bellm et al. 2013; Breton et al. 2013) companions of masses M_c > 0.2 M_⊙, are eclipsed for a large fraction of their orbits, and have orbitally modulated X-ray emission (Gentile et al. 2014; Roberts et al. 2015). Moreover, this X-ray modulation is centered on the pulsar inferior conjunction (Wadiasingh et al. 2017). These are hallmarks of "redback" systems (Roberts 2011).

PSR J1302−3258 also has a short binary period and shows evidence for brief (∼10% of the orbit) eclipses at low frequencies, but does not have excess X-ray emission or a bright optical companion. PSRs J0102+4839, J0307+7443, and J1551−0658 are longer-period (P_B > 1 day) binaries with minimum companion masses typical of helium white dwarf stars orbiting MSPs. This final discovery, PSR J1551−0658, seems to be a chance detection since its position is well outside of the Fermi-LAT error box, and there is no evidence of pulsed emission in the γ-ray data. We show its radio pulse profile and polarimetry in Figure 13.

4.1. Implications of the Number of Spiders

One major surprise of these Fermi searches has been the number of compact (P_b < 1 day) binary systems which show evidence of a pulsar wind interacting with the companion. Such compact, often eclipsing systems come in two varieties, now commonly referred to as "black widows" (see Swihart et al. 2022a, for a census) and "redbacks" (see Strader et al. 2019).

In black widow systems, the pulsar's intense, relativistic wind ablates the companion star, stripping it of mass, and potentially "devouring" it entirely over the course of billions of years. ²⁶ PSR B1957+20, with a degenerate ∼0.02 M_⊙ companion, was the first such system to be found and is considered the prototype of the class. At radio wavelengths, this pulsar is eclipsed for ∼10% of its orbit, likely due to screening by material blown off of the companion star (i.e., not occultation from the companion star itself). It also shows timing irregularities, especially around eclipse ingress/egress (Fruchter et al. 1988), commonly referred to as "eclipse delays," as well as an orbital period derivative due to tidal effects. Whether such systems can fully ablate their companion and turn into the observed isolated MSPs is a matter of ongoing debate. Plausibly, some systems can ablate their companion star entirely, while others may not have sufficiently strong winds to do this within a Hubble time.

Additionally, some eclipsing systems have companions which are significantly more massive (∼0.2–0.4 M_⊙) and probably nondegenerate. These systems were first identified in globular clusters (PSR J1740−5340 in NGC 6397; D'Amico et al. 2001) and have now been termed redbacks ²⁷ from black widows (Roberts 2011). Such systems possibly represent an earlier stage in the recycling process, as argued in the case of PSR J1023+0038, which has displayed the characteristics of both an accreting low-mass X-ray binary and a normal radio MSP in the last decade (Archibald et al. 2009; Stappers et al. 2014). It is also possible that such systems have different progenitors and a different evolution than black widows (e.g., Chen et al. 2013; Benvenuto et al. 2014; Smedley et al. 2015). Prior to Fermi, only four of the known MSPs in the Galactic field were in compact/eclipsing systems.

4.1.1. The New Black Widows

PSR J0023+0923 has a spin period of 3.05 ms and a DM of 14.33 pc cm⁻³, with an implied DM distance of 1.2 kpc using the YMW16 model of Yao et al. (2017). With S₃₅₀ ∼ 2 mJy, the source was easily detected in a search of the first 200 s of the observation. Its orbital period is 3.3 hr, with a minimum companion mass of 0.016 M_⊙. These orbital properties are very similar to those of the black widow PSR J2051−0827. However, it does not show any evidence of an eclipse at either 350 MHz or 2 GHz. The radio profile is quite narrow (W₅₀ ∼ 0.05 ms, see Figure 1), and despite its black-widow-like very-low-mass companion, extensive timing shows no evidence for orbital period variations.

Early optical observations by Breton et al. (2013) showed that the companion is strongly heated but that it appears to have a small radius that is significantly underfilling its Roche lobe. More extensive observations by Mata Sánchez et al. (2023) show the Roche lobe fill factor is ∼50%. Therefore the surface gas may be too strongly gravitationally bound for the pulsar wind to ablate the surface, explaining the lack of eclipses. This may also explain why it does not exhibit orbital variations. The optical emission suggests a moderate inclination ∼ 42°, and so even if the pulsar is as heavy as 2 M_⊙, the companion would still be a very light 0.025 M_⊙, and is therefore consistent with other black widow masses.

The narrow pulse profile and lack of orbital variations allow for high-precision timing and has led to this source being included in the NANOGrav and EPTA pulsar timing arrays (e.g., van Haasteren et al. 2011).

PSR J1124−3653 has a spin period of 2.41 ms and a DM of 44.86 pc cm⁻³, with an implied DM distance of 1.0 kpc. With S₃₅₀ ∼ 0.3 mJy, it is one of the weaker sources presented here and required a search of the full 32 minute discovery observation in order to be found. Its orbital period is 5.45 hr, with a minimum companion mass of 0.028 M_⊙. It is eclipsed around 40% of the time at 350 MHz, 20% of the time at 820 MHz, but apparently not eclipsed at 1400 MHz. However, the 1400 MHz TOAs seem to have excess delays around superior conjunction. As noted above, Gentile et al. (2014) found significant hard X-ray emission suggestive of an intrabinary shock. Optical observations of the system show that the companion is very strongly irradiated by high-energy emission from the pulsar and/or shock (Draghis et al. 2019).

PSR J1810+1744 has a very short ²⁸ spin period of 1.66 ms and a DM of 39.66 pc cm⁻³, which implies a DM distance of 2.4 kpc. Relative to other known MSPs, the source is quite strong at low frequencies, S₃₅₀ = 20 mJy, and it was easily detected in a search of the first 200 s of the observation. At even lower frequencies, it is one of the brightest pulsars in the sky (Kondratiev et al. 2016). Its orbital period is 3.6 hr and the minimum companion mass is 0.045 M_⊙. Like the original black widow PSR B1957+20, PSR J1810+1744 has a very fast spin rate, though its orbital period is significantly shorter. From 350 to 1500 MHz, the radio profile evolves to have a sharp preceding peak (Figure 1), which is reminiscent of the trailing notch seen in the pulse profile of PSR B1937+21 (McKee et al. 2019). The γ-ray and radio peaks are roughly aligned (Figure 10), with a radio precursor becoming clearly detectable at 1500 MHz.

PSR J1810+1744 is eclipsed for roughly 10% of its orbital period, which is a similar eclipse fraction to the original black widow PSR B1957+20 (Fruchter et al. 1988). The eclipse boundaries are very sharp; the observed pulse intensity drops from its full noneclipse value to undetectable within ∼2 min. Eclipse delays, i.e., pulsed emission which has been shifted to a later pulse longitude, are also observed at eclipse egress. Eclipse studies at lower frequencies with LOFAR and WEPSR show the eclipse edges to be highly variable, and that the mass-loss rate is high enough to evaporate the companion within roughly a Hubble time (Polzin et al. 2018). PSR J1810+1744 has one of the highest minimum companion masses (M_c > 0.045 M_⊙) of all known black widows. Extensive observations with Keck indicate a very high neutron star mass of M_ns ∼ 2.1 M_⊙ and a nearly Roche-lobe-filling companion with mass M_c ∼ 0.065 M_⊙ (Romani et al. 2021). This suggests PSR J1810+1744 may be an evolutionary bridge between redbacks and black widows.

4.1.2. The New Redbacks

The longer eclipse fractions (typically 50%) and higher companion masses (0.2–0.4 M_⊙) of these systems identify them as redbacks, a class of system that before the discovery of the X-ray binary/MSP "missing link" PSR J1023+0038 (Archibald et al. 2009) was only observed in globular clusters. The companions of these systems also show strong orbital modulation in the optical (Breton et al. 2013) and X-ray (Archibald et al. 2009; Gentile et al. 2014; Bogdanov et al. 2015). It is unclear whether the redbacks are an evolutionary precursor of black widows.

PSR J2129−0429 has a spin period of 7.61 ms and a DM of 16.88 pc cm⁻³, with an implied DM distance of 1.4 kpc. With S₃₅₀ ∼ 0.5 mJy, it also required the full 32 minute observation to be detected. Its high spin-down rate implies a larger dipole surface magnetic field B ∼ 1.6 × 10⁹ G and a lower characteristic age than other redbacks. PSR J2129−0429 also has a relatively long orbital period (P_B = 15.2 hr) and a large minimum companion mass (P_c ≳ 0.37 M_⊙). The high companion mass and relatively slow rotation period argue that this is a redback that is likely in an early stage of recycling.

The pulsar has large orbital period variations and is eclipsed ∼50% of the time at low frequencies. The double-peaked orbital X-ray light curve shows there is a strong shock which is apparently wrapped around the pulsar, which may indicate that the companion magnetic field is dynamically important (Al Noori et al. 2018). Optical observations show the companion is nearly Roche lobe filled, and also show long-term trends in the brightness and temperature of the companion, suggesting the Roche lobe filling factor may be significantly changing on monthly to yearly timescales (Bellm et al. 2016; Al Noori et al. 2018). The system has a high inclination (i ∼ 77°), which is shown through both optical modeling and eclipses of the pulsar seen in γ-ray and X-ray data (Clark et al. 2023b).

PSR J2215+5135 has a spin period of 2.61 ms and a DM of 69.19 pc cm⁻³, with an implied DM distance of 2.8 kpc. It has the highest spin-down luminosity of any of the sources presented here and is the closest source to the Galactic plane. Its orbital period is 4.1 hr, with a minimum companion mass of 0.21 M_⊙. The pulsar has highly frequency-dependent eclipses, with it being eclipsed for roughly half of the orbit at 350 MHz but only for roughly a quarter of the orbit at 2 GHz. LOFAR studies have shown the eclipse fraction to increase to roughly three-quarters at 54 MHz (Broderick et al. 2016). The orbital and eclipse properties are very similar to those of PSR J1023+0038, and consequently PSR J2215+5135 fits well within the redback category. The radio pulse-profile morphology shows strong evolution with frequency (Figure 1), which makes interpretation of the radio/γ-ray alignment challenging.

The XMM-Newton light curves (Figures 5 and 6) show the double-peaked orbital variations centered at pulsar inferior conjunction that seem to be a common feature of redbacks. The nonthermal emission also has the very hard photon index (Γ ∼ 1) typical of these sources.

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

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

Zoom In Zoom Out Reset image size

This source has garnered significant interest since optical studies have suggested a neutron star mass M_ns > 2.1 M_⊙ (Linares et al. 2018; Kandel & Romani 2020) orbiting the Roche-lobe-filling M_c ∼ 0.33 M_⊙ companion.

PSR J1302–3205, with an orbit of 18.8 hr, a minimum companion mass of 0.15 M_⊙, and evidence for a brief radio eclipse around superior conjunction, is potentially a member of the redback class of systems where the companion has not yet fully evolved and there is the potential for further accretion episodes. We obtained multiple images of PSR J1302−3205 with the Las Cumbres Observatory Global Telescope (LCOGT) network of 1 and 2 m telescopes using various filters. Using the XB-News pipeline analysis (Russell et al. 2019), there were no firm detections of the source. An observation with the Spectral Camera on the LCOGT 2 m telescope at Siding Springs taken near pulsar inferior conjunction, when heating should produce a maximum, puts a conservative upper limit magnitude of ${m}_{{i}^{{\prime} }}\gt 21.5$. An observation with the MuSCAT3 instrument on the LCOGT 2 m telescope at Haleakala taken near pulsar superior conjunction put upper limits on the night side of the companion of ${m}_{{i}^{{\prime} }}\gt 22.3$, ${m}_{{g}^{{\prime} }}\gt 22.3$, ${m}_{{r}^{{\prime} }}\gt 22.6$, and ${m}_{{z}_{s}}\gt 21.9$. With an extinction estimate of A_v ∼ 0.22 from Drimmel et al. (2003) for the nominal distance of 1 kpc, the companion must be no brighter than an M8 subdwarf. The minimum companion mass is consistent with that of a fully evolved helium white dwarf expected for the orbital period–companion mass relationship derived by several authors (e.g., Tauris & Savonije 1999). We further note that PSR J1816+4510 has redback-like eclipses, but optical spectroscopy indicates the companion is more consistent with a proto–white dwarf (Kaplan et al. 2013). Combined with modest, apparently thermal X-ray emission seen from the system, we find no convincing evidence that PSR J1302−3205 should be classified as a redback at this time.

4.2. Implications for γ-ray Pulsar Population Statistics

Our survey contains 22 MSPs (all known as of publication) which are apparently associated with Fermi-LAT sources. Six are black widows, with the four discoveries reported here joined by PSR J2256−1024, discovered in the GBT drift scan survey (Crowter et al. 2015), and PSR J1805+0615, discovered in a survey of Fermi-LAT sources with Arecibo (Cromartie et al. 2016). Four are redbacks and three are redback candidates (Strader et al. 2019). Four of these show evidence of intrabinary shocks in X-rays (Gentile et al. 2014; Roberts et al. 2015), and five of them have been shown to heat their companion significantly (Bellm et al. 2013; Breton et al. 2013). This might suggest that the remarkably high fraction of these systems in this γ-ray-selected sample could be due to excess γ-ray emission from the pulsar wind interacting with the companion. However, it is clear from the estimated background levels in the γ-ray pulse profiles of our sample that essentially all of the observed γ-ray flux can be attributed to pulsed emission, and the γ-ray efficiencies (see Tables 3 and 4) are compatible with those of other MSPs detected by Fermi LAT (see 2PC and 3PC).

We note that all four of our detected pulsars with $\dot{E}\gt 3\times {10}^{34}$ erg s⁻¹ are eclipsing sources, suggesting that such sources are more energetic than the general MSP population. In general, γ-ray luminosity is correlated with $\dot{{\rm{E}}}$, albeit mildly for MSPs in the 10³³–10³⁴ erg s⁻¹ range (see Figure 14 and 3PC). This means that high-$\dot{{\rm{E}}}$ MSPs will be detectable in γ-rays over a larger volume compared to low-$\dot{{\rm{E}}}$ pulsars. Conversely, radio luminosity is only weakly coupled with $\dot{{\rm{E}}}$. Moreover, the eclipsing systems are also harder to detect in radio pulsar surveys since their short orbital periods can require acceleration searches for even short observations, and they can be eclipsed for a significant fraction of the time. Given these considerations, we expect γ-ray searches to be more efficient at detecting "spider" systems than radio surveys.

Zoom In Zoom Out Reset image size

To examine the $\dot{E}$ dependence more concretely, we compared the measurements in 3PC, which provide the best available estimates of γ-ray luminosity and Shklovskii-corrected $\dot{E}$, for isolated, short-period, and long-period binaries, with the results shown in Figure 14. Visually, there seems to be a lack of low-$\dot{E}$ short-period binaries. However, a two-sample Kolmogorov–Smirnov test indicates only modest evidence for a difference between the binary and isolated populations (p-value of 5%) and no evidence for a difference between short- and long-period binaries (p-value of 15%).

Of the 19 sources that we observed that are listed in the 3FGL catalog and remain unidentified or unassociated with possible AGN counterparts, most remain viable pulsar candidates. If they are MSPs, their nondetection is likely due to large accelerations due to orbital motion, large eclipse fractions, or simply radio faintness, as most MSPs appear to have radio and γ-ray beams sweeping similar regions of the sky, with >75% of γ-ray MSPs being radio loud (3PC), although in the case of black widows and redbacks they may be eclipsed more than half of the time at low frequencies. Redbacks may also spend some fraction of their time in a low accretion state with no observable radio pulsations but can have observable unpulsed γ-ray emission in this state (Stappers et al. 2014; Johnson et al. 2015). We therefore estimate that about five of the unidentified sources are MSPs. Some of these may be detectable with reobservations or at higher frequencies (note that PSR J1921+0137, discovered in an 820 MHz survey, is undetectable in our data). In fact, three of our targets, 1FGL J0523.5−2529 (Strader et al. 2014), 1FGL J1119.9−2205 (Swihart et al. 2022b), and 1FGL J1544.5−1127 (Strader et al. 2019), have been suggested to be redbacks through the identification of a plausible optical counterpart. PSR J0955−3947, another redback candidate (Braglia et al. 2020), was recently confirmed as a short-orbit pulsar by TRAPUM ²⁹ (Clark et al. 2023a).

Three of our targets are now known to be nonrecycled, middle-aged pulsars discovered through blind γ-ray searches. These appear to be radio-quiet pulsars, and none are detectable in our radio data. Compared to the 15 sources that have not yet been identified in our sample, these three sources are on average 2.6 times brighter in γ-rays. Assuming the γ-ray blind search detectability is largely flux limited and that there is a roughly spherical distribution of nearby, high Galactic latitude pulsars, this would suggest that roughly a dozen of these fainter sources could be middle-aged (τ∼10⁵ yr), radio-quiet pulsars. We note that nearby γ-ray pulsars are likely to have fairly low spin-down powers and hence are much more likely to be radio faint, compared to younger pulsars deep in the plane (3PC). We therefore find it plausible that essentially all of our remaining unassociated sources are either MSPs or nearby middle-aged pulsars.

Johnson et al. (2014) divide γ-ray pulsars into three classes based on their pulse profiles: γ-ray peaks trail radio (Class I), are aligned (Class II), and γ-ray peaks leading radio (Class III). (See also Espinoza et al. 2013 for a similar classification scheme.) Outer-gap and slot-gap (two-pole caustic) models best fit roughly equal numbers of Class I and III objects, while Class II objects are exclusively fit with pair-starved polar cap models. According to this scheme, PSR J0102+4839 is a Class III pulsar, with the γ-ray peak leading the radio peak. PSR J0307+7443 is a Class II pulsar, with the γ-ray and radio peaks aligned. PSR J0340+4130 is a Class III pulsar, with the primary γ-ray peak leading the primary radio peak, though it should be noted that the primary radio peak is aligned with the secondary γ-ray peak. Johnson et al. (2014) found that Class II MSPs had the shortest spin periods, with five of the six members of this class having spin periods less than 2 ms, while Class I and Class III MSPs had longer spin periods. PSR J0307+7443, with a spin period of 3.2 ms, does not fit that trend.

The close-orbit pulsars are more difficult to classify since nearly all show strong frequency dependence in their radio pulse profiles. At 350 MHz, PSR J1302−3258 is a Class I object, with the main radio peak leading the γ-ray peak. However, we note that an interpulse becomes apparent at 820 MHz, which trails the secondary γ-ray peak and becomes even more prominent at higher radio frequencies. The Fermi-LAT 2PC catalog used the brightest profile peaks (at 2 GHz) of PSRs J0023+0923, J1810+1744, and J2215+5135 when calculating phase lags. However, both PSRs J0023+0923 and J1810+1744 show a narrow "precursor" peak that is not visible at 350 MHz but whose peak brightness begins to approach the main peak's brightness by 2 GHz, and would presumably become the "main peak" at higher frequencies. At low frequencies, PSR J0023+0923 is Class I and PSR J1810+1744 is Class II, but if measured from the flatter-spectrum precursor peak, the classifications would be reversed. The pulse-profile evolution of PSR J2215+5135 is even more pronounced. At 350 MHz, it has two peaks which are fairly well aligned with the γ-ray peaks, implying this would be a Class II source, but at 2 GHz there is a single peak that is neither of the 350 MHz peaks and which trails the γ-ray peaks, and so would be Class III. Pulse profiles which strongly evolve with frequency seem to be a feature of many black widows and redbacks. We note that this has also been noted in the original black widow pulsar PSR B1957+20 (Johnson et al. 2014) as well as in the redbacks PSR J1048+2339 (Deneva et al. 2016) and PSR J1023+0038 (Archibald et al. 2009). For these systems, at least, classification schemes based on relative phases of γ-ray and radio peaks could break down. The pulse-profile evolution common to these also brings up the possibility that the companion is somehow influencing the radio emission of the pulsar if this frequency dependence is truly more common among these systems than in the general MSP population.

4.3. Implications for General Millisecond Pulsar Population Statistics

In the standard theory of MSP formation, a neutron star in a binary system is "recycled" by accreting matter from its companion star (Alpar et al. 1982; Stairs 2004), and consequently, most MSPs are expected to be in binary systems with a white dwarf, the natural end state of the evolution of the light companion following accretion and spinup. However, psrcat currently lists ³⁰ 52 MSPs in the Galactic field without binary companions, about 23% of the total. The formation of these isolated MSPs is still not well understood, though one possibility is that they ablated their companions with their strong relativistic particle winds (Ruderman et al. 1989). Another ∼32 MSPs ³¹ in the Galactic field do not have ordinary white dwarf companions but instead have either undermassive companions or nondegenerate companions. Understanding the various subpopulations of MSPs is critical for improving our view of the pulsar recycling process and whether it proceeds via only one channel or a variety of routes. The rapid rate of MSP discovery is can greatly help in this endeavor by rounding out the population with respect to previous observational biases.

Deep searches targeting γ-ray sources have a very different set of observational biases than shallow, wide-field radio surveys. Due to there being a much stronger correlation between spin-down luminosity and apparent γ-ray luminosity than there is with apparent radio luminosity, a γ-ray targeted survey will be more biased toward nearby, energetic MSPs and less biased toward high radio luminosities. Older low-frequency surveys had much less spectral resolution and were, therefore, more susceptible to DM smearing of the pulse profile, degrading their sensitivity to MSPs. Higher DM and scattering in the Galactic plane would also inhibit the discovery of MSPs. In the past, there were also strong biases against detecting pulsars with short orbital periods since acceleration searches only became common practice after 2000.

Only two out of 18 MSPs (11%) detected in this survey are isolated. 3PC reports 91 MSPs discovered by targeting unidentified Fermi sources, and of these, only 16 are isolated, ³² about 18%. Both are lower than the roughly 25% isolated fraction of the Galactic plane population found before 2009. Furthermore, eight of the pulsars detected in this survey have orbital periods P_b < 1 day. In the same 3PC sample, the fraction of short-period binaries is 54%. In stark contrast, only ∼14% of MSPs found before 2009 have such short binary periods. As a result of the Fermi searches, there are now more known MSPs with short orbits than isolated ones. If we compare the Galactic field population of P_b < 1 day MSPs to the ones in globular clusters, in both cases, the fraction is now ∼25%. This calls into question the importance of binary exchange in the formation of these systems in globular clusters.

Compared to binary MSPs, isolated MSPs could in principle be older or sustain a higher spin-down energy loss ($\dot{E}$) in order to disintegrate their companions fully. To investigate this, we have compared the $\dot{E}$ and characteristic ages of binary and isolated MSPs both in full and restricted to Galactic latitudes ∣b∣ > 5°. In order to obtain a larger sample, we use values of $\dot{P}$ which have not been corrected for the Shklovskii effect to estimate $\dot{E}$ and the characteristic age. (In addition to this caveat, it is important to note that these are derived quantities of P and $\dot{P}$ and are only rough proxies of the true ages and spin-down luminosities.) The distribution function for these two quantities is shown in Figure 15, the inspection of which indicates that the distributions for isolated and binary MSPs are broadly compatible. We quantitatively compared the distributions using a two-sample Kolmogorov–Smirnov test, finding modest evidence (p-value of 5%) that the $\dot{E}$ distributions differ. Further examination indicates that the ∣b∣ > 5 pulsars have only a 2% probability of being drawn from the same distribution, with isolated pulsars generally having lower values of $\dot{E}$, while the ∣b∣ < 5 populations are compatible. There is no substantial difference in any of the age distributions. Further work will be necessary to determine if the observed difference in spin-down luminosities is intrinsic or an observational selection effect against low-$\dot{E}$, high-latitude binaries.

Zoom In Zoom Out Reset image size

This work reports on a groundbreaking survey conducted more than 10 yr ago. Although we wish it had been published sooner, much of the delay was due to the success of the project. The MSP discoveries necessitated timing observations, and the discovery of such an unusually—at the time—large number of "spider" pulsars demanded multiwavelength follow up. It was always tempting to acquire more observational data and better understand this growing source class and its relation to γ-rays. The haul of new MSPs also created its own distractions, because it inspired numerous follow-up surveys, many cited above, of different parts of the sky, at different frequency ranges, and with different Fermi source lists. The success of this method continues through the present day, with powerful new instruments like FAST and MeerKAT and finding ever fainter radio pulsars in Fermi sources (Wang et al. 2021; Clark et al. 2023a).

Despite the long gap between observations and report, this work is the definitive record of the radio observations carried out, the follow-up timing observations, and the phase-aligned γ-ray and radio pulse profiles. Together with a synthesis of multiwavelength observations, particularly X-ray but also benefiting from optical analyses of redback candidates that lack pulsations, we can argue that γ-ray targeted searches correct a long-standing bias in radio surveys, which seem to have overproduced isolated MSPs relative to the general population. And it is only with the benefit of the substantial follow up that we can appreciate how efficient the original target selection was. It is thus fitting that this report finally appear in the literature, and we conclude: better late than never!

This work was partially supported by NASA Fermi Guest Investigator Program grants NNG10PB13P, NNX10AP81G, NNX12AE32G, NNX12AP21G, and NNX13AO96G. The Green Bank Telescope is operated by the National Radio Astronomy Observatory, a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The WSRT is operated by ASTRON/NWO. This work makes use of observations from the Las Cumbres Observatory Global Telescope network, with support from New York University Abu Dhabi. This paper includes observations made with the MuSCAT3 instrument, developed by Astrobiology Center and under financial supports by JSPS KAKENHI (JP18H05439) and JST PRESTO (JPMJPR1775), at Faulkes Telescope North on Maui, HI, operated by the Las Cumbres Observatory. We thank Gyula Jozsa for help scheduling the WSRT observations. J.W.T.H. acknowledges funding from an NWO Vidi fellowship and ERC Starting Grant "DRAGNET" (337062). M.S.E.R. acknowledges funding from the NYU Abu Dhabi faculty research fund. Work at NRL is supported by NASA and by ONR 6.1 basic research funding. E.C.F. is supported by NASA under award number 80GSFC21M0002. S.M.R. is a CIFAR Fellow. M.A.M. is supported by NSF award 2009425. M.E.D., E.C.F., M.K., M.A.M., S.M.R., and P.S.R. are members of the NANOGrav Physics Frontiers Center and are supported by NSF award 2020265.

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden.

Partially based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515.

We thank the NANOGrav collaboration for contributing data on PSR J0340+4130 that were used for this paper. We also thank Duncan Lorimer and Elizabeth Ferrara for maintaining the list of published and unpublished Galactic field MSPs.

Facilities: GBT - , Fermi - , GMRT - , WSRT - , NRT - , Swift - , CXO - , XMM - .