On the Initial Spin Period Distribution of Neutron Stars

We derive the initial spin period distribution of neutron stars by studying the population of young pulsars associated with supernova remnants. Our hierarchical Bayesian approach accounts for the measurement uncertainties of individual observations and selection effects. Without correcting for selection effects, as done in previous studies, we find that pulsar initial spin periods follow a Weibull distribution, peaking at 40 ms, which is favored against a lognormal distribution with a Bayes factor of 200. The known selection effects in radio pulsar surveys, including pulse broadening and period-dependent beaming fraction, have been quantitatively investigated. We show that, based on measurements of pulsar luminosity and spin period from the ATNF Pulsar Catalogue, the impact of pulse broadening on the inference of the pulsar period distribution is likely to be insignificant. Correcting for the beaming selection effect, the Weibull distribution remains the preferred model, while its peak slightly shifts to longer periods at 50 ms. Our method will prove useful in constraining the birth properties of neutron stars in the Square Kilometre Array era.


INTRODUCTION
Pulsars are fast rotating, highly magnetized neutron stars (NSs).Since the first discovery (Hewish et al. 1968), the number of observed pulsars has grown to over 3500 (see the ATNF Pulsar Catalogue1 , Manchester et al. 2005), a majority of which are detected in the radio band.Among them, pulsars observed in association with supernova remnants (SNRs) are of particular interest, since their ages are independently informed by observations of SNRs.The Crab pulsar is the best known example.Firmly established as the remnant star of supernova 1054, its initial spin period is around 20 ms, close to its current spin period of 33 ms.This is widely used as a proxy for pulsar initial spin periods in the community (e.g., Johnston & Karastergiou 2017).
The astrophysical processes that give rise to NS spins are poorly understood.A range of spin periods from milliseconds to seconds are predicted in a variety of processes during supernova explosions.Newborn NSs could inherit the angular momentum of progenitor stars from the collapsing iron cores (Heger et al. 2005;Ott et al. 2006), where the angular momentum transport and mass loss in single stars play a significant role (Fuller et al. 2014(Fuller et al. , 2015(Fuller et al. , 2019;;Ma & Fuller 2019;Eggenberger et al. 2019;Hu et al. 2023).For instance, Ott et al. (2006) found that an NS could be born with periods of tens to hundreds of milliseconds if the spin periods of iron cores are around 50−100 s.NS spins could also stem Du et al.
from the natal kicks caused by the asymmetric mass ejection and anisotropic neutrino emission (e.g., Spruit & Phinney 1998;Ng & Romani 2007;Janka et al. 2022;Coleman & Burrows 2022;Fragione & Loeb 2023;Burrows et al. 2023).In particular, Spruit & Phinney (1998) showed that an off-center kick at a speed of several hundreds km/s could lead to a spin period as short as tens of milliseconds.Other processes that might play a role in producing NS spins include hydrodynamic instabilities during the supernova explosions such as the standing accretion shock instability (e.g., Blondin & Mezzacappa 2007;Guilet & Fernández 2014;Kazeroni et al. 2016), the anisotropic accretion of angular momentum during the pre-explosion phase (Wongwathanarat et al. 2013) and post-explosion phase (e.g., Janka et al. 2022).Therefore, the initial period distribution of NSs is a powerful probe into astrophysical processes during their formation.
Thus far, around a hundred pulsars have been detected in association with SNRs.They are considered to be young (with a typical age less than ∼ 10 5 years), and therefore provide useful insights into the birth properties of NSs, such as the spin period, magnetic field, spatial velocity, and inclination angle between the spin and magnetic axes (Popov & Turolla 2012;Malov 2021;Igoshev et al. 2022).To properly derive the initial distribution of pulsar parameters with the observed pulsar-SNR population, two approaches can be used.First, by performing population syntheses (Emmering & Chevalier 1989;Faucher-Giguère & Kaspi 2006;Bates et al. 2014;Gullón et al. 2014;Cieślar et al. 2020;Dirson et al. 2022), one develops a model pulsar population based on plausible assumptions about their initial properties and time evolution, and then runs the synthetic population through some mock pulsar surveys with an aim to reproduce the observed pulsar sample.The model that best reproduces observations is considered to be a good model.In the second approach, one extrapolates, for example, the initial periods of pulsars back from their measured spins and age estimates with assumptions about spin evolution (e.g., Xu et al. 2023).Then, statistical inference is performed to fit the distribution of initial periods by incorporating selection effects and measurement uncertainties in individual observations.The major challenge in this approach is an appropriate quantification of the possible selection effects that would bias the observed distribution from the true population (e.g., Lorimer et al. 2006Lorimer et al. , 2015)).Igoshev et al. (2022) performed maximum-likelihood estimates for the distributions of initial periods and magnetic fields of NSs using a sample of Galactic pulsars detected in associa-tion with SNRs.However, the selection effects were only qualitatively discussed.
In this work, we revisit the inference of the initial period distribution of NSs.For Galatic pulsars detected within SNRs, we derive the initial period of each pulsar from the measured spins, using the SNR age and assuming a generic spin-down model.We then fit the distribution of initial spin periods through hierarchical Bayesian inference.Our analysis accounts for the uncertainties in age estimates of SNRs, effects of different braking indices, and selection effects in radio pulsar surveys.We also perform Bayesian model selection to determine the best functional form of the initial period distribution.
This paper is organized as follows.In Section 2, we describe the observed sample, our data selection criteria and the spin-down evolution of radio pulsars.Section 3.1 presents the statistical framework.The results and discussions of various effects that might affect our results are described in Section 3.2.We present concluding remarks in Section 4.

DATA AND MODEL
In this section, we first describe the compilation of pulsars observed within SNRs and our data selection criteria.Subsequently, we present the model of spindown evolution of pulsars.

Pulsars in supernova remnants
Recently, Igoshev et al. (2022) compiled 68 pulsars observed in possible association with SNRs by matching the Galactic SNR Catalogue2 (Ferrand & Safi-Harb 2012) with the ATNF Pulsar Catalogue.In this work, we add 20 Galactic pulsar-SNR pairs to the sample adopted by Igoshev et al. (2022).These sources can also be found in the Galactic SNR Catalogue (SNRcat) and have been analyzed in previous works (Popov & Turolla 2012;Malov 2021).In Appendix A, we list these 88 Galactic pulsar-SNR pairs together with their parameters.The association of the Monogem Ring with PSR J0659+1414, which is uncertain in the SNRcat, has recently been confirmed by Yao et al. (2022).
To derive the initial spin period (P 0 ) of pulsars from observed spin parameters and their age estimates by assuming a spin-down model, we select pulsars in our sample based on the following criteria.The pulsars must have their period (P ) and period derivative ( Ṗ ) measured.The values of P and Ṗ are taken from the ATNF Catalogue (version 1.70).We assume that the age of SNR (τ snr ), taken from SNRcat, is the pulsar's true age unless otherwise specified (see Section 2.2 for further discussion).For PSR J1801−2451, the SNR age is unavailable, thus we adopt the kinematic age determined based on the measurements of its proper motion and position (Noutsos et al. 2013).
The pulsar should be uniquely paired with the SNR.Since most of these pulsar-SNR sources are located near the Galactic plane, there might be some overlaps between them by chance.Specifically, there are two issues: (1) an isolated pulsar is located close to multiple SNRs and (2) an isolated SNR is located close to multiple pulsars.We tackle the first issue by only selecting the pairs of which the SNR age is closest to the characteristic age of pulsar.Therefore, the pulsar J1640−4631 will be identified to be associated with SNR G338.3−00.0.To address the second issue, we select the pulsar whose characteristic age is closest to the SNR age.For example, in the case of SNR G035.6−00.4,three pulsars, J1857+0143, J1857+0210, and J1857+0212, were detected close to this remnant; we select PSR J1857+0143 as in association with G035.6−00.4.The included and excluded pairs are labeled as 'Y' and 'N' in the column "Included" of the table in Appendix A, respectively.
Following Igoshev et al. (2022), we exclude magnetars in our sample.Thirteen Galactic pulsars in SNRs are known as magnetars (see Appendix A ), which have been reported in the McGill Online Magnetar Catalogue 3 .Magnetars could represent a distinct population different from normal radio pulsars (e.g., Gullón et al. 2015).The energy source powering magnetars could be more complicated than that of rotation-powered pulsars; the dominant energy loss could be in the forms of, e.g., the decay of magnetar's immense dipole field (Duncan & Thompson 1992) and the significant gravitational radiation due to their large magnetic deformation and (initially) fast rotation (Ioka 2001;Dall'Osso et al. 2009).Regardless of the loss mechanism, magnetars are thought to experience much more efficient spin deceleration before the onset of dipole spin-down, with substantially strong magnetic fields and a wide range of initial periods (e.g., Thompson et al. 2004;Prasanna et al. 2022).
Two pulsars, J0821−4300 and J1210−5226, known as central compact objects (CCOs), are included in our sample.The detection of X-ray pulsations provides strong support for the argument that they are pulsars (Hui & Becker 2006;Gotthelf et al. 2013).The period derivatives of these two objects are measured with phase-coherent X-ray timing by incorporating the measurements of position and proper motion (Gotthelf et al. 2013).Note that Igoshev et al. (2022) did not include the CCOs in their analysis of the initial period.
We also require SNR age less than the characteristic age of pulsars, see the discussions in section 2.2.After applying the aforementioned selection criteria, we are left with 39 pulsar-SNR pairs available for the inference of initial spin period distribution; see Appendix A for details.
We show the pulsar P -Ṗ diagram in Figure 1.The NSs associated with SNRs are denoted by stars: magnetars (cyan) populate the slowest and strongest magnetic filed regions among all species of neutron stars; normal pulsars and two CCOs are highlighted in red and yellow, respectively.In the upper panel of Figure 2, we show the observed spin period distribution of the pulsars (orange histogram) included in our sample.

Spin-down evolution of pulsars
The pulsar spin-down evolution can be modeled with the following torque equation (Manchester & Taylor 1977): where Ω and Ω are the spin angular frequency and its first time derivative, respectively.These two quantities can be accurately measured from pulsar timing observations (Groth 1975).The braking index is defined as  to 3. A value (significantly) differing from n = 3 indicates more complicated and separate braking processes in the rotational evolution, e.g., an evolving dipolar field (e.g., Gao et al. 2017 and references therein), angular momentum loss through particle wind flows (Harding et al. 1999), or gravitational radiation due to multipoles (e.g., Bonazzola & Gourgoulhon 1996;Owen et al. 1998;Miller et al. 2019;Riley et al. 2019).
For a constant n, we can rewrite Equation (1) in terms of P = 2π/Ω and Ṗ = −2π Ω/Ω 2 to describe the spindown evolution of the pulsars on the P -Ṗ diagram, i.e., dP dt ∝ P 2−n . (2) The integration of this equation gives the solution where τ p is the true age of the pulsar, and P 0 is its initial period at τ p = 0. We have assumed in Equation ( 3) that Ṗ is constant over the time-scale of τ p for a young pulsar.In the special case of n = 1, i.e., where the torque is dominated by a stellar particle wind, the solution reduces to with τ c = P/2 Ṗ being the characteristic age of a pulsar under the assumptions of P 0 ≪ P and n = 3. Equations ( 3) and ( 4) allow us to estimate the initial period of each pulsar by utilizing their measured values of P , Ṗ , and n, along with their SNR ages.As P and Ṗ are accurately measured through pulsar timing observations, their uncertainties can be neglected in estimating the initial periods.The long-term braking index, n, has been measured for a few pulsars (Espinoza et al. 2017).
For pulsars without a measured braking index, we first assume a standard value of n = 3 and will examine the effect of n ̸ = 3 in our discussions.Since Ṗ is small and n is generally measured as 0 < n < 3 (see Espinoza et al. 2017, and references therein), the estimated P 0 is not significantly affected by the choice of n, allowing us to disregard its measurement uncertainty.Finally, the main source of uncertainty in estimating the initial period comes from τ snr , which has been determined using different methods for different SNRs (Suzuki et al. 2021).Some SNRs, e.g., G119.5+10.2 and G184.6−5.8, have well-recorded "historical ages" with negligible errors, allowing us to estimate their initial periods with a Dirac δ distribution.For SNRs with ages determined within certain ranges, we estimate the initial period as a uniform distribution over the derived lower and upper bounds.
To avoid imaginary or zero periods (P 2 0 ≤ 0) when taking n = 3 in Equation (3), τ snr < τ c must be satisfied, which is an additional selection criterion for our sample 4 .There are six pulsars with n measurements in this sample.In the upper panel of Figure 2, we show the derived initial periods of the selected pulsars (blue histogram).

Method
With the individual estimates of the initial spin period for N obs = 39 pulsars, we determine the initial spin period distribution through a hierarchical Bayesian inference (Mandel et al. 2019;Thrane & Talbot 2019).We apply the Bayes' Theorem to infer the posterior distribution of the population parameter (Λ) from an ensemble of observations {D} = {P, Ṗ , n, τ snr }, where p(Λ) and p(N ) are the priors for Λ and N respectively, with N being the total number of a population over the observation period.The inclusion of N in Equation (5) allows us to infer the event rate.The probability of detecting the N obs pulsars in association with SNRs can be modeled with an inhomogeneous Poisson process.The first term in the right-hand-side of Equation ( 5) defines a population-level (hierarchical) likelihood (Thrane & Talbot 2019), marginalized over the individual measurements (θ): Here, π(θ|Λ) is the population model; π ϕ (θ) is the prior utilized for individual-source Bayesian analysis; N s represents the number of discrete samples of individual measurements (drawn from a δ function or a uniform distribution, see section 2.2); and ξ(Λ) accounts for the fraction of detectable sources in a population, defined with The probability, p det (θ), denotes the detection probability of observing telescopes.It accounts for the selection bias (see Section 3.2.3).
In this work, we restrict our inference to the shape of the population rather than the rate.Therefore, we marginalize Equation (5) over N by using a log-uniform prior p(N ) ∝ 1/N , which does not affect the inference of Λ, resulting in p(Λ|{D}) ∝ p(Λ) Throughout, we take π ϕ (θ) as a uniform prior.The Bayes factor, which is the ratio of Bayesian evidence, can be used to arbitrate the preference of one population model over another by the data.Here, we assume that the models (H 0 as the null-hypothesis model and H 1 as the alternative model) in question are equally probable before acquiring any knowledge from the data.
We utilize the open-source software BILBY (Ashton et al. 2019) to perform the Bayesian analysis.The posteriors of free parameters are generated by the sampler using the nested sampling techniques implemented in the python package DYNESTY (Speagle 2020), whilst estimating the evidence.We use the Python code ChainConsumer5 to analyze the posterior samples with maximum likelihood statistics by interpolating a Gaussian kernel density function.The results are quoted at the 1-σ credible intervals (i.e., the 68.3% areas below the maximum posterior density) unless noted otherwise.

Results and Discussions
We fit the initial spin period distribution using five parametric functions as population models: the Gaussian (GS), log-normal (LOGN), turn-on-power-law (TOPL), GAMMA, and Weibull distributions.These functions are described in Appendix B. We explore the posteriors of population parameters using uniform priors.The range of the priors and the definition of each parameter are tabulated in Table 1.

Without correcting selection effects
We first assume that the observed pulsars are representative of the entire population without any selection bias, i.e., there is no preference of the detactability for pulsars with different periods.Table 1 provides the 1σ credible intervals of marginalized posterior probability distributions of population parameters corresponding to each model.In Table 2, we present the Bayes factors of different models by taking the GS case as the null hypothesis.We obtain decisive evidence for rejecting the GS model.This result is in agreement with that obtained by Igoshev et al. (2022), who found a log-normal distribution provides a better fit to the measurements than the Gaussian based on the Akaike information criterion.The Weibull distribution emerges as the most favored, being moderately (strongly) preferred over the TOPL (LOGN) by a Bayes factor of 7 (276).Figure 2 plots the posterior predictive distribution (PPD) defined as: The PPD of the Weibull and LOGN model are shown in green and grey, respectively.The shaded areas denotes the 90% credible interval.The Weibull distribution, which is the preferred model, indicates that the majority of the observed pulsars were initially fast spinning, with initial periods peaking at ∼40 ms.The 90% credibility upper limit is constrained to be P 90% 0 ≈ 0.5 s.The cumulative distributions of the Weibull and LOGN are shown in the lower panel of Figure 2, where one can see that the Weibull model provides a better fit than the LOGN.In Appendix C, we compare the initial period distribution inferred with the LOGN model with that obtained by Igoshev et al. (2022).
To additionally provide a quantitative measure of the choices of the distribution functions, we perform a twosample Kolmogorov-Smirnov test implemented by the scipy package (Virtanen et al. 2021)  samples are (1) the initial periods from 39 pulsars (as presented in section 2) and (2) 100 synthetic initial periods drawn from each distribution with the median values of the parameters in Table 1; see the column "no selection (n = 3)".The values of test statistic, defined as the largest deviation between the two compared cumulative distributions, and p-value are presented in Table 3.We observe from Table 3 that the Weibull is the best model; the TOPL and LOGN models come next; and the GAMMA and Gaussian models are disfavoured.A Kolmogorov-Smirnov test repeated with 200 synthetic initial periods gives consistent results.Overall, the Kol-mogorov-Smirnov test results support Bayesian model selection results.

Effects of n ̸ = 3
In the analysis above, we have assumed n = 3 for pulsars without measured brake indices.However, longterm timing observations of young pulsars usually give n ≲ 3 (Espinoza et al. 2017), and n > 3 are also measured for many glitching or non-glitching pulsars (Parthasarathy et al. 2019;Lower et al. 2021).Such measurements hint at various braking mechanisms or unknown noises in the interior of pulsars during their evolution, which inevitably induce uncertainties in the n measurements.
Here, we explore the dependence of the inferred initial period distribution on different values of n.We randomly draw n values from a Gaussian distribution with the mean and standard deviation taken from measured values for pulsars with known n.For pulsars without a measured braking index, we randomly draw n values from a Gaussian distribution centered at 3, with a standard deviation of 1, constrained to 0 ≤ n ≤ 4.
The resultant posteriors, given in Table 1, are consistent with the results in the case of n = 3.The Bayes factors, presented in Table 2, also indicate that the adopted values of n does not impact our results significantly.From the initial period distribution, shown in Figure 3 (in orange), we find that the Weibull distribution is consistent with that obtained in Section 3.2.1 at 90% credibility.Therefore, our results are insensitive to the actual values of n.

Correcting for selection effects
The possible selection effects related to radio pulsars, SNRs, and CCOs have been discussed in detail by Igoshev et al. (2022) and references therein.However, evaluating all the observational effects in our population analysis is challenging, due to different searching strategies and detection methods for different (types of) sources.For simplicity, we focus on the known observational selection effects in radio pulsar surveys, such as pulse broadening and the fraction of radio beam directed towards Earth (Lorimer 2011).We compute the detection fraction, as defined in Equation ( 7), to correct for these selection effects.
The radio selection effects on the spin periods of Galactic pulsar population have been extensively studied using the "current" analysis by analytically estimating the pulsar birthrate as a function of period (Vivekanand & Narayan 1981;Phinney & Blandford 1981;Narayan 1987;Narayan & Ostriker 1990;Lorimer et al. 1993;Vranešević & Melrose 2011).The effects  caused by the limited sensitivity in a particular survey is typically calculated by the scale factors defined as V max /V , where V is the weighted volume in which pulsars are detectable and V max is the weighted volume of the whole Galaxy.With the pulsar-current analysis, Vivekanand & Narayan (1981) proposed an 'injection' of a subpopulation of pulsars with ∼0.5 s (see also Narayan 1987;Narayan & Ostriker 1990).However, the injection has been questioned by Lorimer et al. (1993) and Vranešević & Melrose (2011).In particular, Lorimer et al. (1993) pointed out that the injection could be an artifact induced by including the luminosity selection effects (inverse correlation between luminosity and P ) in the scale factor computations, as which would bias the observed sample towards short periods.
The intrinsic radio-luminosity function remains poorly constrained (see Posselt et al. 2023 and references therein).In Figure 4, we show that there is no significant dependence of the observed radio (pseudo) luminosity on the period for radio pulsars with periods P > 10 ms and B dip ≥ 10 11 G in the ATNF Catalogue.We adopt the observed luminosity distribution to compute the scale factors [see Equation ( D12)].We calculate the scale factors as a function of P based on the sensitivity of the Parkes Multibeam Survey7 .The sensitivity curve of this survey is shown in panel (a) of Figure 5. (The calculation details are described in Appendix D.) Shorter-period and farther pulsars are more susceptible to pulse smearing in the interstellar medium, resulting in lower sensitivities in their detection.However, the sensitivity becomes almost constant for pulsars with longer periods (e.g., P ≳ 0.1 s).The detection probability, estimated from the inverse of the scale factor, is represented by gray dots in the panel (b).The detectability grows slightly with periods for P ≲ 0.1 s, which corresponds to the effect of the sensitivity limit on radio observations shown in panel (a).Overall, the influence of pulse smearing on the shape of the initial period distribution are small, which is consistent with the fact that pulsar luminosity does not depend on spin periods as shown in Figure 4.
Pulsars with shorter periods are found to possess wider radio beams, resulting in higher detectability in pulsar surveys (Lyne & Manchester 1988;Tauris & Manchester 1998).Here, we consider two empirical models, LM88 (Lyne & Manchester 1988) and TM98 (Tauris & Manchester 1998), to quantify the period-beaming fraction correlation (see Appendix E for details), which is plotted in the panel (b) of Figure 5. Incorporating beaming fraction corrections, we obtain similar constraints of population parameters for both models adopted, see Table 1.The Bayes factors given in Table 2 indicate that the Weibull distribution is still the best description of pulsar initial period distribution.In the lower panel of Figure 3, one can see a slight shift towards longer periods, with a peak at around 50 ms.The 90% credibility upper limit is constrained to be P 90% 0 ≈ 0.7 s (0.8 s) for the LM88 (TM98) model.

CONCLUSIONS
We infer the initial spin period distribution of NSs using 39 pulsars detected in association with SNRs.Our hierarchical Bayesian approach accounts for measurement uncertainties and known selection effects in radio pulsar surveys.Using the SNR age as an estimate of the pulsar's true age, we determine the initial spin periods of individual pulsars based on their measured values of P and Ṗ , assuming a generic spin-down model.Assuming the observed pulsar sample is unbiased, we find that their initial spin period distribution is best described as a Weibull distribution (defined in the spin frequency space), peaking at ∼40 ms.The shape of initial period distribution is not affected by the uncertain values of the braking index.We also show that the effect of pulse broadening in pulsar surveys likely has an insignificant impact on the initial period distribution, since the pulsar luminosity appears to be independent of spin periods.After accounting for selection effects due to pulse beaming fraction, the initial spin periods are still best described as a Weibull distribution, with a slight shift towards longer periods and a peak at around 50 ms.
Our analysis represents a first step in uncovering the birth spin distribution of NSs through rigorous Bayesian population inference.There are several caveats or assumptions, which may require further investigations.First, pulsars are assumed to spin down with a constant braking index.Our analysis can be extended to include more complex braking models that allow, e.g., magnetic field decay or decaying magnetic inclination angle.Second, magnetars are assumed to be a distinct population and thus excluded from our calculations.Applying the same analysis to magnetars in SNRs, we obtain initial spin periods P 0 > 2 s, apparently inconsistent with the Weibull distribution found in this work.Third, we exclude pulsars whose characteristic ages are lower than SNR ages because in this case no constraint can be placed on their initial spin periods; either these pulsar are born with O(ms) spins, or a more complicated braking scenario should be considered, or age estimates of SNRs need to be revisited.These problems can be at least partly tackled by including pulsars not in association with SNRs in the analysis as well as using kinematic or thermal age estimates.
We thank the anonymous referee for valuable comments on the manuscript.(V) The Weibull distribution (Patruno et al. 2017;Liu et al. 2024) transformed to the period space: where J and K are known as the shape and scale parameters of the distribution, respectively.The subscript 'ν' means that the corresponding distribution in frequency space is Weibull.The Weibull distribution peaks at K −1 [J/(J + 1)] 1/J .

C. COMPARISON WITH PREVIOUS WORKS
Figure 6 displays the initial periods data along with the log-normal distributions obtained by Igoshev et al. (2022) (in red) and inferred in this work (in green); both without the correction of selection effects.One can see that both sets of results are generally consistent with each other.However, in our work the log-normal model is disfavored.This is due to a different analysis method, and more importantly different data selection criteria, especially the treatment of the observed pulsars with supposedly imaginary initial periods (see Section 2.2).

D. SCALE FACTOR COMPUTATIONS
The scale factor for a pulsar with period, P , is defined as (D9) Here, ρ r represents the space density for the galactocentric radius (r), and ρ z is the density in the height above the Galactic plane (z).We assume that the pulsars distribute uniformly over the galactocentric azimuthal angle (ϕ), which is justified since there is no apparent concentration of pulsars observed in the inner spiral arms of the Galaxy.The selection effects caused by the telescope's limited sensitivity enter in the parameter η(P, L, r, ϕ, z) via the inverse-square law, S ν = Ld −2 , where S ν is the apparent flux density detectable at a center frequency ν center (throughout we take ν center = 1.4 GHz).Given a minimum detectable flux density, S min , in a specific survey, η(P, L, r, ϕ, z) is set to 1 if S ν ≥ S min for a pulsar with period P and luminosity L at coordinates {r, ϕ, z}, or otherwise it is set to 0. The numerator in Equation ( D9) is integrated over the whole volume of the Galaxy, while the denominator is only integrated over the volumes where pulsars can be potentially detectable in any reference surveys.
The integral in Equation (D9) can be numerically computed at a number of values of P utilizing Monte Carlo simulations.For each pulsar with P , S ν and S min are computed at a large number of grid points in the {r, ϕ, z, L} space.Then the numerator, which can be seen as the wighted total volume, is updated by adding the weight ρ(r)ρ(z).The denominator, which can be thought to be the weighted sensitive volume, is updated by adding the weight ρ(r)ρ(z) if S ν ≥ S min .The settings are given below.
Spatial distribution.We use the r-distribution (Yusifov & Küçük 2004) ρ(r) = 37.6 r + 0.55 R ⊙ + 0.55 ) where R ⊙ = 8.5 kpc is the galactocentric distance of the Sun.For the z-distribution, we use the exponential function (Lorimer et al. 2006) Luminosity distribution.We produce the luminosities of model pulsars with the observed luminosity distribution (see Figure 7).We fit the observed radio luminosities (calculated as the product of the mean apparent flux density at 1400 MHz and the square of distance) of 1834 pulsars (P > 10 ms) with a log-normal distribution to obtain where d is in units of pc and n e denotes the free electron density in cm −3 .To account for the errors in estimating the model dispersion measure (DM mod ) by using this electron density model, we take a random error drawn from a Gaussian distribution with a mean of zero and width of 0.2 DM mod .
Instrumental sensitivity and pulse broadening.The sensitivity limits of radio telescopes are determined by the various instrumental parameters, including the system noise temperature (T sys ), antenna gain G, number of polarizations N p , total bandwidth ∆ν, and integration time t int .For a specific survey, the theoretical minimum mean flux density can be characterized by the radiometer equation (Dewey et al. 1985): D14) where ρ s/n = 8 is the signal-to-noise ratio for the detection threshold, and β ≈ 1.5, which includes the loss from one bit digitization (∼ π/2) and other kinds of system losses.We scale the sky temperature measured by Haslam et al. (1982) at 408-MHz to any frequency using the sky background spectrum: where the spectral index −2.6 is taken from Lawson et al. (1987).The sensitivity will drop when the radio photons from the pulsar are located offset about the beam center, which can be accounted for by using a Gaussian pattern to provide Here, G 0 is the antenna gain at the beam center, w is the full width at half maximum of the telescope beam (in arcmin), and R represents the offset from the beam center, randomly drawn from a Gaussian distribution with a mean of zero and standard deviation of w/2.The observed pulse width incorporating the effects of instrumental settings, dispersion, and scattering can be modeled as (Dewey et al. 1985) The first term is the intrinsic pulse width (in second), which is adopted as W int = 0.04P .The second term, τ samp , is the data sampling effects, which account for the details of time resolution of digitized data from the system hardware like antialiasing filters.We take τ samp as the sampling interval (t samp ) of the radio telescope.The third term, τ scatt , is the scatter-broadening time, which arises from the smearing due to multi-path propagation of light in a non-uniform and ionized interstellar medium, given by (Bhat et al. 2004) log τ scatt ms = − 6.46 + 0.154 log DM + 1.07(log DM) 2 − 3.86 log ν center GHz . (D18) The fourth term, τ DM , is the dispersion-broadening time across an individual channel, which is adopted as (Hessels et al.where the channel bandwidth ∆ν chan ≪ ν center is assumed; for the 1.4 GHz surveys, we employ ∆ν chan = 100/256 MHz for DM < 100 cm −3 pc and ∆ν chan = 100/512 MHz for DM > 100 cm −3 pc.The last term, τ ∆DM , is the smearing due to the finite DM step size in the survey, which is neglected here. Pulsar survey.We use the sensitivity limit of the Parkes Multibeam Survey (Manchester et al. 2001) to evaluate the scale factors, its parameters are shown in Table E.The panel (a) of Figure 4 displays the sensitivity curve of this survey.Note that, to obtain the S min in this plot, we have simply adopted a sky temperature of 25 K, the observed pulse width is approximated as [(0.04P ) 2 + t 2 samp + τ 2 DM ] 1/2 , and τ DM ≈ t samp DM DM0 with DM 0 = 28 cm −3 pc.

3Figure 1 .
Figure 1.The P -Ṗ diagram of the observed pulsars.The stars denote the pulsars detected within SNRs, with CCOs and magnetars highlighted in yellow and cyan, respectively.The age and magnetic field strength (B dip ) are calculated in the standard magneto-dipole model.The data set of this work is available at https://github.com/Shen-Shi/NSBirth Population.

Figure 2 .
Figure 2. Upper panel: The distributions of the observed spin periods (orange) and initial periods (blue) of pulsars associated with SNRs.The green (gray) line denotes the posterior predictive distribution of initial period derived from the Weibull (log-normal) distribution (see section 3.2.1),with the shaded area indicating the 90% credible interval.Lower panel: The cumulative distributions of the initial periods and the Weibull and log-normal model are plotted with the same colors as the upper panel.

Figure 3 .
Figure 3.The posterior predictive distributions of NS initial spin periods.Upper panel: The Weilbull distribution derived without correcting for selection effects and by adopting fixed (blue) and random (orange) values of braking index.Lower panel: The Weilbull distribution derived by correcting the beaming fraction using the LM88 (green) and TM98 (gray) models (see Section 3.2.3 for details).The shaded areas indicate the 90% credible intervals.

)Figure 4 .
Figure 4.The observed radio luminosity versus the spin period.The observed radio luminosity is calculated as the product of the mean apparent flux density and the square of distance estimated from dispersion measure.Magenta stars and black dots comprise 583 and 1702 pulsars respectively detected at 400 MHz and 1400 MHz, with P > 10 ms and B dip ≥ 10 11 G in the ATNF Pulsar Catalogue.

Figure 5 .
Figure 5. Panel (a): The sensitivity of Parkes Multibeam Survey.The minimum detectable flux density (Smin) is computed using Equation (D14).Panel (b): The detectability of radio pulsars.The solid and dashed lines account for the beaming fraction of the LM88 and TM98 models, respectively, see Appendix E. The gray dots denote the detectability obtained with the inverse of scale factors (see Appendix D).
Z.-H.Zhu is supported by the National Natural Science Foundation of China under Grants Nos.12021003, 11920101003 and 11633001, and the Strategic Priority Research Program of the Chinese Academy of Sciences, Grant No. XDB23000000.ZQY is supported by the National Natural Science Foundation of China under Grant No. 12305059.XJZ is supported by the National Natural Science Foundation of China (Grant No. 12203004) and by the Fundamental Research Funds for the Central Universities.ZCC is supported by the National Natural Science Foundation of China (Grant No. 12247176 and No. 12247112) and the innovative research group of Hunan Province under Grant No. 2024JJ1006.

Figure 6 .
Figure 6.The log-normal distribution obtained in this work (green) and that in Igoshev et al. (2022) (red).The histograms show the initial periods used by Igoshev et al. (2022) (in blue) along with that analyzed in this work (in orange).

Figure 7 .
Figure 7.The distribution of the observed radio luminosities of 1834 pulsars with P > 10 ms detected at 1400 MHz in the ATNF Catalogue.The red line is obtained by fitting the data to a log-normal distribution (see Equation D12).

Table 1 .
The parameters of each parametric function, their descriptions, priors used in the nested sampling, and posterior credible intervals (1σ) obtained with different considerations of selection effects.We adopt a uniform prior for each parameter.The prior boundaries of 0.0014 s and 24 s are informed by the measured shortest and longest periods in the ATNF Pulsar Catalogue, respectively.

Table 2 .
The Bayes factors against the Gaussian distribution compared to different models with different considerations of selection effects.

Table 3 .
The results of two-sample Kolmogorov-Smirnoff test for different distribution functions.

Table 4 .
Parameters of the candidate pulsars observed in association with supernova remnants in the Galaxy.References for the n

Table 5 .
Parameters of the Parkes multibeam pulsar survey.