© L. Dirson et al. 2022

1. 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 3000¹ (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². 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, $\dot{E}$- 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³ 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. (2002), (2004), (2007a,b). 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₀ = 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, 1996). The slot-gap (along the last open magnetic field 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. 1986, 1986b; Hirotani 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, 2011a). 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 (Viganò et al. 2013) 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 Sect. 2 we outline the model that we used to generate the observed pulsar population, and their detection is addressed in Sect. 3. In Sect. 4, we present the results of our simulations, and discuss their signification in Sect. 5. A summary is proposed in Sect. 6.

2. 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₀ and birth period P₀. We use a Gaussian distribution for the initial spin period and a Gaussian in log space distribution for the magnetic field, such that

$$\begin{array}{c}\hfill p(log(B_0))=\frac{1}{{\sigma }_b\sqrt{2\pi }}{e}^{-{(log{B}_0-log\overline{B})}^2/(2{\sigma }_b^2)}.\end{array}$$(1)

We assume mean values of $\overline{P}=60$ ms and $\overline{B}=2.5\times {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⁹ yr. Hence, a total number of 1.9 × 10⁸ pulsars should normally be simulated. Here we choose to generate only 10⁷ pulsars since a higher number of simulated pulsars, such as 10⁸ pulsars, increases the simulation time significantly. However, we check that in a few instances of simulations using 10⁸ 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 $\dot{E}$ 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₀ (magnetic field at birth), P₀ (period at birth), P and Ṗ (current period and its time derivative), α₀ (inclination angle at birth taken isotropically-uniform in cosα₀), α the current inclination angle, n_Ω (unit vector along the rotation axis), (x₀, y₀, z₀) (birth position in Galactocentric coordinates), v₀ (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.

2.1. 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 parameters that we use in our simulations are summarized in Table 2.

2.1.1. 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

$$\begin{array}{c}\hfill L_{\perp }={\displaystyle \frac{8\pi }{3{\mu }_0{c}^3}}{B}^2{\mathrm{\Omega }}^4{R}^6,\end{array}$$(2)

and for an oblique rotator by

$$\begin{array}{c}\hfill L_{\mathrm{dip}}=L_{\perp }{sin}^2\alpha ,\end{array}$$(3)

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 μ₀ the vacuum permeability, which has the value μ₀ = 4π × 10⁻⁷ 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 $\dot{E}$, such that

$$\begin{array}{c}\hfill \dot{E}=\frac{\mathrm{d}{E}_{\mathrm{rot}}}{\mathrm{d}t}=-I\mathrm{\Omega }\dot{\mathrm{\Omega }}=L_{\mathrm{dip}},\end{array}$$(4)

where $\dot{\mathrm{\Omega }}$ is the spin frequency derivative. The canonical value of the moment of inertia is I = 10³⁸ kg m². Combining Eqs. (3) and (4), the expected evolution of the rotation frequency becomes

$$\begin{array}{c}\hfill \dot{\mathrm{\Omega }}={\displaystyle \frac{8\pi }{3{\mu }_0{c}^3}}{\displaystyle \frac{{(B{R}^{3}sin\alpha )}^2}{I}}{\mathrm{\Omega }}^3.\end{array}$$(5)

The spin-down of a pulsar is generally written in a more concise form:

$$\begin{array}{c}\hfill \dot{\mathrm{\Omega }}=-K_{\mathrm{vac}}{\mathrm{\Omega }}^n,\end{array}$$(6)

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

$$\begin{array}{c}\hfill K_{\mathrm{vac}}={\displaystyle \frac{8\pi }{3{\mu }_0{c}^3}}{\displaystyle \frac{{(B{R}^{3}sin\alpha )}^2}{I}}=K_1{B}^2{sin}^2\alpha ,\end{array}$$(7)

with (K₁ = 8πR⁶)/(3μ₀c³I). For a typical neutron star radius of R = 12 km, as found by recent NICER observations (Riley et al. 2019; Bogdanov et al. 2019), this constant is estimated:

$$\begin{array}{c}\hfill K_{\mathrm{vac}}\approx {10}^{-15}{\mathrm{s}}^{-1}{\left(\frac{B}{4.5\times {10}^8\mathrm{T}}\right)}^2{\left(\frac{R}{12\mathrm{km}}\right)}^6{\left(\frac{I}{{10}^{38}\mathrm{kg}{\mathrm{m}}^2}\right)}^{-1}{sin}^2\alpha .\end{array}$$(8)

From Eq. (5), assuming a constant obliquity α = α₀, we can derive the period P evolution as

$$\begin{array}{c}\hfill P(t)=2\pi \sqrt{2(K_{\mathrm{vac}}t+{\displaystyle \frac{P_0^2}{8{\pi }^2}})}.\end{array}$$(9)

From Eq. (5), we derive P Ṗ = 4 π² K_vac and $\dot{E}=4{\pi }^{2}I\dot{P}{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.

2.2. 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)

$$\begin{array}{c}\hfill L_{\mathrm{ffe}}=\frac{3}{2}{L}_{\perp }(1+{sin}^2\alpha ).\end{array}$$(10)

The constant K for the spin evolution then becomes

$$\begin{array}{c}\hfill K_{\mathrm{ffe}}=\frac{3}{2}{K}_1{B}^2(1+{sin}^2\alpha ).\end{array}$$(11)

The period then also follows Eq. (9) but with K_vac replaced by K_ffe.

2.3. 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

$$\begin{array}{c}\hfill sin\alpha (t)=sin{\alpha }_{0}exp(-t/{\tau }_{\mathrm{align}}^{\mathrm{vac}}),\end{array}$$(12)

where ${\tau }_{\mathrm{align}}^{\mathrm{vac}}=1.5{\tau }_0{cos}^{-2}{\alpha }_0$ is the alignment timescale of a vacuum pulsar, and

$$\begin{array}{c}\hfill {\tau }_0=\frac{I{c}^3}{{\mu }^2{\mathrm{\Omega }}_0^2}\approx {10}^4{\left(\frac{B_0}{{10}^8\mathrm{T}}\right)}^{-2}\left(\frac{P_0}{10\mathrm{ms}}\right)\mathrm{y}{r}^{-1}\end{array}$$(13)

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₀. 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α = Ω₀cosα₀, the asymptotic rotation rate becomes Ω = Ω₀cosα₀, which is not significantly different from Ω₀ 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

$$\begin{array}{c}\hfill \frac{1}{2{sin}^2\alpha (t)}+log(sin\alpha (t))=\frac{t}{{\tau }_{\mathrm{align}}^{\mathrm{MHD}}}+\frac{1}{2{sin}^2{\alpha }_0}+log(sin{\alpha }_0),\end{array}$$(14)

where ${\tau }_{\mathrm{align}}^{\mathrm{MHD}}={\tau }_0{sin}^2{\alpha }_0/{cos}^4{\alpha }_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

$$\begin{array}{c}\hfill \mathrm{\Omega }\frac{{cos}^2\alpha }{sin\alpha }={\mathrm{\Omega }}_0\frac{{cos}^2{\alpha }_0}{sin{\alpha }_0},\end{array}$$(15)

and contrary to the vacuum case, the final asymptotic rotation rate always tends to zero.

2.4. 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:

$$\begin{array}{c}\hfill \frac{\mathrm{d}B}{\mathrm{d}t}=-a{B}^{1+{\alpha }_d},\end{array}$$(16)

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

$$\begin{array}{c}\hfill B(t)=B_0{(1+t/{\tau }_d)}^{-1/{\alpha }_d},\end{array}$$(17)

with the initial magnetic field strength B₀ and the decaying timescale ${\tau }_d=1/({\alpha }_{d}a{B}_0^{{\alpha }_d})$. We notice that, since a α_d is a constant, this decaying time depends on B₀. Thus, for a different initial magnetic field strength B₁, this timescale changes to

$$\begin{array}{c}\hfill {\tau }_1{B}_1^{{\alpha }_d}={\tau }_d{B}_0^{{\alpha }_d}.\end{array}$$(18)

For a decaying field, Eq. (14) is no more valid. It must be replaced by

$$\begin{array}{cc}& ln(sin{\alpha }_0)+\frac{1}{2{sin}^2{\alpha }_0}+K{\mathrm{\Omega }}_0^2\frac{{cos}^4{\alpha }_0}{{sin}^2{\alpha }_0}\frac{{\alpha }_d{\tau }_d{B}_0^2}{{\alpha }_d-2}[{(1+\frac{t}{{\tau }_d})}^{1-2/{\alpha }_d}-1]\hfill \\ \hfill & =ln(sin\alpha )+\frac{1}{2{sin}^2\alpha }.\hfill \end{array}$$(19)

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⁸ T is 4.6 × 10⁵ 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.

2.5. 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 & 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 $d=\sqrt{{(x-x_{\odot })}^2+{(\mathit{y}-{\mathit{y}}_{\odot })}^2+{(z-z_{\odot })}^2}$. 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:

$$\begin{array}{c}\hfill \{\begin{array}{c}x=Rcos\theta \hfill \\ y=Rsin\theta .\hfill \end{array}\end{array}$$(20)

2.5.1. Initial distribution

The initial location of the pulsars is given by the distributions found by Paczynski (1990) for the radial direction:

$$\begin{array}{c}\hfill {\rho }_R(R)=\frac{a_R{e}^{-R/R_{\exp }}R}{R_{\exp }^2},\end{array}$$(21)

and for the altitude by

$$\begin{array}{c}\hfill {\rho }_z(z)=\frac{1}{2z_{\exp }}{e}^{-|z|/z_{\exp }},\end{array}$$(22)

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.0683, with R_max = 20 kpc. Following van der Kruit (1995), we take z_exp = 75 pc.

2.5.2. 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 σ_vyoung = 265 km s⁻¹ for pulsars with age < 3 Myr and σ_vold = 75 km s⁻¹ for older ones. Therefore, one can write the element of velocity space as d³v = dv_xdv_ydv_z, for velocities in a standard Cartesian coordinate system. In this case, the distribution for a single direction is a normal distribution with mean μ_vz = μ_vx = μ_vy = 0 and standard deviation σ_vx = σ_vy = σ_vx. The mean velocity in three dimensions and the velocity dispersion in one dimension are connected by $\overline{\mathit{v}}\sim {\sigma }_{\mathit{v}}\sqrt{8/\pi }\simeq 1.6{\sigma }_{\mathit{v}}$.

3. 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.

3.1. 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 ${\mathit{n}}_{\mathbf{\Omega }}\mathbf{=}{\displaystyle \frac{\mathbf{\Omega }}{\parallel \mathbf{\Omega }\parallel }}$ 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:

$$\begin{array}{c}\hfill ({\displaystyle \frac{x-x_{\odot }}{d}},{\displaystyle \frac{y-y_{\odot }}{d}},{\displaystyle \frac{z-z_{\odot }}{d}}).\end{array}$$(23)

When n, μ, and n_Ω are in the same plane, we can define the line of sight inclination angle as $\xi \mathbf{=}\mathbf{(}\widehat{{\mathit{n}}_{\mathbf{\Omega }}\mathbf{,}\mathit{n}}\mathbf{)}$, and the magnetic axis inclination angle as $\alpha \mathbf{=}\mathbf{(}\widehat{{\mathit{n}}_{\mathbf{\Omega }}\mathbf{,}\mu }\mathbf{)}$. From these angles, we deduce the usual minimum angle between the magnetic axis and line of sight as $\beta \mathbf{=}\mathbf{(}\widehat{\mu \mathbf{,}\mathit{n}}\mathbf{)}$. The ‘impact angle’, β represents the closest path of the line of sight to the magnetic axis.

3.2. 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.

3.2.1. 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

$$\begin{array}{c}\hfill \rho =3\sqrt{{\displaystyle \frac{\pi h_{\mathrm{em}}}{2Pc}}}\mathrm{radians}\end{array}$$(24)

Equation (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⁵ 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)

$$\begin{array}{c}\hfill cos\rho =cos\alpha cos\xi +sin\alpha sin\xi cos(w_{\mathrm{r}}/2).\end{array}$$(25)

This relation serves to compute the observed pulse width w_r of our sample, knowing all other parameters.

3.2.2. 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

$$\begin{array}{c}\hfill F_{\mathrm{r}}=9\text{mJy}{\left(\frac{d}{1\text{kpc}}\right)}^{-2}{\left(\frac{\dot{E}}{{10}^{29}\text{W}}\right)}^{1/4}{10}^{F_j},\end{array}$$(26)

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 $\dot{E}={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_{\mathrm{survey}}^{\min }$, then the detection signal-to-noise ratio (S/N) is

$$\begin{array}{c}\hfill S/N=\frac{F_{\mathrm{r}}}{S_{\mathrm{survey}}^{\min }}.\end{array}$$(27)

The minimum flux that is related to its period P and width of radio emission w_r is

$$\begin{array}{c}\hfill S_{\mathrm{survey}}^{\min }=S_0\sqrt{\frac{\stackrel{\sim }{w_{\mathrm{r}}}}{P-\stackrel{\sim }{w_{\mathrm{r}}}}}.\end{array}$$(28)

Here w̃_r = w_rP/2π. The scaling factor S₀ 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 width w̃_r = 0.1P, the estimated S₀ ∼ 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, and w̃_r and if the S/N >  10, then we identify the pulsar as detected.

3.3. Gamma-ray emission model

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

3.3.1. Gamma-ray geometry

As shown by Pétri & Mitra (2021), 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:

$$\begin{array}{c}\hfill cos(\pi \mathrm{\Delta })=|cot\xi cot\alpha |.\end{array}$$(29)

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.

3.3.2. 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_{\gamma }=f({\epsilon }_{\mathrm{cut}},B,\dot{E}),$ 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:

$$\begin{array}{c}\hfill L_{\gamma (2D)}={10}^{26.15\pm 2.6}\text{W}{\left(\frac{B}{{10}^8\text{T}}\right)}^{0.11\pm 0.05}{\left(\frac{\dot{E}}{{10}^{26}\text{W}}\right)}^{0.51\pm 0.09}.\end{array}$$(30)

This expression remains very similar to the one given by Pétri (2011a) where he showed that

$$\begin{array}{c}\hfill L_{\gamma }=\dot{E}\sqrt{\frac{{10}^{26}\text{W}}{\dot{E}}},\end{array}$$(31)

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

$$\begin{array}{c}\hfill F_{\gamma }=\frac{L_{\gamma (2D)}}{4\pi f_{\mathrm{\Omega }}{d}^2}.\end{array}$$(32)

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⁻¹⁵ W m⁻² and 16 × 10⁻¹⁵ W m⁻² for blind searches.

4. 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 Sect. 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 $\dot{E}>{10}^{31}$ W and $\dot{E}>{10}^{28}$ W). 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.

4.1. 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⁸ T and ${10}^{23}W<\dot{E}<{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.

Fig. 1. P−Ṗ diagram of the simulated population, along with the observations, including both gamma-ray and radio pulsars. (a) Decaying magnetic field case. (b) Constant magnetic field case.

Fig. 2. Distribution of the observed period taken from the ATNF catalogue, along with the simulations.

Fig. 3. Distribution of the observed period derivative Ṗ taken from the ATNF catalogue, along with the simulations.

4.2. 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 $\dot{E}={10}^{28}$ W and $\dot{E}={10}^{31}$ W. The resulting numbers are within the observations, and are better constrained by the decaying magnetic field model (see Tables 3 and 4). We represent the P−Ṗ diagram for N_r, N_g and N_rg in Figs. 4–6 for both the simulations and the observations. We find good agreement between the observed and simulated population separately for each sub-class.

Fig. 4. P−Ṗ diagram of the radio-only pulsars for both the simulations and the observations.

Fig. 5. P−Ṗ diagram of the gamma-only pulsars for both the simulations and the observations.

Fig. 6. P−Ṗ diagram of the radio-loud gamma ray pulsars for both the simulations and the observations.

4.3. 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. 8, which represents the Galactic coordinate y as a function of x. As was mentioned earlier, since we detect more pulsars with low $\dot{E}$, 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. 9. 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.

Fig. 7. Distribution of the distance to Earth for the simulated and observed populations.

Fig. 8. Detected and observed spatial distribution of pulsars projected onto the Galactic plane.

Fig. 9. Distribution of the observed Galactic latitude, b, in degrees, taken from the ATNF catalogue, along with the simulations.

4.4. 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–10. 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/50 yr, 1/100 yr, and 1/200 yr (see Table 5).

The mean magnetic field also drastically impacts N_tot. The number of detected pulsars decreases with increasing field strength (Table 6). This is mainly due to the fact that high magnetic pulsars align faster than low magnetic field pulsars. Furthermore, the magnetic field decay is also faster for higher magnetic fields.

The initial mean period moderately impacts the detected pulsar number, as shown in Table 7. 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 8).

Increasing the spread in magnetic field strength allows for lower B field values and, as seen above, a larger number of pulsars will be detected (Table 9). 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 10 allows us to better understand how the detection changes with the S/N.

Decreasing the S/N naturally increases the number of detected pulsars. For instance, switching from S/N = 10 to S/N = 5 multiplies the number of detected pulsars by a factor of 1.7.

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₁₀ = (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.

Fig. 10. log(wr) as a function of log(P). The quantity wr is expressed in degrees and P is expressed in s. The data are the width at the 10% 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.

4.6. Inclination angle

Figure 11 shows the initial obliquity cos(α₀) 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. 12 highlights that most of the detected pulsars in our sample are older than 10⁶ yr. This is expected since the birth rate is constant up to Gyr. Figure 13 shows the distributions of initial cos(α₀) and current cos(α) obliquity for the detected pulsars. The current obliquity distribution cosα is almost uniform, whereas the distribution at birth cos(α₀) 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α₀ (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α₀.

Fig. 11. Distribution of cosα and cosα0 for a sub-sample of 106 simulated population (‘before’ the detection).

Fig. 12. Age distribution of the detected and observed pulsars.

Fig. 13. Distribution of cos(α) and cosα0 for the ‘detected’ pulsars.

Fig. 14. Alignment timescale, ${\tau }_{\mathrm{align}}^{\mathrm{MHD}}$ as a function of cos(α0). Simulation for 107 pulsars.

5. 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, Pétri & Mitra (2021) 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.

6. 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 gamma-ray 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 Sect. 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.

Acknowledgments

This work has been supported by the CEFIPRA grant IFC/F5904-B/2018 and ANR-20-CE31-0010. We acknowledge the High Performance Computing Centre of the University of Strasbourg for supporting this work by providing scientific support and access to computing resources. We thank David Smith for stimulating discussions. DM acknowledges the support of the Department of Atomic Energy, Government of India, under project no. 12-R&D-TFR-5.02-0700.

References

All Tables

All Figures

	Fig. 1. P−Ṗ diagram of the simulated population, along with the observations, including both gamma-ray and radio pulsars. (a) Decaying magnetic field case. (b) Constant magnetic field case.
In the text

	Fig. 2. Distribution of the observed period taken from the ATNF catalogue, along with the simulations.
In the text

	Fig. 3. Distribution of the observed period derivative Ṗ taken from the ATNF catalogue, along with the simulations.
In the text

	Fig. 4. P−Ṗ diagram of the radio-only pulsars for both the simulations and the observations.
In the text

	Fig. 5. P−Ṗ diagram of the gamma-only pulsars for both the simulations and the observations.
In the text

	Fig. 6. P−Ṗ diagram of the radio-loud gamma ray pulsars for both the simulations and the observations.
In the text

	Fig. 7. Distribution of the distance to Earth for the simulated and observed populations.
In the text

	Fig. 8. Detected and observed spatial distribution of pulsars projected onto the Galactic plane.
In the text

	Fig. 9. Distribution of the observed Galactic latitude, b, in degrees, taken from the ATNF catalogue, along with the simulations.
In the text

Fig. 10. log(wr) as a function of log(P). The quantity wr is expressed in degrees and P is expressed in s. The data are the width at the 10% 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.

In the text

	Fig. 11. Distribution of cosα and cosα0 for a sub-sample of 106 simulated population (‘before’ the detection).
In the text

	Fig. 12. Age distribution of the detected and observed pulsars.
In the text

	Fig. 13. Distribution of cos(α) and cosα0 for the ‘detected’ pulsars.
In the text

	Fig. 14. Alignment timescale, ${\tau }_{\mathrm{align}}^{\mathrm{MHD}}$ as a function of cos(α0). Simulation for 107 pulsars.
In the text