The Galactic population of canonical pulsars

Context. Current wisdom suggests that the observed population of neutron stars are manifestations of their birth scenarios and their thermal and magnetic ﬁeld evolution. Neutron stars can be observed at various wavebands as pulsars, and radio pulsars represent by far the largest population of neutron stars. Aims. In this paper, we aim to constrain the observed population of the canonical neutron star period, its magnetic ﬁeld, and its spatial distribution at birth in order to understand the radio and high-energy emission processes in a pulsar magnetosphere. For this purpose we design a population synthesis method, self-consistently taking into account the secular evolution of a force-free magnetosphere and the magnetic ﬁeld decay. Methods. We generated a population of pulsars and evolved them from their birth to the present time, using the force-free approximation. We assumed a given initial distribution for the spin period, surface magnetic ﬁeld, and spatial Galactic location. Radio emission properties were accounted for by the polar cap geometry, whereas the gamma-ray emission was assumed to be produced within the striped wind model. Results. We ﬁnd that a decaying magnetic ﬁeld gives better agreement with observations compared to a constant magnetic ﬁeld model. Starting from an initial mean magnetic ﬁeld strength of B = 2 . 5 × 10 8 T with a characteristic decay timescale of 4 . 6 × 10 5 yr, a neutron star birth rate of 1 / 70yr and a mean initial spin period of 60ms, we ﬁnd that the force-free model satisfactorily reproduces the distribution of pulsars in the P − ˙ P diagram with simulated populations of radio-loud, radio-only, and radio quiet gamma-ray pulsars similar to the observed populations.


Introduction
Pulsars are rapidly rotating, highly magnetized neutron stars surrounded by a plasma-filled magnetosphere that emits regular pulses of radiation at their spin frequency. The term 'pulsar' comes from 'pulsating radio source' since they were first observed at radio wavelengths. Currently, the number of observed radio pulsars is roughly 3,000 1 (Manchester et al. 2005). Subsequently, however, it was found that pulsars were also bright in X-rays, optical, and gamma-rays. In particular, the Large Area Telescope (LAT) on board the Fermi satellite has forcefully changed our understanding of gamma-ray pulsars by discovering dozens of radio-quiet gamma-ray pulsars as well as millisecond pulsars (MSPs). After more than ten years of operation, the LAT has detected over 250 pulsars (Abdo et al. (2013); Smith et al. (2019)), establishing them as the dominant class of GeV sources in the Milky Way. This signified a new era for gamma-ray pulsar statistics, first started by EGRET with the detection of only seven pulsars (Thompson et al. 1995). These recent discoveries have raised fundamental questions, one of the most challenging amongst them being to find a physical consistent theory of pulsar multi-wavelength emission mechanisms. With the advancement of instrumental techniques, future surveys with greatly improved sensitivities are being envisaged. Amongst them, the Square Kilometre Array (SKA) is expected to be 50 times more sensitive, and is predicted to survey the sky 10,000 times faster than any existing imaging radio telescope array 2 . In this context, it is of utmost importance to be able to anticipate the discovery rate of new pulsars, since it enables us to validate our current understanding of pulsar emission mechanisms, as well as their evolution track and the factors impacting their detection.
Pulsar population synthesis (hereafter PPS) studies are powerful tools aiming to predict the detectability of pulsars. Inferring the underlying properties responsible for the observed pulsar population is an arduous task. Indeed, both the impact of the detection procedure and the inhomogeneous properties of the interstellar medium (ISM) have to be accurately estimated to account for possible observational biases.
Much of the progress in PPS has come from thorough Monte Carlo simulations that generate pulsars and test whether they fulfil the criteria for detection according to geometrical factors and sensitivity issues. It is then possible to develop and optimize a model for the underlying pulsar population, which informs us about the important intrinsic neutron star parameters and distribution, enabling predictions for future surveys. Usually, PPS studies follow two simple approaches. The first one is to take a 'snapshot' of the Galaxy as it appears today, where no assumptions are made regarding the prior evolution of the pulsar population. Instead, this population is generated assuming various distribution functions (typically spatial distribution, spin period P and,Ė-orṖ), which are optimized to match the observations. Inspired by earlier studies from Taylor & Manchester (1977) and Lyne et al. (1985), Lorimer et al. (2006) applied the snapshot approach to the canonical 3 pulsar population to determine best-fitting probability density functions in Galactocentric radius (R), luminosity (L), height with respect to the Galactic plane (z), and the period P for the currently observed population of pulsars. Alternatively, one may consider 'evolution' strategies where the pulsars are evolved from birth up to the present era, starting from an initial spatial distribution, and an initial period and magnetic field distribution. A fine example of the latter genre is the comprehensive study of Faucher-Giguère & Kaspi (2006), which quite successfully reproduced the properties of the main part of the radio pulsar population using a model in which the luminosity has a power-law dependence on P andṖ. PPS studies can be also extended to the population of neutron stars observed in other bands, such as X-rays (see for instance Popov et al. (2010)) and γ−rays. Indeed, with the broad increase in γ-ray pulsar numbers, a statistical treatment of the γ-ray population in combination with deep radio surveys of the Galactic plane is now feasible. Early works on radio-loud gamma-ray pulsar populations carried out before the Fermi era include Gonthier et al. (2002Gonthier et al. ( ), (2004Gonthier et al. ( ), (2007aGonthier et al. ( ), (2007b. With the advent of Fermi/LAT, new studies emerged (see for instance Gonthier et al. (2018); Ravi et al. (2010); Takata et al. (2011);Watters & Romani (2011);Pierbattista et al. (2012)), trying to constrain the geometry and the location of the gamma-ray emission sites. Watters & Romani (2011) showed that an initial spin period of P 0 = 50 ms and a birth rate of one neutron star per 59 yr were required to reproduce the observed γ-ray population. They made the prediction that after ten years of operations, Fermi should detect ∼120 young γ-ray pulsars, of which about one half would be radio-quiet. Gonthier et al. (2004) included an exponentially decaying magnetic field with a 2.8 Myr timescale and displaying a satisfactory agreement with the P −Ṗ distribution at that time. Later, more accurate magnetic field decay models were elaborated, especially for magnetars, as presented in Viganò et al. (2013).
Gamma-ray emission modelling also drastically benefited from the second Fermi pulsar catalogue. Based on fluid simulations in a dissipative regime by extension of the force-free regime, Kalapotharakos et al. (2014) constrained the magnetic axis inclination with respect to the rotation axis and computed the associated light curves for curvature radiation. Later, Kalapotharakos et al. (2017) also constrained the dissipation rate depending on the pulsar spin down, differentiating millisecond pulsars from young pulsars. Eventually, Kalapotharakos et al. (2018) gave a simple analytical fit for the gamma-ray luminosity depending on the cut-off energy, spin down, and surface magnetic field. We used these results to predict the gamma-ray flux in our population synthesis.
The appropriate emission mechanisms of γ-rays emission from pulsars are still under active investigation. Various models have been proposed to explain the particle acceleration sites as the origin of pulsed γ-rays emission. The first model suggested is the polar cap region, which is confined in the open magnetosphere at low altitudes (Sturrock (1971); Ruderman & Sutherland (1975); Daugherty & Harding (1982); Daugherty & Harding (1996)). The slot-gap (along the last open magnetic field 3 Defined here as the pulsars that are non-binaries, with P > 20 ms and that are not magnetars. 28  197  101  35 61   total  2665 2553 63 84   Table 1: Number of known pulsars withĖ above and below 10 28 W, and above 10 31 W. 4 4 The quantities N r , N g , and N rg are the number of radio only, gamma only, and radio-loud gamma-ray pulsars, respectively. It should be noted that we excluded the binary pulsars as well as pulsars with P < 20 ms. The data have been taken from the ATNF catalogue and from https://confluence.slac.stanford.edu/display/ GLAMCOG/Public+List+of+LAT-Detected+Gamma-Ray+Pulsars lines, Arons (1983); Muslimov & Harding (2004); Harding et al. (2008); Harding & Muslimov (2011)) and the outer-gap model (extending to the edge of the light cylinder, Cheng et al. (1986Cheng et al. ( ), 1986bHirotani (2008), Takata et al. (2011)) are both located at high altitudes in the outer magnetosphere. Polar cap models predict a sharp, super-exponential cut-off at several GeV, which is steeper than those of the γ -ray spectra of the Crab and Vela pulsars measured by Fermi (Abdo et al. (2013)). In the first few years of the 2000s, a new picture emerged for the origin site of gamma-ray pulsars. It is located well outside the light cylinder in the striped wind (Kirk et al. 2002) based on the structure discussed by Coroniti (1990) and by Michel (1994). Relativistic beaming and the spiral structure of the emitting current sheet produces a pulsed radiation (Pétri (2009), Pétri (2011. No PPS studies have been carried out so far by considering the latter emission site. Furthermore, although pulsar magnetospheres are filled with a relativistic electron-positron plasma (Goldreich & Julian (1969)), most of the work on pulsar population studies are carried out assuming a spherical neutron star rotating in a vacuum. However, the radiation from neutron stars is produced by charged particles flowing within their magnetosphere at ultra-relativistic speeds. Therefore, it is necessary to take the plasma back reaction into account for the pulsar period and magnetic inclination angle evolution. State-of-the-art simulations of the magneto-thermal evolution ) and revised magnetospheric models (Philippov et al. 2014) allow for a more accurate prediction of the long-term behaviour of magnetic field strength, angular momentum loss rate, and the inclination angle (angle between the magnetic dipole moment and the rotational axis). Following these studies, Gullón et al. (2014) revisited the population synthesis of isolated radio-pulsars incorporating a magnetic field decay in the framework of the vacuum approximation and realistic magnetospheric models to be compared with each other. Following the work of Gullón et al. (2014), we developed an evolution model for the population study of canonical pulsars in both radio and gamma-ray bandwidths by taking into account both the evolution of the dipole magnetic inclination angle and the decay of the magnetic field in the force-free approximation. The striped wind model was used for the first time in a PPS study to describe the observed γ-ray pulsar population. We used a polar cap model for radio emission. We compare the results of our simulations with the observations by using the shape of the P −Ṗ diagram, as well as the number of radio-only, gamma-only, and gamma-ray radio loud pulsars (see Table 1). In order to better constrain the radio emission sites from the polarization data following the rotating vector model, we focus only on young (non-recycled) pulsars for which the radio emission height is reasonably well constrained.
The paper is organised as follows. In Section 2 we outline the model that we used to generate the observed pulsar population, and their detection is addressed in Section 3. In Section 4, we present the results of our simulations, and discuss their signification in Section 5. A summary is proposed in Section 6.

P −Ṗ evolution model
We start our population synthesis analysis by describing the underlying model of the period and luminosity evolution, starting from an initial sample of neutron stars with magnetic field strength B 0 and birth period P 0 . We use a Gaussian distribution for the initial spin period and a Gaussian in log space distribution for the magnetic field, such that We assume mean values ofP = 60 ms andB = 2.5 × 10 8 T, and standard deviation σ p = 10 ms and σ b = 0.5. These distributions are similar to the ones used by Gullón et al. (2014), Johnston et al. (2020), Faucher-Giguère & Kaspi (2006, Yadigaroglu & Romani (1995), and Watters & Romani (2011). We also assume an isotropic distribution of magnetic inclination angles α with respect to the rotation axis, meaning that the variable cos α is uniformly distributed between zero and one. Moreover, in our model, we generate a population of pulsars with a constant birth rate of 1/70 yr. It is important to note that the age of each pulsar should be taken uniformly from age zero up to the age of the Milky Way, that is about the age of the Universe 13 · 10 9 yr. Hence, a total number of 1,9 · 10 8 pulsars should normally be simulated. Here we choose to generate only 10 7 pulsars since a higher number of simulated pulsars, such as 10 8 pulsars, increases the simulation time significantly. However, we check that in a few instances of simulations using 10 8 pulsars the final results only change by a few percent. We include in our model a magnetic field decay with a decay rate τ d = k τ d · τ v with τ v the decay rate given by Viganò et al. (2013) by extrapolation to a lower field strength compared to the magnetars they studied. The pulsar, then evolves with a decreasing total spin down powerĖ which serves as a proxy for the radio and gamma-ray luminosity, as shown later.
In order to derive the pulsar's detectability, we need to know its distance d to the Earth. To summarize, each pulsar has intrinsic characteristics such as: B 0 (magnetic field at birth), P 0 (period at birth), P andṖ (current period and its time derivative), α 0 (inclination angle at birth taken isotropically-uniform in cos α 0 ), α the current inclination angle, n Ω (unit vector along the rotation axis), (x 0 , y 0 , z 0 ) (birth position in Galactocentric coordinates), v 0 (birth kick velocity), and d (distance to Earth).
Pulsars are spread all over the Galaxy with a radial and height distribution above the Galactic plane deduced from observations. In a Cartesian coordinate system attached to the galaxy, they are located at individual positions denoted by (x, y, z). We try to retrieve the pulsar current spatial distribution from an initial distribution where they are concentrated in the Galactic plane.

Birth and evolution of the pulsars
The current period P and magnetic moment inclination angle α are evolved from their initial values at birth, according to some spin evolution models. We consider a force-free evolution scenario, as well as the evolution of the inclination angle α. The τ birth (1/yr) P mean (ms) B mean (T) σ P (ms) σ B α d k τ d 70 60 2.5 × 10 8 10 0.5 1.5 5 Table 2: Parameters used in our simulations.
parameters that we use in our simulations are summarized in Table 2.

Vacuum dipole
The most studied rotator is a magnetic dipole radiating in vacuum. Exact solutions have been given by Deutsch (1955) but the point dipole formula is sufficient to accurately evolve the period (Jackson 2001). The spin down luminosity is given for the orthogonal rotator by and for an oblique rotator by where α is the magnetic obliquity, B the magnetic field strength at the magnetic equator, Ω = 2π/P the rotation speed, R the neutron star radius, and µ 0 the vacuum permeability, which has the value µ 0 = 4π × 10 −7 H/m. This spin down removes energy and angular momentum from the star, leading to a braking with an increase in the period at a rateṖ related to the rate of rotational kinetic energy lossĖ, such thaṫ whereΩ is the spin frequency derivative. The canonical value of the moment of inertia is I = 10 38 kg m 2 . Combining Eq.
(3) and Eq. (4), the expected evolution of the rotation frequency becomeṡ The spin-down of a pulsar is generally written in a more concise form: where n is the braking index. From Eq. (5), we deduce its value to be n = 3, which also holds for a force-free magnetosphere, (see for instance Pétri (2016)), and with (K 1 = 8 .  From Eq. (5), assuming a constant obliquity α = α 0 , we can derive the period P evolution as From Eq. (5), we derive PṖ = 4π 2 K vac andĖ = 4π 2 IṖP −3 from Eq. (4). These last equations plot with straight lines in the P −Ṗ diagram and are usually shown in the approximation of a vacuum field with no inclination angle evolution. The real path of a single pulsar can significantly deviate from the line of constant B.

Force-free
The vacuum model can be adapted to the force-free model by replacing the vacuum spin down by its force-free counterpart, as given by Spitkovsky (2006) and Pétri (2012) The constant K for the spin evolution then becomes The period then also follows Eq. (9) but with K vac replaced by K ffe .

Evolution of the inclination angle
Concomitantly, the magnetic axis tends to align with the rotation axis, with approximately the same timescale as the spin down. Therefore, the stellar braking and spin alignment must be evolved consistently in accordance with the electromagnetic torque exerted on its surface. A simple analytic solution for the time evolution of α for a vacuum rotator was found by Michel & Goldwire (1970), who determined an exponentially fast alignment according to where τ vac align = 1.5τ 0 cos −2 α 0 is the alignment timescale of a vacuum pulsar, and is the characteristic spin-down time of a pulsar, with µ the norm of the magnetic moment. The subscript '0' refers to the values of pulsar variables at birth, t = t 0 . Thus, the vacuum pulsar evolves to the aligned configuration exponentially fast, without a significant slow-down of its rotation. When alignment is almost reached, due to the integral of motion Ω cos α = Ω 0 cos α 0 , the asymptotic rotation rate becomes Ω = Ω 0 cos α 0 , which is not significantly different from Ω 0 for high initial inclination. For a force-free magnetosphere, the situation is radically different. Philippov et al. (2014) indeed investigated the evolution of non-spherical pulsars with a plasma filled magnetosphere with the help of magneto-hydrodynamic (MHD) simulations. An analytic solution to the evolution of pulsar obliquity α(t) derived by Philippov et al. (2014) is where τ MHD align = τ 0 sin 2 α 0 / cos 4 α 0 is the MHD pulsar alignment timescale. The inclination angle asymptotes at later times to a power law, α(t) ∝ t −1/2 . The alignment timescale is therefore much longer than for a vacuum model. The actualṖ values are drastically different. Whereas in the vacuum case the spin down vanishes for an aligned rotator, for an MHD model, even the aligned rotator spins down. Consequently the period derivativeṖ decreases much slower for the latter model. The integral of motion is now and contrary to the vacuum case, the final asymptotic rotation rate always tends to zero.

Magnetic field decay
The neutron star magnetic field is known to decay with time, depending on the temperature and strength of the field. In order to take this effect into account in our population synthesis, for simplicity, we assume that the magnetic field decays according to a power law: where a and α d are constant parameters controlling the speed of the magnetic field decay. They are assumed to be independent of the neutron star model. Integrating in time, the magnetic field evolves as with the initial magnetic field strength B 0 and the decaying timescale τ d = 1/(α d a B α d 0 ). We notice that, since a α d is a constant, this decaying time depends on B 0 . Thus, for a different initial magnetic field strength B 1 , this timescale changes to For a decaying field, Eq. (14) is no more valid. It must be replaced by ln(sin α 0 )+ 1 2 sin 2 α 0 +KΩ 2 0 cos 4 α 0 sin 2 α 0 If the decaying is very slow, τ d t and the equation simplifies into Eq. (14). The magnetic field decay timescale is affected by the initial field strength and the current age of each pulsar, as given by Eq. (18). The typical decay timescale for a mean magnetic field of 2.5 × 10 8 T is 4.6 × 10 5 yr, which is consistent with the estimate of Viganò (2013) for the same field strength. However, a pulsar with a magnetic field lower than the mean value will have longer decaying timescales and vice versa. Moreover, since we use Eq. (18), a value for a is not needed.

Galactic distribution
It is usually agreed that the z-distribution above the Galactic plane of a population I object is approximately exponential (Binney & Merrifield 1998). The scale height of the z-distribution of pulsar progenitors is of the order 50 − 100 pc (Mdzinarishvili & Article number, page 4 of 12 Ludmilla Dirson , Jérôme Pétri, and Dipanjan Mitra : The Galactic population of canonical pulsars Melikidze (2004)). Today, this disparity is explained by the fact that pulsars are moving away from their initial birthplace due to a large kick velocity.
To describe the position of neutron stars within the galaxy, we work in the right-handed Galactocentric coordinate system (x, y, z) that has the Galactic centre at its origin, with y increasing in the disc plane towards the location of the Sun, and z increasing towards the direction of the north Galactic pole. In order to calculate the distance d between the Earth and the pulsar, we need to know the Sun's position with respect to the Galactic centre and with respect to the Galactic plane. Currently, neither of these quantities are absolutely known. The Sun is thought to be located at about z = 15 pc ((Siegert 2019)) above the Galactic plane and about y = 8.5 kpc from the Galactic centre. Therefore, the distance of a pulsar from us is We also use the cylindrical coordinate system (R, θ, z), with R the distance to an axis passing through the Galactic centre in the plane of constant z. The Cartesian coordinates (x, y, z) are related to the cylindrical coordinates (R, θ, z) in the usual way: (20)

Initial distribution
The initial location of the pulsars is given by the distributions found by Paczynski (1990) for the radial direction: and for the altitude by where R is the axial distance from the z−axis, and z is the distance from the Galactic disc. Numerical values are taken to be, R exp = 4.5 kpc and a R ≡ [1 − e 1−R max /R exp (1 + R max /R exp )] −1 = 1.0683, with R max = 20 kpc. Following van der Kruit (1987), we take z exp = 75 pc.

Birth kick velocity
To describe the supernova kick velocity, we proceed similarly to Hobbs et al. (2005), and use a Maxwellian distribution with a characteristic width of σ v young = 265 km s −1 for pulsars with age < 3 Myr and σ v old = 75 km s −1 for older ones. Therefore, one can write the element of velocity space as d 3 v = dv x dv y dv z , for velocities in a standard Cartesian coordinate system. In this case, the distribution for a single direction is a normal distribution with mean µ v z = µ v x = µ v y = 0 and standard devia- The mean velocity in three dimensions and the velocity dispersion in one dimension are connected bȳ v ∼ σ v √ 8/π 1.6σ v .

Detection
In order to see if our generated set of pulsars is responsible for the observed set, we need to determine if they are actually detected. The pulsar detection is governed by three main factors. The first is the beaming fraction, which indicates the fraction of the sky covered by the radiation beam (in a particular wavelength, here we suppose radio or gamma-rays), and thus the fraction of the population that is possibly detectable from Earth. The beaming fractions of radio and γ-ray pulsars differ from each other and depend on the pulsars' spin, geometry, and location of the emission regions. In this paper, we use the force-free model to calculate the γ−ray beaming fraction. The second and third factors are the pulsars' luminosities and the sensitivity of a given radio or γ−ray survey. Indeed, the detection of pulsars depends on their brightness, and location, and on the sensitivity of the instruments. This requires an emission model for a pulsar. The most basic requirement of such a model is the presence of a dense plasma around the strongly magnetized neutron star. The plasma streams along the open magnetic field line regions and the radio emission rises close to the neutron star's surface, while the gamma ray emission is constrained to originate close to the light cylinder and preferably outside it, within the striped wind. When combined, these two factors indicate if a given pulsar can be actually detected by a given instrument based on the geometry. To calculate the luminosity of each pulsar, we need its distance d from us. These three factors together provide us with the actual number of detected pulsars.

Beaming fraction
We assume isotropic distributions for the Earth viewing angle ξ and a conventional width of the radio beam denoted by ρ. The orientation of the unitary rotation vector n Ω = Ω Ω is assumed to be isotropic, meaning it is uniform in φ and in cos θ.
The Cartesian coordinates of the unit rotation vector n Ω are sin θ n Ω cos φ n Ω , sin θ n Ω sin φ n Ω , cos θ n Ω . n is the unit vector along the line of sight, with coordinates: When n, µ, and n Ω are in the same plane, we can define the line of sight inclination angle as ξ = ( n Ω , n), and the magnetic axis inclination angle as α = ( n Ω , µ). From these angles, we deduce the usual minimum angle between the magnetic axis and line of sight as β = ( µ, n). The 'impact angle', β represents the closest path of the line of sight to the magnetic axis.

Radio emission model
The polar cap geometry can account for most of the observed radio emission properties. This generally accepted model suggested that the radio emission is enclosed in a cone-shaped beam centred on the magnetic axis. The opening angle of the conical envelope is usually defined as the half-opening angle ρ, which depends upon the width of the open field region at the emission height, h em . The observed pulse width w r , depends on geometrical factors, namely how the observer's line of sight cuts the emission cone. The detailed geometry and efficiency of the radio emission is now presented.

Radio beam geometry
The radio beam geometry is explained in depth, for instance, in Chapter 3 of Lorimer & Kramer (2004). We summarize useful results in this section. A simple geometrical relationship between the half opening angle of the emission cone (ρ), the emission height (h em ), and the spin period is given via ρ = 3 π h em 2 P c radians (24) Eq. (24) is well verified for slow pulsars (Lorimer & Kramer 2004). The emission height h em is roughly constant, with a mean value of h em ∼ 3.10 5 m across the population. This can be inferred from observations via Eq. (24), as a number of studies found consistently ρ ∝ P −0.5 (Mitra 2017). With knowledge of the height h em , the period distribution P, and α for each pulsar, the radio beaming fraction can then be computed for all simulated pulsars.
The pulsar is detectable in radio if β = |ξ − α| ≤ ρ (corresponding to the cone located in one hemisphere), or if |ξ − (π − α)| ≤ ρ (if it is located in the other hemisphere). We must also satisfy α ≥ ρ and α ≤ π − ρ in order to effectively see radio pulsation, because the line of sight must go in and out from the emission cone to observe pulsations.
Equation (24) represents the opening angle of a 'fully filled' open-field line region. Consequently, the width of the radio profile, w r , can be calculated from ρ, and the geometry depicted by the angles α and ξ via (see Lorimer & Kramer (2004)) cos ρ = cos α cos ξ + sin α sin ξ cos(w r /2).
This relation serves to compute the observed pulse width w r of our sample, knowing all other parameters.

Radio luminosity
To estimate if our simulated pulsars can be detected in current pulsar surveys, we exactly follow the procedure suggested by Johnston et al. (2020), where the radio flux density F r of a pulsar at 1.4 GHz is written as where d is the distance in kpc and F j is the scatter term, which is modelled as a Gaussian with a mean of µ = 0.0 and a variance of σ = 0.2. Then, for example, a pulsar withĖ = 10 29 W at a distance of 1 kpc has a mean flux density of 9 mJy. Thus, if in a radio survey the minimum flux is S min survey , then the detection signal-to-noise ratio (S/N) is The minimum flux that is related to its period P and width of radio emission w r is Herew r = w r P/2π. The scaling factor S 0 reflects the survey parameters. Johnston et al. (2020) noted that currently the two deepest pulsar surveys that detected a large number of pulsars are the Parkes multibeam survey (Manchester et al. (2001)) and the Arecibo survey (Cordes et al. 2006), where both these surveys for normal pulsars have a sensitivity of ∼0.15 mJy. Hence, for detecting a pulsar with S/N∼10, and using F r ≈ 0.15 mJy and pulsar widthw r = 0.1P, the estimated S 0 ∼ 0.05 mJy. We use the same criteria as those of Johnston et al. (2020) for pulsar detection, and in our simulation, we obtain F r , P, andw r and if the S/N> 10, then we identify the pulsar as detected.

Gamma-ray emission model
Our gamma-ray photon production relies on the emission emanating from the current sheet within the striped wind. The ultrarelativistic plasma velocity directed along the radial direction radiates along the velocity vector within a small cone of halfopening 1/Γ, where Γ is the wind Lorentz factor. The detailed geometry and efficiency of our model is now described.

Gamma-ray geometry
As shown by , the striped wind geometry in the split monopole force-free magnetosphere resembles the actual dipole force-free magnetosphere, connecting the current sheet wobbling angle around the equator to the pulsar magnetic obliquity and the gamma-ray peak separation ∆. The analytical relation reads: A distant observer detects only gamma-ray photons when the current sheet crosses his or her line of sight. This constrains the line of sight to remain along the equator, with an inclination angle bounded by |ξ − π/2| ≤ α. In our population synthesis, we do not investigate the precise pulse profile, discriminating between a single or double peak gamma-ray profile. To complete the gamma-ray emission part, we also need to prescribe the associated luminosity and its possible detection by current instrumentation.

Gamma-ray luminosity
The gamma-ray luminosity function L γ is extracted from the recent study of Kalapotharakos et al. (2019), where they show that it is described by a fundamental plane such that L γ = f (ε cut , B,Ė), where cut is the cut-off energy and B the surface magnetic field strength. Since cut is not accessible in our model, we choose their 2D model, not requiring information about the cut-off energy: This expression remains very similar to the one given by (Pétri 2011a) where he showed that except for the weak dependence on the magnetic field introduced by Kalapotharakos et al. (2019). We use Expression (30) for the gamma-ray luminosity prescription.
The γ-ray flux as detected on Earth, F γ , is then simply L γ corrected for the line-of-sight cut with beaming fraction f Ω and divided by the square of the distance d to a given pulsar. Explicitly we get For the striped wind model, the f Ω factor has been computed by Pétri (2011b), and varies between 0.22 and 1.90. We use a rough approximation for the beaming fraction term f Ω such that if α < −ξ +0.6109 then f Ω = 1.9 and f Ω = 1 otherwise. From the sensitivity expectation of the Fermi/LAT instrument, we assume that the sensitivity to pulsars at Galactic latitudes < 2 • is F min = 4 × 10 −15 W m −2 and 16 × 10 −15 W m −2 for blind searches.

Simulations
In this paper, we mainly focus our attention on reproducing the form of the spin period versus the period derivative P −Ṗ diagram and the number of pulsar types. To compare our simulations with the observations, we exclude the binary pulsars and the pulsars with P < 20 ms from the ATNF catalogue. There is no need to exclude the pulsars with P < 20 ms in our simulations, since there are virtually none. Given that the intrinsic parameters provided in Section 2 are poorly constrained, they have been varied within their likely minimum and maximum values until the model satisfactorily reproduces both the P −Ṗ diagram and the number of radio-only, radio-loud gamma-ray, and gamma-ray pulsars (forĖ > 10 31 W andĖ > 10 28 W (see Table. 4)). This enables us to constrain the allowed regions in the parameter space of the pulsars' intrinsic parameters. We get a reasonable estimate of the observations with the parameters already given in Table 2. It should be noted that the parameters are changed manually and, to save time, we do not use a multidimensional root finder algorithm.

P −Ṗ diagram
The resulting P −Ṗ diagram extracted from our simulations is shown in Fig. 1 panel (a) for the decaying magnetic field and panel (b) for the constant magnetic field, along with the pulsar data from the ATNF catalogue. Those diagrams are obtained with the parameters given in Table 2. As it can be seen, the model including the magnetic field decay, better matches the observations than the model of the constant magnetic field. We notice that with a decaying magnetic field, more pulsars are lying in the lower part of the diagram (meaning a smallerṖ), whereas for a constant magnetic field, more pulsars are lying in the right part of the diagram (meaning longer periods). This is due to the effect of magnetic field evolution, which clearly causes the longer period pulsars to slow down less rapidly than with a constant magnetic field. This mismatch could be improved by refining the distributions of periods and magnetic fields at birth, but other parameters could also impact the current P −Ṗ distribution, for instance the detection limits. The radio emission mechanism can also affect this discrepancy. The exact nature of these dependencies that accurately represents the observed pulsar population needs a detailed investigation, which is beyond the scope of this work. In the P −Ṗ diagram, the region with a magnetic field greater than ∼ 10 8 T and 10 23 W <Ė < 10 28 W is also depleted. Indeed, the high magnetic field pulsars are not taken into account in our study because they probably belong to a distinct class of magnetized neutron stars. For a better representation of the P −Ṗ diagram, we also show the histograms of periods P in Fig. 2 and period derivativesṖ in Fig. 3, along with the data.

Pulsar types
We compute the total number N r , N g , and N rg , which are the number of radio-only, gamma-only, and radio-loud gamma-ray pulsars, respectively, and within the energy bandsĖ = 10 28 W andĖ = 10 31 W. The resulting numbers are within the observations, and are better constrained by the decaying magnetic field model (see Table 3 and Table 4). We represent the P −Ṗ diagram for N r , N g and N rg in Fig. 4, Fig. 5, and Fig. 6 for both the simulations and the observations. We find good agreement between the observed and simulated population separately for each subclass.   Table 4: Quantities N r , N g , and N rg are the number of radioonly, gamma-only, and radio-loud gamma-ray pulsars, respectively, obtained from our simulations with a decaying magnetic field.

Spatial distribution
The distances d from Earth to the detected pulsars are shown in Fig. 7 for both the simulations and the observations. We find that more pulsars are detected at closer distances, and less at further distances. The projected 2D distribution is shown in Fig. 9, which represents the Galactic coordinate y as a function of x. As was mentioned earlier, since we detect more pulsars with lowĖ, we detect more pulsars that are closer to us (see Eq. 26). The distribution of the latitude for the simulated and observed population is shown in Fig. 8. In our simulations, we detect fewer pulsars at a lower Galactic latitude. This is due to the spatial distribution of pulsars, which shows pulsars closer to the Sun compared to observations, and therefore artificially increases their latitude. This may indicate that a more refined treatment of the sensibility depending on the latitude, and a better description of the initial spatial distribution as well as the birth kick velocity, should be implemented. Furthermore, our current simulation does not take into account the Galactic potential. The simulated pulsars' latitude spreads would shift more towards the Galactic plane, thus increasing the number of pulsars in lower Galactic latitudes, and perhaps better representing the data.

Influence of the parameters
Here we want to show the influence of the parameters on the quantities N r , N g , and N rg , which are the number of radio only, gamma only, and radio-loud gamma-ray pulsars, respectively, from our simulations with the parameters given in Table 2. This is summarized in Tables 5, 6, ??, 7, 8, and 9. The most sensitive parameters are the birth rate, the initial mean magnetic field, and σ b . First the number of detected pulsars N tot scales linearly with the birth rate τ birth . Increasing the birth rate by a factor of two decreases N tot by the same factor. We indeed detect statistically more pulsars if more are present within our Galaxy. This trend is clearly seen for the rates 1/50yr, 1/100yr, and 1/200yr (see Table 5).
The mean magnetic field also drastically impacts N tot . The number of detected pulsars decreases with increasing field strength. This is mainly due to the fact that high magnetic pul-Article number, page 7 of 12 A&A proofs: manuscript no. 43305corr     sars align faster than low magnetic field pulsars. Furthermore, the magnetic field decay is also faster for higher magnetic fields.  Table 5: Influence of the birth rate on the quantities N r , N g , and N rg .
The initial mean period moderately impacts the detected pulsar number, as shown in Table ??. Increasing this period by a factor of ten only increases N tot by a factor of two. Moreover, N tot remains insensitive to the spread in this birth period (see Table 7).
Increasing the spread in magnetic field strength allows for lower B field values and, as seen above, a larger number of pul-    Table 6: Influence of the initial magnetic field on the quantities N r , N g , and N rg , which are the number of radio only, gamma only, and radio-loud gamma-ray pulsars, respectively, from our simulations with the parameters given in Table 2. sars will be detected. In this paper, we use S /N = 10 to get the detected pulsars, but in reality pulsars have been detected with telescopes having different sensitivities. Table 9 allows us to better understand how the detection changes with the S/N.    Table 7: Influence of the σ P on the N r , N g , and N rg .   Table 9: Influence of the S/N on the quantities N r , N g , and N rg .
4.5. Radio width w r as a function of the period Figure 10 represents the logarithm of the width of the radio pulse profile log(w r ) as a function of the logarithm of the period log(P), to compare our simulations with the data. We retrieve the trend shown in the observations, namely a decrease of the width with an increase in the period. The slope in log-log scale is −0.53 ± 0.02, meaning w r = (1.04 ± 0.01)P −0.53±0.02 as expected by the emission cone model. A linear fit to the measured pulsar width data from Posselt et al. (2021) gives w 10 = (1.17 ± 0.01)P −0.34±0.03 , which, although indicative of the expected trend, still does not agree with the index of ∼ −0.5. We believe that a more robust analysis and a larger data set will be more appropriate to verify the width period relation, which is beyond the scope of this work. However, although expected from geometrical considerations, no small values for w r have been obtained, particularly for faster periods. This may be explained by the fact that fast rotating pulsars usually show a larger opening angle for the emission beam, naturally leading to a broadening of the pulse profile. Figure 11 shows the initial obliquity cos(α 0 ) at birth and the current obliquity cos(α) distribution of the pulsars, irrespective of their detection. Most of the pulsars tend to the aligned configuration when aging. This is due to the fact that older pulsars are in larger number and are more likely to be aligned. Indeed, Fig. 13 highlights that most of the detected pulsars in our sample are older than 10 6 yr. This is expected since the birth rate is constant up to Gyr. Figure 12 shows the distributions of initial cos(α 0 ) and current cos(α) obliquity for the detected pulsars. The cur- Fig. 10: log(w r ) as a function of log(P). The quantity w r is expressed in degrees and P is expressed in s. The data are the width at the10% level taken from Posselt et al. (2021). The solid blue line and the shaded area represent the result of a linear fit to the data and its uncertainties, and the solid red line and the shaded area represent a linear fit to the simulations, with its uncertainties. rent obliquity distribution cos α is almost uniform, whereas the distribution at birth cos(α 0 ) peaks for values less than 0.4. This means that mostly only pulsars with larger initial inclination angles, close to orthogonal rotators, will be detected. Moreover, the alignment timescale is longer for high obliquities, meaning small cos α 0 (see Fig. 14). Since the pulsed emission becomes harder to detect when the pulsar approaches an aligned rotator, pulsars with a longer alignment timescale are favoured, and therefore there is a bias towards small cos α 0 .

Discussion
Population synthesis of neutron stars represents a valuable and reliable tool to constrain the basic physics of neutron star electrodynamics, from its surface up to large distances well outside the light cylinder. Following the work of Gullón et al. (2014), we have shown that a realistic pulsar model such as the force-free evolution scenario can self-consistently account for the observed  population of canonical pulsars, and thus verified the findings of Gullón et al. (2014). Additionally, we have considered the sub-population of gamma-ray pulsars associated with the striped wind emission model. We find a birth rate of 1/70 yr, which is consistent with previous studies (Gonthier et al. 2004), (Watters & Romani 2011). In the case of a constant magnetic field, we find that a birth rate of 1/300 yr improves the matching to the observations. However, this value is not consistent with expectations from previous works. The best results are obtained with a decaying magnetic field model.
In our study, we focused only on canonical pulsars for several reasons. First, the radio emission height and geometry were well constrained only for pulsars with periods P > 20 ms. Second, the binary pulsars were mostly revived through an accretion phase that was not included in our model. Third, magnetars with super-critical magnetic field strengths require a special initial magnetic field distribution, deviating from the normal canonical pulsar population we used. These magnetars are probably produced by a dynamo effect right at their birth during the proto-neutron star stage, increasing their internal magnetic field by several orders of magnitude through turbulent motion or magneto-rotational instabilities. Although achievable, we did not include such processes that would lead to a second distribution of magnetic field strengths designed especially for magnetars.
The spatial distribution of neutron stars at birth was postulated and adapted from the literature to fit as well as possible the current observed Galactic distribution, taking into account selection effects related to the sensitivity of radio and high energy telescopes. We found a reasonable agreement between the current spatial distribution and the simulated distribution.
The radio emission geometry constrained by the dipolar region combined with the striped wind emission model for gamma-rays are able to reproduce the entire population of young pulsars summarized in the P −Ṗ diagram. The number of radio only, gamma-ray only, or radio-loud gamma-ray pulsars fits the observed population as given by the ATNF catalogue. In a previous study,  demonstrated that the radio-loud gamma-ray pulsar emission model accurately constrained the geometry of the magnetic axis and line of sight axis with respect to the rotation axis. In this paper, we went one step further, showing that our model is consistent with the entire canonical pulsar population.
There are few other selection effects that can affect the detection of a neutron star in the same way as a radio pulsar, and we have not included these in our analysis. One such selection effect is the propagation of the radio signal in the ISM, which is subject to a scattering effect that smears out the signal, rendering the detection of some pulsars very difficult because the pulsation disappears at low frequencies. Interstellar scattering is associated with the fluctuation of the electron distribution along the observer's line of sight, where the pulsar signal suffers a multi-path scattering as the signal traverses the ISM, (see Scheuer (1968)). The scattering effect, therefore, can be described as the pulse width to be convolved with the response function of the ISM, and as a result the observed pulse width appears to be smeared. If the smearing exceeds the pulse period, then the pulsar cannot be detected as a pulsed signal. The effect of scattering is known to increase with pulsar dispersion measure (see, e.g. Krishnakumar et al. (2015)) and hence, also the distance. Thus, the farther the pulsars, the more difficult they becomes to detect. The scattering timescale is a strong function of the observing frequency ν with the timescale decreasing as ∼ ν −4.4 (see, e.g. (Löhmer et al. 2004)), hence sensitive pulsar searches at higher frequencies can in principle be used to detect more pulsars. However, due to steep radio pulsar spectra, pulsars are weaker at higher frequencies, and thus it may not be viable to detect a highly scattered pulsar.
The other selection effect is related to the coherent radio emission mechanism for radio pulsars. Our analysis here considers the geometrical model of a star centred magnetic dipole for the radio pulsars. However, the condition under which the radio emission mechanism can operate in pulsars are not included. Various studies suggest that the presence of a surface multi-polar field facilitates the generation of pair plasma, which in turn can assist mechanisms of coherent radio emission. Thus, the nature of the multi-polar surface field can affect the number of detectable pulsars, and several studies suggest a death-line in the P−Ṗ diagram (see e.g. Chen & Ruderman (1993); Mitra et al. (2020)). In the current study, however, the pulsars detected via simulations are similar to the actual pulsars, and thus the above selection effects do not appear to have any major impact on the number of detected pulsars.

Summary
In this paper, we have shown that the secular evolution of the neutron star spin rate and its braking is accurately captured by the force-free magnetosphere model, in which radio photons emanate from regions near the polar cap, whereas the gammaray photons are produced within the current sheet of the striped wind. Following simple prescriptions for the birth kick, the birth period, and the birth magnetic field strength, assuming an isotropic obliquity and an isotropic viewing angle, we were able to pin down the essential parameters of the canonical pulsar population. Almost orthogonal rotators at birth represent the vast majority of pulsars currently detected due to the alignment effect, the geometry of the radio, and gamma-ray beams.
Our investigation could be improved in several ways, as exposed in the previous Section 5. Moreover, knowing the geometry of each individual pulsar, we could model the thermal X-ray emission from the polar caps, for instance for PSR J1136+1551 (Pétri & Mitra 2020). The upcoming advent of the SKA telescope will increase by one order of magnitude the number of known pulsars and constrain even further the evolution scenario and emission physics of these stellar remnants.