On the Spin Period Distribution of Millisecond Pulsars

Spin period distribution provides important clues to understand the formation of millisecond pulsars (MSPs). To uncover the intrinsic period distribution, we analyze three samples of radio MSPs in the Galactic field and globular clusters. The selection bias due to pulse broadening has been corrected but turns out to be negligible. We find that all the samples can be well described by a Weibull distribution of spin frequencies. Considering MSPs in the Galactic field or globular clusters and in isolation or binary systems, we find no significant difference in the spin distribution among these subpopulations. Based on the current known population of MSPs, we find that submillisecond pulsars are unlikely to be discovered by the Square Kilometre Array, although up to ∼10 discoveries of pulsars that spin faster than the current record holder of P = 1.4 ms are expected.


INTRODUCTION
Pulsars are fast-spinning neutron stars whose electromagnetic radiation regularly sweeps over the Earth.Their strong gravitation, dense matter and diverse observational features at various wavelengths have made them favourable objects to study for a range of topics, including the search for nanohertz gravitational waves (Agazie et al. 2023;Antoniadis et al. 2023;Reardon et al. 2023;Xu et al. 2023), theories of gravity (Weisberg & Huang 2016;Kramer et al. 2021), neutron star equation of state (Haensel et al. 2007), ionised interstellar medium (Yao et al. 2017) and the evolution of binary systems (Tauris et al. 2012).
Currently, more than 3000 pulsars have been discovered, among which ∼ 560 are millisecond pulsars (MSPs) (see PSRCAT 1 , Manchester et al. 2005).Both numbers will increase significantly in the next few years when new telescopes like the Five-hundred-meter Aper-Corresponding author: Xing-Jiang Zhu zhuxj@bnu.edu.cn 1 https://www.atnf.csiro.au/research/pulsar/psrcat/ture Spherical radio Telescope (FAST) and MeerKAT complete their first rounds of pulsar surveys (Smits et al. 2009b;Li et al. 2018;Han et al. 2021;Padmanabh et al. 2023).More importantly, when the square kilometer array (SKA) is used for pulsar search, the numbers are expected to further grow several folds (Smits et al. 2009a), possibly finding a large fraction of all Milky-Way pulsars that are beaming toward the Earth (Keane et al. 2015).
With such a large number of pulsars, statistical uncertainties due to small sample size become less significant and population analysis has become an effective way to study various population characteristics, including luminosity, birth rate, spatial (Narayan 1987;Cordes & Chernoff 1997;Lorimer et al. 2006) and kinematic distribution (Hobbs et al. 2005;Verbunt et al. 2017;Igoshev 2020).It is also an important method to constrain braking index and study the decay of surface magnetic field, which plays a key role in pulsar evolution (Phinney & Blandford 1981;Johnston & Karastergiou 2017;Dirson et al. 2022).
The intrinsic period distribution is a long-standing problem in pulsar population studies.A well-established period distribution is an essential input in population simulations (Bates et al. 2014), which are used to pre-dict survey yields (e.g.Keith et al. 2010;Ng et al. 2015), to estimate the total number of pulsars in the Galaxy, and, by comparing the simulated sample with survey results, to study other population characteristics (Faucher-Giguère & Kaspi 2006).
For MSPs, an accurate period distribution is critical to infer the value of the minimum spin period, which may impose significant constraints on theories of dense matter, since the maximum spin frequency (i.e., the Keplerian frequency) of pulsars depends sensitively on the equation of state governing their internal composition (Haensel et al. 2009;Watts et al. 2015).It is also important in answering whether there are sub-millisecond pulsars (Possenti et al. 1998;Burderi et al. 2001;Levin et al. 2013;Lorimer et al. 2015), which are expected theoretically (e.g.Friedman et al. 1986;Du et al. 2009) but have not been discovered despite numerous search attempts (e.g., Edwards et al. 2001;Han et al. 2004;Davoust et al. 2011).In addition, joint analysis of period and period derivative distributions were used to study the spin-up history of MSPs (Pringle & Rees 1972;Ghosh & Lamb 1979;Tauris et al. 2012;Liu et al. 2022) and the accretion physics during their low-mass X-ray phases (Chakrabarty 2008;Papitto et al. 2014;Patruno et al. 2017;C ¸ıkıntoglu & Ekşi 2023).
The period distribution of MSPs was firstly studied by Cordes & Chernoff (1997).Using a sample of 22 MSPs, they studied a power-law model with a lower cut-off and found the minimum period > 0.65 ms with 99 percent confidence.However, the power-law model was challenged by Lorimer et al. (2015), which studied a sample of 56 MSPs discovered by various surveys that used the Parkes (Murriyang) radio telescope.Using maximum likelihood analysis, Lorimer et al. (2015) studied four two-parameter models and found that a log-normal distribution was favoured over the power-law model.Using a sample of 337 radio MSPs, Patruno et al. (2017) found that, in the spin-frequency space, a Weibull model performed better than the log-normal distribution.
The studies of pulsar population are plagued with selection effects (Lorimer 2011), and so is the study of period distribution.Several factors can introduce period dependency into the pulsar detection rate.Notable ones include pulse dispersion and scattering, which increase the effective pulse width and thus the flux threshold to make a detection (e.g.Dewey et al. 1985;Cordes & Chernoff 1997;Lorimer 2008).The selection effects thus should be taken into account in the studies of period distribution.
In this paper, we study the period distribution of MSPs using the largest available sample from known pulsar surveys, which is more than four times the sam-ple size of Lorimer et al. (2015).We have also developed a self-consistent Bayesian framework that accounts for selection effects in exploring the underlying period distribution.Compared with the grid-searching method used in Lorimer et al. (2015), the Bayesian framework is more efficient in sampling high-dimensional likelihoods, enabling us to study more complex models and reduce the possibility of model misspecification.More importantly, the Bayesian method provides a natural way to incorporate selection effects into the analysis (e.g.Hinton et al. 2017;Vitale et al. 2020) and allows robust model selection and parameter estimation.
The structure of the paper is as follows.We firstly describe the pulsar data and corresponding survey parameters in Section 2, then describe the Bayesian framework for parameter estimation and model comparison in Section 3.After that, Section 4 outlines the computation of the selection effects, while Section 5 gives the period distribution models in consideration and Section 6 presents our results.Finally, the results are discussed and concluded in Section 7.

THE DATA
Here, MSPs are defined as pulsars with a spin period, P , shorter than 20 ms.Our fiducial sample is Sample A, whose selection criteria are specified below.
We only use radio pulsars and excluded pulsars discovered at other wavelengths, as the selection effects may vary significantly.We also exclude pulsars in globular clusters, as their spatial distribution differs from that of field pulsars (see Appendix A.1), making them unfit for the treatment of selection effects considered here.Furthermore, we focus on pulsars discovered by surveys around 1400 MHz (L-band), avoiding the uncertainty introduced by extrapolating flux density measured at different frequencies.
Since the computation of selection effects depends sensitively on survey parameters, we could only use surveys with well-defined parameters, i.e. the parameters should not change too many times over the survey lifetime.Using the Galactic MSP catalogue 2 , which collects both MSP parameters and corresponding discovery surveys, 46 different surveys were found.To make an efficient analysis, we set a yield threshold of 10 pulsars and 13 surveys satisfies the requirement.Since most radio MSPs were discovered in blind surveys, we thus exclude targeting searches of Fermi sources to maintain sample homogeneity.We are left with five surveys: the Parkes Multibeam Pulsar Survey (PMPS), High Time Resolution Universe Southern survey (HTRUS), Pulsar Arecibo L-band Feed Array (PALFA) survey, Commensal Radio Astronomy FasT survey (CRAFTS) and Galactic Plane Pulsar Snapshot (GPPS) survey.The five surveys and corresponding parameters are listed in Table 1.We note that both CRAFTS 3 and GPPS 4 are on-going projects and have a number of MSPs that are not yet included in either the Galactic MSP catalogue or PSRCAT.We thus query the corresponding project webpages and add these new discoveries.
In total, there are 260 MSPs in Sample A, which is more than four times the sample size of Lorimer et al. (2015).The normalized period histogram of Sample A and that of Lorimer et al. (2015) are shown in Fig 1 .It is clear that Sample A has a higher portion of faster MSPs than Lorimer et al. (2015).Since dispersion measure (DM) also affects the inference of period distribution (see Section 4 for details), Fig. 1 shows the DM value of both samples.Compared with the Lorimer et al. (2015) sample, Sample A has a higher portion of large DM pulsars.
For a more complete picture, we also compile a list of all the radio MSPs in the Galactic field (Sample B, 473 MSPs) and a list of all MSPs in globular clusters 5 (GCs, Sample C, 273 MSPs).Compared with Sample A, Sample B has the fastest Galactic field pulsar (J0952−0607, 1.41 ms, Bassa et al. 2017) and slightly more MSPs below 4 ms, while Sample C has the fastest pulsar J1748−2446ad (1.396 ms, Hessels et al. 2006).Both samples are also shown in Fig. 1.

THE BAYESIAN METHOD OF PARAMETER ESTIMATION AND MODEL COMPARISON
Given a sample of MSP period {d} and an underlying period distribution model M which has a set of parameters, Λ, Bayesian inference (e.g.Thrane & Talbot 2019) constrains the model parameters through P (Λ|d, M) = π(Λ|M)L(d|Λ)/Z, where π(Λ|M) is the prior probability density function of the parameters, L(d|Λ) is likelihood and Z is model evidence.The likelihood is given by where P denotes the period predicted by the underlying distribution and π(P |Λ) is hyper-prior.Note that in the second equality of Eq. 1, we have used L(d|P ) = δ(d − P ), which is an excellent approximation in practice, as the measurement of pulse period is usually very precise.Given selection effects C, however, parts of expected data are filtered out by selection effects, leading to bias in both parameter estimation and model selection.Deriving the biased likelihood, L(d|Λ, C), is thus a key step in Bayesian analysis.To obtain L(d|Λ, C), further using Bayes theorem and Eq. 1, then (e.g.Hinton et al. 2019;Mandel et al. 2019) where p(C|d, Λ) is the selection effect for data {d} and p(C|P, Λ) is the selection effect for a general data set.For the observed data set {d}, p(C|d, Λ) = 1 according to definition, while for a general set of data, a detailed analysis of the selection effect is required and given in Section 4. Note that the selection effect usually has no direct dependence on population parameters, thus p(C|P, Λ) = p(C|P ).When the data are from different surveys, the total likelihood becomes where i is the index of i−th survey.By sampling over the parameter space of the likelihood, one can estimate both model parameters and evidence.For two models, M 1 and M 2 , the naturallog-based Bayes factor is usually defined as ln BF M1 M2 = Table 1.The survey parameters and corresponding number of MSPs (NMSP) used in Sample A. The meaning of parameters are explained in Section A.2. Notes: (a) Two different receivers were used in both PALFA and CRAFTS, but most PALFA MSPs were discovered using Mock spectrometer system, while most CRAFTS MSPs were discovered by 19-beam receiver.We picked out these MSPs and used the survey settings in corresponding receiver phases.(b) As G of multibeam receivers decreased from the core to outer beams, typical values were used here for simplicity.G is in units of Jy

THE SELECTION EFFECT
For pulsar surveys, the major selection effect is caused by the limited survey sensitivity.The sensitivity depends on the survey parameters and scales as (W/(P − W )) 1/2 , where W is effective pulse width and P is the pulse period.The effect is thus period-dependent and affects the inference of underlying period distribution.Following Lorimer et al. (2015), the effect can be modelled in two steps: first deriving the detection rate of a given MSP, D smear , then estimating the selection effect, p smear , by fitting D smear with a function of P and DM.Specifically, D smear is given by where p(S) is the probability density function of flux density S, S min is survey sensitivity, and S 0 is the limit of S min at large P and zero DM.The effect is obtained by fitting D smear with Since S min varies from survey to survey, each survey has a set α and β.For readability, we give the fitting results of α and β for the five surveys in Table 2, and show the computation details of D smear in Appendix A.

MODELLING THE PERIOD DISTRIBUTION
Since the underlying period distribution is unknown, we consider ten different models to minimize the possibility of model misspecification.The models are mainly motivated by the observed distribution and can be classified into three categories as described below.

Single-component models
Following Lorimer et al. (2015), we use normal, lognormal, and Γ distribution to model the period distribution of MSPs.The normal distribution is denoted by N (µ, σ), where µ is the mean and σ is the width.The log-normal distribution is defined as where µ and σ are the mean and width, respectively.The Γ distribution is defined as where Γ(k) is the Gamma function and k is shape parameter and θ is scale.We also consider a Weibull distribution, which reads as where both scale θ > 0 and shape k > 0. A Weibull distribution of spin frequency as considered in Patruno et al. (2017) is equivalent to the following where the subscript indicates that the corresponding distribution in the frequency (ν) space is Weibull.The priors of all the parameters are listed in Table 3.The models above are normalized before being applied for analysis.For the normalization, we use 20 ms (the sample selection criterion) as the upper boundary and use a fixed lower boundary of 0.6 ms, which is the Keplerian period limit computed from (Haensel et al. 2009) using M = 2 M ⊙ and R = 10 km.

Lower cut-off models
Theoretically, rotating neutron stars should have a maximum spin frequency, or equivalently, a minimum spin period.However, the exact value of the minimum period is unclear and it is interesting to search for evidence of a minimum spin period in the spin distribution of MSPs.Therefore, we consider models with a lower cut-off, P min , which is set as a free parameter.
The first cut-off model is modified from the normal distribution and given by N cut = H(P − P min )N (µ, σ), where H(P − P min ) is the Heaviside function.Similarly, the second cut-off model is a modified version of the log-normal distribution and reads as logN cut = H(P − P min )logN (µ, σ), while the third model is an adaption of Γ(k, θ): Γ cut = H(P − P min )Γ(K, θ).Finally, we also consider the once favoured power-law model, which also has a lower cut-off and reads as L Pow = H(P − P min )P γ .
In all the cut-off models, we use a uniform prior, U (0.1, 2) ms, for P min .The priors of other parameters are also uniform and given in Table 3. Slightly different from the single-component models, all the cut-off models are normalized in the range of [P min , 20 ms].

Two-component models
The population of MSPs can be divided into different groups, e.g., isolated and binary MSPs, it is reasonable to speculate if they follow different period distributions.To study this possibility, we consider three two-component models.
The first model consists of two normal distributions: , where w is weight and we set µ 1 < µ 2 to avoid confusion between the two components.Similarly, the second model is a mixture of normal and log-normal model, L (N +logN ) = wN (µ 1 , σ 1 ) + (1 − w)logN (µ 2 , σ 2 ), while the third model is a mixture of the normal and Gamma models: L (N +Γ) = wN (µ, σ) + (1 − w)Γ(k, θ).For all the three models above, a uniform prior, U (0, 1), is given to w.The priors for other parameters are also uniform and given in Table 3.The models are normalized in the range of [0.6, 20] ms, the same as that of singlecomponent models.

RESULTS
Using the Bayesian framework and the ptemcee sampler6 (Vousden et al. 2021), we obtain the posterior distributions of the parameters and corresponding model Bayesian evidence, which are shown in Table 3.In the analysis of Sample A, we consider both the case with and without selection effect modelling.We also analyze Samples B and C, which contain all the radio MSPs in the Galactic field and globular clusters, respectively.Due to the complex discovery history in Sample B and the significantly different spatial distribution in Sample C, we only considered the case without selection effect and show the results in Table 3.
Table 3 shows that, no matter the selection effect is modelled or not, both the W ν and logN cut model are strongly favoured by Sample A, although logN cut is marginally (ln BF ≤ 0.9) better.While using Sample B, W ν is strongly favoured over logN cut with ln BF = 8.3.In Sample C, W ν remains the best model, but the logN cut model cannot be ruled out as the relative ln BF is only 1.4.In all cases, two-component models are strongly disfavoured.
For posteriors, there are only slight changes for the parameters of W ν for all the cases, as exemplified in Fig. 2 for Samples A and B. The value of both parameters (shape k ∼ 2 and scale θ ∼ 0.29 ms −1 = 290 Hz) are consistent with that of Patruno et al. (2017).The fitness of the W ν and logN cut model to Samples A and B are shown by the posterior predictive distributions in Fig. 3.
Using the posterior predictive distribution of the W ν model, we computed the fraction of sub-millisecond pulsars (f sub ) and the fraction of fast pulsars, i.e. pulsars with P ≤ 1.396 ms, f ≤1.396 .Both fractions are shown in Table 4.
We also tested the robustness of our method by simulating mock samples using the W ν model (see Appendix B for details).Our analysis successfully recover the injected parameters and the ln BF decisively favour the true model over others.

Robustness of the results
The results can be affected by several factors, including selection effects and the choice of MSP samples.We discuss these factors below.
First, using Sample A, one sees that for both logN cut and W ν model, the selection effect due to pulse broadening has very limited impacts on both the posteriors (Fig. 2) and the log Bayes factors (Table 3).The small impacts of selection effects, especially for samples discovered with modern instruments, were also noticed by Lorimer et al. (2015).The small selection effects thus warrant us to extend the analysis from Sample A to Samples B and C.
Second, in addition to pulse broadening, beaming also introduces period-dependent bias to the observed MSP sample (e.g.Lyne & Manchester 1988;Emmering & Chevalier 1989;Lorimer 2008), as it affects the probability of a pulsar beaming toward the Earth hence its detectability.For MSPs, there is no empirical relation of beaming fraction; however, the fraction is believed to be between 0.4 and 1, and possibly close to unity (Kramer et al. 1998;Levin et al. 2013).More importantly, faster pulsars usually have a higher beaming fraction and thus are less affected by beaming.Therefore, the impact of beaming on the constraint of P min should be small.
Additionally, binary motion modulates the observed spin frequency and leads to period-dependent loss of signal-to-noise ratios (Lorimer 2008).This impact has been largely overcome in recent surveys, especially for wide-orbit binaries, by the extensive use of acceleration search.For the impact on compact binaries with a small spin period, it remains to be investigated.However, no significant difference was found in the period distribution between the isolated and binary systems, for both Lorimer et al. (2015) and our sample (see Appendix C for further details).
Finally, the parameter posteriors are relatively stable with respect to the sample in use, but the Bayes factors are clearly affected by the sample size.To see this more clearly, we also analyze the Lorimer et al. (2015) sample for the case with and without considering selection effects in their paper.In both cases, the logN cut model is marginally better (ln BF ∼ 2) than the W ν model.Taking this into account and considering the large size of Sample B, we thus tend to believe, for radio MSPs, W ν is more likely to be the underlying period distribution7 than the other models considered here.

Sub-millisecond pulsars and equation of state
Irrespective of the sample and selection effect, the fraction of sub-millisecond pulsars obtained for the W ν model remains < 10 −5 .The small f sub is consistent with the lack of fast pulsars close to the break-up limit, a puzzle noticed in both nuclear-and accretion-powered MSPs (Chakrabarty et al. 2003;Patruno 2010), which are free from the selection effect (Hessels et al. 2006) but have a much smaller sample size.
Since the SKA will likely find thousands of new MSPs (Smits et al. 2009a;Keane et al. 2015), we also compute the expected number of sub-millisecond pulsars and fast pulsars.Given the yield uncertainty, we assume a low case of 1500 MSPs discoveries and a high case of 6000 discoveries by the SKA.In either case, the discovery of sub-millisecond pulsars is unlikely, with an expectation < 0.1, while the chance to discover fast pulsars looks higher, with expectations ranging from ∼ 1 to around 10 discoveries.The lack of pulsars faster than ∼ 1 ms has attracted great attention and several models have been proposed, including spin-down via gravitational radiation at extremely small spin periods (e.g.Chakrabarty 2008), and the failure to spin-up to the break-up limit during the recycling phase (e.g.Bhattacharyya & Chakrabarty 2017;Ertan & Alpar 2021;C ¸ıkıntoglu & Ekşi 2023).
The lack of sub-millisecond pulsars will have an important impact on the study of dense matter equations of state: a practical equation of state must yield Kepler frequencies not less than the observed spin frequencies (Watts et al. 2015).The possible discovery of sub-millisecond pulsars has long been seen as a strong evidence of the existence of quark stars (Glendenning 1990).However, as shown in Fig. 4, the result of P min ≥ 1 ms, i.e. f max ≤ 1000 Hz, does not allow for discrimination between different models of equation of state, or prove the existence of quark stars in this way.

Do different groups of MSPs follow the same period distribution?
Usually, population studies treat GC and field pulsars separately, partly due to the apparent difference in spatial distribution and formation history or environments.Their spin periods, however, seem to be well described  (Zhou et al. 2018;Zhu et al. 2018;Li et al. 2020, and references therein).
The observational constrains are from GW170817 (Abbott et al. 2018) and PSRs J0030+0451 (Miller et al. 2019;Riley et al. 2019) and J0740+6620 (Miller et al. 2021;Riley et al. 2021).The black lines are the lower boundaries of the mass-radius space, using the current maximum spin rate (716 Hz, Hessels et al. 2006) and an assumed maximum rate of 1000 Hz.The range for Pmin = 1.43 +0.03 −0.09 from the logNcut model is also indicated near the 716-Hz curve.
by the W ν model with comparable parameters (k ∼ 2 and θ ∼ 0.3 ms −1 ), although for GC pulsars, Sample C  has a weaker preference for W ν over logN cut , probably due to the smaller sample size.Another interesting issue is whether isolated and binary pulsars follow the same period distribution.We test the hypothesis with the isolated (N = 148) and binary (N = 200) radio MSPs in the Galactic field (taken from PSRCAT V1.70).Using the best-fit model, W ν , both samples give consistent posteriors within 1σ uncertainty, see Fig. 5.

Conclusions
Analyzing the spin period distribution of known MSPs, we find that a Weibull model best describes both MSPs from the selected surveys and all the radio MSPs in the Galactic field, with the best-fitting parameters of k = 2.09 +0.09 −0.09 and θ = 0.30 +0.01 −0.01 ms −1 , in agreement with Patruno et al. (2017).Moreover, we find that the spin periods of MSPs in globular clusters also follow a Weibull distribution with consistent parameters.We find no statistical support for a minimum spin period cutoff in the spin distribution of MSPs.For the MSPs in the Galactic field, our analysis shows that both the isolated and binary pulsars follow the same parent period distribution.We also incorporate selection effects due to pulse broadening into the Bayesian analysis, and find that they are negligible.
Based on the Weibull model as inferred for the known MSP population, we find that the SKA will likely find a few pulsars faster than 1.396 ms (PSR J1748−2446ad).However, the fraction of sub-millisecond pulsars is < 10 −5 , which means an extremely low chance for the SKA to discover such pulsars.The lack of sub-millisecond pulsars suggests that observed pulsar spins might not be an effective way to discriminate unusual equations of state of neutron stars.) where S and D are in units of mJy and kpc, respectively.Given Galactic longitude (l) and latitude (b), the term p(D) is (Lorimer et al. 2006;Verbiest et al. 2012)

A.2. The survey sensitivity and detection rate
Another key ingredient in computing the selection effect is S min .Following Dewey et al. (1985); Lorimer & Kramer (2012), S min can be estimated by where S/N is detection threshold, β is a factor describting digitization performance (e.g.Jenet & Anderson 1998), T sys is system temperature, G is telescope gain, n p is the number of polarization, t int is integration time, B is receiver bandwidth, W e is the effective pulse width and P is pulse period.The value of all these parameters (except T sys and W e ) is listed in Table 1.Below, we give the prescriptions to obtain T sys and W e .
The system temperature is given by T sys = T rec + T sky , where T rec and T sky are receiver and sky temperature, respectively.We computed T sky using the Haslam sky model (Haslam et al. 1981) implemented in the PyGDSM package8 (Price 2021), which uses for extrapolation a spectral index of −2.85 between 400 MHz and a few GHz (Dickinson et al. 2019).
The effective pulse width is given by where W i is the intrinsic pulse width, t samp is sampling time, DM 0 is the DM value that smears the pulse in one bin to t samp and t scat is scattering time.For W i , we use an 8 per cent duty cycle, i.e.W i /P = 0.08, which is the median value of the MSPs in our sample.DM 0 is given by where N chan is the number of frequency channel.For t scat , we use the result from Cordes (2002)  Note that Eqn.A4 is prone to overestimating the flux density at small period for large DM values (Cordes & Chernoff 1997;Lazarus et al. 2015), as S min increases rapidly when W e approaches P .This feature will bias the estimation of selection effect, as it will reduce the value of the numerator in Eqn. 4, hence D smear .A more appropriate, though complicated, way to estimate S min was given by Cordes & Chernoff (1997).Since most MSPs in Sample A have relatively small DM thus small degree of overestimation, we use Eqn.A4 for clarity.
Using the results above and Eqn. 4, we computed D smear for all the five surveys and show the results in Fig.

B. ROBUSTNESS OF THE BAYESIAN METHOD
To test the robustness of our method, we firstly generate 473 spin periods with the model W ν (2, 0.3), simulating both the sample size and posteriors of Sample B. For simplicity, no selection effect was considered.The model and mock data are shown in Fig. 7.We then analyze the mock data with the W ν and logN cut model in the same way as the analysis of real data.The log Bayes factor is 33.6, giving a decisive preference for W ν over logN cut , which is clearly shown in the reconstructed results in Fig. 7.The posteriors are also consistent with the injected model parameters.The simulations thus prove that our method is robust in both parameter estimation and model selection.

C. OTHER SELECTION EFFECTS
Binary motion may also induce a selection bias, as the motion modulates the apparent spin frequency and reduces the chance to discover binary pulsars.Specifically, the frequency shift is aT /(P c), where a is the acceleration, T is observation length, P is the spin period and c is the speed of light (Lorimer 2008).For the impact on period analysis, most recent pulsar surveys have done extensive acceleration searches and have found a lot of binary pulsars.The efforts thus have significantly reduced the impact, especially for binaries with wide orbits.Actually, our analysis find no significant difference in the period distribution between isolated and binary systems (see Section 7.3).While for compact binaries and for short spin periods, the effect may still be severe.An adequate treatment of this effect requires further investigation but is beyond the scope of this paper.

Figure 3 .
Figure 3.Comparison of the posterior predictive distribution between the logNcut and Wν model.The shaded area indicates the 90% confidence regions, while the dashed lines show the histograms of the samples.Left: The distributions for Sample A considering the selection effect due to pulse broadening.Right: Results for Sample B without modelling the selection effect.

Figure 4 .
Figure4.Mass-radius relations of neutron stars (with or without phase transition in cores; grey solid curves) and quark stars (in non-superfluid or superfluid states; grey dashed curves) for different equations of state(Zhou et al. 2018;Zhu et al. 2018; Li et al. 2020, and references therein).The observational constrains are from GW170817(Abbott et al. 2018) and PSRs J0030+0451(Miller et al. 2019;Riley et al. 2019) and J0740+6620(Miller et al. 2021;Riley et al. 2021).The black lines are the lower boundaries of the mass-radius space, using the current maximum spin rate (716 Hz,Hessels et al. 2006) and an assumed maximum rate of 1000 Hz.The range for Pmin = 1.43 +0.03 −0.09 from the logNcut model is also indicated near the 716-Hz curve.

Figure 5 .
Figure 5.The posteriors of the isolated and binary radio MSPs in the Galactic field using the Wν model.
A3) where R = R 2 0 + (D cos b) 2 − 2R 0 D cos b cos l, and R 0 is the Solar distance to the Galactic centre.We use R 0 = 8.2 kpc (Gravity Collaboration et al. 2020).Numerically integrating Eqn.A1 gives the p(S) for a given set of (l, b) or line-of-sight.Note that the adoption of Eqn.A2 restricts us to use MSPs discovered by surveys at ∼1.4 GHz.

Figure 6 .
Figure6.The detectability rate of Sample A pulsars from the five surveys and the fitting results to psmear.The blue dots are the computed Dsmear and the yellow mesh is the fitting result with the fitting error indicated by vertical red lines.

W
Figure 7. Left: The mock data, the injected and reconstructed Wν model, and the reconstructed logNcut model.Right: Comparision between the posteriors of the Wν model and the injected parameters used in generating the mock data.
Jeffreys 1998, where Z M1 and Z M2 are the evidence for each model.A positive log Bayes factor indicates preference of M 1 over M 2 (e.g.Jeffreys 1998).

Table 2 .
The coefficients of psmear for the surveys in Sample A.

Table 3 .
Summary of period distribution models used in this study.Columns are model names, parameters, priors, maximum a posteriori values (with 1-σ error) and natural-log Bayes factors (with respect to that of the normal model).

Table 4 .
The fraction of sub-millisecond and fast pulsars (P ≤ 1.396 ms) for different cases (considering different samples and the inclusion of selection effects or not) using the Wν model.Also listed are the corresponding expected discoveries of sub-millisecond and fast pulsars, assuming a low and high yield of the SKA pulsar survey.Posterior distributions of the Wν model.For Sample A, the results with selection effect modelling are in blue, while those without the effect are in yellow.The results for Sample B without the effect are in green.