11institutetext: Université de Strasbourg, CNRS, Observatoire astronomique de
Strasbourg, UMR 7550, 67000 Strasbourg, France 22institutetext: National
Centre for Radio Astrophysics, Tata Institute for Fundamental Research, Post
Bag 3, Ganeshkhind, Pune 411007, India 33institutetext: Janusz Gil Institute
of Astronomy, University of Zielona Góra, ul. Szafrana 2, 65-516 Zielona Góra,
Poland

# The Galactic population of canonical pulsars II

Mattéo Sautron 11 Jérôme Pétri 11 Dipanjan Mitra 2233 Ludmilla Dirson 11

(Received / Accepted)

###### Abstract

Context.  Pulsars are highly magnetized rotating neutron stars, emitting in a
broad electromagnetic energy range. These objects have been discovered more
than 50 years ago and are astrophysical laboratories for studying physics at
extreme conditions. Reproducing the observed population of pulsars refines our
understanding of their formation, evolution as well as radiation processes and
geometry.

Aims. In this paper, we improve our previous pulsar population synthesis by
studying the impact of the Galactic gravitational potential and of the death
line putting emphasis on the $\gamma$-ray pulsar population. In order to
elucidate the necessity of a death line, refined initial distributions for
spin periods and radial positions at birth were implemented, elevating the
sophistication of our previous simulations to the most recent state-of-the-
art.

Methods. The individual motion of pulsars in the Galactic potential is evolved
by a fourth order symplectic integration scheme. Our pulsar population
synthesis takes into account the secular evolution of a force-free
magnetosphere and the magnetic field decay. Each pulsar is evolved from its
birth up to the present time. The radio and $\gamma$-ray emissions are modeled
respectively by the polar cap geometry and the striped wind model.

Results.  When simulating a cohort of ten million pulsars akin to previous
endeavors, the integration of a death line proves essential for a better
agreement with observational trends. Yet, discrepancies emerges when modeling
a reduced ensemble of only one million pulsars, endowed with a higher birth
rate. This divergence yields a remarkably refined $P-\dot{P}$ diagram,
conspicuously with or without any delineating death line. Intriguingly, this
observation suggests an overestimation of the ages of detected pulsars in
existent data sets, thereby challenging the necessity of death lines in pulsar
population synthesis. Validation through Kolmogorov-Smirnov testing underlines
the statistical congruence between the spin period and spin period derivative
distributions of both observations and simulations. Moreover simulations with
a higher $\gamma$-ray instrumental sensitivity show that the GeV excess in the
Galactic centre could be of pulsar origin.

###### Key Words.:

pulsars: general - Gravitation - radio continuum: stars - Gamma rays: stars -
methods: statistical

##  1 Introduction

Pulsars, discovered by Jocelyn Bell Burnell in 1967 (Hewish et al., 1968) are
rotating neutron stars born in a core-collapsed supernova. Pulsars are highly
magnetized and surrounded by a plasma-filled magnetosphere emitting regular
pulses of radiation at their spin frequency. Due to the magneto-dipole losses,
they lose energy as their spin period increases. Thanks to radio surveys,
performed for instance by the Parkes and Arecibo radio-telescopes, nearly 3000
pulsars have been discovered (Manchester et al., 2005) (see Table 1, with the
data on the canonical pulsars). Although pulsars were first observed in radio,
it was found later that they are also bright in X-rays, optical and gamma-
rays. The Large Area Telescope (LAT) on board of the Fermi satellite, has
discovered dozens of radio-quiet gamma-ray pulsars as well as millisecond
pulsars (MSPs). Since its launch, the LAT has detected about 300 pulsars (Abdo
et al., 2013; Smith et al., 2019, 2023). The LAT has unveiled a new
perspective, expanding our understanding about pulsars by presenting a growing
collection of neutron stars that are detectable through their high energy
emission. Moreover, the distinct emission process in high-energy and its
unique beaming properties offer an alternative viewpoint on the global pulsar
population.

Table 1: Number of known pulsars with spindown luminosity $\dot{E}$ above and below $10^{28}$ W, and above $10^{31}$ W. 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 log($\dot{E}$) (W) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
¿31 | 2 | 0 | 0 | 2
¿28 | 197 | 101 | 35 | 61
Total | 2665 | 2553 | 63 | 84

With the advent of future surveys such as the Square Kilometre Array (SKA
111https://www.skao.int/en/science-users/118/ska-telescope-specifications) or
the Cherenkov Telescope Array (CTA 222https://www.cta-
observatory.org/about/how-ctao-works/), much more pulsars are expected to be
detected. Those instruments are scheduled to start collecting data by the end
of the decade. SKA is going to be 50 times more sensitive and is predicted to
survey the sky 10 000 times faster than any existing imaging radio telescope
(McMullin et al., 2020). CTA is aimed at detecting gamma rays in the energy
range from a few tens of GeV up to hundreds of TeV (Hofmann & Zanin, 2023)
while Fermi/LAT operates in the energy range from 20 MeV to several hundred
GeV 333https://fermi.gsfc.nasa.gov/science/instruments/table1-1.html.

Pulsar population synthesis (PPS) is a powerful tool aiming at predicting the
discovery rate of new pulsars as well as better understanding their emission
processes. PPS studies also constrains the overall properties of pulsars.
Pulsars are generated from their birth and evolved up to the present. Once the
set of pulsars has been generated we check whether they fulfill some
prescribed detection criteria. The detectability of an individual pulsar in a
given energy band (radio or $\gamma$) depends on whether the associated
emission beam intersects our line of sight and on the sensitivity of the
corresponding instrument.

Most of PPS studies assume that neutron stars are rotating in vacuum (Faucher-
Giguère & Kaspi, 2006; Popov et al., 2010; Johnston et al., 2020), or take
only radio or gamma-ray emission into account (Watters & Romani, 2011; Gullon
et al., 2014), or assume a constant magnetic field (Gonthier et al., 2002;
Faucher-Giguère & Kaspi, 2006; Johnston & Karastergiou, 2017). For the first
time in a PPS study, Dirson et al. (2022) have modeled the population of both
the gamma-ray and radio pulsars by taking into account the force-free
magnetosphere model in conjunction with the magnetic field decay. The neutron
star radiation is produced by relativistic charged particles flowing within
their magnetosphere. Therefore it is necessary to take the plasma back
reaction into account for the pulsar period and dipole magnetic inclination
angle evolution.Whereas the gamma-ray emission site from pulsars is still
under debate, different locations have been suggested. These are for instance
the polar cap region (Sturrock, 1971; Ruderman & Sutherland, 1975; Daugherty &
Harding, 1982, 1996) at low altitudes, the slot gap along the last open
magnetic field lines (Arons, 1983; Muslimov & Harding, 2004; Harding et al.,
2008; Harding & Muslimov, 2011), or the outer gaps at high altitudes in the
outer magnetosphere (Cheng et al., 1986; Hirotani, 2008; Takata et al., 2011).
The striped wind model, which has for emission site the pulsar wind outside
the light cylinder, (Pétri, 2009, 2011) is used for the first time by Dirson
et al. (2022) in a PPS study. The radio emission is described as usual by the
polar cap model.

In this work, and the work of Dirson et al. (2022) only the canonical pulsar
population in our Galaxy is reproduced. The canonical pulsar population is
defined by isolated pulsars which are not magnetars, and with a spin period
greater than 20 ms.

While the study yielded satisfactory results when simulating a population
closely resembling the observed pulsar population in the $P-\dot{P}$ diagram,
there was room for improvement regarding the pulsar motion. Specifically, the
study is limited by the assumption that pulsars move in straight lines at a
constant speed from their birth to the present time. A noticeable discrepancy
arises when comparing the observed position of pulsars to the simulation
expectations, revealing a significant disparity in their spatial
distributions.

Since the results obtained by Dirson et al. (2022) did not require the
modeling of a death line, it was not included in their PPS. The death line
separates in the $P-\dot{P}$ diagram the neutron stars that generate particles
through pair cascading from those which cannot radiate anymore and are hence
called dead pulsars. This contrasts with the PPS performed by Graber et al.
(2023) who did not require a death line. Therefore in this paper we examine
the necessity of a death line to explain the observed canonical pulsars
population.

This paper improves the PPS of Dirson et al. (2022), questions the necessity
of a death line in PPS work and discusses in depth the $\gamma$-ray pulsar
population. The paper is organised as follows. In Sect. 2 we detail our model,
recalling the processes for generation and evolution of the pulsars in the
PPS, a description of the Galactic potential, the death line and the multi-
wavelength detection criteria. Then the results are shown in Sect. 3. The
paper ends with a discussion in Sect. 4 and a summary in Sect. 5.

##  2 Description of the PPS model

Our PPS study is based on a self-consistent state-of-the-art evolution of the
pulsar geometry and proper motion within the Milky Way. In this section we
detail the model by exposing first the generation of individual pulsars, their
initial position, magnetic field and spin period, then how these quantities
should evolve in time, giving explicit expressions for the Galactic potential
used in our simulations, before introducing the death line and finally the
radio and $\gamma$-ray detection criteria.

###  2.1 Generation of pulsars

We generate a population of pulsars whose ages are generated in ascending
order, every X years for example if the birth rate is 1/X yr-1 between 0 yr
and the age of the Milky Way at maximum. The age of the Milky Way is almost
equal to the age of the Universe, which is of about $t_{H}\approx 13.8\times
10^{9}$ yr, therefore, a total number of $t_{H}/X$ pulsars should be
simulated. However as we will see, we chose to emulate $10^{7}$ pulsars in a
first run and $10^{6}$ in another run. This choice was made for two reasons :
firstly, the computation time is significantly decreased while pulsars older
than $10^{8}$ yr are mostly not detected, secondly, the simulation that
emulates $10^{7}$ pulsars allows to emulate older pulsars and require a death
line while the other works well without. This approach was also chosen because
it allows to check whether our birth rate is relevant by comparing our results
from the simulation to the observations. It is different from other works like
Johnston & Karastergiou (2017) and Gullon et al. (2014), where they stop the
simulation when the number of pulsars detected in the simulation are the same
as the number of observed pulsars. The value chosen for the birth rate, which
is 1/41 yr-1, is taken as best educated guess from the interval [1/150,1/33]
yr-1 because these studies (Faucher-Giguère & Kaspi, 2006; Gullon et al.,
2014; Johnston & Karastergiou, 2017) evaluated the birth rate of pulsars in
this interval.

To describe the position of the pulsars in the Galaxy, we use the right-handed
Galacto-centric coordinate system $(x,y,z)$ with the Galactic centre at its
origin, $y$ increasing in the disc plane towards the location of the Sun, and
$z$ increasing towards the direction of the north Galactic pole. The initial
spatial distribution of the pulsars is given by Paczynski (1990). He found two
distributions, one describing the radial spread and one the spread in
altitude, however only his distribution in altitude is used in our work
because the radial distribution in the galactic plane is better depict in
Ahlers et al. (2016) by the Milky Way’s pulsar surface density defined as

$\displaystyle\rho(R)$
$\displaystyle=A\left(\frac{R+R_{1}}{R_{\odot}+R_{1}}\right)^{a}\exp\left(-b\left(\frac{R-R_{\odot}}{R_{\odot}+R_{1}}\right)\right)$
(1a) $\displaystyle\rho_{z}(z)$ $\displaystyle=\frac{e^{-|z|/h_{c}}}{h_{c}}$
(1b)

where $R$ is the axial distance from the $z$-axis, and $z$ is the distance
from the Galactic disc. The numerical values for the constants are, $A=37.6$
kpc-2, $a=1.93$, $b=5.06$, $R_{1}=0.55$ kpc, $h_{c}=180$ pc and
$R_{\odot}=8.5$ kpc and in agreement with the distribution of young massive
stars in our galaxy (Li et al., 2019). The purpose of choosing equation (1a)
to describe the radial spread is to put neutron stars birth positions within
the spiral arms of the Galaxy. More precisely the Galactic spiral structure
contains four arms with a logarithmic shape function allowing to obtain the
azimuthal coordinate $\phi$ as a function of the distance from the Galactic
center

$\phi(R)=k\ln\left(\frac{R}{r_{0}}\right)+\phi_{0}$ (2)

The values of the models describing each arm are given in Table 2. These
values were used by Ronchi et al. (2021) and taken from Yao et al. (2017) to
match the shape of the arms of the Galaxy. The local arm is not modeled,
because of its very low density, much smaller than the four other arms. Each
star has an equal probability to be in one of the four arms, its angular
coordinate $\phi$ for a given $R$ being deduced from (1a). Furthermore, the
Galaxy is not static, and its arms are moving with an approximated period
$T=250$ Myr (Skowron et al., 2019). Since we also know that the Galaxy is
rotating in the clockwise direction, by knowing the age of a pulsar we can
know its angular position at birth. Following the same procedures as Ronchi et
al. (2021), in order to avoid artificial features near the Galactic center,
noise is added to both coordinates R and $\phi$. For instance
$\phi_{corr}=\phi_{rand}\,\exp(-0.35\,R)$ where $\phi_{rand}$ is randomly
drawn from a uniform distribution between $0$ and $2\,\pi$. $r_{corr}$ will be
taken from a normal distribution with a mean of $0$ and a standard deviation
$\sigma_{corr}=0.07\,R$ and we add this two corrections to $\phi$ and $R$
respectively. Therefore $R_{birth}=R+r_{corr}$ and
$\phi_{birth}=\phi+\phi_{corr}+\frac{2\,\pi\,t_{age}}{T}$, where $t_{age}$ is
the age of the pulsar. As soon as $R_{birth}$ and $\phi_{birth}$ are known we
can convert these coordinates in $x$ and $y$ Galacto-centric coordinates for
each pulsar.

Table 2: Parameters of the Milky Way Spiral Arm structure. Arm Number | Name | $k$ | $r_{0}$ | $\phi_{0}$
---|---|---|---|---
|  | (rad) | (kpc) | (rad)
1 | Norma | 4.95 | 3.35 | 0.77
2 | Carina-Sagittarius | 5.46 | 3.56 | 3.82
3 | Perseus | 5.77 | 3.71 | 2.09
4 | Crux-Scutum | 5.37 | 3.67 | 5.76

Next, we need to describe the individual properties of each pulsar, namely,
its inclination angle $\alpha$, which is the angle between its rotation axis
and its magnetic axis, its initial spin period and initial magnetic field at
birth. The inclination angle $\alpha$ is assumed to follow an isotropic
distribution generated from a uniform distribution $U\in[0,1]$ and given by
$\alpha=\arccos(2\,U-1)$. In our work, the initial spin period and magnetic
field follow both a log-normal distribution, as suggested from the results of
a study of 56 young neutron stars by Igoshev et al. (2022) while other works
simply use a Gaussian distribution (Faucher-Giguère & Kaspi, 2006; Gullon et
al., 2014; Johnston et al., 2020). Explicitly, the probabilities to born with
an initial magnetic field $B_{0}$ and period $P_{0}$ are given by

$\displaystyle p(\log(B_{0}))=\frac{1}{\sigma_{b}\sqrt{2\pi}}e^{-(\log
B_{0}-\log\bar{B})^{2}/(2\sigma_{b}^{2})}$ (3) $\displaystyle
p(\log(P_{0}))=\frac{1}{\sigma_{p}\sqrt{2\pi}}e^{-(\log
P_{0}-\log\bar{P})^{2}/(2\sigma_{p}^{2})}.$ (4)

The distribution for the magnetic field is identical to the one used in other
pulsar population synthesis like for instance Faucher-Giguère & Kaspi (2006),
Gullon et al. (2014), Johnston et al. (2020), Yadigaroglu & Romani (1995) and
Watters & Romani (2011).

Compared to the study of Dirson et al. (2022) we introduce a novelty, namely
the alignment between the kick velocities at birth and the rotation axis of
the pulsar, as suggested by Rankin (2007). Moreover according to Hobbs et al.
(2005) a Maxwellian distribution best replicates the observations, with a
standard deviation $\sigma_{\rm v}=265$ km/s. The mean velocity of the
distribution is related to the standard deviation by $\bar{v}=\sigma_{\rm
v}\sqrt{{8}/{\pi}}$. The velocity vector is then distributed along the unit
vector of the rotation axis generated as detailed in the subsection 2.5.

###  2.2 Pulsar evolution

The pulsar initial period, inclination angle and magnetic field must be
evolved in time self-consistently due to electromagnetic radiation and
internal dissipation within the neutron star. To this end we employ the pulsar
force-free magnetosphere model taking into account the plasma current and
charge within the magnetosphere. The spin down luminosity $\dot{E}$ is

$\displaystyle\dot{E}=\frac{dE_{\rm rot}}{dt}=-I\Omega\dot{\Omega}=L_{\rm
ffe}$ (5) $\displaystyle L_{\rm ffe}=\frac{4\pi
R^{6}B^{2}\Omega^{4}(1+\sin^{2}\alpha)}{\mu_{0}c^{3}}$ (6)

where $I$ is the neutron star moment of inertia its value being approximately
$I\approx 10^{38}$ kg.m2, $\Omega=2\pi/P$ is the rotation frequency of the
pulsar ($P$ is the spin period) and $\dot{\Omega}$ is its time derivative,
$R=12$ km is the typical radius of a neutron star as found by recent NICER
observations (Riley et al., 2019; Bogdanov et al., 2019), $B$ is the magnetic
field at the magnetic equator, $\alpha$ is the inclination angle, $\mu_{0}$ is
the vacuum permeability constant, $\mu_{0}=4\pi\times 10^{-7}$ H/m and $c$ is
the speed of light. The term $L_{\rm ffe}$ for the spin down luminosity was
given by Spitkovsky (2006) and Pétri (2012). Equation (5) clearly shows the
correlation between the obliquity, the magnetic field and the rotation
frequency. Combining equation (5) and (6) leads to

$\displaystyle\dot{\Omega}=-K_{\rm ffe}\ \Omega^{n}$ (7) $\displaystyle K_{\rm
ffe}=\frac{4\pi R^{6}B^{2}(1+\sin^{2}\alpha)}{\mu_{0}c^{3}I}$ (8)

which is written in a more concise form where $n$ is the braking index, its
value being $n=3$ for magnetic dipole radiation.

The integral of motion between $\Omega$ and $\alpha$ is also very useful and
reads

$\Omega\frac{\cos^{2}\alpha}{\sin\alpha}=\Omega_{0}\frac{\cos^{2}\alpha_{0}}{\sin\alpha_{0}}$
(9)

where quantities with index 0 indicate their initial value. The quantities
without the index are their current time value.

As shown by Gullon et al. (2014) a decaying magnetic field better accounts for
comparison between pulsar population synthesis and observations. Therefore, we
prescribe a magnetic field decay according to a power law

$B(t)=B_{0}(1+t/\tau_{d})^{-1/\alpha_{d}}$ (10)

where $\alpha_{d}$ is a constant parameter controlling the speed of the
magnetic field decay and $\tau_{d}=k_{\tau_{d}}\tau_{v}$ is the magnetic field
decay time scale. The decay rate $\tau_{v}$ was chosen randomly between 3
values : $5\times 10^{4}$ yr or $7\times 10^{4}$ yr or $5\times 10^{5}$ yr,
with a probability of 0.23, 0.46 and 0.31 respectively (these probabilities
were chosen as best educated guesses after several run with different values)
for each pulsar, corresponding to the $\tau_{v}$ value for canonical pulsars
in Viganò et al. (2013). The values of $\alpha_{d}$ and $k_{\tau_{d}}$ can be
found in Table 5 of Sect. 3.

In line with the spin down luminosity, the inclination angle satisfies another
evolution equation that after integration was found by Phillipov et al. (2014)
for a spherically symmetric neutron star with a constant magnetic field. With
our decaying prescription the inclination angle $\alpha$ is found by solving
for the root of

$\displaystyle\begin{split}\ln(\sin\alpha_{0})+\frac{1}{2\sin^{2}\alpha_{0}}+K\Omega_{0}^{2}\frac{\cos^{4}\alpha_{0}}{\sin^{2}\alpha_{0}}\frac{\alpha_{d}\tau_{d}B_{0}^{2}}{\alpha_{d}-2}\left[\left(1+\frac{t}{\tau_{d}}\right)^{1-2/\alpha_{d}}-1\right]\\\
\phantom{{}={}}=\ln(\sin\alpha)+\frac{1}{2\sin^{2}\alpha}\end{split}$ (11)

where $t$ is the time representing the age of the pulsar and
$K={R^{6}}/{I\,c^{3}}$. The typical decay timescale for a mean magnetic field
of $2.5\times 10^{8}$ T is $4.6\times 10^{5}$ yr in Viganò (2013), however in
this study the typical decay timescale associated to this magnetic field is
$3.5\times 10^{5}$ yr which is consistent with the estimate of Viganò (2013)
for the same field strength.

Finally, we need to follow the particle motion within the Galactic potential.
Each generated pulsar evolves in the gravitational potential $\Phi$ subject to
an acceleration $\@vec{\ddot{x}}$ according to

$\@vec{\ddot{x}}=-\@vec{\nabla}\Phi\ .$ (12)

Equation (12) is integrated numerically thanks to a Position Extended Forest
Ruth-Like (PEFRL) algorithm (Peter Young, 2014), a fourth order integration
scheme, presented in Appendix A with a convergence and accuracy study. The
galactic potential model used in discussed in the next section.

###  2.3 Description of the Galactic potential

The Galaxy is divided in four distinct regions with different mass
distributions and associated gravitational potentials. The four potentials
are: the bulge $\Phi_{b}$, the disk $\Phi_{d}$, the dark matter halo
$\Phi_{h}$ and the nucleus $\Phi_{n}$. The total potential of the Milky Way
$\Phi_{tot}$ is the sum of these potentials

$\Phi_{tot}=\Phi_{b}+\Phi_{d}+\Phi_{h}+\Phi_{n}.$ (13)

The expressions for these potentials are taken from Bajkova & Bobylev (2021)
with parameters given in Table 3. The nucleus mass was found in Bovy (2015).

Table 3: Constants values used for the different potentials, with the solar mass $M_{\odot}=1.99\times 10^{30}$ kg. Parameters | Values
---|---
$M_{h}$ | $2.9\times 10^{11}\pm 7.6\times 10^{10}M_{\odot}$
$a_{h}$ | 7.7 $\pm$ 2.1 kpc
$a_{d}$ | 4.4 $\pm$ 0.73 kpc
$b_{d}$ | 0.308 $\pm$ 0.005 kpc
$M_{d}$ | $6.50\times 10^{10}\pm 1.9\times 10^{9}M_{\odot}$
$a_{b}$ | 0.0 kpc
$b_{b}$ | 0.267 $\pm$ 0.009 kpc
$M_{b}$ | $1.02\times 10^{10}\pm 6.3\times 10^{8}M_{\odot}$
$M_{n}$ | $4.10^{6}\pm 0.42\times 10^{6}M_{\odot}$

The potential for the bulge $\Phi_{b}$ and for the disk $\Phi_{d}$ were both
chosen to have the form proposed by Miyamoto & Nagai (1975) which are
typically used for models of the gravitational potential of the Milky Way.
They read

$\Phi_{i}(R,z)=-\frac{GM_{i}}{\left[R^{2}+\left(a_{i}+\sqrt{z^{2}+b_{i}^{2}}\right)^{2}\right]^{1/2}}$
(14)

where $R^{2}=x^{2}+y^{2}$, and $i=b$ is for the bulge and $i=d$ is for the
disk. In those formula $R$ depends on the coordinates $x$ and $y$, $a_{i}$ and
$b_{i}$ are the scale parameters of the components in kpc, $M_{i}$ is the mass
(for the disk or the bulge) and $G$ is the gravitational constant. The
constant values can be found in Table 3.

Solving the pulsar equation of motion requires to compute the gradient of the
potentials, the derivatives of the disk and the bulge potentials are

$\frac{\partial\Phi_{i}}{\partial
x_{j}}=\frac{GM_{i}x_{j}}{\left[R^{2}+\left(a_{i}+\sqrt{z^{2}+b_{i}^{2}}\right)^{2}\right]^{3/2}}$
(15) $\frac{\partial\Phi_{i}}{\partial
z}=\frac{GM_{i}z\left(a_{i}+\left(z^{2}+b_{i}^{2}\right)^{\frac{1}{2}}\right)}{\left[R^{2}+\left(a_{i}+\sqrt{z^{2}+b_{i}^{2}}\right)^{2}\right]^{3/2}\
\left(z^{2}+b_{i}^{2}\right)^{1/2}}$ (16)

$x_{j}$ being either $x$ or $y$. Equation (15) is the derivative for the
potential of Miyamoto Nagai, where $i=d$ or $i=b$ depending on disk or bulge
derivative. Equation (16) is the derivative of the Miyamoto Nagai potential
with respect to the $z$ coordinate.

The potential for the dark matter halo $\Phi_{h}$, according to the frequently
used potential of Navarro et al. (1997), it is expressed as

$\Phi_{h}(r)=-\frac{GM_{h}}{r}\ \ln\left(1+\frac{r}{a_{h}}\right)$ (17)

where $M_{h}$ is the mass of the halo, $r=x^{2}+y^{2}+z^{2}$, and $a_{h}$ is
the length scale whose values are found in Table 3. The derivative of equation
(17) with respect to any coordinate $x_{j}$ $(x,y,z)$ is

$\frac{\partial\Phi_{h}}{\partial x_{j}}=\frac{GM_{h}x_{j}}{r^{2}}\
\left[\frac{1}{r}\
\ln\left(1+\frac{r}{a_{h}}\right)-\frac{1}{r+a_{h}}\right].$ (18)

Finally for the last part of the Galaxy, the nucleus of the Milky Way is
simply represented by a Keplerian potential $\Phi_{n}$ such as

$\Phi_{n}(r)=-\frac{GM_{n}}{r}$ (19)

where $M_{n}$ is the mass of the nucleus. The derivative of equation (19) with
respect to any coordinate $x_{j}\ (x,y,z)$ is given by

$\frac{\partial\Phi_{n}}{\partial x_{j}}=\frac{GM_{n}\,x_{j}}{r^{3}}.$ (20)

###  2.4 The death line

Several equations for the death line have been proposed in the literature
(Chen & Ruderman, 1993; Zhang et al., 2000; Gil & Mitra, 2001; Mitra et al.,
2019). The most suitable one in our study is that of Mitra et al. (2019) where
they found an expression for the death line

$\dot{P}_{\rm line}=\frac{3.16\times 10^{-19}\ T_{6}^{4}\
P^{2}}{\eta^{2}b\cos^{2}\alpha_{l}}\ .$ (21)

Here $T_{6}=T/10^{6}$ K, where $T$ corresponds to the surface temperature of
the polar cap. The parameter $\eta=1-\rho_{i}/\rho_{GJ}$ is the electric
potential screening factor due to the ion flow, where $\rho_{i}$ correspond to
the ion charge density and $\rho_{GJ}$ to the Goldreich-Julian charge density
above the polar cap. The quantity $b$ is the ratio of the actual surface
magnetic field to the dipolar surface magnetic field and $\alpha_{l}$ is the
angle between the local magnetic field and the rotation axis. The parameters
we used for this death line are : $\alpha_{l}=45\degr$, $b=40$, $\eta=0.15$,
$T_{6}=2$. Spread in the parameter values of the model causes significant
variations in the death line, and thus a death valley rather than a single
death line describes the condition for pulsar extinction in the $P-\dot{P}$
diagram.

To decide whether a pulsar is still active in radio, its critical spin down
rate $\dot{P}_{\rm line}$ is computed knowing its spin period $P$. If the
pulsar lies in the death valley defined by the inequality
$10^{-0.55}<\dot{P}_{\rm line}/\dot{P}<10^{1.15}$ then a value for
$\alpha_{l}$ is taken from a uniform distribution between $0\degr$ and
$65\degr$. Moreover a value for $T_{6}$ is also taken from a uniform
distribution between 1.9 and 2.8. The parameter $b$ is taken from a uniform
distribution between 30 and 60. The boundary parameter values for the death
valley were chosen to traverse the $P-\dot{P}$ diagram from the region in the
bottom right quadrant where the density starts to decrease, to the detectable
pulsar with the lowest $\dot{P}$ and highest $P$ in the same quadrant, as
shown in Fig. 1. The edges of the death valley given by the values
$10^{-0.55}$ and $10^{1.15}$ were found by taking the extreme values for the
parameters that are drawn from the different uniform distributions. Then we
compare its actual spin down rate $\dot{P}$, computed from the evolution
model, to the $\dot{P}^{\rm new}_{\rm line}$ computed with equation (21) with
the parameters which were drawn from uniform distributions. For pulsars that
are not in the death valley, their $\dot{P}$ is only compared with the spin
down rate $\dot{P}_{\rm line}$ obtained with the default parameters. If the
pulsar has its $\dot{P}$, which was computed from the evolution model, below
the one computed with the death line ($\dot{P}_{\rm line}$ or $\dot{P}^{\rm
new}_{\rm line}$), then the pulsar is considered as dead, otherwise it is
considered as still emitting photons. These choices were made in order to have
more like a death valley near the death line than a strict condition of death,
because the parameters $T_{6}$, $\alpha_{l}$ and $b$ may be different from one
pulsar to another.

Figure 1: $P-\dot{P}$ diagram of the canonical pulsars along with the death
line, green solid line, and death valley, shaded green area.

###  2.5 Detection

For each generated pulsar, we check whether it fulfills the criteria for
detection depending on three factors. First, the beaming fraction indicates
the fraction of the sky covered by the radiation beam and depends on the
considered wavelength, here radio or $\gamma$-ray. The beaming fraction varies
also with the pulsar spin rate, geometry and location of the emission regions.
Before explaining how the beaming fraction in radio and in $\gamma$-ray are
computed, it is important to define several angles.

The angle between the line of sight and the rotation axis is denoted by
$\xi=(\widehat{\@vec{n_{\Omega}},\@vec{n}})$ where
$\@vec{n}_{\Omega}={\@vec{\Omega}}/{\|\@vec{\Omega}\|}$ is a unit vector along
the rotation axis and $\@vec{n}$ the unit vector along the line of sight. The
inclination angle $\alpha=(\widehat{\@vec{n_{\Omega}},\@vec{\mu}})$ is the
angle between the rotation axis and the magnetic moment, $\@vec{\mu}$ being
the unit vector along the magnetic moment. Traditionally the impact angle is
also introduced as the angle between the magnetic moment and the line of sight
$\beta=(\widehat{\@vec{\mu},\@vec{n}})$. Moreover, it is related to the
previous angles by $\alpha+\beta=\xi$.

We choose a isotropic distribution for the Earth viewing angle $\xi$ as well
as for the orientation of the unitary rotation vector. The Cartesian
coordinates of the unit rotation vector $\@vec{n_{\Omega}}$ are
$(\sin\theta_{n_{\Omega}}\cos\phi_{n_{\Omega}},\sin\theta_{n_{\Omega}}\sin\phi_{n_{\Omega}},\cos\theta_{n_{\Omega}})$.
We set the Sun position at $(x_{\odot},y_{\odot},z_{\odot})=$ (0 kpc, 8.5 kpc,
15 pc) (Siegert, 2019). The coordinates for $\@vec{n}$ are

$\@vec{n}=\left(\frac{x-x_{\odot}}{d},\frac{y-y_{\odot}}{d},\frac{z-z_{\odot}}{d}\right).$
(22)

To compute the pulsar distance from Earth, we use the formula for the distance
which is

$d=\sqrt{(x-x_{\odot})^{2}+(y-y_{\odot})^{2}+(z-z_{\odot})^{2}}.$ (23)

####  2.5.1 Radio detection model

The beaming fraction in radio depends on the half opening angle of the radio
emission cone $\rho$ computed according to

$\rho=3\sqrt{\frac{\pi\,h_{\rm em}}{2Pc}}$ (24)

where $h_{\rm em}$ is the emission height, $P$ is the spin period of the
pulsar and $c$ is the speed of light. The emission height is taken constant
with an average value of $h_{\rm em}=3.10^{5}$ m, estimated from observations
of a sample of pulsars by (Weltevrede & Johnston, 2008; Mitra, 2017; Johnston
& Karastergiou, 2019; Johnston et al., 2023). The cone half-opening angle
$\rho$ estimated in eq. (24) holds only for the last open field lines of a
magnetic dipole and can only be applied for slow pulsars where the radio
emission altitude is high enough for the multipolar components to decrease
significantly and become negligible.

The pulsar is detected in radio if $\beta=|\xi-\alpha|\leq\rho$ or if
$\beta=|\xi-(\pi-\alpha)|\leq\rho$ corresponding to the north and south
hemisphere respectively. It must also satisfy the condition $\alpha\geq\rho$
and $\alpha\leq\pi-\rho$ in order to effectively see radio pulsation, because
the line of sight must cross the emission cone to observe pulsations. Another
useful quantity is the observed width of the radio profile $w_{r}$ which is
computed as (see Lorimer & Kramer (2004))

$\cos(\rho)=\cos(\alpha)\cos(\xi)+\sin(\alpha)\sin(\xi)\cos(w_{r}/2).$ (25)

The second factor for detection is the luminosity, which is also different
between radio and gamma-rays. First, concerning the radio flux density, the
formula used is the same as in Johnston et al. (2020) for a pulsar at 1.4 GHz,
in order to model the detection carried out by the Parkes radio telescope in
the Southern Galactic plane (Kramer et al., 2003; Lorimer et al., 2006;
Cameron et al., 2020) and the Arecibo telescope in the northern plane (Cordes
et al., 2006)

$F_{r}=9\ \text{mJy}\left(\frac{d}{1\
\text{kpc}}\right)^{-2}\left(\frac{\dot{E}}{10^{29}\
\text{W}}\right)^{1/4}\times 10^{F_{\rm j}}$ (26)

where $d$ is the distance in kpc and $F_{\rm j}$ is the scatter term which is
modeled as a Gaussian with a mean of 0.0 and a variance of $\sigma$ = 0.2. The
detection threshold in radio is set by the signal to noise ratio defined by

$S/N=\frac{F_{r}}{S_{\rm survey}^{\rm min}}.$ (27)

The pulsar is detected if the signal to noise ratio $S/N$ is greater or equal
to 10. $S_{\rm survey}^{\rm min}$ is the minimum flux which is related to the
period of rotation $P$, the width profile of radio emission $w_{r}$ and the
sensitivity of the survey which is the last factor for detection, and will be
detailed later. We directly compute the radio flux, without computing the
luminosity in radio which overlooks the fact that the luminosity received will
depend on the geometry of the beam.

The third and last factor is the sensitivity, depending on the survey,
therefore on the instrument used and on the pulse profile observed. A pulsar
is detected in radio if the signal to noise ratio is greater than 10, but this
ratio depends on the minimum flux, another function of the instrumental
sensitivity. With the aim of computing the observed width pulse profile, we
use the same formula as Cordes & McLaughlin (2003)

$\displaystyle\tilde{w}_{r}=\sqrt{\left(w_{r}\frac{P}{2\pi}\right)^{2}+\tau_{samp}^{2}+\tau_{DM}^{2}+\tau_{scat}^{2}}$
(28a) $\displaystyle
S_{\text{survey}}^{\text{min}}=S_{0}\sqrt{\frac{\tilde{w}_{r}}{P-\tilde{w}_{r}}}$
(28b)

Equation (28a) takes into account density fluctuations in the interstellar
medium (ISM), the dispersion and scattering during the propagation of the
radio pulse when interacting with the free electrons. Moreover, the
instrumental effect is also taken into account with $\tau_{samp}$, which is
the sampling time of the instrument. All this mechanisms lead to a broader
pulse profile. In order to determine the influence of the ISM through
$\tau_{DM}$, the formula of Bates et al. (2014) is used

$\tau_{DM}=\frac{e^{2}\Delta f_{ch}\,DM}{\pi\,m_{e}\,c\,f^{3}}$ (29)

$e$ is the electron charge, $m_{e}$ its mass, $\Delta f_{ch}$ the width of
frequency of the instrument channel, $f$ the observing frequency and $DM$ the
dispersion measure. Furthermore, in order to determine the influence of
scattering by an inhomogeneous and turbulent ISM through $\tau_{scat}$, the
empirical fit relationship from Krishnakumar et al. (2015) is used

$\tau_{scat}=3.6\times 10^{-9}DM^{2.2}(1+1.94\times 10^{-3}DM^{2})\ .$ (30)

Both $\tau_{DM}$ and $\tau_{scat}$ depend on the dispersion measure, which is
found for each pulsar by running the code of Yao et al. (2017) converting the
distance of a pulsar into the dispersion measure thanks to their state-of-the
art model of the Galactic electron density distribution. The radio survey
parameters used are the ones from the Parkes Multibeam Pulsar Survey (PMPS),
see Table 4.

Table 4: Survey parameters of the Parks Multibeam Pulsar Survey (PMPS). $\Delta f_{ch}$ | $f$ | $\tau_{samp}$
---|---|---
(kHz) | (GHz) | (µs)
3000 | 1.374 | 250

Concerning equation (28b), $S_{0}$ represents the survey parameters. Johnston
et al. (2020) estimated that $S_{0}$ should approximately be equal to 0.05 mJy
to have a signal to noise ratio greater than 10 for normal pulsars in this
survey. We chose equation (28b) to compute the minimum flux, to decide whether
a pulsar will be detected in radio or not, however it does not reproduce the
whole complexity of this quantity which could be computed more precisely, (see
equation (24) of Faucher-Giguère & Kaspi (2006) for instance, where $S_{0}$ is
computed more precisely) even though it remains a good approximation. In
addition we only use the parameters of PMPS to compute $\tau_{DM}$, we could
have also used the parameters of The Pulsar Arecibo L-band Feed Array (PALFA)
Survey, however the differences in the results in the end are not large if we
use these parameters, consequently, it persists as a favorable approximation
to exclusively utilize the parameters delineated by PMPS.

####  2.5.2 Gamma-ray detection model

The $\gamma$-ray emission model relies on the striped wind model, describing
gamma-ray photons production emanating from the current sheet within the
striped wind. In order to detect the gamma-ray, the line of sight of the
observer must remain around the equator plane with an inclination angle
constrained by $|\xi-\pi/2|\leq\alpha$. The gamma-ray luminosity is extracted
from a study of Kalapotharakos et al. (2019), where they showed that the
luminosity is described by a fundamental plane. The 3D model depends on the
magnetic field $B$, the spin down luminosity $\dot{E}$ and the cut-off energy
$\epsilon_{cut}$. However, in our PPS the cut-off energy is not computed
therefore we use their two dimensional version

$L_{\gamma(2D)}=10^{26.15\pm 2.6}\ 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}$ (31)

The spindown $\dot{E}$ is computed with equation (5) and the associated
$\gamma$-ray flux detected on earth is computed with

$F_{\gamma}=\frac{L_{\gamma(2D)}}{4\,\pi\,f_{\Omega}\,d^{2}}$ (32)

where $f_{\Omega}$ is a factor depending on the emission model reflecting the
anisotropy. For the striped wind model Pétri (2011) showed that this factor
varies between 0.22 and 1.90. Nevertheless, an approximation is made : if
$\alpha<-\xi+0.6109$ then $f_{\Omega}=1.9$ otherwise $f_{\Omega}=1$. The
pulsar is detected in gamma depending on the instrumental sensitivity as
described below.

The sensitivity in gamma-ray is based on the expectation of the Fermi/LAT
instrument
444https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone
/Cicerone_LAT_IRFs/LAT_sensitivity.html. Two conditions are required to count
the detection, first we need a sufficient value for the flux. If the galactic
latitudes of the pulsar is ¡ 2°, then $F_{\rm min}=4\times 10^{-15}$ W.m-2,
and if blind searches are assumed we set $F_{\rm min}=16\times 10^{-15}$
W.m-2. Hence, if the gamma-ray flux $F_{\gamma}$ is greater than $F_{\rm
min}$, the pulsar is detected in gamma-rays if the second condition is
fulfilled. Moreover, in order to take the ISM into account, we set a threshold
on the dispersion measure DM (a similar approach was done in Gonthier et al.
(2018) for millisecond gamma-ray pulsars) because gamma-ray canonical pulsars
typically originate from supernova explosions and are associated with regions
of recent star formation, such as supernova remnants or star-forming regions.
These environments can have higher densities of free electrons in the
interstellar medium (ISM), leading to higher DM values for pulsars in these
regions. Since gamma-rays pulsars are often young, they are likely still
remaining in these regions. Adjusting the DM threshold can help optimize the
balance between maximizing sensitivity to pulsar signals and minimizing
contamination from background sources. Therefore, the DM threshold chosen is
15 cm${}^{-3}.$pc, following the lowest DM in the ATNF catalogue for high-
energy pulsars.

##  3 Results

Let us now move to the PPS simulation results. First we show runs simulating
10 million pulsars with and without a death line, allowing to simulate old
pulsars (with the oldest aged $4.1\times 10^{8}$ yr). Then we show runs
simulating 1 million pulsars (with the oldest aged $4.1\times 10^{7}$ yr) with
and without considering a death line. We explore a possible spin-velocity
misalignment effect due to the galactic potential and the impact of the
initial spin period distribution, normal vs log-normal. Finally, the gamma-ray
simulations and detection are discussed in details.

###  3.1 Simulating 10 millions pulsars

####  3.1.1 Without a death line

The first simulations performed with the Galactic potential implemented
discarded the death line. The parameters used for these runs are given in
Table 5. The total number of pulsars is summarized in Table 6.

Table 5: Parameters used in the simulations. $\tau_{\text{birth}}$ (1/yr) |  $P_{\text{mean}}$ (ms) |  $B_{\text{mean}}$ (T) | $\sigma_{p}$ | $\sigma_{b}$ | $\alpha_{d}$ | $k_{\tau_{d}}$
---|---|---|---|---|---|---
41 | 129 | $2.75\times 10^{8}$ | 0.45 | 0.5 | 1.5 | 5
Table 6: Number of pulsars detected without the implementation of the death line. log($\dot{E}$) (W) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
¿31 | 11 | 0 | 8 | 3
¿28 | 161 | 9 | 75 | 77
Total | 2653 | 2336 | 139 | 178

The quantities $N_{r}$, $N_{g}$, and $N_{rg}$ are the number of radio-only,
gamma-only, and radio-loud gamma-ray pulsars, respectively, obtained from our
simulation.

Compared to the results obtained by Dirson et al. (2022), at very high spin
down luminosity, $\dot{E}>10^{31}$ W, the detection is almost identical but at
$\dot{E}>10^{28}$ W there are much less pulsars detected in radio and
radio/gamma in this work. There are probably two causes for this discrepancy
in the number of detected pulsars at $\dot{E}>10^{28}$ W : firstly the
difference in the spin period distribution at birth. Pulsars have a higher
probability to start with a higher $P_{0}$ since the mean of the distribution
is higher (in the previous study the mean was 60 ms while it is 129 ms now),
therefore the slowing down of the pulsars can only decrease their rotation
frequency $\Omega$, and equation (5) shows that $\dot{E}$ depends on $\Omega$.
However, it does not affect the detection in gamma-rays because usually
canonical gamma pulsars have a high $\dot{P}$ and a low $P$, being young,
therefore their $\dot{E}$ is higher than for radio pulsars. Secondly, the ISM
dispersion and scattering reduces this number of radio pulsars detected (the
results without the ISM effect are not shown but without it there is much more
radio pulsars detected). Nonetheless, in terms of total number of detected
radio pulsars, the Galactic potential allows to have an increase compared to
the number found in Dirson et al. (2022), a result expected because of the
attractive nature of the potential of the Galaxy. In its previous version the
PPS moved the pulsars in random spatial directions at constant speed, as a
result no pulsar was bound to the Milky Way in closed orbits this is
especially true for old pulsars which could have time to leave the Galaxy.

Figure 2: Distribution of the observed period taken from the ATNF catalogue,
along with the simulations without the implementation of the death line. In
green, the period of the simulated pulsars with age greater than $10^{7.5}$ .

Fig. 2 shows the distribution of the spin period of the detected pulsars in
the simulation, in red and in observations, in blue. However the comparison is
not satisfactory because there is an excess of pulsars detected in the
simulation. Furthermore, this excess is also present in Fig. 3 where the age
of the detected pulsars from the observations and simulation are plotted.

In order to verify if there is a correlation between old pulsars and the ones
with a high spin period, we plotted in Fig. 2, in green the pulsars with an
age greater than $10^{7.5}$ yr, corresponding to the oldest pulsars of the
simulation and it corresponds to the excess found in Fig. 2. Due to the
Galactic potential, more pulsars remain bound to the Galaxy, especially older
pulsars which were escaping the Galaxy more easily in the work of Dirson et
al. (2022) because of the absence of potential. Fig. 4 shows that many pulsars
are in what is commonly called the graveyard of pulsars at the bottom right of
the $P-\dot{P}$ diagram.

Figure 3: Distribution of the observed age taken from the ATNF catalogue,
along with the simulations without the implementation of the death line.

Figure 4: $P-\dot{P}$ diagram of the simulated population, along with the
observations, without the implementation of the death line.

Given this excess of long period pulsars without a death line, we decided to
implement one, in order to stick closer to reality.

####  3.1.2 With the death line

Adding a death line decreases the number of pulsars detected in radio compared
to the situation without a death line, see Table 7.

Table 7: Number of pulsars detected with the death line implementation. log($\dot{E}$) (W) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
¿31 | 10 | 0 | 5 | 5
¿28 | 118 | 7 | 69 | 42
Total | 2077 | 1807 | 116 | 154

Figure 5: Distribution of the observed period taken from the ATNF catalogue,
along with the simulations with the death line implementation.

Figure 6: $P-\dot{P}$ diagram of the simulated population, along with the
observations, with the death line implementation.

Figure 7: Density plot of the $P-\dot{P}$ diagram in comparison with
observations

The death line improves the fit as seen in Fig. 5 for the distribution of spin
period and compared with Fig. 2, even though there is a slight excess of
pulsars which have spin period between $10^{-0.5}$ s and $1$ s. In Fig 6,
these pulsars in excess lie at $\dot{P}$ between $10^{-18}$ and $10^{-17}$,
where there is no data and which correspond to the oldest pulsars in the
simulation. Therefore the oldest pulsars simulated here do not correspond to
any real data.

In addition, there are too many pulsars close to the death line compared to
the observations, there is a pile up effect created by the death line which is
not in the observations.

The goodness-of-fit can be checked using a density plot of the $P-\dot{P}$
diagram as in Johnston & Karastergiou (2017) instead of putting points for
individual pulsars. In order to construct such diagrams we binned the data
evenly in a log scale in $P$ and $\dot{P}$ and count the number of pulsars in
each bin. This leads to a 2D histogram to be compared with the observed
density plot 2D histogram. Because the total number of simulated pulsars
$N_{totsim}$ is different from the total number of observed pulsars
$N_{totobs}$, we normalise the number of pulsars in each bin in order to
compare proportions. We therefore introduce the quantity $R$ as follows

$R=\frac{N_{sim}/N_{totsim}-N_{obs}/N_{totobs}}{N_{sim}/N_{totsim}+N_{obs}/N_{totobs}}.$
(33)

With this definition the number $R$ remains in the interval $[-1,+1]$.
$N_{sim}$ is the number of pulsar count in a bin of the simulation, $N_{obs}$
the number of pulsar detected in a bin of the observations. A good fit in a
bin correspond to a value of $R$ close to $0$, meaning that the simulation
results are close to the observations. Inspecting fig. 7, we note a stronger
deviation from the observations at the boundaries of the density plot, for
instance it confirms too many pulsars are piling up close to the death line.
The high values of $R$ at these boundaries originates from the weak number of
pulsars counted in this bins leading to bad statistics and large variations in
the proportions, a side effect of the density plot approach. The fact that the
colors are very dark, shows that the proportion of pulsars in most of the bins
are not good.

A Kolmogorov Smirnov (KS) test (Smirnov, 1948) was conducted on both $P$ and
$\dot{P}$ which gave p-values substantially below 0.05 respectively when no
death line is implemented. While we obtained p-values of 0.05 and another one
well below 0.05 respectively for $P$ and $\dot{P}$ with a death line
implementation. The null hypothesis that observations and simulations
originate from the same distribution for the spin period can not be rejected
(p-value $\geq$ 0.05) only when a death line is considered, but is clearly
rejected for $\dot{P}$ in both case.

###  3.2 Simulating 1 millions pulsars

####  3.2.1 With the death line

The parameters used for these runs are shown in Table 5. This new set of
simulations is intended to verify if generating pulsars younger than
$4.1\times 10^{7}$ yr is sufficient, since we observe a small excess of very
old pulsars. The death line is also taken into account.

Table 8: Number of pulsars detected with the death line implementation for 1 million pulsars simulated. log($\dot{E}$) (W) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
¿31 | 2 | 0 | 2 | 0
¿28 | 124 | 7 | 73 | 44
Total | 1694 | 1414 | 154 | 126

Figure 8: Distribution of the observed period taken from the ATNF catalogue,
along with the simulations with the death line implementation for 1 million
pulsars simulated.

Figure 9: Distribution of the observed period derivative taken from the ATNF
catalogue, along with the simulations with the death line implementation for 1
million pulsars simulated.

Figure 10: $P-\dot{P}$ diagram of the simulated population, along with the
observations, with the death line implementation for 1 million pulsars
simulated.

Table 8 shows that a few less pulsars were detected in this simulation
compared to the one with the death line and 10 millions pulsars (see Table 7).
Moreover, inspecting the histograms for spin period and spin period derivative
respectively in Fig. 8 and 9, we conclude that observations and simulations
are very similar in both histograms. A KS test was conducted for both
distributions $P$ and $\dot{P}$ and yielded p-values of $0.24$ and $0.07$
respectively. Therefore, it is impossible to reject the null hypothesis that
the simulated distributions are similar to the observed distribution.

Figure 11: Density plot of the $P-\dot{P}$ diagram in comparison with
observations for 1 million pulsars simulated.

Fig 10 shows the good agreement between simulations and observations.
Furthermore, we can compare now Fig. 11 and Fig. 7. It is obvious that Fig. 11
is better, because globally the colors are lighter which signifies the
proportion of pulsars in each 2D bin are almost alike. In addition, there is
not a pile up effect anymore next to the death line, which corresponds better
to the data.

Figure 12: Distance distribution to Earth for observed and simulated
populations (1 million pulsars), including the death line implementation.

The distribution of the distance to Earth is shown in Fig. 12. It has a
slightly better similarity to the observations than the work of Dirson et al.
(2022). Nevertheless the improvement is marginal compared to the improvements
of $P$ and $\dot{P}$ histograms. However, we noticed that removing the effect
of dispersion and scattering of the ISM on the radio pulse profile, we would
obtain very similar distribution of distances between observations and
simulations. Therefore, a refinement in the description of the effect of the
ISM on the pulse profile might be necessary.

Indeed, Fig. 13 confirms that the positions of the pulsars in the X-Y plane of
the Galaxy are slightly more scattered than before, but the improvement is
marginal.

Figure 13: Spatial distribution of the simulated pulsars, along with the
observations, projected onto the Galactic plane for 1 million pulsars
simulated.

Figure 14: Distribution of the observed latitude taken from the ATNF
catalogue, along with the simulations with the death line implementation for 1
million pulsars simulated.

Finally fig. 14 compares the Galactic latitudes from the previous study (see
Fig. 9 of Dirson et al. (2022)), now the simulated distribution resembles more
the observations, especially around all the latitudes that are different from
0°. Therefore, the Galactic potential helped tracking the pulsar trajectories
more realistically.

####  3.2.2 Without a death line

Because the simulation where 1 million pulsars are simulated do not allow to
obtain pulsars older than $4.1\times 10^{7}$ yr, we undertook a thorough
examination of the usefulness of a death line in modeling the canonical pulsar
population. Upon conducting simulations with and without a death line (see
Table 5 for the parameter values used in the simulations), it emerged that the
resultant number of detected pulsars and the corresponding plots exhibited
remarkable similarities. In addition, when we conduct the KS test on $P$ and
$\dot{P}$ between 2 simulations with and without death line, we get p-values
of 0.71 and 0.16 respectively, which confirms the similarity between the two
simulations. Consequently, our attention turned towards exploring the impact
of lower birth rates to ascertain if they provide a more congruous fit for
both spin period and spin period derivative distributions.

Figure 15: $P-\dot{P}$ diagram of the simulated population, along with the
observations, without death line implementation, birth rate of 1/70 yr-1 for 1
million pulsars simulated.

Figure 16: $P-\dot{P}$ diagram of the simulated population, along with the
observations, without death line implementation, birth rate of 1/150 yr-1 for
1 million pulsars simulated.

Running simulations with a birth rate of $1/70$ yr-1 or a birth rate of
$1/150$ yr-1 (as a consequence, increasing the oldest simulated pulsar),
conspicuously revealed an insufficient number of detected pulsars in both Fig.
15 and Fig. 16. Intriguingly, a discernible paucity of pulsars congregating in
the bottom right quadrant of the diagram, colloquially known as the pulsar
graveyard, persisted even under these conditions, therefore regardless of the
birth rates, thus casting doubt on the necessity of a death line in this
configuration. Ultimately, our analysis underscores the pivotal importance of
delineating the maximum age of simulated pulsars, as elucidated in Subsection
3.1, while the incorporation of a death line do not seem indispensable for
attaining realistic results when simulating 1 million pulsars.

###  3.3 Spin-velocity misalignment effect

With the aim of observing the effect of the Galactic potential, we looked at
the alignment of the rotation axis and the proper motion velocities of the
pulsars. Results are taken from the simulation containing 1 million pulsars,
see Table 5 for the parameters used. Noutsos et al. (2013) showed by studying
a catalog of 58 pulsars that pulsars younger than 10 Myr keep their rotation
axis and their velocities vector still aligned or anti-aligned, because they
were aligned at birth according to Rankin (2007). Meanwhile, pulsars older
than 10 Myr do not show this trend anymore. This progressive misalignment
effect is due to the movement through the Galactic potential.

Figure 17: Distribution of the spin-velocity angles for all the detected
pulsars in the simulation for 1 million pulsars simulated.

Figure 18: Distribution of the spin-velocity angles for the detected pulsars
younger than 10 Myr (blue) and pulsars older than 10 Myr (red) in the
simulation.

Fig. 17 demonstrates that most pulsars have their rotation axis aligned with
their velocities, but for spin-velocity angles between 25° and 180° the
numbers of pulsars is considerably lower than around 0° and 180°. On the one
hand, for pulsars younger than 10 Myr, shown in blue in Fig. 18, only a small
fraction of pulsars have their spin-velocity angle greater than 25° and below
150°. On the other hand, for pulsars older than 10 Myr, shown in red in Fig.
18, their spin-velocity angle becomes greater than 25° or below 150°.

Figure 19: The spin-velocity angles for the detected pulsars in function of
the number of orbits made by the pulsars.

Furthermore, in Fig. 19, we can confirm that, the older the pulsar is, the
better is its chance to orbit more in the Galaxy, allowing it to deviate its
spin-velocity angle from alignment (or anti-alignment). In addition, even if
the pulsar is not that old, if it made many orbits in the Galaxy, its spin-
velocity angle will be between 25° and 150°, which confirms the fact that its
the Galactic potential (and not only age) which makes this spin-velocity angle
increase (or decrease if its is anti-aligned at birth). In the contrary, there
are pulsars which are old which did not made many orbits in the Galaxy, as a
consequence they have their spin-velocity angle near 0° or 180°. These plots
clearly shows the influence of the Galactic potential making the pulsars lose
their alignment between their rotation axis and their velocities through time
because they evolve in the Galactic potential.

###  3.4 Log normal vs normal distribution for spin period at birth

Figure 20: Distribution of the observed period taken from the ATNF catalogue, along with the simulations with the death line implementation and the Gaussian distribution at birth for the spin period. Table 9: Number of pulsars detected with the death line implementation and the Gaussian distribution at birth for the spin period. log($\dot{E}$) (W) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
¿31 | 7 | 0 | 5 | 2
¿28 | 275 | 9 | 138 | 128
Total | 1642 | 1188 | 223 | 231

In order to estimate the impact of a new distribution of spin period at birth,
we compare the results with the previous distribution of spin period

$p(P_{0})=\frac{1}{\sigma_{p}\sqrt{2\pi}}e^{-(P_{0}-\bar{P})^{2}/(2\sigma_{p})^{2}}$
(34)

The simulation in this subsection uses $\bar{P}=60$ ms and $\sigma_{p}=10$ ms
as parameters of equation (34). The Galactic potential and the death line were
both taken into account to compare the influence of the initial spin period
distribution and 1 million pulsars were simulated, see Table 5 for the
parameters used.

Table 9 highlights the fact that much more pulsars are detected above
$\dot{E}>10^{28}$ W compared to the simulation with the Gaussian in log space
for the spin period, see Table 7. It confirms the suspicion that the
distribution of initial spin period influences the number of pulsars at higher
spin down luminosity.

The spin period distribution in Fig. 20 shows some simulated pulsars with a
low spin period, in excess compared to the observations. This artifact
disappears when the Gaussian distribution in log space for the spin period is
used. The log-normal distribution enables to remove those pulsars in excess
and make the distribution of the spin period better fit the observations. The
KS test was also conducted on both $P$ and $\dot{P}$ and we find p-values,
once again, considerably below 0.05 for both quantity, therefore rejecting the
null hypothesis. We conclude that the initial distribution for the spin period
found by Igoshev et al. (2022) allows to better reproduce the canonical
pulsars population.

###  3.5 Gamma-ray detection

Let us now focus on the $\gamma$-ray band of the pulsar population. First, we
note that the total number of pulsars detected simultaneously in gamma-ray and
in radio is significantly larger in our simulation compared to the
observations, when using the Fermi/LAT instrument sensitivity. However, the
simulated gamma-ray pulsars are in the expected area of Fig. 21, meaning the
striped wind model allow to correctly reproduce these pulsars. The explanation
of this higher number of detection would be that there are many unidentified
sources in the Fourth Fermi-LAT Catalog of Gamma-ray sources (4FGL). Since
many unidentified Galactic sources are listed in these catalogs, many gamma-
ray pulsars might have been detected without being identified as such yet.
Furthermore, we can compare our results with the PPS work of Pierbattista et
al. (2012), where they simulated the young gamma-ray pulsars population with 4
different models : Outer Gap (OG), Slot Gap (SG), Polar Cap (PC) and Outer
Polar Cap (OPC) models. In summary their best model was the OPC model, which
shows almost identical results as ours (the $P-\dot{P}$ shape and the number
of detection), except the fact that we detect gamma-ray pulsars a little
further away from Earth.

Figure 21: $P-\dot{P}$ diagram of the gamma-only pulsars for both the
simulations and the observations, with 1 million pulsars simulated.

What if we run another simulation (with the same parameters as in Table 5), by
only changing the threshold of detection for gamma-ray ? We did this in order
to see how many more pulsars we would detect in gamma-ray and where they would
be in the $P-\dot{P}$ diagram, and also to deal with the problematic of the
GeV excess in the Galactic centre. Therefore we increased $F_{\rm min}$ by 10
: if the galactic latitudes of the pulsar is ¡ 2°, then $F_{\rm min}=0.4\times
10^{-15}$ W.m-2, and concerning blind searches $F_{\rm min}=1.6\times
10^{-15}$ W.m-2.

Figure 22: Spatial distribution of the simulated gamma-ray detected pulsars,
projected onto the Galactic plane for 1 million pulsars simulated, with a
lower $F_{min}$.

The gamma-ray detected pulsar positions in the X-Y plane of the Galaxy are
shown in Fig. 22. In this simulation, 1832 gamma-ray pulsars are detected,
among them only 168, in red in Fig. 22, would have been detected if the same
sensitivity as before was kept. Therefore we detect 11 times more gamma-ray
pulsars compared to the number found with the previous sensitivity.
Furthermore, the other 1664 pulsars, in green in Fig. 22, detected with this
increased sensitivity highlight several features. We obviously detect more
pulsars close to us, and also many more in the centre of the Galaxy. This
result indicates that the GeV excess in the Galactic centre could be linked to
pulsars not yet identified, while usually the hypothesis is that the GeV
excess corresponds to self-annihilating dark matter particles (Hooper, 2022).
In addition, millisecond pulsars, not modeled in this work, could however also
contribute to this GeV excess.

##  4 Discussion

For the initial spin period distribution Igoshev et al. (2022) suggested
values for the mean period $P_{\text{mean}}$ in the range from approximately
$57$ ms to $129$ ms and for the spread $\sigma_{p}$ the range from $0.45$ to
$0.65$. The values $P_{\text{mean}}=129$ ms and $\sigma_{p}=0.45$ are best
educated guesses. However it is obvious that further investigations could
benefit from an automated method to find the best set of values. Hence, we are
currently working on an optimization algorithm to find the optimal initial
parameters for the population synthesis, not only for the spin period, but for
the whole parameter space presented in Table 5, except for $\alpha_{d}$ and
$k_{\tau_{d}}$. Our optimization algorithm will help to constrain the
properties of the Milky Way pulsar population by constraining parameters like
the birth rate which is supposed to be between $1/33$ yr-1 and $1/150$ yr-1,
see Faucher-Giguère & Kaspi (2006), Gullon et al. (2014) or Johnston &
Karastergiou (2017).

We also added the effect of the ISM on the radio pulse profile, rendering more
difficult to observe pulsars farther away than $10$ kpc. Furthermore, because
of this plasma propagation effect fewer radio pulsars are detected compared to
the observations. If the simulations discard the ISM effect, much more pulsars
are detected at distances larger than $10$ kpc and with higher $\dot{E}$ as
well, improving the match between the simulation and the observation. Perhaps,
the effect of the ISM is overestimated when computing the scattering and
dispersion effect converted in equation (28a) with $\tau_{scat}$ and
$\tau_{DM}$ is too important and should have less influence on the observed
width pulse profile, therefore with formulas that would give lower $\tau_{DM}$
and $\tau_{scat}$, we might get more realistic results.

In the main part of the paper, we abstain from showing outcomes stemming from
simulations featuring a constant magnetic field. However we ran such
simulations that clearly highlight the inferior performance of a simulation
without magnetic field evolution as can be checked in appendix B. Moreover
appendix C explores different parameter spaces for the PPS in order to
meticulously scrutinizing whether the parameters showcased in the main part of
the paper optimally encapsulate observational data. This demonstrates that
even if a substantial number of parameters are required for the PPS, they can
be meaningfully constrained by such studies with reasonable accuracy.

Furthermore, in Subsection 3.5 we explore the results when the sensitivity in
gamma-ray is increased by a factor of 10. However we did not show what would
happen if we tried to make a prediction for the radio detection by increasing
the sensitivity by a factor of 50 to match approximately the properties of
SKA. A simulation was run with this increased sensitivity and 52276 pulsars
were detected (when simulating 1 million pulsars, with a death line),
therefore approximately 31 times more than with the sensitivity of Parkes and
Arecibo telescopes.

##  5 Summary

Without the implementation of the Galactic potential, in the PPS of Dirson et
al. (2022), the pulsars were moving in straight lines along one direction
chosen randomly for each pulsar. This evolution of the proper motion was not
realistic, even though the evolution model had already replicated a
significant portion of the canonical pulsar population. However all the new
improvements put into this PPS, meaning : the Galactic potential, the effect
of the ISM on the radio pulse profiles and the fact that the decay timescale
of the magnetic field is taken randomly within a certain range, allow a more
accurate reproduction of the population in the $P-\dot{P}$ diagram.

The first results where we emulate 10 millions pulsars showed an excess of old
pulsars detected in the simulation. In the $P-\dot{P}$ diagram we saw that
this excess was associated to pulsars which lie below the death line,
therefore, not creating pairs anymore. With this finding we realized that it
was imperative to implement the death line. The results obtained with the
death line were then in better agreement with observations when plotting the
most relevant distributions such as period and period derivative.

Nonetheless, simulating 1 million pulsars allowed to show that with or without
death line the results were alike and therefore, in the case where the oldest
pulsar simulated is $4.1\times 10^{7}$ yr, a death line is not really
necessary to explain the observations. In this paper we showed two approach to
reproduce the observed canonical pulsars population, however when we consider
the approach where the oldest pulsar simulated is only $4.1\times 10^{7}$ yr
we state that no older pulsars might exist in the Galaxy, while the Milky Way
is very much older than this. There is apparently no reason to not simulate
pulsars that are older than this, and if there is one we do not know it (they
become harder to detect because of their wide pulse profile with their age,
but some are still detectable). If this situation is the reality then it means
two things : the age of the observed pulsars are overestimated (because we
observe pulsars older than $4.1\times 10^{7}$ yr in the canonical population,
even if they are not numerous) and the death line is not necessary in pulsar
population synthesis work to explain what we observe. When the simulation is
ran with 10 millions pulsars, the death line becomes a necessity. However this
simulation gives results less similar to the observations, because of a few
simulated data which do not match any observational data.

It was also shown that the old pulsars with age larger than 10 Myr have a
tendency to lose their alignment between rotation axis and velocity, a direct
effect of the Galactic potential as shown by the observations of Noutsos et
al. (2013), but it is the first time it is done in a simulation.

The striped wind model seems to reproduce the gamma-ray population of
canonical pulsars very well, however in the simulation a higher number of
these pulsars are detected compared to reality. It probably indicates that a
part of unidentified sources from the 4FGL catalog are pulsars. Testing an
increased sensitivity for the detection of gamma-ray pulsars lead us to think
that a part of the origin of the GeV excess in the Galactic centre could be
young gamma-ray pulsars, and when the millisecond pulsars will be modeled in
our simulation we would also be able to check that these pulsars might be
responsible for this GeV excess. The usual picture is to suppose that only
millisecond pulsars are responsible, while young gamma-ray pulsars could also
play a role.

In summary, the best parameters found for this pulsar population synthesis are
for $P_{\text{mean}}=129$ ms, $\sigma_{p}=0.45$ for the spin period
distribution at birth, $B_{\text{mean}}=2.75\times 10^{8}$ T, $\sigma_{b}=0.5$
for the magnetic field distribution at birth and for a birth rate of 1/41
yr-1. Ultimately, this work provides a simulated population of pulsars within
the Milky Way more similar to the observations compared to Johnston et al.
(2020) and Dirson et al. (2022), the two most recent population synthesis
considering both radio and gamma-ray emission. Moreover, such investigations
will be useful to make more predictions for future radio surveys such as SKA,
or NenuFAR.

The PPS method used in this study could be applied to model populations of
pulsars in the Milky Way other than the canonical ones. The magnetars, for
instance, would need : another distribution at birth for their magnetic field
to be modeled correctly, taking into account the emission in X-rays and having
a faster decay for the magnetic field. The magnetars are the pulsars just
above the canonical pulsars in the $P-\dot{P}$ diagram, at the top right of
the diagram. A more complex population is the recycled or millisecond pulsars,
that were ”brought back to life” by accreting matter from a companion star.
This accretion mechanism spins up the neutron stars moving them from the
graveyard of pulsars to the bottom left of the $P-\dot{P}$ diagram, crossing
again the death line in the opposite sense and allowing them to create
electron-positron pair cascades again. However, they are more complex to model
because this evolution scenario should take into account the accretion process
and spin-up. We are also currently working on both populations to reproduce
them in a PPS based on this work.

###### Acknowledgements.

This work has been supported by the grant 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. D.M. acknowledges the support of the Department of Atomic Energy,
Government of India, under project No. 12-R&DTFR-5.02-0700.

## References

  * Abdo et al. (2013) Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17
  * Ahlers et al. (2016) Ahlers, Markus ; Bai, Yang ; Barger, Vernon ; Lu, Ran, Physical Review D, Volume 93, Issue 1, id.013009
  * Arons (1983) Arons, J. 1983, ApJ, 266, 215
  * Bajkova & Bobylev (2021) Bajkova, A. T., Bobylev, V. V. Astronomical & Astrophysical Transactions, 2021, Vol. 32, Issue 3, pp. 177-206
  * Bates et al. (2014) Bates, S. D., Lorimer, D. R., Rane, A., & Swiggum, J. 2014, MNRAS, 439, 2893
  * Bogdanov et al. (2019) Bogdanov, S., Lamb, F. K., Mahmoodifar, S., et al. 2019, ApJ, 887, L26
  * Bovy (2015) Bovy, J. The Astrophysical Journal Supplement Series, Volume 216, Issue 2, article id. 29, 27 pp. (2015)
  * Cameron et al. (2020) Cameron A. D., et al., 2020, MNRAS, 493, 1063
  * Chen & Ruderman (1993) Chen, Kaiyou ; Ruderman, Malvin, January 1993, Astrophysical Journal v.402, p.264
  * Cheng et al. (1986) Cheng, K. S., Ho, C., & Ruderman, M. 1986, ApJ, 300, 500
  * Cordes & McLaughlin (2003) Cordes, J. M., & McLaughlin, M. A. 2003, ApJ, 596, 1142
  * Cordes et al. (2006) Cordes J. M., et al., 2006, ApJ, 637, 446
  * Daugherty & Harding (1982) Daugherty, J. K., & Harding, A. K. 1982, ApJ, 252, 337
  * Daugherty & Harding (1996) Daugherty, J. K., & Harding, A. K. 1996, ApJ, 458, 278
  * Dirson et al. (2022) Dirson, Pétri, and Mitra 2022, Astronomy & Astrophysics, Volume 667, id.A82, 12 pp.
  * Faucher-Giguère & Kaspi (2006) Faucher-Giguère, C.-A., Kaspi, V. M. 2006, ApJ, 643, 332
  * Gil & Mitra (2001) Gil, Janusz ; Mitra, Dipanjan, The Astrophysical Journal, Volume 550, Issue 1, pp. 383-391, March 2001
  * Gonthier et al. (2018) Gonthier, Peter L.; Harding, Alice K.; Ferrara, Elizabeth C. ; Frederick, Sara E.; Mohr, Victoria E. ; Koh, Yew-Meng, The Astrophysical Journal, Volume 863, Issue 2, article id. 199, 27 pp. (2018)
  * Gonthier et al. (2002) Gonthier, Peter L. ; Ouellette, Michelle S. ; Berrier, Joel ; O’Brien, Shawn ; Harding, Alice K. The Astrophysical Journal, Volume 565, Issue 1, pp. 482-499, January 2002
  * Graber et al. (2023) Graber, Vanessa ; Ronchi, Michele ; Pardo-Araujo, Celsa ; Rea, Nanda, eprint arXiv:2312.14848
  * Gullon et al. (2014) Gullón, M., Miralles, J. A., Viganò, D., Pons, J. A. 2014, MNRAS, 443, 1891
  * Harding et al. (2008) Harding, A. K., Stern, J. V., Dyks, J., Frackowiak, M. 2008, ApJ, 680, 1378
  * Harding & Muslimov (2011) Harding, A. K., Muslimov, A. G. 2011, ApJ, 743, 181
  * Hewish et al. (1968) Hewish, A. ; Bell, S. J. ; Pilkington, J. D. H. ; Scott, P. F. ; Collins, R. A. Nature, Volume 224, Issue 5218, pp. 472 (1969)
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D. R., Lyne, A. G., Kramer, M. 2005, MNRAS, 360, 974
  * Hofmann & Zanin (2023) Hofmann & Zanin, May 2023, eprint arXiv:2305.12888
  * Hooper (2022) Hooper, Dan, September 2022, eprint arXiv:2209.14370
  * Hirotani (2008) Hirotani, K. 2008, ApJ, 688, L25
  * Igoshev et al. (2022) Igoshev, Andrei P., Frantsuzova, Anastasia, Gourgouliatos, Konstantinos, N. Monthly Notices of the Royal Astronomical Society, Volume 514, Issue 3, pp. 4606-4619
  * Jawor & Tauris (2022) Jawor, J.A.,Tauris, T.M. Monthly Notices of the Royal Astronomical Society, Volume 509, Issue 1, pp.634-657, 2022
  * Johnston & Karastergiou (2017) Johnston, S., Karastergiou, A. Monthly Notices of the Royal Astronomical Society, Volume 467, Issue 3, p.3493-3499, 2017
  * Johnston & Karastergiou (2019) Johnston, Simon ; Karastergiou, A., Monthly Notices of the Royal Astronomical Society, Volume 485, Issue 1, p.640-647
  * Johnston et al. (2020) Johnston, S., Smith, D. A., Karastergiou, A., Kramer, M. 2020, MNRAS, 497, 1957
  * Johnston et al. (2023) Johnston et al., Monthly Notices of the Royal Astronomical Society, Volume 520, Issue 4, pp.4801-4814
  * Kalapotharakos et al. (2019) Kalapotharakos, C., Harding, A. K., Kazanas, D., Wadiasingh, Z. 2019, ApJ, 883, L4
  * Kramer et al. (2003) Kramer M., et al., 2003, MNRAS, 342, 1299
  * Krishnakumar et al. (2015) Krishnakumar, M. A., Mitra, D., Naidu, A., Joshi, B. C., & Manoharan, P. K. 2015, ApJ, 804, 23
  * Li et al. (2019) Li, C., Zhao, G., Jia, Y., et al. 2019, ApJ, 871, 208
  * Lorimer & Kramer (2004) Lorimer, D. R., Kramer, M. 2004, Handbook of Pulsar Astronomy
  * Lorimer et al. (2006) Lorimer D. R., et al., 2006, ApJ, 640, 428
  * Manchester et al. (2005) Manchester, R. N., Hobbs, G. B., Teoh, A., Hobbs, M. 2005, AJ, 129, 1993
  * McMullin et al. (2020) McMullin et al., Proceedings of the SPIE, Volume 11445, id. 1144512 18 pp. (2020)
  * Mitra (2017) Mitra, D. 2017, J. Astrophys. Astron., 38, 52
  * Mitra et al. (2019) Mitra, D., Basu, R., Melikidze, G. I., Arjunwadkar, M. Monthly Notices of the Royal Astronomical Society, Volume 492, Issue 2, p. 2468-2480
  * Miyamoto & Nagai (1975) Miyamoto & Nagai, Publications of the Astronomical Society of Japan, Vol. 27, 533-543 (1975)
  * Muslimov & Harding (2004) Muslimov, A. G., Harding, A. K. 2004, ApJ, 617, 471
  * Navarro et al. (1997) Navarro, Frenk, White, The Astrophysical Journal, Volume 490, Issue 2, pp. 493-508
  * Noutsos et al. (2013) Noutsos, A., Schnitzeler, D. H. F. M., Keane, E. F., Kramer, M., Johnston, S. Monthly Notices of the Royal Astronomical Society, Volume 430, Issue 3, p. 2281-2301
  * Paczynski (1990) Paczynski, Bohdan, Astrophysical Journal v.348, p. 485
  * Peter Young (2014) Peter Young. The leapfrog method and other “symplectic” algorithms for integrating Newton’s laws of motion. 2014.
  * Pétri (2009) Pétri, J. 2009, A&A, 503, 13
  * Pétri (2011) Pétri, J. 2011, MNRAS, 412, 1870
  * Pétri (2012) Pétri, J. 2012, MNRAS, 424, 605
  * Phillipov et al. (2014) Philippov, A., Tchekhovskoy, A., Li, J. G. 2014, MNRAS, 441, 1879
  * Pierbattista et al. (2012) Pierbattista, M. ; Grenier, I. A.; Harding, A. K. ; Gonthier, P. L., Astronomy & Astrophysics, Volume 545, id.A42, 23 pp., September 2012
  * Popov et al. (2010) Popov, S. B. ; Pons, J. A. ; Miralles, J. A. ; Boldin, P. A. ; Posselt, B, MNRAS, Volume 401, Issue 4, pp. 2675-2686, February 2010
  * Rankin (2007) Rankin, Joanna M, The Astrophysical Journal, Volume 664, Issue 1, pp. 443-447
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, ApJ, 887, L21
  * Ronchi et al. (2021) Ronchi, M. ; Graber, V.; Garcia-Garcia, A.; Rea, N.; Pons, J. A., The Astrophysical Journal, Volume 916, Issue 2, id.100, 20 pp.
  * Ruderman & Sutherland (1975) Ruderman, M. A., Sutherland, P. G. 1975, ApJ, 196, 51
  * Siegert (2019) Siegert, T. 2019, A&A, 632, L1
  * Skowron et al. (2019) Skowron, D. M., Skowron, J., Mróz, P., et al. 2019, Sci, 365, 478
  * Smirnov (1948) Smirnov, N. (1948). Table for estimating the goodness of fit of empirical distributions. Annals of Mathematical Statistics, 19(2), 279-281.
  * Smith et al. (2019) Smith, D. A., Bruel, P., Cognard, I., et al. 2019, ApJ, 871, 78
  * Smith et al. (2023) Smith et al., The Astrophysical Journal, Volume 958, Issue 2, id.191, 72 pp.
  * Spitkovsky (2006) Spitkovsky, Anatoly, The Astrophysical Journal, Volume 648, Issue 1, pp. L51-L54
  * Sturrock (1971) Sturrock, P. A. 1971, ApJ, 164, 529
  * Takata et al. (2011) Takata, J., Wang, Y., Cheng, K. S. 2011, ApJ, 726, 44
  * van der Kruit (1995) van der Kruit, P. C. in Comparison of the Galaxy with External Spiral Galaxies, eds. G. Gilmore, & B. Carswell (Dordrecht: Springer, Netherlands), 27
  * Viganò (2013) Viganò, D. 2013, ArXiv e-prints [arXiv:1310.1243]
  * Viganò et al. (2013) Viganò, D., Rea, N., Pons, J. A., et al. 2013, MNRAS, 434, 123
  * Watters & Romani (2011) Watters, K. P., Romani, R. W. 2011, ApJ, 727, 123
  * Weltevrede & Johnston (2008) Weltevrede, Patrick ; Johnston, Simon, Monthly Notices of the Royal Astronomical Society, Volume 391, Issue 3, pp. 1210-1226
  * Yadigaroglu & Romani (1995) Yadigaroglu, I.-A., Romani, R. W. 1995, ApJ, 449, 211
  * Yao et al. (2017) Yao, J. M. ; Manchester, R. N. ; Wang, N, The Astrophysical Journal, Volume 835, Issue 1, article id. 29, 32 pp. (2017).
  * Zhang et al. (2000) Zhang, Bing ; Harding, Alice K. ; Muslimov, Alexander G., The Astrophysical Journal, Volume 531, Issue 2, pp. L135-L138, March 2000

##  Appendix A The integration scheme used

###  A.1 The scheme

One of the new aspect in this work is to solve the equation of motion (12) in
order to follow more realistically the motion of our sample of pulsars. To
achieve this objective we implemented a fourth order integration scheme called
Position Extended Forest Ruth-Like (PEFRL). This algorithm was chosen because
of its high order allowing a high precision on the trajectory of the pulsar.
In addition, it is a symplectic algorithm, meaning that it conserves the total
energy of the pulsar (kinetic plus potential), minimizing the error even if
the integration is done for a long time with any orbits. It works in the same
way as a Verlet scheme of second order, the velocities and coordinates are
shifted separately but in the end they get synchronized. The different steps
of the algorithm are as follows

$\displaystyle\@vec{r}^{n+\frac{1}{5}}$ $\displaystyle=\@vec{r}^{n}+\xi
h\@vec{v}^{n}$ (35a) $\displaystyle\@vec{v}^{n+\frac{1}{4}}$
$\displaystyle=\@vec{v}^{n}-(1-2\lambda)\frac{h}{2}\@vec{\nabla}\Phi(\@vec{r}^{n+\frac{1}{5}})$
(35b) $\displaystyle\@vec{r}^{n+\frac{2}{5}}$
$\displaystyle=\@vec{r}^{n+\frac{1}{5}}+\chi h\@vec{v}^{n+\frac{1}{4}}$ (35c)
$\displaystyle\@vec{v}^{n+\frac{2}{4}}$
$\displaystyle=\@vec{v}^{n+\frac{1}{4}}-\lambda
h\@vec{\nabla}\Phi(\@vec{r}^{n+\frac{2}{5}})$ (35d)
$\displaystyle\@vec{r}^{n+\frac{3}{5}}$
$\displaystyle=\@vec{r}^{n+\frac{2}{5}}+(1-2(\chi+\xi))h\@vec{v}^{n+\frac{2}{4}}$
(35e) $\displaystyle\@vec{v}^{n+\frac{3}{4}}$
$\displaystyle=\@vec{v}^{n+\frac{2}{4}}-\lambda
h\@vec{\nabla}\Phi(\@vec{r}^{n+\frac{3}{5}})$ (35f)
$\displaystyle\@vec{r}^{n+\frac{4}{5}}$
$\displaystyle=\@vec{r}^{n+\frac{3}{5}}+\chi h\@vec{v}^{n+\frac{3}{4}}$ (35g)
$\displaystyle\@vec{v}^{n+1}$
$\displaystyle=\@vec{v}^{n+\frac{3}{4}}-(1-2\lambda)\frac{h}{2}\@vec{\nabla}\Phi(\@vec{r}^{n+\frac{4}{5}})$
(35h) $\displaystyle\@vec{r}^{n+1}$
$\displaystyle=\@vec{r}^{n+\frac{4}{5}}+\xi h\@vec{v}^{n+1}$ (35i)

here $\@vec{r}$ is the position vector, $\@vec{v}$ is the velocity vector, $h$
is the time step, $\@vec{\nabla}\Phi$ is the gradient of the potential and $n$
is the time at which we compute the position or the velocity. $\lambda$, $\xi$
and $\chi$ are constants that are tabulated in Table 10, see for instance
Peter Young (2014).

Table 10: Value of the constants used for the PEFRL scheme. Parameters | Values
---|---
$\lambda$ | -0.212341831062605
$\xi$ | 0.178617895844809
$\chi$ | -0.06626458266981849

###  A.2 Precision of the PEFRL algorithm

Before implementing the gravitational potential and the integration scheme in
the whole synthesis, the precision of the trajectory of a single pulsar in the
Milky Way was checked. The initial conditions for the pulsar are the following
: $x=5.1$ kpc, $y=-1.9$ kpc, $z=0.05$ kpc, $v_{x}=-47.8$ km/s $v_{y}=97.0$
km/s, and $v_{z}$ = 22.6 km/s. Those initial conditions were chosen randomly,
the only constrain being to have the pulsar in the Galaxy with realistic
coordinates on the one hand and with realistic speeds for a pulsar on the
other hand. Furthermore, the total integration time is 1.25 Gyr and the time
step used is $10^{5}$ yr. In order to check the precision of the scheme, the
total energy of the pulsar, sum of the kinetic plus potential energy, that
should be conserved is computed as

$E_{tot}=\frac{v_{x}^{2}}{2}+\frac{v_{y}^{2}}{2}+\frac{v_{z}^{2}}{2}+\Phi_{tot}.$
(36)

The trajectory obtained after the integration with the PEFRL algorithm is
shown in Fig. 23. This trajectory shows that the pulsar, with these initial
conditions, remains bound to the Galaxy. The pulsar here takes approximately
201.3 Myr to complete an orbit in the Galaxy. The time step chosen is
significantly smaller than the orbital period ($h/P_{\rm orb}=0.0005$), and we
obtain approximately 6 orbits by integrating the trajectory during 1.25 Gyr.
Bajkova & Bobylev (2021) showed several trajectories in their paper, the one
obtained here seems coherent compared to the different trajectories that they
got.

Figure 23: Trajectory of a single pulsar in the Galactocentric frame.

Fig. 24 shows the relative error on the total energy. It is indeed almost
constant along the whole integration of the trajectory. The maximum relative
error on the energy is of the order of $10^{-9}$, a very good accuracy for
this trajectory. Globally, when a single pulsar has its trajectory integrated
in the PPS, the relative error is at maximum 0.001 and the minimum seen was
$10^{-14}$.

Figure 24: Relative error on the energy of a single pulsar in the
Galactocentric frame.

##  Appendix B Results for a constant magnetic field

Depending on the PPS, the magnetic field will be considered constant (Johnston
& Karastergiou 2017; Faucher-Giguère & Kaspi 2006) or decaying (Gullon et al.
2014; Dirson et al. 2022; Jawor & Tauris 2022). The goal is here is to show
the results with a constant magnetic field in order to see which model is
better between the decaying magnetic field and the constant magnetic field.

Hence, this simulation was run with 1 million pulsar simulated, without a
death line and without a decaying of the magnetic field. The KS test gave
p-values of $6.3\times 10^{-22}$ and $9.8\times 10^{-77}$, therefore the null
hypothesis is rejected for both $P$ and $\dot{P}$. The parameters used are in
Table 5. It is less good than with a decaying magnetic field. As can be seen
in Fig. 25, too many pulsars from the simulation appear on the top right of
the diagram and almost none at the bottom. In Table 11, we note a total of
pulsars which is smaller from what we found for all of our previous results in
Sect. 3.

Figure 25: $P-\dot{P}$ diagram of the simulated population, along with the observations, with the death line implementation and a constant magnetic field for each pulsar. Table 11: Number of pulsars detected for the simulation with constant magnetic field for each pulsar and the death line implementation. log($\dot{E}$) (W) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
¿31 | 15 | 0 | 12 | 3
¿28 | 168 | 4 | 94 | 70
Total | 1356 | 1034 | 162 | 160

##  Appendix C Influence of the different parameters

In order to show that the parameters chosen in the main part of this paper are
meticulously picked and that the uncertainties brought by the large number of
parameters is manageable, several simulations were run with different values
for the main parameters of the initial distributions, showing a deviation from
the best values gives bad results. The simulations were run with in the same
conditions as Subsection 3.2, with only one parameter changing between the
different sets of simulation.

Table 12: Influence of the birth rate. Birth rate (in yr-1) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
1/30 | 2388 | 1964 | 204 | 220
1/80 | 963 | 799 | 79 | 85
1/120 | 696 | 597 | 50 | 49
Table 13: Influence of the mean initial period. $P_{\rm mean}$ (in ms) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
57 | 1499 | 1134 | 193 | 172
70 | 1623 | 1253 | 215 | 155
100 | 1652 | 1348 | 151 | 153
Table 14: Influence of the mean initial magnetic field. $B_{mean}$ (in T) | $N_{tot}$ | $N_{r}$ | $N_{g}$ | $N_{rg}$
---|---|---|---|---
$10^{7}$ | 8407 | 6943 | 406 | 1058
$5\times 10^{7}$ | 4447 | 3646 | 306 | 495
$5\times 10^{8}$ | 1111 | 933 | 109 | 69

The $P-\dot{P}$ diagram of each simulation is not shown here because almost
every simulation finds a too high (or too low) number of pulsars, therefore
the $P-\dot{P}$ diagrams are showing large differences between observations
and simulations. Nonetheless, the number of pulsars detected for two of the
simulations with different $P_{\rm mean}$ (with $P_{\rm mean}=70$ ms and
$P_{\rm mean}=100$ ms) gives number closer to reality. However their p-values
for $P$ and $\dot{P}$ when using the KS test were significantly below $0.05$,
which allowed us to reject the null hypothesis for the distributions of $P$
and $\dot{P}$, for the results shown in Table 12, 13 and 14.

Moreover in order to get an idea on a good range of values for the parameters,
the d-value obtained thanks to the KS test are shown in Table 15, 16, 17 and
18. For each parameter, 10 simulations were made in order to get a mean and a
standard deviation for the d-value associated. These simulations helped
finding the best values for the parameters by considering the parameter is
well constrained when its d-values for both distribution of $P$ and $\dot{P}$
are below, or really close to 0.05 (meaning 5 % similarity between the
observations and simulation distribution). Even though it is not an
optimisation (and therefore the results can be improved on that point), it
allows us to conclude that $P_{\rm mean}$ should be between 120 ms and 135 ms
(with 129 ms as best value). $B_{\rm mean}$ should be between $2.25\times
10^{8}$ T and $2.8\times 10^{8}$ T (with $2.75\times 10^{8}$ T as best value).
The birth rate should be between 1/37 yr-1 and 1/45 yr-1 (with 1/41 yr-1 as
best value). Finally, $\sigma_{b}$ is well constrained between 0.49 and 0.51
(with 0.5 as best value). Concerning $\sigma_{p}$, it always gives d-values
between 0.04 and 0.057 when $\sigma_{p}$ is below 0.6, therefore the value
0.45 was chosen as suggested by Igoshev et al. (2022).

Table 15: D-value for different $P_{\rm mean}$. $P_{\rm mean}$ (in ms) | d-value ($\dot{P}$) | d-value ($P$)
---|---|---
90 | $0.072\pm 0.0073$ | $0.0577\pm 0.0080$
100 | $0.0651\pm 0.00692$ | $0.0503\pm 0.00636$
110 | $0.0652\pm 0.00688$ | $0.0381\pm 0.0074$
120 | $0.0554\pm 0.00714$ | $0.0302\pm 0.00468$
129 | $0.0557\pm 0.00740$ | $0.0284\pm 0.00583$
135 | $0.0544\pm 0.00768$ | $0.0329\pm 0.00767$
150 | $0.0572\pm 0.00779$ | $0.0388\pm 0.00935$
Table 16: D-value for different $B_{\rm mean}$. $B_{\rm mean}$ (in T) | d-value ($\dot{P}$) | d-value ($P$)
---|---|---
$1\times 10^{8}$ | $0.113\pm 0.00823$ | $0.203\pm 0.00675$
$2\times 10^{8}$ | $0.0418\pm 0.00459$ | $0.0658\pm 0.00705$
$2.25\times 10^{8}$ | $0.0456\pm 0.00455$ | $0.0455\pm 0.0061$
$2.75\times 10^{8}$ | $0.0557\pm 0.00740$ | $0.0284\pm 0.00583$
$2.8\times 10^{8}$ | $0.0576\pm 0.00907$ | $0.0298\pm 0.00728$
$3\times 10^{8}$ | $0.0648\pm 0.01092$ | $0.0426\pm 0.00775$
$3.25\times 10^{8}$ | $0.0761\pm 0.00930$ | $0.0480\pm 0.00997$
Table 17: D-value for different $\sigma_{b}$. $\sigma_{b}$ | d-value ($\dot{P}$) | d-value ($P$)
---|---|---
0.45 | $0.0767\pm 0.00668$ | $0.0513\pm 0.01258$
0.48 | $0.0621\pm 0.00880$ | $0.0385\pm 0.01234$
0.49 | $0.0646\pm 0.00974$ | $0.0304\pm 0.00501$
0.5 | $0.0557\pm 0.00740$ | $0.0284\pm 0.00583$
0.51 | $0.056\pm 0.00735$ | $0.0306\pm 0.00701$
0.52 | $0.0547\pm 0.00786$ | $0.0358\pm 0.00579$
Table 18: D-value for different birth rate (BR). BR (in yr -1) | d-value ($\dot{P}$) | d-value ($P$)
---|---|---
1/33 | $0.0611\pm 0.01025$ | $0.0313\pm 0.00636$
1/37 | $0.0565\pm 0.00851$ | $0.0295\pm 0.00645$
1/41 | $0.0557\pm 0.00740$ | $0.0284\pm 0.00583$
1/45 | $0.0585\pm 0.00510$ | $0.0302\pm 0.00565$
1/70 | $0.0691\pm 0.00830$ | $0.0342\pm 0.01165$
1/150 | $0.1165\pm 0.01564$ | $0.0408\pm 0.00671$