Discovery and Timing of Pulsars in the Globular Cluster M13 with FAST

The voltage data from the UWB were first channelized, down-sampled, summed in polarization, and written out in PSRFITS search-mode format⁹ (Hotan et al. 2004). During the conversion, the data were split into two subbands in order to decrease computing power and DM smearing within a channel. One frequency band spanned 250–400 MHz¹⁰ and the other spanned 400–800 MHz. The channel numbers are 21847 (with a channel bandwidth of 6.9 kHz) and 29128 (with a channel bandwidth of 13.7 kHz) for the two subbands, respectively, which were chosen by both minimizing the dispersion smearing within a channel for the cluster DM of 30 pc cm⁻³ and equalizing the DM smearing and sampling time. This number of channels resulted in sampling times in the two subbands of 145.64 μs and 72.82 μs, respectively. Figure 1 shows DM smearing within a channel as a function of frequency; the upper and lower dashed lines represent the sample times in the two subbands. We also estimated the effects of scattering using the empirical formula (t_scatt is in ms) (Bhat et al. 2004)

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{lg}({t}_{\mathrm{scatt}}) & = & -6.46+0.154\,\mathrm{lg}(\mathrm{DM})+1.07{\left(\mathrm{lg}(\mathrm{DM}\right)}^{2}\\ & & -3.86\,\mathrm{lg}\left(\displaystyle \frac{\nu }{\mathrm{GHz}}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

Even at 250 MHz, the scattering effect is estimated to be 27 μs for M13, which is much smaller than the sampling time. The data from the 19-beam receiver were recorded in the search-mode PSRFITS format. As described in Section 2, it has 4096 channels and a sampling time of 50 μs. The data cover a frequency range of 1000–1500 MHz.

All the data were processed using the PRESTO software package (Ransom et al. 2002). We used the rfifind package to generate a Radio Frequency Interference (RFI) mask in a similar way to that described in Hessels et al. (2007). The mask was used with prepsubband to generate 100 dedispersed time series at a central DM = 30 pc cm⁻³ and steps of 0.1 pc cm⁻³. The DM step size was chosen to have small DM smearing caused by an incorrect DM and sensible requirements on processing. DM smearing caused by an incorrect DM in channels is given by

$$\begin{eqnarray}&&{t}_{{\rm{\Delta }}\mathrm{DM}}=4.1\left[{\left(\displaystyle \frac{{\nu }_{\mathrm{low}}}{\mathrm{GHz}}\right)}^{-2}-{\left(\displaystyle \frac{{\nu }_{\mathrm{high}}}{\mathrm{GHz}}\right)}^{-2}\right]\left(\displaystyle \frac{{\rm{\Delta }}\mathrm{DM}}{\mathrm{pc}\,{\mathrm{cm}}^{-3}}\right)\mathrm{ms},\end{eqnarray} \tag{ 2 }$$

where ν_low and ν_high are the lowest and highest edges of the frequency bandwidth. ΔDM is the DM step size that we have chosen. In our cases, it introduces a maximum DM smearing within a channel of 0.4 μs, which is much smaller than the sampling time. The DM range was determined by considering the DM difference between the pulsar with the largest and smallest DM in GCs with more than one pulsar. The largest difference is 9.93 pc cm⁻³, thus we use a DM range of 10 pc cm⁻³ centered around the average DM of the known pulsars in M13.

A fast Fourier transform was applied to each of these dedispersed time series to obtain the power spectra and red noise was removed from each power spectra using the default settings for rednoise. After dedispersion, the effects from the rotation of the Earth and its motion around the Sun were removed.

Finally, an acceleration search was performed using accelsearch on the spectra with an acceleration range specified with the parameter zmax = 300. This parameter is the largest drift of the pulse frequency across frequency bins caused by acceleration, and an explanation can be found in Ransom et al. (2002). We also searched each data set after dividing them into 20 minute subintegrations in order to be sensitive to systems with minimum orbital periods of 3 hr. The initial candidates were sifted using the code sift, which is part of the pulsarTool¹¹ software package. Each candidate was searched for periods that have harmonic ratios from 1 to 16. The initial parameters resulting from the search for the candidates were used to refold the raw search-mode data and the results were visually inspected.

The new pulsar, with a spin period of approximately 3.00 ms, was first detected in the 400–800 MHz UWB data and a summary of the detection can be found in Figure 2. Here the signal was seen across the whole frequency band, and detected at a DM of 30.4 pc cm⁻³, which is close to the average DM value of the known pulsars in M13. The P–Pdot diagram in the bottom right panel shows that the signal has a significant period derivative, indicating the existence of acceleration to the pulsar most likely due to a binary companion. The 250–400 MHz UWB data were also searched following the same approach and also with respect to the detection parameters from the discovery in the upper band, but the signal-to-noise (S/N) ratio was too low to confirm the presence of the signal at this frequency. The discovery was later confirmed with follow-up observations using the 19-beam at the L band (see Figure 3). Of the 24 follow-up observations, 19 detected M13F over the duration of approximately one year. The nondetections are most likely attributable to interstellar scintillation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to derive the orbital parameters of the system, for each epoch the raw search-mode data were first folded into subintegrations of 1 minute each, using the DSPSR package (van Straten & Bailes 2011) and the ephemerides derived from the search results of the same data. From each sub-integration we generated its time-of-arrival (ToA) using the pat routine from the PSRCHIVE package. A plugin, stridefit2, from the TEMPO2 package (Hobbs et al. 2006) was used to generate a series of measurements of spin period, which was used in the fitorbit¹² software package to fit for the binary parameters. We used the measurements from fitorbit as an initial guess and carried out a timing analysis of the ToAs with tempo2 to obtain more accurate measurements of the timing parameters. These procedures were iterated a few times until the final coherent timing solution was converged to. DM was measured with tempo2 using TOAs from multifrequency subbands. Since low eccentricity was shown in the initial timing solution, we thus use binary model ELL1 (Lange et al. 2001) for M13F and the other binary pulsars in M13 in order to break degeneracy between orbital parameters. In Table 2 we present the position, rotational, and binary parameters of M13F. In Figures 4 and 5, the timing residuals are shown as functions of time and orbital phase.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2.  Timing Parameters for the Six Pulsars in M13, as Obtained from Fitting the Observed ToAs with TEMPO2

Pulsar	M13A	M13B	M13C
R.A., α (J2000)	16:41:40.87019(5)	16:41:40.39144(7)	16:41:41.00748(3)
Decl., δ (J2000)	+36:27:14.9788(4)	+36:25:58.4880(5)	+36:27:02.7438(2)
Spin frequency, f (Hz)	96.362234567(3)	283.44094318360(2)	268.66686619507(6)
First spin frequency derivative, $\dot{f}$ (10−15 Hz s−1)	0.675(3)	0.012(10)	−0.089(5)
Start of timing data (MJD)	58396.28	58396.28	58396.28
End of timing data (MJD)	58741.46	58741.46	58741.46
Dispersion measure, DM (pc cm−3)	30.4386(5)	29.4456(6)	30.1320(2)
Number of ToAs	760	750	758
Residuals rms (μs)	8.506	8.735	3.307
Binary Parameters
Binary model	⋯	ELL1	⋯
Projected semimajor axis, xp (lt-s)	⋯	1.388545(1)	⋯
Epoch of passage at ascending node, Tasc (MJD)	⋯	58561.7806094(1)	⋯
EPS1	⋯	2(1) × 10−6	⋯
EPS2	⋯	−2(1) × 10−6	⋯
Orbital period, Pb (days)	⋯	1.259112631(1)	⋯
Derived Parameters
Angular offset from center in α, θα (arcsec)	−4.4527	−10.237	−2.8051
Angular offset from center in δ, θδ (arcsec)	−20.5212	−97.012	−32.7562
Spin period, P (ms)	10.3775094516(3)	3.5280718048989(3)	3.7220816029988(8)
First spin period derivative, $\dot{P}$ (10−21 s s−1)	−7.3(3)	−0.1(1)	1.23(6)
Orbital eccentricity, e	⋯	2(1) × 10−6	⋯
Minimum companion mass, Mc,min (M⊙)	⋯	0.1605	⋯
Median companion mass, ${M}_{{\rm{c}},\mathrm{med}}$ (${M}_{\odot }$)	⋯	0.1876	⋯
Pulsar	M13D	M13E	M13F
R.A., α (J2000)	16:41:42.395232(7)	16:41:42.0221(1)	16:41:44.6058(3)
Decl., δ (J2000)	+36:27:28.2021(3)	+36:27:34.9676(9)	+36:28:16.0034
Spin frequency, f (Hz)	320.6886956910(2)	402.0937023406(4)	332.9448049492(9)
First spin frequency derivative, $\dot{f}$ (10−15 Hz s−1)	2.42(1)	−2.82(2)	−1.55(6)
Start of timing data (MJD)	58397.29	58396.29	58397.29
End of timing data (MJD)	58740.46	58741.45	58740.46
DM (pc cm−3)	30.451(3)	30.54(2)	30.366(4)
Number of ToAs	750	147	156
Residuals rms (μs)	3.345	7.46	15.07
Binary Parameters
Binary model	ELL1	ELL1	ELL1
Projected semimajor axis, xp (lt-s)	0.9243183(4)	0.035869(1)	1.251702(3)
Epoch of passage at ascending node, Tasc (MJD)	58397.38429971(5)	58397.27621667(1)	58398.0011780(7)
EPS1	0.0001578(9)	6(6) × 10−5	5(4) × 10−6
EPS2	0.0005508(5)	−0.00011(5)	−2(4) × 10−6
Orbital period, Pb (days)	0.5914408990(1)	0.11261745834(4)	1.378005120(6)
Derived Parameters
Angular offset from center in α, θα (arcsec)	+13.9368	+9.4353	+40.6053
Angular offset from center in δ, θδ (arcsec)	−7.2979	−0.532	40.5030
Spin period, P (ms)	3.118288899599(2)	2.486982497311(2)	3.003500835979(8)
First spin period derivative, $\dot{P}$ (10−20 s s−1)	−2.36(1)	1.75(2)	1.40(6)
Orbital eccentricity, e	0.0005730(6)	0.00012(5)	5(4) × 10−6
Minimum companion mass, ${M}_{{\rm{c}},\min }$ (M⊙)	0.1782	0.01942	0.1347
Median companion mass, ${M}_{{\rm{c}},\mathrm{med}}$ (M⊙)	0.2085	0.0225	0.1571

Note. M13A and M13C are isolated pulsars. M13B, M13D, M13E, and M13F are in binary systems. The companion masses are calculated assuming a pulsar mass of 1.4 M_⊙.

Download table as:  ASCIITypeset image

Five pulsars have been discovered in previous observations of M13. The timing solution of M13A was published in Kulkarni et al. (1991), but has not been updated. There are no published timing solutions for M13B to E. To detect those pulsars, we first attempted to fold our data using the ephemerides from the ATNF pulsar catalog¹³ (Manchester et al. 2005). However, we found those ephemeris are not accurate enough to keep phase coherence in our data, except for M13A. Thus, to detect M13B, M13C, M13D, and M13E for each epoch observation, we performed a search for their periodic signal and then folded the data using the ephemerides from the search. Once the detections were made, we then followed the procedures as described in the previous subsection to obtain an accurate timing solution for all five pulsars.

As shown in Table 1, M13A to C were detected in all 19-beam observations. M13D and E were detected in 92% and 83% of the 19-beam observations, respectively. Only M13A and M13B were detected in the two subbands (250–400 MHz and 400–800 MHz) of the UWB data. Timing solutions of these pulsars are summarized in Table 2. M13A and M13C are the only isolated pulsars known in this cluster, the periods of which are 10.38 ms and 3.72 ms, respectively. M13B is a 3.53 ms MSP in a binary system with an orbital period of 1.3 days. M13B has the largest position offset among the six pulsars in this cluster so far, which places it at the edge of the half-light radius of M13 (see Figure 6). M13B has the smallest DM (see Table 2) among the six pulsars in this cluster, which suggests that this pulsar is located on the nearside of the cluster, or the contribution to the gas internal to the cluster is lower here. M13D is also a binary MSP with spin and orbital periods of 3.12 ms and 0.6 days, respectively. M13E has the fastest rotational period in the cluster and is closest to the center of the cluster. M13E displays eclipses near the superior conjunction (at orbital phase 0.25). The eclipse event lasts for approximately 15 minutes, which is about 10% of the orbital period. Figure 7 shows the timing residuals and DM variation during the two eclipsing events of M13E. The DM variation is measured with 1 minute subintegrations using the frequency-resolved template-matching method developed in Liu et al. (2014). The peaks in the plot, i.e., the delay to the arrival time of the pulsar signal, are caused by increased DM, which is contributed by the material outflowing from the companion star. It can be seen that during the eclipses the DM increased by at least 0.05 pc cm⁻³.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The profiles of all the known pulsars at different frequency bands are illustrated in Figure 3. We only detected M13A, M13B, and M13F in the frequency band 400–800 MHz. We divided this band into four subbands, 400–500 MHz, 500–600 MHz, 600–700 MHz, and 700–800 MHz respectively. M13A and M13B were seen in the four subbands and showed profile variation across the frequency bands (see Section 5 below), while M13F was only seen in the subband 400–500 MHz.

M13F has the second shortest spin period and the longest orbital period among the six pulsars that have been discovered in M13 (see Table 2). The spin period derivative of M13F is 1.4 × 10⁻²⁰, which is a typical value of MSPs in GCs.¹⁴ The position offset of M13F (Figure 6) provides evidence that M13F is located at the edge of the cluster core. The median companion mass of M13F is 0.16 M_⊙, which indicates that the companion is likely a white dwarf. The eccentricity of the system, 5 × 10⁻⁶, is small compared to the "normal" MSP–WD systems in 47 Tuc (Freire et al. 2017). This might be due to the low stellar density in this cluster, which usually results in less encounter interactions. Encounter interactions were shown to be the main cause for eccentricities of low-mass binary MSPs in GCs, which are supposed to be born with $e\approx {10}^{-6}\mbox{--}{10}^{-3}$ (Rasio & Heggie 1995; Heggie & Rasio 1996).

The profile of M13A shows two components; the leading component decreases in amplitude as the observing frequency increases compared to the trailing component , indicating these two components have different spectral indices. Similar phenomena have been seen in other pulsars (e.g., Dai et al. 2015). M13B shows less apparent profile evolution in frequency, although the component at the leading edge fades out at frequencies above 500 MHz. M13F also has multiple components, as seen in Figure 3. It is hard to discern its profile evolution due to the low S/N detection in the frequency range 400–800 MHz.

Of all the observations we have made, the longest one is 228 minutes, which corresponds to a minimum flux density of 0.4 μJy to a candidate with S/N = 7. The sensitivity is calculated using the modified radiometer equation (Lorimer & Kramer 2004)

$$\begin{eqnarray}&&{S}_{\min }=\displaystyle \frac{\sigma \xi {T}_{\mathrm{sys}}}{G\sqrt{n{\rm{\Delta }}\nu {T}_{\mathrm{obs}}}}\sqrt{\displaystyle \frac{W}{P-W}},\end{eqnarray} \tag{ 3 }$$

where σ represents the minimum S/N of candidates. Imperfections arise due to digitization of the signal, and other effects are incorporated in ξ (Lorimer & Kramer 2004); we use ξ = 1.0, as our data are sampled in 8 bit. T_sys is the equivalent temperature of the observing system and the sky, which is 35 K for the FAST 19-beam receiver. G is the gain of the telescope and we use the value of 16 K Jy⁻¹ (Jiang et al. 2019). We use two polarizations, so n = 2. Δν is the bandwidth and is effectively 400 MHz in our L-band data. T_obs is the integration time, W and P are pulse width and period, respectively, and we take 10% for W/P. If we simply consider luminosity L = Sd², where S is flux density and d is the distance of the cluster, we obtain a limit on luminosity with FAST of 19.9 μJy kpc². The weakest pulsar we have detected in GC M13, M13F, has an L-band luminosity of L = 3.1 mJy kpc², and we did not find any pulsar with luminosity lower than this value. We also probed different orbital phases of binaries and thus covered epochs where accelerations may have been small. This may indicate that searches with higher sensitivity are needed if we are to find more pulsars in M13 and therefore require the SKA.

Our results show that the eccentricities of pulsars in M13 are very small compared to what is found in other GCs. As plotted in Figure 8, most GC binary pulsars have relatively large eccentricities, while the eccentricities of pulsars in M13 are close to those in the Galactic plane. Excluding those pulsars in M13 only a small fraction of pulsars in GCs have eccentricities below  10⁻³ . The majority of these fractions are from the non-core-collapse cluster 47 Tuc (Howell et al. 2000) and one is located in the cluster M62, which has an intermediate interaction rate per binary (Verbunt & Freire 2014). This is in agreement with the prediction that more low-eccentricity pulsars are found in clusters with very low interaction rates per binary (Rasio & Heggie 1995; Heggie & Rasio 1996). The eccentricities and positions in Table 2 and Figure 6 seem to confirm that eccentricity is related to the stellar density around the pulsar. The pulsars in M13 close to the core have relatively large eccentricities, like M13D and M13E. M13B has low eccentricity and is located farther from the core. The other pulsar with a large offset from the core is M13F, and it also has low eccentricity.

Zoom In Zoom Out Reset image size

In Verbunt & Freire (2014) the encounter rate for a single binary, γ, is used to characterize the difference between pulsar populations in GCs. In clusters with a low γ, a binary will have a relatively long life. In such low-density environments, binaries will evolve undisturbed, and result in an MSP population similar to that observed in the Galaxy, where binary MSPs with low-mass WD and low-eccentricity significantly outnumber isolated pulsars. M13 has the smallest γ among the 14 GCs listed in Verbunt & Freire (2014). The number of discovered binary pulsars is twice more than that of isolated pulsars in M13, thus the result is well explained by the above assumption. A low γ GC allows X-ray binaries to live long enough that the neutron star can be spun up to short periods (Sigurdsson & Phinney 1995). The fact that the periods of all six discovered pulsars in M13 are shorter than 0.015 s is consistent with this expectation.

The interstellar medium (ISM) along the line of sight affects the pulsar signal, which manifests as density variation (Lorimer & Kramer 2004). One of our initial thoughts is that if the pulsars are spatially close enough, the turbulence of the ISM along their signal propagation might be turbulent in a similar way, which might lead to correlations between flux density variations of these pulsars. M13 has a DM value of approximately 30 cm⁻³ pc and its pulsars have been seen to exhibit significant variations in flux density due to interstellar scintillation. In Figure 9, for each pulsar we plot the S/Ns per unit integration time from each observation, and calculate the correlation coefficient of each pair. No significant correlation is found between any pulsar pair. This indicates that the distance of each pulsar pair is not small enough that the pulsars can scintillate simultaneously.

Zoom In Zoom Out Reset image size

The difference between the largest and smallest DM values (δDM) in M13 pulsars is 1.1 cm⁻³ pc, which is the sixth smallest δDM among all the clusters that have been detected with more than one pulsar. The other five clusters with δDM smaller than M13 are 47 Tuc, M30, M3, M5, and NGC 6749, which are all clusters with relatively low stellar density. The ratio of δDM and average DM in M13 is 3.6%, which is in agreement with the value around 3% in 47 Tuc (Freire et al. 2001a). The contribution from the turbulence of Galactic electron column density to variations in DM is at the level of 0.05 cm⁻³ pc (Nordgren et al. 1992). The difference we observed is 20 times larger than this value. Such a DM difference could be caused by the gas distribution within the cluster. If a considerable population of pulsars can be discovered in M13, the gravitational field can be constrained by the DMs, positions, period derivatives, and proper motions of pulsars in three dimensions (Freire et al. 2001a, 2001b, 2017). This could possibly reveal the existence of an intermediate-mass black hole (e.g., Kızıltan et al. 2017; Prager et al. 2017; Perera et al. 2017; Abbate et al. 2019a, 2019b) in the center of M13.

Two of the four binaries in M13 have measured eccentricities, namely, they both have noncircular orbits. This might allow measurement of the rate of periastron advance $(\dot{\omega })$ and a method in Özel & Freire (2016) can be used to measure the mass of pulsars.

The companion star in a black widow system can be strongly heated on the side facing the pulsar, and this heating pattern can be observed in the optical light curve (e.g., Callanan et al. 1995; Romani et al. 2015). Considering the accurate position we derived for the black widow system M13E, we expected the system to follow this template. However, no evidence shows there is a variable star at the same position of M13E (Kopacki et al. 2003; Pietrukowicz & Kaluzny 2004).