The Parkes multibeam pulsar survey – VII. Timing of four millisecond pulsars and the underlying spin-period distribution of the Galactic millisecond pulsar population

Abstract

We present timing observations of 4-ms pulsars discovered in the Parkes 20-cm multibeam pulsar survey of the Galactic plane. PSRs J1552−4937 and J1843−1448 are isolated objects with spin periods of 6.28 and 5.47 ms, respectively. PSR J1727−2946 is in a 40-d binary orbit and has a spin period of 27 ms. The 4.43-ms pulsar J1813−2621 is in a circular 8.16-d binary orbit around a low-mass companion star with a minimum companion mass of 0.2 M_⊙. Combining these results with detections from five other Parkes multibeam surveys, gives a well-defined sample of 56 pulsars with spin periods below 20 ms. We develop a likelihood analysis to constrain the functional form which best describes the underlying distribution of spin periods for millisecond pulsars. The best results were obtained with a lognormal distribution. A gamma distribution is less favoured, but still compatible with the observations. Uniform, power-law and Gaussian distributions are found to be inconsistent with the data. Galactic millisecond pulsars being found by current surveys appear to be in agreement with a lognormal distribution which allows for the existence of pulsars with periods below 1.5 ms.

methods: statistical, stars: neutron, pulsars: general

1 INTRODUCTION

Millisecond radio pulsars (MSPs) are fascinating objects to study. Their phenomenal rotational stability allows them to be used for a wide variety of fundamental physics experiments including as a Galactic-scale observatory to search for low-frequency gravitational waves (see e.g. Hobbs et al. 2010). Ever since the discovery of the first MSP (Backer et al. 1982) it has been clear that the difficulties in detection imply that the Galactic population of MSPs is substantial. Early studies showed that the population of MSPs is comparable to that of normal pulsars (see e.g. Kulkarni & Narayan 1988; Johnston & Bailes 1991).

The continued improvement of data acquisition systems over the past 20 yr has led to a dramatic increase in survey sensitivity to MSPs. The number of known MSPs in the Galactic disc (i.e. those not associated with globular clusters) is now 228 (Manchester et al. 2005).¹ Because of the great success of blind surveys of the Galactic field (for a review, see Stovall, Lorimer & Lynch 2013), and targeted searches of Fermi sources (Ray et al. 2012), for the first time in a decade, Galactic MSPs outnumber their counterparts in globular clusters.

The Parkes multibeam pulsar survey (PMPS) of the Galactic plane is the most successful large-scale search for pulsars so far undertaken. Six previous papers in this series have presented timing parameters for 742 newly discovered pulsars and have discussed various aspects of the survey results (Manchester et al. 2001; Morris et al. 2002; Kramer et al. 2003; Faulkner et al. 2004; Hobbs et al. 2004; Lorimer et al. 2006; Crawford et al. 2013). Over the past 5 yr, several re-analyses of the survey data have been carried out. Keith et al. (2009) discovered a further 28 pulsars by applying new candidate sorting algorithms to the data processed earlier by Faulkner et al. (2004). Eatough, Keane & Lyne (2009) applied a new interference removal technique to a small portion of the data and discovered a further four pulsars. In a further reanalysis, Keane et al. found one fast radio burst (Keane et al. 2012) and 10 rotating radio transients (Keane et al. 2010) in addition to the 11 found originally by McLaughlin et al. (2006). Mickaliger et al. (2012) reported the discovery of the 34.5-ms binary pulsar J1725−3853 as well as four other millisecond pulsars. One other binary MSP, J1753−2814, has also been discovered as a result of this processing effort (Mickaliger et al., in preparation). Following earlier discoveries using ‘stack-slide’ acceleration searches by Faulkner et al. (2004), Eatough et al. (2013) report the discovery of 16 pulsars in a coherent acceleration search of the data. Ongoing processing by Einstein@Home volunteers (Knispel et al. 2013) has resulted in the discovery of a further 23 pulsars.

In this paper, we present timing solutions for four MSPs discovered in the PMPS. Preliminary discovery and confirmation observations of these pulsars were previously published by Faulkner et al. (2004). The total number of pulsars found in the survey so far stands at 833. Since extensive population studies of the normal pulsar population as revealed by the PMPS have already been carried out (Faucher-Giguère & Kaspi 2006; Lorimer et al. 2006), in this paper we focus our discussion on the spin-period distribution of the MSP population. The plan for this paper is as follows. In Section 2, we present the basic timing parameters, pulse widths, mean profiles and flux densities for the four new MSPs. In Section 3, we compile a sample of MSPs and use it to carry out a likelihood analysis to constrain the underlying distribution of spin periods. The main conclusions are summarized in Section 4.

2 FOUR MILLISECOND PULSARS

The pulsars were discovered using the processing schemes described by Faulkner et al. (2004). Following the confirmation and positional refinement procedures described by Morris et al. (2002), each pulsar was observed regularly at Parkes using initially the 512 × 0.5 MHz analogue filterbank system (Manchester et al. 2001) and subsequently the digital filterbank systems (Manchester et al. 2013). For each pulsar, pulse times of arrival were determined from the individual observations using standard pulsar timing techniques (see e.g. Lorimer & Kramer 2005) implemented in the psrchive software package (Hotan, van Straten & Manchester 2004).² A model containing the spin, astrometric and (if necessary) any binary parameters was fitted to the arrival times using the tempo2 timing package (Hobbs, Edwards & Manchester 2006). Arrival times were referred to TT(TAI) and the DE421 planetary ephemeris (Folkner, Williams & Boggs 2008) was used. Timing parameters are expressed in ‘TCB’ units native to tempo2 (see Hobbs et al. 2006, for the definition of TCB). Integrated pulse profiles are shown in Fig. 1.

Figure 1.

Integrated pulse profiles each showing 360° of rotational phase at 20 cm wavelength for the four MSPs described in this paper. Data were obtained with the Parkes digital filterbank systems.

Timing parameters from these analyses along with various derived quantities are presented in Tables 1 and 2. For PSR J1727−2947, time-of-arrival uncertainties were multiplied by a factor ranging between 0.85–1.5 for different backend systems to maintain a reduced χ² value close to unity. Also listed in Table 1 is the post-fit root-mean-square residual. The values obtained from our timing so far are relatively large (17–78 μs) and indicate that these pulsars are unlikely to be useful additions to MSP timing arrays. Although proper motions in right ascension have been measured for PSRs J1813−2621 and J1843−1448, we are unable to measure a significant proper motion in declination because of the low ecliptic latitude of these pulsars. Flux densities at 1400 MHz and pulse widths at 50 per cent of the peak level based on the profiles shown in Fig. 1 are listed in Table 1.

Table 1.

Spin, astrometric and derived parameters from the timing analysis of four MSPs. All astrometric parameters are given in the J2000 coordinate system. The reduced χ² values from each fit is listed as |$\chi ^2_r$|⁠. Figures in parentheses represent the uncertainties in the least significant digit and are the nominal 1σ tempo2 uncertainties. Distance estimates are based on the pulsar DM using the Cordes & Lazio (2002) NE2001 electron-density model. Pseudo-luminosities are computed by multiplying flux density by distance squared. Characteristic ages, magnetic fields and spin-down luminosities are based on the spin period and period derivative (see e.g. Lorimer & Kramer 2005) and account for the contributions due to the Shklovskii effect and Galactic acceleration.

Parameter	PSR J1552−4937	PSR J1727−2946	PSR J1813−2621	PSR J1843−1448
RA (hh:mm:ss.s)	15:52:13.2709(4)	17:27:15.09493(17)	18:13:40.59165(10)	18:43:01.3750(3)
Dec. (dd:mm:ss.s)	−49:37:49.744(11)	−29:46:36.797(17)	−26:21:57.055(18)	−14:48:12.61(3)
Proper motion in RA (mas y−1)	−3(3)	0.6(9)	−7.3(9)	10.5(19)
Proper motion in Dec. (mas y−1)	−13(8)	0(8)	−22(16)	12(15)
Epoch of position	J2000	J2000	J2000	J2000
Spin period (ms)	6.2843113814174(12)	27.0831832440066(12)	4.4300116286341(18)	5.4713308755095(6)
First derivative of spin period	1.900(4) × 10−20	2.4632(3) × 10−19	1.2466(6) × 10−20	6.209(18) × 10−21
Dispersion measure (cm−3 pc)	114.19(8)	60.74(3)	112.524(9)	114.51(7)
Rotation measure (rad m−2)	–	–61(32)	136(8)	–
Epoch of period (MJD)	54 033	54 723	54 058	53 934
Data span (MJD)	52 860–55 206	52 666–56 781	52 696–55 419	52 696–55 152
|$\chi ^2_r$| / degrees of freedom	1.00/132	1.45/208	0.90/134	1.16/115
Post-fit rms residual (μs)	78	43	17.5	49
Flux density at 1.4 GHz (mJy)	0.14	0.25	0.65	0.57
Pulse width at 50 per cent of peak (ms)	0.9	1.8	0.66	1.0
Distance (kpc)	4.8	1.4	2.9	2.9
Pseudo-luminosity (mJy kpc2)	3.2	0.49	2.1	4.8
Intrinsic period derivative	0.7(1.1) × 10−21	2.426(7) × 10−19	−0.6(1.7) × 10−21	−0.5(1.1) × 10−21
Characteristic age (Gy)	>15	1.77	>5.6	>14
Surface magnetic field (108 G)	<2.1	25.9	<2.4	<1.9
Spin-down luminosity (1033 ergs s−1)	<1.1	0.48	<5.7	<1.5

Table 2.

Measured and derived orbital parameters for PSRs J1727−2946 and J1813−2621 which use the ‘BT’ binary model (Blandford & Teukolsky 1976) and the ‘ELL1’ binary model (Lange et al. 2001), respectively. Parameters listed are the binary period (P_b), projected semimajor axis (a sin i), orbital eccentricity (e), first and second Laplace–Lagrange parameters (ϵ₁ and ϵ₂), longitude and epoch of periastron (ω and T₀) and epoch of ascending node (T_asc). Figures in parentheses represent 1σ tempo2 uncertainty in the least significant digits. For PSR J1813−2621, we also list the corresponding values of e, ω and T₀ computed from the Laplace–Lagrange parameters. The Keplerian mass function (⁠|$4\pi ^2\,G (a \sin i)^3/P_{\rm b}^2$|⁠, where G is Newton's constant), and the minimum companion mass (calculated assuming a 1.4 M_⊙ pulsar and setting, i = 90°) are listed.

PSR	J1727−2946	J1813−2621
Binary model	BT	ELL1
Pb (d)	40.307 710 94(3)	8.159 760 702(10)
a1sin i (light second)	56.532 497(5)	5.592 583(3)
e	0.045 629 43(16)	|$2.7^{+1.2}_{-1.1}\times 10^{-6}$|
ω (°)	320.396 25(20)	|$289_{-18}^{+28}$|
T0 (MJD)	54 711.471 69(3)	|$54\,061.5_{-4.0}^{+6.0}$|
ϵ1	–	−2.5(10) × 10−6
ϵ2	–	9(8) × 10−7
Tasc (MJD)	–	54 054.932 8319(6)
Mass function (M⊙)	0.1194	0.002 82
Min. comp. mass (M⊙)	0.827	0.188

For two of the MSPs, J1727–2947 and J1813–2621, significant levels of polarized emission was measured and these are shown in the integrated pulsed profiles in Fig. 2. Rotation measures for both these pulsars were determined using the rmfit tool within psrchive with conservative estimates of the uncertainties.

Figure 2.

Integrated polarization profiles for the two MSPs for which significant RMs were detected. The top panels of each plot show the polarization position angle as a function of pulse phase. The bottom panels show total intensity (bold line), linear polarization (solid line) and circular polarization (dotted line).

PSRs J1552−4937 and J1843−1448 bring the total number of isolated MSPs known in the Galactic disc to 37. When compared to the sample of 172 MSPs for which an orbiting companion has been confirmed, the fraction of observed isolated MSPs currently stands at 18 per cent. An outstanding issue in our understanding of MSP population is to explain this population in a self-consistent fashion. In particular, an open question is whether isolated MSPs formed in a different way from binary MSPs. We discuss this issue further in Section 3.6.

PSR J1727−2947 is a relatively long-period MSP (P ∼ 27 ms) in a mildly eccentric (e ∼ 0.04) 40-d binary system. With a minimum companion mass of ∼0.8 M_⊙, the system is most likely a member of the so-called ‘intermediate-mass binary pulsar’ (IMBP) class (Camilo et al. 2001) with a relatively massive CO white dwarf companion. The parameters for PSR J1813−2621 imply that it is very representative of the low-mass binary MSP population. Interpreting the orbital parameters in the standard way (see e.g. Lorimer & Kramer 2005), we infer a companion mass of at least 0.2 M_⊙, typical of a low-mass white dwarf. Optical studies of these companions may provide further insights into the nature of these two binary systems.

For nearby MSPs, it is well known (see e.g. Damour & Taylor 1991) that two significant contributions to the observed period derivative are the effects of secular acceleration (sometimes referred to as the ‘Shklovskii effect’, Shklovskii 1970) and Galactic acceleration. For a pulsar of period P, transverse speed V, acceleration |${\vec{a}}$|⁠, distance D with an intrinsic period derivative |$\dot{P}_{\rm int}$|⁠, the observed period derivative

\begin{equation} \dot{P}_{\rm obs} = \dot{P}_{\rm int} + P \left( \frac{{\boldsymbol a}\cdot {\boldsymbol {\hat{n}}}}{c} + \frac{V^2}{cD} \right), \end{equation}

(1)

where |${\boldsymbol {\hat{n}}}$| is a unit vector along the line of sight to the pulsar and c is the speed of light. Following the discussion in section 3.1 of Nice & Taylor (1995) to compute these effects, and assuming a 25 per cent uncertainty on the distances, computed using the NE2001 electron density model (Cordes & Lazio 2002), we calculated or placed limits on |$\dot{P}_{\rm int}$| for each pulsar and list our results in Table 1. The resulting characteristic age and magnetic field strength estimates for these pulsars are also indicated in Table 1. As can be seen, for all but PSR J1727−2946, these effects account for most of the observed period derivative.

3 THE SPIN-PERIOD DISTRIBUTION OF GALACTIC MSPs

The large sample of over 1000 normal pulsars detected in the various Parkes multibeam surveys has provided significant advances in our understanding of the normal pulsar population (see e.g. Faucher-Giguère & Kaspi 2006; Lorimer et al. 2006). The observed sample of MSPs is substantially less numerous because of their generally lower luminosity and observational selection effects. Nevertheless, the large sky coverage and uniformity of the observing systems of the multibeam surveys provides an excellent sample to begin characterizing their population. Recent work by Lorimer (2013) models this population via Monte Carlo realizations of synthetic pulsars drawn from distribution functions. These synthetic populations are subsequently ‘observed’ with realistic models of the surveys to produce samples that can be compared with the observed data. In this paper, we focus on constraining the spin-period distribution of MSPs.

For this study, we define an MSP as a pulsar with P < 20 ms. Our final sample of 56 MSPs is drawn from detections by this survey (PMPS), the Swinburne Intermediate Latitude Survey (SWIL; Edwards et al. 2001), the Swinburne High Latitude Survey (SWHL; Jacoby et al. 2007), the Parkes High Latitude Survey (PH; Burgay et al. 2006), the Perseus Arm Survey (PA; Burgay et al. 2013), and the Deep Multibeam Survey (DMB; Lorimer, Camilo & McLaughlin 2013). The basic parameters of this sample of pulsars are summarized in Table 3. For the purposes of population analyses, this sample of pulsars is a very natural one to analyse, since the pulsars were found using the same telescope, receiver and data acquisition system. Since the sensitivity of this system is well understood (Manchester et al. 2001; Lorimer et al. 2006), our survey models are reliable. In addition, since the surveys were all carried out at 20 cm wavelength, we need not make assumptions about MSP flux density spectra in order to extrapolate results from surveys carried out at other frequencies.

Table 3.

The 56 MSPs used in the population study. For each pulsar, we list the spin period (P), dispersion measure (DM), Galactic longitude (l), Galactic latitude (b), whether this is a binary pulsar as well as the survey which detected the pulsar. The surveys considered were the PMPS, PH, PA, DMB, SWIL and SWHL.

PSRJ	P	DM	l	b	Bin	Survey
	(ms)	(pc cc−1)	(°)	(°)
0437−4715	5.76	2.6	253.4	−42.0	Y	PH
0610−2100	3.86	60.7	227.7	−18.2	Y	PH
0711−6830	5.49	18.4	279.5	−23.3	N	SWHL
0721−2038	15.54	76.1	234.7	−2.9	Y	PA
0900−3144	11.11	75.7	256.2	9.5	Y	PH
0922−52	9.68	122.4	273.8	−1.4	Y	PMPS
1022+1001	16.45	10.2	231.8	51.1	Y	PH
1024−0719	5.16	6.5	251.7	40.5	N	PH
1045−4509	7.47	58.2	280.9	12.3	Y	SWIL
1125−6014	2.63	53.0	292.5	0.9	Y	PMPS
1147−66	3.72	133.5	296.5	−4.0	Y	PMPS
1216−6410	3.54	47.4	299.1	−1.6	Y	PMPS
1435−6100	9.35	113.7	315.2	−0.6	Y	PMPS
1546−59	7.79	168.2	323.5	−3.8	Y	PMPS
1552−4937	6.28	114.6	330.0	3.5	N	PMPS
1600−3053	3.60	52.3	344.1	16.5	Y	SWIL
1603−7202	14.84	38.0	316.6	−14.5	Y	SWIL
1618−39	11.99	117.5	340.8	7.9	Y	SWIL
1629−6902	6.00	29.5	320.4	−13.9	N	SWIL
1643−1224	4.62	62.4	5.7	21.2	Y	SWHL
1652−48	3.78	187.8	337.9	−2.9	Y	PMPS
1708−3506	4.50	146.8	350.5	3.1	Y	PMPS
1713+0747	4.57	16.0	28.8	25.2	Y	SWHL
1721−2457	3.50	47.8	0.4	6.8	N	SWIL
1723−2837	1.86	19.9	357.3	4.2	Y	PMPS
1725−38	4.79	158.4	349.4	−1.8	Y	PMPS
1730−2304	8.12	9.6	3.1	6.0	N	SWIL
1732−5049	5.31	56.8	340.0	−9.5	Y	SWIL
1738+0333	5.85	33.8	27.7	17.7	Y	SWHL
1741+1351	3.75	24.0	37.9	21.6	Y	SWHL
1744−1134	4.07	3.1	14.8	9.2	N	SWIL
1745−0952	19.38	64.5	16.4	9.9	Y	SWIL
1748−30	9.68	420.2	359.2	−1.1	Y	PMPS
1751−2857	3.92	42.8	0.6	−1.1	Y	PMPS
1753−2814	18.62	298.4	1.4	−1.2	Y	PMPS
1757−5322	8.87	30.8	339.6	−14.0	Y	SWIL
1801−1417	3.62	57.2	14.5	4.2	N	PMPS
1801−3210	7.45	176.7	358.9	−4.6	Y	PMPS
1802−2124	12.65	149.6	8.4	0.6	Y	PMPS
1804−2717	9.34	24.7	3.5	−2.7	Y	PMPS
1813−2621	4.43	122.5	5.3	−3.9	Y	PMPS
1826−24	4.70	81.9	8.3	−5.7	Y	PMPS
1835−0115	5.12	98	29.9	3.0	Y	PMPS
1843−1113	1.85	60.0	22.1	−3.4	N	PMPS
1843−1448	5.47	114.6	18.9	−4.8	N	PMPS
1853+1303	4.09	30.6	44.9	5.4	Y	PMPS
1857+0943	5.36	13.3	42.3	3.1	Y	PMPS
1905+0400	3.78	25.7	38.1	−1.3	N	PMPS
1909−3744	2.95	10.4	359.7	−19.6	Y	SWHL
1910+1256	4.98	38.1	46.6	1.8	Y	PMPS
1911+1347	4.63	31.0	47.5	1.8	N	PMPS
1918−0642	7.65	26.6	30.0	−9.1	Y	SWIL
1933−6211	3.54	11.5	334.4	−28.6	Y	SWHL
1934+1726	4.20	62.0	53.2	−1.1	Y	DMB
1939+2134	1.56	71.0	57.5	−0.3	N	DMB
2010−1323	5.22	22.2	29.4	−23.5	Y	SWHL

3.1 Likelihood analysis description

Following earlier work (e.g. Cordes & Chernoff 1997), we adopt a likelihood analysis to constrain the period distribution. In our approach to this problem, the probability p_i of detecting pulsar i in the sample with period P_i, and dispersion measure DM_i can be written as follows:

\begin{equation} p_i = f(P|\,a,b) {\cal D}(P_i,{\rm DM}_i). \end{equation}

(2)

In this expression, as usual in statistical parlance, the ‘|’ symbol denotes a conditional probability. The quantity f(P| a, b) represents the probability density function (PDF) of the period which we seek to constrain. All the models considered in this work can be described by two parameters which we refer to here generally as a and b. Specific parameters will be defined below. The detectability function, |${\cal D}$|⁠, reflects the probability of detecting the pulsar in one of the six surveys mentioned above and therefore appropriately accounts for their non-uniform period sensitivity. Note that in this analysis, we assume that any loss of sensitivity due to binary motion is negligible given the significant number of acceleration searches that have been carried out on these data (Faulkner et al. 2004; Eatough et al. 2013; Knispel et al. 2013).

Given a model period PDF and a detectability model which we describe in detail below, the likelihood function |${\cal L}(a,b)$| for a given combination of a and b is simply the product of all the 56 individual p_i values. The optimal set of model parameters |$\hat{a}$| and |$\hat{b}$| are those which maximize |${\cal L}$| in a grid of parameter space over a and b. Once the maximum likelihood |${\cal L}_{\rm max}$| has been found, the marginalized PDF for a is then simply the distribution of |${\cal L}(a,\hat{b})$| (and vice versa for the PDF of b). This approach provides PDFs for a and b as well as a means for evaluating different period PDFs from the ratio of the maximum likelihoods (i.e. the Bayes factor, K). For example, given two model PDFs ‘x’ and ‘y’, |$K={\cal L}_{\rm max,x}/{\cal L}_{\rm max,y}>1$| if model x better describes the sample than model y. According to Jeffreys (1961), a Bayes factor of between 1 and 3 is deemed to be essentially indistinguishable from the best model while Bayes factors in the range 3–10 begin to favour x over y. Bayes factors higher than 100 decisively favour x.

3.2 Detectability model

The detectability of a given pulsar, |${\cal D}$|⁠, reflects how likely it is to be found in the sample of 56 MSPs we will ultimately be applying this analysis to. Calculating |${\cal D}$| therefore requires accurately accounting for the difficulties in detecting each pulsar. Two approaches that can be brought to bear on this problem are to (i) run large numbers of Monte Carlo simulations which model the detectability; (ii) develop a simple analytical model. After initial experimentation with the first approach, it became clear to us that the Monte Carlo simulations require a large number of assumptions and significant computational resources to carry out a sufficient number of realizations necessary to estimate |${\cal D}$|⁠. We therefore followed the second approach and calculate |${\cal D}$| for each MSP in our sample. The essence of our approach, described in detail below, is to find for any line of sight the flux density distribution of pulsars: p(S). The detectability of a pulsar along this line of sight is then the fraction of such pulsars visible by a survey. For the ith pulsar, we may therefore write

\begin{equation} {\cal D}_i = \frac{\int _{S_{{\rm min},i}}^{\infty } p(S) {\rm d}S}{\int _{S_0}^{\infty } p(S) {\rm d}S}. \end{equation}

(3)

Here, S_{min, i} is the minimum detectable survey flux density of this pulsar which we compute from its pulse period, DM and pulse width (an intrinsic pulse duty cycle of 10 per cent was assumed for each pulsar). The term S₀ represents the lowest flux density detectable in the survey, i.e. the limit of S_{min, i} as P becomes large and DM tends to zero. In paper VI of this series (Lorimer et al. 2006), we gave expressions for computing S_{min, i} and refer the reader to this work for details. The advantage of this approach to the problem is that it does not depend on knowledge of individual distances to MSPs, which are uncertain. The analysis also does not depend on detected signal to noise or flux densities for any individual pulsars. Instead, it takes advantage of prior information about the pulsar population to calculate p(S) rigorously. As we show below, the final determination of the detectability function for this sample of pulsars can be given in terms of only two parameters which are robust to uncertainties inherent in the assumptions.

The calculation of p(S) is most readily achieved from an application of Bayes’ theorem which implies for a given line of sight the following relationship between PDFs in flux density S and distance D,

\begin{equation} p(D|S) \propto p(S|D) p(D). \end{equation}

(4)

Because the distribution we seek, p(S), is simply p(S|D) marginalized over distance, we can use equation (4) to show that

\begin{equation} p(S) = \int _0^{\infty } p(S|D) {\rm d}D \propto \int _0^\infty \frac{p(D|S)}{p(D)} {\rm d}D. \end{equation}

(5)

To get expressions for p(D|S) and p(D), we use the results described in section 3.3 of Verbiest et al. (2012). In this work, assuming a lognormal pulsar luminosity function (see e.g. Faucher-Giguère & Kaspi 2006) with mean μ and standard deviation σ, it is shown that

\begin{equation} p(D|S) \propto \frac{1}{D} \exp \left[ -\frac{1}{2} \left( \frac{\log S + 2 \log D - \mu }{\sigma } \right)^2 \right]. \end{equation}

(6)

Verbiest et al. (2012) also show that, along a line of sight defined by Galactic longitude l and latitude b, an axisymmetric distribution of pulsars with the radial density profile found in paper VI leads to the result

\begin{equation} p(D)\propto R^{1.9} D^2 \exp \left[-\frac{|z|}{h}\right] \exp \left[-\rho \frac{|R-R_0|}{R_0}\right]. \end{equation}

(7)

In this expression, z = D sin b is the vertical height off the Galactic plane, h is the scaleheight of pulsars, ρ is a dimensionless parameter used to scale the population over the Galactic disc, R₀ = 8.5 kpc is the Galactocentric radius of the Sun and the pulsar Galactocentric radius

\begin{equation} R = \sqrt{R_0^2 + (D \cos b)^2 - 2 R_0 D \cos b \cos l}. \end{equation}

(8)

With these analytical results, it is then straightforward using numerical integration of equation (5) to find the appropriate form of p(S) for each l and b. This PDF is then numerically integrated according to the limits in equation (3) to find |${\cal D}_i$|⁠.

Following the results of Bagchi, Lorimer & Chennamangalam (2011) and Lorimer (2013), we adopted nominal parameter values of h = 500 pc, ρ = 5, μ = −1.1 and σ = 0.9. Fig. 3 shows the results of this calculation on our sample of 56 pulsars as scatter plots of |${\cal D}$| as a function of P and DM. As expected, |${\cal D}$| is lower for shorter period and/or higher DM pulsars. To approximate this trend in the likelihood analysis, we set

\begin{equation} {\cal D} = [1 - \exp (-P/\alpha )] \, \exp (-{\rm DM}^2/2 \beta ^2), \end{equation}

(9)

where simple fits to the data show that α ≃ 10 ms and β ≃ 110 cm⁻³ pc provide a good description of these trends, as shown by the smooth surface in Fig. 3. The root-mean-square deviation of the data from this surface is 0.1. As described in Section 3.5, while the choice of duty cycle or population parameters assumed in the detectability analysis impacts the values of α and β somewhat, the particular values of α and β do not significantly affect our conclusions.

Figure 3.

Detectability as a function of period (in ms) and DM (in cm⁻³ pc). The data points show our estimates of detectability for each pulsar which we compute as described in the text. The grid shows our best approximation to this behaviour using the two-dimensional detectability function defined in equation (9) assuming the parameters α = 10 ms and β = 110 cm⁻³ pc. For each data point, a vertical line is drawn showing its distance from the best-fitting surface.

3.3 Period distributions investigated

We considered a variety of analytical functions to find the PDF which best describes the spin period of the MSP population. The simplest case of a uniform distribution is clearly not favoured by the data. In preliminary investigations, we found Bayes factors relative to other models of the order of 10⁻¹² and disregarded it from further analysis. Better approximations to the true period PDF can be found by considering functions with some well-defined peak. All four functional forms we investigate henceforth (i.e. Gaussian, gamma, lognormal and power-law distributions) require two parameters. In a similar way to paper VI, we refer to these parameters using capital letters. The Gaussian distribution has a mean A and standard deviation B,

\begin{equation} f(P)_{\rm gauss} \propto \exp \left[\frac{-(P-A)^2}{2B^2}\right]. \end{equation}

(10)

The gamma distribution is parametrized by C and D,

\begin{equation} f(P)_{\rm gamma} \propto \exp (-P/C) (P/C)^{D-1}. \end{equation}

(11)

The lognormal distribution is parametrized by E and F,

\begin{equation} f(P)_{\rm lnorm} \propto \frac{1}{P} \exp \left[\frac{-(\ln (P)-E)^2}{2F^2}\right]. \end{equation}

(12)

We also considered a power-law distribution parametrized by an exponent G and a minimum period H as follows:

\begin{eqnarray} f(P)_{\rm power} & = & 0 \quad {\rm for} \,P \le H, \end{eqnarray}

(13)

\begin{eqnarray} f(P)_{\rm power} & \propto & P^G \quad {\rm for} \,P>H \,\, {\rm and} \,P < 20\ {\rm ms}, \end{eqnarray}

(14)

\begin{eqnarray} f(P)_{\rm power} & = & 0 \quad {\rm for} \,P \ge 20\ {\rm ms}. \end{eqnarray}

(15)

Note that the last boundary condition simply reflects our definition of an MSP as a pulsar with P < 20 ms.

3.4 Application to the observed sample

Using the above method, we maximize |${\cal L}$| for each of the period distribution models. A program was written³ to implement the analysis and derive marginalized PDFs of the resulting model parameters. For each model, we normalized the detection probability and period distribution such that

\begin{equation} \int _0^{\infty } {\cal D}(P,{\rm DM}) \, f(P|\,a,b) \,\,{\rm d}P = 1. \end{equation}

(16)

This normalization ensures that the resulting likelihood values can be compared with one another to compute Bayes factors. In the results below, we give the Bayes factors for the best model relative to each model under consideration.

The results of our analysis when applied to the observed sample of 56 pulsars are summarized in Fig. 4 and Table 4. Fig. 4 shows the marginalized posterior PDFs for each of the model parameters. Table 4 lists the 95 per cent credible intervals for all the model parameters. The highest likelihood values were obtained for the lognormal model. The Bayes factors of the other models relative to this one are also given in Table 4. These results indicate that the lognormal and gamma distributions give by far the most plausible descriptions of the MSP spin-period distribution.

Figure 4.

Left and centre columns: marginalized posterior PDFs for each of the model parameters A–H (see equations 10–15 for definitions) obtained from the likelihood analysis described in the text. Right column: corresponding PDF for f(P) for each distribution when the nominal parameter values given in Table 4 are adopted.

Table 4.

Results of the likelihood analysis for the period distribution models considered. For each model, we list the median and 95 per cent confidence interval on the parameters defined in equations (10)–(15) along with the Bayes factor (K) computed by dividing the lognormal model likelihood by the likelihood of that model.

Model	First parameter	Second parameter	K
Gaussian	|$A=0.7_{-0.6}^{+1.7}$| ms	|$B=5.8_{-0.8}^{+1.0}$| ms	738
Gamma	C = 2.3 ± 0.4 ms	|$E=2.2_{-0.3}^{+0.4}$|	13
Lognormal	E = 1.5 ± 0.2	|$F=0.58_{-0.09}^{+0.12}$|	1
Power law	G = −1.7 ± 0.4	|$H=1.51_{-0.20}^{+0.05}$| ms	182

3.5 Testing the validity of the analysis

Before discussing the impact of our results, it is important to demonstrate the reliability of the parameter estimation approach and its sensitivity to assumptions. To do this, we generated fake samples of detectable pulsars with known period distributions and passed these as input to the likelihood analysis. We used the psrpop software package⁴ introduced in paper VI (see also Lorimer 2013) to generate synthetic populations of MSPs for this purpose. As a starting point, we distributed the model pulsars with model parameter values of h = 500 pc, ρ = 5, μ = −1.1 and σ = 0.9. For the period distribution, we then chose each of the four distributions in turn and set the parameters A–H to be the notional values given in Table 4 from our analysis of the real data. In each simulation, we generated enough synthetic pulsars such that a total of 56 of them were detectable by models of the PMPS, PH, PA, DMB, SWIL and SWHL surveys available in psrpop. With each synthetic sample, we first carried out a detectability analysis as described in Section 3.2 to determine values of the detectability-model parameters α and β and then applied these in our likelihood analysis. We found that the returned parameter values A–H from the likelihood analysis were entirely consistent with the input values of the period distribution of the parent population. In addition, we found that the method consistently favoured the correct form of the input distribution by assigning it the maximum likelihood. For example, when we generated synthetic populations assuming a lognormal distribution, we consistently found the Bayes factors for the lognormal likelihood model to be lower than the other distributions, as is seen for the actual sample of MSPs. Similar results were found when other underlying period distributions were assumed.

While the above results are very encouraging, they represent idealized conditions in which we input the actual values of h, ρ, μ and σ into the detectability analysis to determine α and β. In reality, of course, these numbers are not known and are only approximations to the true distribution of MSPs. To examine how robust the analysis is to changes in the assumed duty cycle, h, ρ, μ and σ, we repeated the above procedure over a range of values to determine α and β. The ranges we explored were 5–30 per cent duty cycles, 300 < h < 900 pc, 4 < ρ < 6, −2.5 < μ < −1.5 and 0.3 < σ < 1.5. Although these led to variations in the detectability parameters in the ranges 2 < α < 15 ms and 100 < β < 300 cm⁻³ pc, we still found that the input parameter distributions were recovered and that the correct distribution was favoured. An example of this is shown in Fig. 5 in which we see the inferred PDFs from an analysis of a fake population with a lognormal period distribution. These results give us confidence that our analysis on the observed sample of 56 MSPs is providing reliable insights into their underlying spin-period distribution, f(P).

Figure 5.

Marginalized PDFs for the lognormal parameters E and F deduced from a fake population in which the true values were E = 1.4 and F = 0.46 shown as dashed vertical lines. In this case, the assumed population parameters for the detectability analysis were intentionally biased to be h = 900 pc, ρ = 4, μ = −2.1 and σ = 0.5 and lead to a detectability model with α = 2.2 and β = 200. Even with such a bias, the PDFs successfully encompass the true population values and favour the lognormal model by an order of magnitude over the three other models.

3.6 Discussion

Based on the analysis presented in this paper, we have found evidence favouring the underlying spin-period distribution of Galactic MSPs to be lognormal in form. While a gamma distribution is compatible with the data, it is less favoured than the lognormal. Uniform, power-law and Gaussian distributions are decisively ruled out in our likelihood analysis as being good descriptions to f(P). We note that the strong preference for a lognormal model found here is in contrast to the power-law model proposed by Cordes & Chernoff (1997) based on a much smaller sample of MSPs. While the exponent of our power-law model tested here (–1.7) is consistent with theirs, the likelihood analysis strongly favours the lognormal model.

While our likelihood analysis weighs the different distributions we tested against each other, some measure of the absolute agreement between the lognormal model and the observed sample of 56 MSPs can be found by comparing the sample with the predicted observed period distribution for this model. Combining our detectability model and lognormal period distribution, the observed period distribution takes the form

\begin{equation} f_{\rm obs}(P) \propto \frac{1}{P} \exp \left[\frac{-(\ln (P)-E)^2}{2F^2}\right] \left(1 - \exp \left[ \frac{-P}{\alpha } \right] \right), \end{equation}

(17)

where the lognormal parameters E = 1.5 and F = 0.58 and the detectability parameter α = 10 ms. As can be seen from the comparison of this function with the binned data from the 56-MSP sample in Fig. 6, the agreement is excellent, with the reduced χ² value being 1.1.

Figure 6.

A comparison of the sample of 56 MSPs considered in this paper with our best-fitting period distribution from equation (17) (solid line).

Since the sample of MSPs used in this analysis is based on surveys carried out a decade ago, it is useful to confront the distribution we obtained with the present sample of objects. This is shown in cumulative form in Fig. 7 where it is seen that the 95 per cent credible region of lognormal functions we derive is broadly compatible with the present sample of 228 MSPs which have been detected in the Parkes High Time Resolution Universe Surveys (Keith et al. 2010; Barr et al. 2013), targeted searches of Fermi sources (Ray et al. 2012) and also in surveys at lower frequencies with Arecibo and Green Bank (Deneva et al. 2013; Stovall et al. 2014). We note that the observed sample lies to the upper end of the 95 per cent credible region shown in Fig. 7. Future studies of this newer larger sample of MSPs should, therefore, provide more stringent constraints on the period distribution.

Figure 7.

Cumulative distribution functions showing the observed sample of 56 MSPs (rightmost dashed curve), the 95 per cent credible region of our best-fitting lognormal model (shaded band) and the current sample of 206 MSPs with P < 20 ms (leftmost dashed curve). The deviation between the shaded band and the rightmost dashed curve highlights the observational selection effects at short periods in our sample of 56 MSPs.

The general agreement with our lognormal model and the present sample of MSPs suggests that the period-dependent selection effects on these ‘first generation’ Parkes multibeam surveys (i.e. PMPS, PM, PA, SWIL, SWHL and DMB) which we model in our detectability function are much less severe in the present generation of MSP surveys.

4 CONCLUSIONS

We have presented timing models for four MSPs found as part of the PMPS of the Galactic plane. From a likelihood analysis of the sample of 56 MSPs detected with this earlier generation of PMPS, we demonstrate that the underlying population of spin periods for MSPs is compatible with a lognormal distribution. When this distribution is confronted with more recent discoveries from other surveys, we see that it is broadly consistent with the new results. It is important to note that the distribution we have derived here applies to the present-day MSP population. Although to first order, because of the very low spin-down rates of MSPs, the birth spin-period distribution may not be significantly different (see e.g. Camilo, Thorsett & Kulkarni 1994), further investigations are necessary to confirm this conjecture.

Although the true period distribution for MSPs may not be as simple as our analysis might initially suggest, it is clear that the distributions considered here all allow for the existence of a small fraction of pulsars with P < 1.5 ms. Based on the smooth curves shown in Fig. 7, the fraction of such pulsars in the population is around 3 per cent. The true fraction could even be higher than this if we have overestimated the detectability of such rapidly spinning pulsars. Given our estimate of the analytic form of the observed period distribution given in equation (17), we find that the probability of not detecting a pulsar in our sample of 56 MSPs with period P < 1.5 ms in the current sample is 99.2 per cent. This is entirely consistent with the lack of such pulsars in the sample so far. Further discussion about the possibility of submillisecond pulsars can be found in Levin et al. (2013) and references therein.

An inspection of the current sample of MSPs shows no statistically significant difference between the spin-period distributions of isolated objects versus binary systems. A useful approach which is currently being pursued (Lazarus, private communication) is to artificially add short period signals to existing search data sets and directly test the effectiveness of pulsar search codes in recovering these signals. The modelling techniques presented here may be useful in further analyses of the MSP population which need to take account of the different selection biases and observing frequencies that have taken place since the completion of the initial Parkes multibeam surveys. The techniques may also be applied to other population parameters in which selection effects may be apparent, for example the |$P-\dot{P}$| distribution.

The Parkes radio telescope is part of the Australia Telescope which is funded by the Commonwealth of Australia for operation as a National Facility managed by CSIRO. DRL acknowledges support from the Royal Society as a University Research Fellow during the early phases of this project. DRL and MAM acknowledge support from Oxford Astrophysics while on sabbatical leave in 2013. PE acknowledges a Fulbright Research Scholar grant administered by the US–Italy Fulbright Commission and is grateful to the Harvard–Smithsonian Center for Astrophysics for hosting him during his Fulbright exchange. The Fulbright Scholar Program is sponsored by the US Department of State and administered by CIES, a division of IIE. Pulsar Research at UBC is supported by an NSERC Discovery Grant. Current support to DRL and MAM is provided by the National Science Foundation PIRE award 0968296. DRL thanks Joris Verbiest for useful discussions and to Simon Johnston and the CSIRO staff for their hospitality during the final stages of this work. We thank the referee for useful comments on the originally submitted version of this paper.

1

http://astro.phys.wvu.edu/GalacticMSPs

2

http://psrchive.sourceforge.net

3

http://psrpop.phys.wvu.edu/pdist

4

http://psrpop.sourceforge.net