Pulsar Population Synthesis with Magnetorotational Evolution: Constraining the Decay of the Magnetic Field

We present a population synthesis model for normal radio pulsars in the Galaxy incorporating the latest developments in the ﬁ eld and the magnetorotational evolution processes. Our model considers spin-down with a force-free magnetosphere and the decay of the magnetic ﬁ eld strength and its inclination angle. The simulated pulsar population is ﬁ t to a large observation sample that covers the majority of radio surveys using the Markov Chain Monte Carlo technique. We compare the distributions of four major observables — spin period ( P ) , spin-down rate (  P ) , dispersion measure, and radio ﬂ ux density — using accurate high-dimensional Kolmogorov – Smirnov statistics. We test two B - ﬁ eld decay scenarios, an exponential model motivated by ohmic dissipation and a power-law model motivated by the Hall effect. The former clearly provides a better ﬁ t, and it can successfully reproduce the observed pulsar distributions with a decay timescale of -+ 8.3 3.03.9 Myr. The result suggests that signi ﬁ cant B - ﬁ eld decay in aged pulsars and ohmic dissipation could be the dominant process.


INTRODUCTION
Since the first pulsar was detected in 1967 (Hewish et al. 1969), to date there are over 3000 known pulsars in the sky according to the Australia Telescope National Facility (ATNF) catalog (Manchester et al. 2005).They are very interesting objects and powerful tools for probing many different aspects of physics, from relativity to gravity.For example, understanding their initial spin period (P 0 ) and magnetic field (B 0 ) distributions could reveal the angular momentum transport process and the origin of B-field during core collapse.Modeling the pulsar population could provide tests to the emission mechanisms and constrain the event rates for the related high-energy astrophysical phenomena, including supernovae, fast radio bursts, gravitational waves, and Gamma-ray bursts.
Normal pulsars lose their rotational energy in the form of particle wind and electromagnetic radiation, which results in spin-down.The spin parameters (P and Ṗ ) thus depend on the B-field strength and configuration.The evolution of the magnetic field could leave observable signatures on the pulsar rotational properties.This is the so-called magneto-rotational evolution scenario (see Gullón et al. 2014).
That being said, the detailed evolution of pulsar Bfield remains an open question.There is indirect evidence that suggests magnetic field decay during the pulsar lifetime.For instance, the strongest B-fields are found in young pulsars and weak fields are found in old millisecond pulsars (MSPs), which is generally consistent with the idea of field decay expected from theories (e.g.Goldreich & Reisenegger 1992;Igoshev et al. 2021).Moreover, some pulsars are found to have braking index deviated from 3, which could be the result of decreasing magnetic field strength (e.g., Viganò et al. 2013;Igoshev 2019) due to Ohmic or Hall effects (Flowers & Ruderman 1977;Goldreich & Reisenegger 1992;Igoshev & Popov 2015;Igoshev 2019) or evolution of the magnetic alignment angle (e.g., Weltevrede & Johnston 2008) by plasma effect (Philippov et al. 2014).
Different theories suggest a wide range of decay time scales from a few 10 5 yr (Igoshev & Popov 2015) to 150 Myr (Bransgrove et al. 2018).Observationally, while the decay rate is difficult to measure directly, there are previous attempts to obtain constraints from popula-tion modeling.Most of the reported values are in the range of a few Myr (e.g.Gonthier et al. 2004;Cieślar et al. 2020), but some found no evidence of decay (e.g., Faucher-Giguere & Kaspi 2006).In this study, we revisit this problem by developing a state-of-the-art population model that includes the latest developments in the field as realistic physics inputs, and fit it to a large pulsar sample to constrain the model parameters.Our model adopts the spin-down formula in a plasma-filled magnetosphere (Spitkovsky 2006;Philippov et al. 2014) with the explicit dependence of α, the updated kick velocity distribution (Verbunt et al. 2017;Igoshev 2020), sophisticated radio beam model (Gonthier et al. 2018), an updated free electron density distribution of the Galaxy (Yao et al. 2017), and the latest Gamma-ray luminosity law and beam geometry (Watters et al. 2009;Kalapotharakos et al. 2017Kalapotharakos et al. , 2019Kalapotharakos et al. , 2022)).In previous works, model comparison is often only limited to pulsars from the Parkes and Swinburne Multibeam (PMB) surveys (Manchester et al. 2001;Edwards et al. 2001) and the High Time Resolution Universe (HTRU) survey (Keith et al. 2010).We attempt to fit to a much larger observation sample that covers all major normal radio surveys.Finally, we employ, for the first time, a correct formula for the KS statistic at high dimension for model comparison (Fasano & Franceschini 1987;Justel et al. 1997).
This paper is organized as follows: in Section 2, we introduce the simulation methods including five parts: birth properties, evolution, emission, detection, and comparison with observations.In Section 3, we report our simulation results for the optimal parameters and the corresponding population properties.The discussion is in Section 4 and the conclusion is drawn in Section 5.

Overview of the simulation and fitting procedure
We adopt the evolution approach for the population synthesis.We follow the basic procedures outlined in Faucher-Giguere & Kaspi (2006) but with more realistic physics and model inputs from recent studies.First, we assign birth properties to newborn pulsars, i.e.P 0 , B 0 , α, location, kick velocity and direction, and viewing angle.Then each pulsar starts two independent evolution processes: magneto-rotational evolution and dynamical motion under the Galactic gravitational potential.The evolution time of each pulsar is randomly assigned from a uniform distribution between 0 and 100 Myr.Since it is computationally more intensive to trace the pulsar motion in the Galaxy, we first calculate the evolution of P , Ṗ , and α, and use the final values to determine if a pulsar is visible.We assume a distance of 1 pc and set a threshold of 10 −6 mJy MHz, which is much lower than the sensitivity of any radio surveys used in this study.If a pulsar is beamed away from the Earth or has radio flux below the above limit, we can safely assume that it cannot be detected and move on to simulate another pulsar.Only for pulsars pass this initial screening we calculate their orbits.We stress that this procedure has no effect on the simulated pulsar population since the two evolutions are independent.At the end of the simulation, we model the radio and Gamma-ray emission based on their final spin property, location, and viewing geometry.A pulsar is only counted if it is detectable by the radio surveys in our list.
In each simulation run, we generate 10 times more detectable pulsars than the observation sample (i.e. a total of 18590 pulsars).This number was chosen to minimize the statistical fluctuation while keeping the simulation efficient.See Appendix B for more details.After each run, we then determine the goodness of fit for a given set of model parameters, by comparing the simulated pulsar distribution with the observations using a high-dimensional Kolmogorov-Smirnov (KS) statistic.Finally, based on the KS statistic, we fit the model parameters using the Markov chain Monte Carlo (MCMC) technique.

Spatial distribution
Normal pulsars are believed to be born in core-collapse supernova of massive stars.Due to short lifetime of the massive stars, the progenitors are not expected to travel far before their birthplace.We can therefore assume that pulsars are born very close to the spiral arms of the Galaxy.We adopt a logarithmic function to depict the four spiral arms of the Milky Way where k and θ 0 are constants from Wainscoat et al. (1992), r xy = x 2 + y 2 , θ = arctan(y/x), where x, y, z are coordinate values in a left-handed Galactocentric Cartesian coordinate system such that the Sun locates at (x, y, z) = (8.5, 0, 0.025) kpc.We then adopt the radial distribution model of surface density (ρ(r xy )) for newborn pulsars (2) where R ⊙ is the distance from the Sun to the Galactic center, and a, b, and R 1 are constants from Yusifov & Küçük (2004).Next, we apply a polar angle correction θ corr and distance adjustment r corr to each pulsar Pulsar population with magneto-rotational evolution 3 following Faucher-Giguere & Kaspi (2006) to blur the pulsar distribution, so that they are not all piled up at the centroids of arms.Finally, a two-sided exponential function is applied to describe the distribution of vertical distance above and below the Galactic plane with a scale height z 0 of 50 pc (3)

Kick velocity
Several direct observation evidence supports that neutron stars receive natal kicks during supernova explosions.For instance, young pulsars can reach a velocity range from tens to over 1000 km s −1 while their progenitors, i.e.O-B stars, have a low velocity of only ∼10-30 km s −1 .Moreover, the scale height of spatial distribution of pulsars in the Galaxy is ∼300 pc (Lorimer et al. 2006), significantly larger than that of the progenitors (∼50 pc).Previous studies try different ways to parameterize the velocity distribution, including bimodal (Arzoumanian et al. 2002), single Maxwellian (Hobbs et al. 2005), and bimodal Maxwellian (Verbunt et al. 2017;Igoshev 2020).In our simulation, we adopt the twocomponent isotropic Maxwellian distribution suggested by a recent work (Igoshev 2020) 2σ 2 dv and ( 4) (5) We fix w = 0.2, σ 1 = 56 km s −1 , and σ 2 = 336 km s −1 , obtained by fitting a sample of young pulsars with characteristic age younger than 3 Myr old (Igoshev 2020).The velocity is assumed to have a random direction drawn from an isotropic distribution.
Note that we do not attempt to fit the pulsar spatial nor velocity distributions, as they are insensitive to the magneto-rotational evolution of pulsars, which is the focus of this work.

Initial spin and magnetic field
We regard the initial period P 0 and magnetic field B 0 as independent from other properties although they could be potential connection (e.g., between P 0 and v; see Ng & Romani 2007).
In the previous population synthesis study, the pulsar initial magnetic field strength B 0 and period P 0 are assumed to draw from log-normal distributions (e.g., Graber et al. 2023).This is supported by recent studies, which show that log-normal is preferred over normal distribution for both B 0 and P 0 (Igoshev et al. 2022;Xu et al. 2023).This is adopted in our simulation where µ P , σ P , µ B , and σ B , are free parameters determined by fitting.Finally, for each simulated pulsar, we assigned a random initial inclination angle α (between the magnetic and rotational axes) and a random viewing angle ξ (between the spin axis and line of sight), both drawn from isotropic distribution in 3D, i.e. cos α, cos ξ ∼ uniform(0,1).

Magneto-rotational evolution
The spin-down process of pulsars converts its rotational kinetic energy into electromagnetic radiation and particle outflow, resulting in an increasing period at a rate of Ṗ .For the simplest model of a rotating dipolar magnetic field in vacuum, Ṗ is related to field strength B at the magnetic equator and the magnetic inclination angle α where I is the pulsar moment of inertia, c is the speed of light, and R is the pulsar radius.Taking the canonical values of R = 10 6 cm and I = 10 45 g cm 2 gives While this standard spin-down formula was commonly adopted in previous population synthesis studies (e.g., Faucher-Giguere & Kaspi 2006;Gonthier et al. 2007;Popov et al. 2010).Recent numerical simulations of pulsar magnetosphere found a more realistic force-free spin-down model of i.e., sin 2 α dependence becomes (k 0 + k 1 sin 2 α), such that even an aligned rotator (α = 0) can spin down (Spitkovsky 2006;Kalapotharakos & Contopoulos 2009;Pétri 2012;Tchekhovskoy et al. 2013;Philippov et al. 2014).During spin-down, the electromagnetic torque exerted by the plasma-filled magnetosphere also makes the magnetic axis to progressively align with the rotational axis (Beskin et al. 1988;Philippov et al. 2014;Igoshev & Popov 2020).We follow Philippov et al. (2014) to model this with where the numerical factor k 0 , k 1 , k 2 have a week dependence on the ratio between neutron star radius and light cylinder radius.For typical pulsar, k 0 ≈ k 2 ≈ 1.0, k 1 ≈ 1.2 (Philippov et al. 2014).As mentioned in the introduction, there is observational evidence suggesting the decay of pulsar's magnetic field over its lifetime.Taking this into account in the population synthesis affects the evolution of Ṗ and α, and hence the pulsar evolutionary tracks in the P -Ṗ diagram.We consider in this work the B-field decay in the crust due to Ohmic loss and Hall effect.The decay in the core is ignored since the core can be considered as a nearly neutral superconductor, which has a very long decay timescale.Ohmic decay is the result of scattering between phonons and impurities of the crystalline lattice (e.g., Pethick & Sahrling 1995;Igoshev & Popov 2020).It can be parameterized by an exponential function where timescale τ Ohm is related to the impurities, conductivity, and spatial scale of the crust (Igoshev et al. 2021).On the other hand, Hall effect redistributes the magnetic energy, leading to non-linear decay as with timescale τ Hall that depends on the initial magnetic field B 0 These two effects can be described by a more general form of B-field decay as suggested by numerical calculations (Colpi et al. 1999;Beniamini et al. 2019;Jawor & Tauris 2022).Here the power-law index α B indicates the decay contribution of Ohmic dissipation and Hall effect.Integrating the equation above, the magnetic field evolution thus follows: It is obvious that when α B = 0, the equation above reduces to Eq. 12 with τ Ohm = 1/f B .For the special case of α B = 1, it reduces to Eq. 13 with τ Hall = 1/B 0 f B .
In our simulations, we implement the two functions in Eq. 15 separately to avoid numerical instability when α B is close to 0. We fit τ Ohm for the exponential decay case and fit f B and α B for the power-law decay case.Finally, we also introduce a minimum magnetic field strength B min = 10 8 G (Zhang & Kojima 2006) so that it does not decrease to unphysical values.We will show in Section 3 that our result is insensitive to B min , since almost no pulsar B-fields can reach this limit at the end of the simulation.
At each timestep of integration, the value of B is first calculated, then used to determine changes of P and α according to Eqs. 10 and 11.The simulation of each pulsar is terminated at a maximum age randomly chosen within 100 Myr.

Dynamical evolution
After a pulsar is born, we trace its orbit in the Galactic gravitational potential during its lifetime as in previous studies (e.g., Cieślar et al. 2020).The gravitational potential mainly consists of contributions from the Galactic disk, bulge, and dark matter halo.We use the model given by Miyamoto & Nagai (1975) to characterize the gravitational potential of the bulge and the disc respectively, where the bulge mass is M b = 1.12 × 10 10 M ⊙ and disc mass is 277 kpc, a d = 4.2 kpc, and b d = 0.198 kpc.R and r xy are the distance and radial distance to the Galactic Center, respectively.For the halo potential, we follow the model by Paczynski (1990): (18) where the halo mass is M h = 5×10 10 M ⊙ and r c = 6 kpc.Since the mass of halo is divergent, we applied a cutoff radius of r cut =100 kpc, beyond which the halo density drops to zero and Φ halo ∝ r −1 , following Belczynski et al. (2010).

Pulsar Emission
After simulation, we calculate the radio and Gammaray emission of a pulsar using their final properties (e.g., P , B, viewing angle, distance, DM, etc.) to determine if they are detectable by the surveys.

Radio Emission
The observed pulsar radio emission geometry can be modeled by core and cone components (Rankin 1993).This empirical model has been widely used in previous population synthesis studies (e.g., Arzoumanian et al. 2002;Gonthier et al. 2004;Harding et al. 2007;Gonthier et al. 2018).In this work, we follow Gonthier et al. (2018) to describe the beam geometry and radio flux density S ν .
First, we adopt an empirical power law function to parameterize the total radio luminosity where the total radio luminosity L ν , the cone luminosity L cone and the core luminosity L core are measured in mJy kpc 2 MHz, P −3 is the spin period in 10 −3 s, Ṗ−21 is the spin-down rate in 10 −21 s/s, f ν , α ν and β ν are free parameters determined by fitting.The ratio R between the core and cone luminosities depends on P , Ṗ with a broken power law, The simulated pulse profile of radio emission is axisymmetric and depends on the survey frequency ν and polar angle θ as where ρ core and w e are the effective width of core and cone beam, θ is the hollow cone beam annulus size.The index i refers to the cone or core component.The spectral index α i is given by Harding et al. (2007) and D is the distance in kpc.The polar angle θ is related to the magnetic inclination angle α, the viewing angle ξ, and the rotational phase ϕ by cos θ = sin ξ sin α cos ϕ + cos α cos ξ.
The normalization factors are given by The core beam is a Gaussian function that scales with spin period as The cone beam has a hollow Gaussian form with the annulus size and width as respectively, where δ w = 0.18 (Gonthier et al. 2006) and the opening angle depends on the pulsar emission altitude r KG and spin period P (Kijak & Gil 1998) as with r KG estimated using the pulse width method (Kijak & Gil 2003) where ν GHz is the observation frequency in GHz.Eq. 22 shows how radio flux density drops as the polar angle θ increases.Given ξ and α, the observed pulse intensity variation with the rotational phase ϕ can be calculated according to Eq. 23.This gives the radio pulse profile.
We then calculate the phase-averaged radio flux density to determine if the pulse is brighter than the minimum detectable threshold of different radio surveys.
We follow some previous studies not to assume a death line for the radio emission (e.g., Gullón et al. 2014;Graber et al. 2023), such that the intrinsic radio emission of old pulsars does not completely shut down.On the other hand, the assumed luminosity law (Eq.19) will naturally make old pulsars become too faint to be detected.

Gamma-ray luminosity and beaming geometry
One of the applications of our population synthesis study is making predictions for the number of detectable Gamma-ray pulsars.We model the Gammaray luminosity with the fundamental plane relation (Kalapotharakos et al. 2022) where ε cut is the cutoff energy in units of MeV obtained from the empirical ε cut -Ė relations (Kalapotharakos et al. 2017), B P is the field strength at the magnetic pole, i.e. twice the value of B in our simulation, and L γ and Ė are in units of erg s −1 .The observed phaseaveraged Gamma-ray flux is given by where the beaming factor f Ω depends on the viewing angle ξ, α, and Gamma-ray efficiency (w =)η (∝ Ė1/2 ).We adopt the results from the outer gap (OG) and two-pole caustic (TPC) models to estimate f Ω (Watters et al. 2009).For the OG model, where β = ξ −α and ξ I relates to the boundary of hollow cone above the null charge surface We do not attempt to fit any of the parameters here due to the relatively small sample size of Gamma-ray pulsars.Therefore, this has no impact on other parameters of our population model and B-field evolution.

Radio Detection
To determine the visibility of a simulated pulsar for a specific radio survey, we first check if the pulsar's final location lies inside the survey region.If so, we calculate the dispersion measure (DM) from its final location with the latest free electron density model of the Galaxy (Yao et al. 2017).Then we employ the standard radiometer formula (Dewey et al. 1985) to compare the pulsar flux density and the detection threshold S min of the survey where C is the detection signal-to-noise threshold, β is the degradation factor caused by system loss, T rec is the receiver system temperature in K, T sky (l, b) is the sky background temperature in K obtained from the all-sky atlas (Haslam et al. 1982), G is the effective gain of the telescope antenna in units of K Jy −1 , N p is the number of polarization, t obs is the observation time in s, ∆F is the total bandwidth in MHz, P is the pulsar period in s, and W e is the effective pulse width in ms.We estimate W e based on the formulas given by Bhat et al. (2004) and we follow previous studies (e.g., Lorimer et al. 1993;Bates et al. 2014) to estimate G from the original telescope gain G 0 listed in Table 6 by This accounts for the decrease in sensitivity when a pulsar has an offset r from the telescope pointing center.
Here w is the full width at half maximum (FWHM) of the telescope beam and r is randomly chosen assuming uniform distribution of r 2 between 0 and w 2 /4 (Bates et al. 2014).
If a simulated pulsar lies within the survey area with radio flux density above the survey sensitivity, then it is considered as detected and its parameters will be compared with the observed pulsar sample using the procedure described in Section 2.6 below.Regardless of whether a simulated pulsar is detected or not, it is always counted for the birth rate calculation.

Pulsar sample
We select our pulsar sample from major surveys listed in the ATNF catalog 1 (Manchester et al. 2005) with the following criteria: 1. P > 0.03 s and Ṗ > 0 to rule out millisecond pulsars.
2. Exclude rotating radio transients (RRATs) and binary pulsars, as these could have different formation channels and evolution paths than normal pulsars.
3. Exclude radio surveys with fewer than 25 normal pulsars to speed up the computing time.These are listed in Table 5 in Appendix A. This only leads to 38 fewer pulsars, not significant.
4. Exclude surveys that are fully overlapped with others (see Table 5 in Appendix A).
6. Exclude 87 pulsars that have DM-estimated distance larger than 25 kpc, which are likely to be extra-galactic.
The above selection criteria give a total of 1859 radio pulsars from 14 surveys.covering the vast majority of normal radio pulsars.Table 1 shows the number of pulsars from each survey and the survey parameters are presented in Table 6.
Our simulations assume that the radio flux density is measured at 1400 MHz (S 1400 ).However, there are 235 pulsars in our sample that only have measurements in other bands.We scale their flux densities to 1400 MHz assuming a typical spectral index of −1.60 (Jankowski et al. 2018).We show in Appendix A that this does not change the flux density distribution.

Fitting of the model parameters
To evaluate the goodness of the model fit given a set of parameters, we compare the distribution of simulated and observed pulsar properties using multivariate KS statistic (Fasano & Franceschini 1987;Justel et al. 1997).The procedure is outlined in Appendix B. We follow previous works (e.g., Gonthier et al. 2018) to focus on four observables: P , Ṗ , DM, and S 1400 as other observables are not sensitive to the magneto-rotational evolution model that we are interested in.The fitting parameters are µ P , σ P , µ B , σ B , α B , f B , f ν , α ν , and β ν for the power-law decay model.For the Ohmic decay model, f B is replaced with τ Ohm and there is one fewer free parameter as α B = 0. We employed the MCMC technique to search for the optimal model parameters.Given a KS statistic D ks , the p-value can be calculated as where n 1 and n 2 are the sample size of simulation and observation (Fasano & Franceschini 1987).This is then is used to construct the likelihood function in MCMC.
There are no truncated values for the free parameters.We assume log-uniform priors for τ Ohm , f ν and uniform priors for all other fitting parameters, and employed the emcee package (Foreman-Mackey et al. 2013) to run the MCMC chains.

Best-fit model parameters
The best-fit power-law model suggests a very small value of α B close to 0, which reduces the model to exponential decay.As we will show in Sect.3.4 below, the model with large α B produce pulsars narrowly distributed in P -Ṗ space, which does not fit the observed pulsar population.
Figure 1 shows the corner plot for the exponential model parameters, obtained from MCMC with 4 independent chains of 16 walkers and 2500 steps each, a total of 160,000 simulations.The optimal parameters determined from the most probable values (MPVs), i.e. at the maximum marginal probability density, with 68% confidence intervals, are listed in Table 2.As indicated in the figure, the model parameters µ B , σ B , τ Ohm and β ν are well constrained, but f ν , the lower bound of µ P and upper bound of σ P are not very well determined.

Best-fit pulsar population
Table 1 shows a comparison between the observed and simulated number of pulsars in different surveys.The latter is obtained from the average of 100 simulations using the best-fit model parameter.These numbers are generally consistent and are very well matched for the palfa and pkshl surveys.However, some surveys, including ar4, gb3, and pksmb, show a significant discrepancy of up to 50%.The possible cause of this will be discussed in Section 4.4.In Figure 2 we compare the distributions of observables for the simulation and pulsar sample.We plot the four observables in the first row of Figure 2 (P , Ṗ , DM, and S 1400 ) that are used to determine the goodness-of-fit, and also the pulsar Galactic coordinates, the B-field, and the characteristic age τ c ≡ P/2 Ṗ .The results show that our model can generally reproduce the observed pulsar population.The only exception may be the Galactic longitude, which has obvious excess at l ∼ 30 • that our model fails to capture.This problem was first noted by Faucher-Giguere & Kaspi (2006) and is due to the spiral arm structure that cannot be accurately reproduced by the simple spatial model in Eq. 1.
We also show in the figure a comparison between the true B-field in the simulation and the B-field estimated from spin down using Eq. 9 with sin α = 1, i.e.
which is commonly adopted in the literature.They have very similar distributions and the true B-field is slightly smaller.This is expected as a direct comparison between Eqs. 10 and 37 gives B/B SD = 3/(2 + 2.4 sin 2 α), which varies between 0.55 and 0.82.In other words, Eq. 37 is indeed a good approximation to the realistic force-free spin-down model.Finally, we show the age distribution of the simulated pulsars and its comparison with τ c .The mean true age is 4.8 Myr and over 40% of them are older than τ Ohm .τ c shows a good agreement with the true age for young pulsars up to ∼ 1 Myr, and then becomes much larger beyond that.This can be attributed to the decay of B-field, which makes Ṗ to drop rapidly (see discussion in Sect.3.4 below), thus, resulting in large τ c .
Figure 3 shows the P -Ṗ diagram of the observed pulsar sample and a random realization of the best-fit model with the same number of pulsars.They generally show a good match.We notice that our model fails to reproduce pulsars in the lower left corner of the diagram (i.e. with fast spin but small Ṗ ).These could be mildly recycled pulsars and their formation is not modeled in our simulation.

Dynamical evolution
We plot in Figure 4 the location of the pulsar sample in Galactic coordinates and a Monte Carlo realization of  our model for comparison.They show a good agreement in general.Results of pulsar dynamical evolution are plotted in Figure 5.These include distributions of pulsar displacements and velocities.The pulsar displacement plot suggests that most pulsars do not travel far over their lifetime.Over 78% of them are found within 3 kpc of their birthsite at the end of the simulation.Similarly, the velocity distribution plot shows that the final velocities (v final ) do not change much from the initial values (v 0 ).They follow the same double Mawellian distribution with identical peak locations, but slightly fewer (∼1.3% drop) pulsars in the low-velocity component (< 130 km s −1 ) and more (∼1.1% increase) with 200−300 km s −1 .All these results suggest that acceleration due to Galactic gravitational potential is not significant for most of pulsars.

Magneto-rotational evolution
To investigate the effect of magnetic field decay on pulsar magneto-rotational evolution, we plot in Figure 6 the pulsar evolutionary tracks in the P -Ṗ diagram and the time evolution of B-field strength and inclination angle.These are calculated using the best-fit exponential de-  6.The red open circles and triangles show radio-loud and radio-quiet Gamma-ray pulsars, respectively, from the Third Fermi-LAT catalog (Smith et al. 2023).
cay model with τ Ohm = 8.3 Myr.We set P 0 = 0.03 s and initial α = 60 • , and try different initial B-fields from B 0 = 10 11.5 G to 10 13.5 G.For illustration purposes, we also plot the power-law decay model with α B = 1.0 and log f B = −26.8.
As shown in the top panel of the figure, the tracks in both scenarios show similar behavior at the beginning.They all follow straight lines of constant B-field, since the field decay is not significant at early times.For the exponential decay model, the tracks turn vertically downward when the pulsars are older than τ Ohm , which is due to B-field decay, such that Ṗ becomes very small ( Ṗ ∝ B 2 ; see Eq. 10) and hence P stays nearly constant.
The power-law decay model shows a similar overall trend but the detailed evolutionary tracks after the turning point depend on B 0 .The large B 0 tracks start to turn downward earlier.This is because the B-field decay timescale depends on 1/B 0 according to Eq. 15, i.e. larger initial B-field decays faster.For this reason, all pulsars end up with similar B-fields around 10 11 G at 100 Myr, much higher than that in the exponential model.This is clearly seen in the B-field evolution plot in Figure 6.As a consequence, pulsars show a narrow range of P and Ṗ after evolution, which is not preferable by the observations.We also note that pulsars evolve much slower in this model.
In addition to the B-field strength, the inclination angle α between the magnetic and spin axes is another key factor that determines Ṗ .Its evolution is plotted in Fig. 6.It is clear that the evolution depends critically on B 0 .α drops rapidly for large B 0 case but it takes a long time to evolve if B 0 is low.This is expected as α ∝ −B 2 (see Eq. 11).At later times, α stops evolving when the B-field decays to very small values.In the exponential decay scenario, this happens when the evolution time exceeds τ Ohm .The final value of α therefore remains relatively large for pulsars started with low B 0 .This effect, on the other hand, is less prominent for the power-law decay model.It can give significantly smaller α for old pulsars, since the B-field stays comparatively larger till the end of the evolution.In Figure 7 we plot the change of α for the simulated pulsars that are detected by the surveys.The pulsars are assumed to have isotropic α distribution at birth, but those with larger initial α (i.e.orthogonal rotators) have a higher chance to be detected.This is mainly due to the geometrical effect as the emission beams of these pulsars are more likely to sweep across the Earth.Moreover, the plot also clearly illustrates the effect of α decay, which makes most of the detected pulsars to become nearly aligned rotators (∼18% with α < 10 • ) at the end of the evolution.
In Figure 8 we plot the observed α versus characteristic age for 80 pulsars reported by Malov & Nikitina (2011) and compare with our best-fit model.The observations are generally consistent our model prediction, although not many pulsars are found in the most likely parameter space suggested by the model.We attribute the discrepancy to a small sample size and also the different beam geometry we adopt (core and cone beam) than the one (cone beam) assumed by Malov & Nikitina (2011) to deduce α.

Gamma-ray pulsar
After a pulsar is simulated, we calculated its phaseaveraged Gamma-ray flux according to the equations in Section 2.4.2.The result is then compared with the 14-year Fermi all-sky sensitivity map2 (Smith et al. 2023) to determine if it is visible.Table 3 lists the number of detectable radio-quiet and radio-loud Gamma-ray pulsars based on the OG and TPC models, averaged from 100 simulations.We find that the OG performs better than the TPC.It can produce the number of radio-loud Gamma-ray pulsars as observed in the Third Fermi-LAT  Catalog (Smith et al. 2023), while the TPC model predicts too many.Figure 3 shows a random simulation of Gamma-ray pulsars in the P -Ṗ diagram using the preferred OG model.They occupy the same parameter space as the observed pulsars, both are in the upper left of the diagram where Ė is high.

DISCUSSION
Table 3. Number of Gamma-ray pulsars in the third Fermi-LAT pulsar catalog (Smith et al. 2023) compared with our simulation results assuming the outer gap (OG) and the twopole caustic (TPC) models (Watters et al. 2009

Magnetic field decay model
In this paper, we present a detailed population study of normal pulsars using realistic physics inputs, including an updated spin-down formula, evolution of magnetic field strength and inclination angle, and sophisticated radio beam geometry.We fit the model to a large sample of pulsars using a rigorous statistical test and determine the optimal model parameters using MCMC.We parameterize the B-field evolution with the general form dB/dt ∝ B (1+α B ) and the best-fit result indicates that α B = 0 is preferred.This implies that the exponential decay model can better reproduce the observed pulsar population than the power law.We are able to constrain the decay timescale to be 8.3 Myr (1σ confidence interval of 5.2−12.3Myr).This is comparable to the mean age of 4.8 Myr among the detected pulsars in the simulation, implying that B-field decay is nonnegligible in most pulsars.
It was suggested that Ohmic dissipation in the crust due to phonon scattering is the primary cause of the magnetic field decay in normal pulsars (Igoshev & Popov 2015) and the Hall effect is only important for B-fields stronger than 10 14 G (Viganò et al. 2013), which is much higher than that in our simulated pulsars (mean log B 0 = 12.3 G).In this picture, the B-field would show a rapid decay at first driven by the Hall cascade.The Ohmic decay then takes over at later times.This could help explain why the exponential model is preferred as our result shows.The Ohmic decay timescale τ Ohm depends sensitively on a few physical parameters, including the temperature, impurity factor, and electron fraction (Igoshev et al. 2021).Using T = 3.0 × 10 8 K and the typical values for other parameters, we can obtain τ Ohm =8.3 Myr as our best-fit value.

Comparison with previous studies
Magnetic field decay has been considered in a number of pulsar population synthesis studies.These works adopt different models with various decay timescales ranging from 1 kyr to 10 Myr (e.g., Dirson et al. 2022;Xu et al. 2023;Graber et al. 2023).Only one attempt to constrain the timescale directly from fitting, and a timescale τ = 4.3 ± 0.4 Myr was suggested for the exponential decay model (Cieślar et al. 2020).While this is compatible with our result, we caution that they employ a simple spin-down formula without α evolution and their pulsar sample is much smaller than ours.More importantly, they smoothed the observed pulsar distribution for model comparison, which may not accurately capture the statistic (see Graber et al. 2023).
The different decay timescales in previous studies could have an impact on the initial B-field required.Our best-fit model prefers slightly lower B 0 with µ B = 12.3.This could be attributed to the strong correlation between µ B and τ Ohm as Fig. 1 shows.Our relatively long time scale suggests a late onset of the B-field decay, such that the initial field could be lower.The situation could be more complicated for the power-law model, since the decay timescale would depend on B 0 (see Dirson et al. 2022;Xu et al. 2023) instead of fixed as assumed in some works (e.g., Graber et al. 2023).
Overall, all the population models for normal pulsars do not require strong B 0 (µ B > 13.0) even with the consideration of field decay.Our best-fit model produces only 0.2% of pulsars with B 0 > 10 13 G and none above 10 14 G.This is not sufficient to explain magnetars, which have field strengths of 10 14 − 10 15 G (Kaspi & Beloborodov 2017).This suggests that this class of objects could be a distinct population born with another channel.
The evolution of the magnetic inclination angle was studied in a recent work (Dirson et al. 2022), but a slightly different model was adopted for which the decay rate is independent of the B-field strength.Similar to our result, it was also found that the effect of α decay is significant, making the majority of pulsars nearly aligned rotators at the end of the evolution.However, their detected pulsars consist of mainly orthogonal rotators, which is contrary to our result.The discrepancy could be due to different radiation beam models.In particular, our core and cone emission components with Gaussian profiles allow detection even with very small alignment angles.
For the birth spin period distribution, our fit requires similar µ P (µ P = −1.40)as in other studies even using the updated sample (see Fig. 3).However, we note that µ P is not very well constrained in the fitting and the large σ P value suggests a broad distribution of P 0 .Such a result is consistent with the finding that the observed population is insensitive to the P 0 distribution and any values in the range of P 0 < 0.5 s could fit well (Gonthier et al. 2004;Gullón et al. 2014).
We parameterize the radio luminosity with the general form L ν ∝ P αν Ṗ βν as in previous studies.The bestfit result we found (α ν , β ν ) = (−1.76,0.67) is close to the typically reported value of (−1.5, 0.5) (e.g., Faucher-Giguere & Kaspi 2006;Bates et al. 2014).Our slightly steeper dependence on P and Ṗ could be attributed to the B-field decay, which results in more aged pulsars in the lower part of the P -Ṗ diagram than in previous works.A larger β ν is therefore needed to make them non-detectable, in order to match the observations.This also leads to a more negative value of α ν due to the strong correlation between α ν and β ν (see Fig. 1).We note that a very shallow dependence (−0.39, 0.14) is recently reported based on the MeerKAT Survey (Posselt et al. 2023).This is however not seen in other surveys (e.g., Anumarlapudi et al. 2023) and more observations are needed to confirm this.
Alternatively, some studies investigate the correlations between L ν and Ė (e.g., Szary et al. 2014;Wu et al. 2020).We did not attempt to express L ν in terms of Ė in our modeling, since this requires a ratio R αβ ≡ −α ν /β ν = 3, not supported by our results.The values of R αβ reported in the literature always lie between 2 and 3 (e.g., Faucher-Giguere & Kaspi 2006;Ridley & Lorimer 2010;Bates et al. 2014;Johnston & Karastergiou 2017;Posselt et al. 2023), and the value R αβ = 2.6 we obtain is consistent with this range.On the P -Ṗ diagram, this ratio corresponds to the slope of constant L ν lines, which can be regarded as observationlimit lines (see Wu et al. 2020).It is interesting that the pulsar emission death lines expected from theories also share similar slopes (see Abolmasov et al. 2024, and references therein).This could imply a deeper physical connection underneath.We however stress that it is highly non-trivial to infer the pulsar intrinsic luminosity from observations, as this requires detailed knowledge of the radio beam geometry, the magnetic inclination angle and the viewing angle.All of these are difficult to obtain and any uncertainties will blur the correlation between L ν , P , and Ṗ .
Our study shows that it is possible to reproduce the pulsar P -Ṗ distribution without assuming a death line, same conclusion as Gullón et al. (2014) and Graber et al. (2023).Our luminosity law makes aged pulsars fade away and naturally become undetected.This results in an observation-limit line, which has the same effect as a death line from an observational point of view.On the other hand, we are unable to rule out the death line either.We tried fitting with the classical death line (Bhattacharya et al. 1992) but there is not much change in the model parameters.This is because in our model most pulsars are already too faint to detect before they Pulsar population with magneto-rotational evolution reach the death line.We note that the details of the death line could depend on many physical parameters, including the B-field strength, the magnetic inclination α (Beskin & Litvinov 2022), and even the equation of state (Zhou et al. 2017).Testing these ideas is beyond the scope of this work.

Pulsar birth rate
Our best-fit population model gives a normal pulsar birth rate of ∼0.68 per century in the Galaxy.This value is not affected by the age cutoff of 100 Myr assumed in the simulation, since pulsars that old cannot be detected.The birth rate we found is lower than those in previous pulsar population studies, e.g., 2.8 per century from Faucher-Giguere & Kaspi (2006), 1.4 per century from Lorimer et al. (2006), and 1.6−2 per century from Graber et al. (2023), but still much higher than 0.16 per century from Cieślar et al. (2020).Our result indeed better aligns with the recent Galactic core-collapse supernova rate estimate of 1.63±0.46(Rozwadowska et al. 2021).However, we stress that the birth rates inferred from population modeling depend critically on the radio luminosity law, for which many of the details (such as the beam geometry) is not fully understood and could be subject to large uncertainties as discussed above.Additionally, including different surveys in the observation sample could give different birth rates (see Graber et al. 2023).This could be caused by unreliable survey parameters (see Section 4.4 below).Cieślar et al. (2020) and Chakraborty & Bagchi (2020), the MeerKAT TRAPUM survey (Ridolfi et al. 2021), and the Canadian Hydrogen Intensity Mapping Experiment Pulsar Survey (CHIME; Good et al. 2021).The expected number of detections is listed in Table 4.We note that these could be very different from the actual numbers, possibly due to inaccurate survey parameters.For instance, the Giant Metrewave Radio Telescope High-Resolution Southern Sky (ghrss) and the Green Bank Northern Celestial Cap (gbncc) surveys were expected to discover 80 new normal pulsars and 130 new MSPs, respectively (Bhattacharyya et al. 2016;Stovall et al. 2014), but at the end only 22 and 57, respectively, were found.

CONCLUSION
We carry out a population synthesis study for normal radio pulsars employing the most updated model inputs, including log-normal distributions for P 0 and B 0 (Igoshev et al. 2022), the spin-down formula from force-free magnetohydrodynamic simulations with exact dependence on the magnetic inclination angle (Philippov et al. 2014), double Maxwellian distribution for the kick velocity (Verbunt et al. 2017;Igoshev 2020), free electron density model by Yao et al. (2017), core and cone radio emission model (Gonthier et al. 2018), and latest empirical relation for Gamma-ray luminosity (Kalapotharakos et al. 2017(Kalapotharakos et al. , 2019(Kalapotharakos et al. , 2022)).We compared our simulation results with a large pulsar sample that covers all major surveys and performed a fitting utilizing the MCMC approach with 4D KS statistic to constrain the model parameters.Our main findings are: 1.The observed pulsar population can be well reproduced with the exponential B-field decay model.This model also performs substantially better than the power-law decay in terms of goodness of fit.
The power-law model can only provide a good fit when α B = 0, which reduces into the exponential form.The results suggest that the Ohmic dissipation could play an important part in the B-field decay process in neutron stars.2. Our fitting is able to constrain µ B to be 12.3 G and the magnetic field decay timescale to be 8.3 +3.9 −3.0 Myr for the exponential model.This implies that the decay is non-negligible for most pulsars.3. Due to the evolution of the magnetic inclination angle, we expect a large fraction of pulsars to have a good alignment between the magnetic and spin axes.
4. Our best-fit radio luminosity law L ∝ P −1.76 Ṗ 0.67 can reproduce the pulsar P -Ṗ distribution without the need for a death line.However, we note that the exact values of the indices are modeldependent and the largest uncertainty is from the emission beam geometry.

5.
Employing the OG emission model, our population model can well reproduce the number of gammaray pulsars observed in the Galaxy.
Our study shows that pulsar population synthesis can provide a powerful tool to probe the physics of neutron stars.In future works, more microphysics could be included, such as the evolution of B-field geometry and its dependence on temperature and even the equation of state.Also, comparing the simulations with a large observational sample from upcoming surveys will provide better constraints on the model parameters.
We thank Wynn C. G. Ho for the discussion and useful suggestions.We thank the referee for the constructive advice.We thank Di Li for providing us with the FAST GPPS survey information.Z.S. and C.-Y.N. are supported by a GRF grant of the Hong Kong Government under HKU 17303221.Among the 1859 pulsars in our sample, only 1624 have flux density measurements at 1400 MHz in the ATNF catalog.In order to keep those 235 pulsars, we estimate S 1400 from the measurement at the closest frequency available assuming a power law spectrum with a spectral index of −1.60.Figure 9 shows a comparison of the flux density distributions with and without these 235 pulsars.We apply the KS test and obtain a statistic of 0.0256 with p-value of 0.599, suggesting that using this S 1400 estimate method has no significant effect on the distribution and hence should not affect our analysis results.
Table 5 lists the radio surveys that are not used in this study, either fully overlapped by other surveys (the left column), or detected fewer than 25 normal pulsars according to our selection criteria mentioned in Section 2.6.1.Including all these surveys would only add 38 new pulsars to our sample.After each simulation run, we need to compare the distributions of M simulated pulsars (M = 10N in this study) and N observed pulsars (N = 1859 in our sample) for k observables (A, B, C, D, . ..) (k = 4 in our case).We evaluate the goodness of fit using the multivariate KS statistic (Fasano & Franceschini 1987;Justel et al. 1997), which is calculated as the following: The value of D ks indicates how well the simulated pulsar population matches that of the observed sample, in terms of the observables i.e.P , Ṗ , DM, and S 1400 in this work.The smaller the value the better the match.
Since our simulation involves random sampling processes, this leads to fluctuation of the KS statistic, i.e. slightly different D ks values are obtained every time even with identical model parameters.To investigate this effect, we run simulations with the exponential decay model and the best-fit parameters, and generate 1, 5, 10, 20, 50 times of simulated pulsar as the observation sample.The whole process is repeated for 100 times to determine the fluctuation in D ks .The result is shown in Figure A10 and the 1-σ (i.e.68%) range and standard deviations are listed in Table A7.We find that simulating ten times of pulsars as the observation sample (i.e.18590 simulations for 1859 pulsars observed) gives sufficiently small fluctuation of D ks and more simulations provide no significant improvement.Therefore this is adopted in our simulation.Table 7. 1σ range and standard deviation (SD) of KS statistic obtained from 100 simulations, each with 1, 3, 5, 10, 20, 50 times of pulsars as the observation sample.All simulations are done with the exponential model using the best-fit parameters.

Figure 1 .
Figure 1.The MCMC marginal parameter space for Ohmic dissipation assumption.The 2D marginal plot is smoothed using a Gaussian kernel.

Figure 2 .
Figure 2. Distributions of observables for the simulated (blue) and observed (red) pulsar populations.The simulation is based on the exponential B-field decay model with the best-fit parameters.The initial spin period (P0), initial spin-down rate ( Ṗ0) and birth magnetic field strength (B0) distributions (yellow histograms) are also shown for comparison.The magnetic field distributions shown by the red and blue lines are BSD from Eq. 37, and the true B-field from the simulations is plotted in grey color.The characteristic age τc is calculated with P/2 Ṗ , and the grey line shows the true age from simulations.The vertical dashed line in the plot indicates the best-fit Ohmic decay timescale of 8.3 Myr.

Figure 3 .
Figure 3. P -Ṗ diagram for the observed (red) and simulated (blue) pulsars.The simulation is based on the exponential B-field decay model with the best-fit parameters.The red dots show the pulsar sample from radio surveys listed in Table6.The red open circles and triangles show radio-loud and radio-quiet Gamma-ray pulsars, respectively, from the Third Fermi-LAT catalog(Smith et al. 2023).

Figure 4 .
Figure 4. Galactic location distributions for the simulated (blue) and observed (red) pulsar populations.The simulation is based on the exponential B-field decay model with the best-fit parameters.

Figure 5 .
Figure 5. Dynamical evolution of simulated pulsars that are detected by the surveys.Left: distribution of pulsar displacement between their birth sites and final locations.Right: distribution of initial and final pulsar velocities at the end of the simulation.

Figure 6 .
Figure 6.Pulsar evolution tracks in the P -Ṗ diagram and the evolution of the magnetic field and the magnetic inclination angle (α) for the best-fit exponential B-field decay model.All pulsars start with P0 = 0.03 s and initial α = 60 • .Different colors lines represent different B0 from 10 11.5 G to 10 13.5 G, uniform in logarithmic space.

Figure 7 .
Figure 7. Distributions of the initial and final values (yellow and red lines, respectively) of the magnetic inclination angle α for the simulated pulsars detected with the surveys.

Figure 8 .
Figure 8. Magnetic inclination angle α versus characteristic age for observed and simulated pulsars.The black dots show 80 measurements of α reported by Malov & Nikitina (2011).The underlying density plot is obtained by simulations of 50 times the pulsar sample (i.e.92950 pulsars) using our bestfit population model.

Figure 9 .
Figure 9. Distribution of radio flux density (S1400) for 1624 pulsars with measurements from the ATNF pulsar catalog (red line) compared with the entire sample of 1859 pulsars (blue line).For pulsars that have measurements at other frequencies only, we scale their flux densities to 1400 MHz assuming a spectral index of −1.60.
1.For each simulated or observed pulsar, call it the i-th pulsar (i = 1, . . ., N from observation or i = 1, . . ., M from simulation), which has the properties: (A i , B i , C i , D i , . ..).We first calculate the cumulative fraction of the observation sample compared with this i-th pulsar (O i ): if n observation pulsars that satisfy the criterion(A < A i , B < B i , C < C i , D < D i , . ..), then O i,1 ≡ n N (B1)There are totally 2 k different values ofO i : O i,1 (A < A i , B < B i , C < C i , D < D i , . ..), O i,2 (A < A i , B > B i , C < C i , D < D i , . ..), . .., O i,2 k (A > A i , B > B i , C > C i , D > D i , .. .)accounting for all the different definitions of O i , and 0 ≤ O i ≤ 1.2.Repeat similar step calculating the cumulative fraction of the simulation sample compared with the i-th pulsar: if m simulation pulsars that satisfy the criterion (A< A i , B < B i , C < C i , D < D i , . ..), then S i,1 ≡ m M (B2)Again there are totally 2 k different values of S i :S i,1 (A < A i , B < B i , C < C i , D < D i , . . .), S i,2 (A < A i , B > B i , C < C i , D < D i , . ..), . . ., S i,2 k (A > A i , B > B i , C > C i , D > D i , .. .), and 0 ≤ S i ≤ 1.3.For this i-th pulsar, calculate the maximum distance between 2 k O i and corresponding 2 k S i asD i ≡ max (j=1,2,...,2 k ) |O i,j − S i,j |.If this i-th pulsar is from observation, label the distance as D obs,i .If this i-th pulsar is from simulation, label the distance as D sim,i .4. Repeat (1)-(3) above to obtain the distance for each observed pulsar D obs,i (i = 1, 2, . . ., N ) and each simulated pulsar D sim,i (i = 1, 2, . . ., M ).Then calculate the maximum distance among the observed sample and the simulated sample respectively: D obs,max ≡ max(D obs,1 , D obs,2 , . . ., D obs,N ) and D sim,max ≡ max(D sim,1 , D sim,2 , . . ., D sim,M ).5.The KS statistic value is then D ks = (D obs,max + D sim,max )/2, with 0 ≤ D ks ≤ 1.

Figure 10 .
Figure 10.Fluctuation of KS statistic values for 100 simulations, all performed using the best-fit exponential decay model but with a different number of simulated pulsars, from 1 to 50 times as the observation sample.

Table 1 .
Number of normal pulsars detected from radio surveys compared with our best-fit model simulation results.See Table6in Appendix A for the full name of the surveys.

Table 2 .
Best-fit parameters with 1σ confidence intervals for the exponential and power-law B-field decay models. ).

Table 4 .
Prediction about the detectable number of future radio surveys.N det and N dis are the number of detectable and newly discovered Galactic normal pulsars by the surveys.

Table 5 .
Radio surveys not used in this work.