# Binary population synthesis of the Galactic canonical pulsar population

Yuzhe Song  Centre for Astrophysics and Supercomputing, Swinburne University
of Technology, Hawthorn, VIC 3122, Australia  [ Simon Stevenson  Centre for
Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn,
VIC 3122, Australia  Debatri Chattopadhyay  Gravity Exploration Institute,
School of Physics and Astronomy, Cardiff University, Cardiff, CF24 3AA, United
Kingdom

###### Abstract

Pulsars are rapidly rotating neutron stars that emit radiation across the
electromagnetic spectrum, from radio to $\gamma$-rays. We use the rapid binary
population synthesis suite COMPAS to model the Galactic population of
canonical pulsars. We account for both radio and $\gamma$-ray selection
effects, as well as the motion of pulsars in the Galactic potential due to
natal kicks. We compare our models to the catalogs of pulsars detected in the
radio, and those detected in $\gamma$-rays by Fermi, and find broad agreement
with both populations. We reproduce the observed ratio of radio-loud to radio-
quiet $\gamma$-ray pulsars. We further examine the possibility of low spin-
down luminosity ($\dot{E}$ ) pulsars emitting weak, unpulsed $\gamma$-ray
emission and attempt to match this with recent stacking results. We
demonstrate that the apparent correlation between the latitude of a pulsar and
its $\dot{E}$ arises due to natal kicks imparted to pulsars at birth, assuming
that all pulsars are born in the Galactic disk.

###### keywords:

Pulsars: general, stars: massive, binaries: general, gamma-rays: general

OzGrav, ARC Centre for Excellence of Gravitational Wave Discovery, Hawthorn,
VIC 3122, Australia Y<EMAIL_ADDRESS>OzGrav, ARC Centre for
Excellence of Gravitational Wave Discovery, Hawthorn, VIC 3122, Australia

##  1 Introduction

Pulsars are rotating neutron stars most commonly observed by their periodic
radio emission. Around 3500 pulsars are known to date Manchester:2004bp.
Modern radio facilities such as MeerKAT, ASKAP and FAST are uncovering new
pulsars (Han:2021RAA; Padmanabh:2023arXiv, e.g.,), whilst modern search
techniques are still discovering pulsars in decades-old archival surveys
(Sengar:2023MNRAS, e.g.,). The ever-growing pulsar population provides an
invaluable insight into the final remnants of massive stellar evolution,
whilst relativistic binary pulsars allow for exquisitely precise tests of
theories of gravity (Kramer:2021jcw, e.g.,).

Whilst most pulsars are discovered and studied through radio surveys, some are
also observed at other wavelengths, including optical Cocke:1969Natur;
Nather:1969Natur, X-rays Bradt:1969Natur; Fritz:1969Sci; Oegelman:1993Natur
and $\gamma$-rays Kniffen:1974Natur; Tumer:1984Natur; Thompson:1992Natur. A
growing number of pulsars are being detected in $\gamma$-rays by the Large
Area Telescope (LAT) onboard the Fermi space telescope Fermi-LAT:2023zzt,
including some pulsars which were first detected in $\gamma$-rays and then
later detected as radio pulsars, a situation exemplified by the Geminga pulsar
Halpern:1992Natur; Bertsch:1992Natur; Malofeev:1997Natur.

Typically, young pulsars with high spindown luminosity ($\dot{E}$) are found
at low Galactic latitudes, and millisecond pulsars (MSPs) are observed in the
$\gamma$-rays (Abdo:2009Sci; Fermi-LAT:2013svs; Pletsch:2013iva; Fermi-
LAT:2023zzt, e.g.,), where $\dot{E}$ is the spin down luminosity given by

$\dot{E}=4\pi^{2}I\dot{P}P^{-3}.$ (1)

where $P$ is the spin period of the pulsar, $\dot{P}$ is its spin derivative
and $I$ is its moment of inertia. Comparison of the observed $\gamma$-ray
luminosity $L_{\gamma}$ to the spindown luminosity $\dot{E}$ indicates that
$\gamma$-ray emission accounts for a large fraction of the energy radiated
away by many pulsars (Fermi-LAT:2013svs, with this fraction $\eta\gtrsim
10\%$;). Around half of the pulsars detected in $\gamma$-rays are also
observed in radio, whilst the other half (typically referred to as radio quiet
pulsars) are not Fermi-LAT:2008wkm; Abdo:2009Sci; Fermi-LAT:2013svs; Fermi-
LAT:2023zzt. The distinction between these two populations likely depends on
whether or not their radio emission is beamed towards the Earth Abdo:2009Sci.

The current sample of $\gamma$-ray detected pulsars is around 300 Fermi-
LAT:2013svs; Smith:2019ApJ; Fermi-LAT:2023zzt111A list of Fermi detected
pulsars is available at
https://confluence.slac.stanford.edu/display/GLAMCOG/Public+List+of+LAT-
Detected+Gamma-Ray+Pulsars and is limited by the sensitivity of Fermi-LAT.
There is an apparent ‘death-line’ around a spindown luminosity of $\dot{E}\sim
10^{33}$ erg s-1, which may be a result of the limited sensitivity Fermi-
LAT:2013svs. An extensive search of $\gamma$-ray pulsations from a thousand
pulsars Smith:2019ApJ provides a similar result that the lack of detection of
pulsars with lower $\dot{E}$ might be due to the limitation of the current
instruments, and suggests population synthesis work to explore this specific
subpopulation of pulsars.

Recently, Song2023 performed a stacking analysis of 12 years of Fermi data,
targeting the positions of 362 known radio pulsars with no previous
$\gamma$-ray detections. In order to reduce the impact of the Galactic
$\gamma$-ray background, they selected pulsars that lie at high galactic
latitudes ($|b|>20^{\circ}$). This probed a population of relatively old, low
$\dot{E}$ pulsars, different to the usual population of $\gamma$-ray loud
pulsars. Song2023 report the significant detection of $\gamma$-rays from this
set of pulsars compared to the expected background. Song2023 argue that this
suggests that the usual $\gamma$-ray death line is predominantly a result of
the sensitivity of Fermi. In this paper we attempt to explore the implications
of this detection for massive binary and pulsar evolution.

Theoretical population synthesis studies can be used to predict the population
of observable pulsars under a set of assumptions. There have been several
previous population synthesis studies of the $\gamma$-ray pulsar population.
In a seminal work, Faucher-Giguere:2006ApJ performed a population synthesis of
the Galactic radio pulsar population. It accounted for almost all of the
relevant effects, including radio selection effects, magnetic field decay and
pulsar orbits in the Galaxy. Assumptions were made on distributions of pulsar
birth spin period and birth magnetic field strength. Since then, there have
been a number of similar studies building on this work (Takata2011;
Dirson2022, e.g.,) that include radio and $\gamma$-ray selection effect and
pulsar evolution theories. Typically, existing pulsar population synthesis
works, such as those done by using PsrPopPy Bates:2014MNRAS, do not model the
details of the full evolution of the binary stellar systems. The study of
Kiel:2008MNRAS is closest to ours in terms of methodology, but they focus on
radio pulsars and do not explore the $\gamma$-ray pulsar population. In this
work, we attempt to explore the canonical radio and $\gamma$-ray pulsar
populations produced through massive binary stellar evolution. In a follow-up
work, we will attempt to understand the full effects on pulsar populations
from the stellar evolution physics.

Neutron stars—and hence pulsars—originate from massive stars, and the majority
of massive stars are born in binaries (Sana:2012Sci, e.g.,). Hence, whilst the
majority of “normal”, young pulsars (with $P\gtrsim 0.1$ s) are isolated, most
of these objects are expected to have been born in binaries, which are
subsequently disrupted by the supernova that formed the neutron star. Hence,
the binary fraction and spatial distributions (including scale height) of
pulsars encode information about their birth locations and natal kicks. It is
then important to account for a full picture of stellar and binary evolution
for a better understanding of the Galactic pulsar population. Binary
population synthesis code suites, such as COMPAS COMPAS, are ideally suited
for this task, as they are able to rapidly simulate the evolution of large
populations of massive binary stars, incorporating all relevant processes such
as natal kicks and pulsar evolution. Understanding the properties of the
Galactic pulsar population across the radio and $\gamma$-ray wavelengths will
aid in constraining pulsar emission physics and the formation and evolution of
massive binary stars.

As a first study of the Galactic $\gamma$-ray pulsar population using COMPAS,
one of the main purposes of this work is to verify the binary and pulsar
physics included in COMPAS. This work and its follow-ups can also be used to
make predictions for the population of $\gamma$-ray pulsars that will be
observed in the future with continued operation of Fermi-LAT. Other
interesting binary pulsar systems include neutron star X-ray binaries
1998Natur.394..344W, spider pulsars 1988Natur.333..237F and double neutron
star binaries (Hulse:1975ApJL; PhysRevLett.119.161101, e.g.,), and these
events can only be studied through the inclusion of binary physics. In follow-
up studies, we will investigate known peculiar binary systems that contain at
least one neutron star using COMPAS and other suitable tools.

The remainder of this paper is structured as follows. In Section 2 we describe
our methodology. We present our results in Section 3. Section 4 focuses on
follow-up studies from this work, and we conclude in Section 5.

##  2 Methods

In this section we introduce the methods and theory used in this study,
including evolution theory of isolated pulsars, radio observation and
selection criteria, and $\gamma$-ray observation and selection criteria. We
also introduce the basics of COMPAS relevant to synthesising pulsar
populations.

###  2.1 COMPAS

We simulate a population of massive binary stars using the rapid binary
population synthesis suite COMPAS Stevenson:2017tfq; COMPAS. COMPAS is based
on the Single Star Evolution (SSE, SSE,) code and the Binary Star Evolution
(BSE, BSE,) code with updates and modifications made. Binary interactions such
as mass transfer and common envelope evolution are included, but we neglect
the effects of tides COMPAS.

We use default binary property parameters for our COMPAS simulations COMPAS.
These parameters and their values/ranges are based on past theoretical and
observational studies. The mass of the initially more massive star in the
binary is drawn from an initial mass function (Kroupa:2000iv, IMF,) within the
range of 5 to 150 M⊙ as we need stars massive enough to leave behind neutron
stars at the end of their evolution. The mass ratio ($q=m_{2}/m_{1}$,
$m_{2}<m_{1}$) of the binary is drawn from a uniform distribution between 0.01
and 1, with a minimum secondary mass of $0.1$ M⊙ Stevenson:2019ApJ. We assume
initially circular orbits BSE. We draw the orbital separation from a uniform
distribution in log space with a minimum of 0.01 and maximum of 1000 AU
Sana:2012Sci. Since we are interested in pulsars formed recently in the Milky
Way, we assume a representative value of Galactic metallicity
$Z_{\mathrm{Gal}}=0.014$ Asplund:2009ARA&A. A list of default COMPAS option
values, including mass transfer rate and common envelope evolution, can be
found in the link in footnote222https://github.com/FloorBroekgaarden/Double-
Compact-Object-
Mergers/blob/main/otherFiles/Table_with_detailed_binary__population_synthesis_simulation_settings.pdf.

In this work, neutron stars formed through accretion-induced supernovae, and
neutron stars undergoing any recycling/rejuvenation processes are omitted.
Only neutron stars formed through core collapse supernovae and electron-
capture supernovae that are not in wide, non-interacting binaries
Willcox:2021kbg are included in the following analysis. The reasoning for this
decision is discussed in § 2.2.

###  2.2 Neutron star kicks

Neutron stars receive natal kicks at birth Lyne:1994Nature; Hobbs:2005MNRAS;
Verbunt:2017A&A; Igoshev:2020MNRAS; Igoshev:2021MNRAS. We use the physical
remnant prescription from Mandel:2020qwb to relate the remnant mass and kick
of a neutron star to the properties of its progenitor. This prescription is
based on semi-analytic supernova models, and is calibrated to the Galactic
pulsar population Kapil:2022blf; specifically, we use $v_{\mathrm{NS}}=520$
km/s for the scaling pre-factor for neutron star kicks, and
$\sigma_{\mathrm{NS}}=0.3$ for the fractional stochastic scatter in the kick
velocities. Kapil:2022blf find the distribution of kicks in single and binary
neutron stars to be similar, so we opt to use a single distribution for both.
As discussed by Mandel:2022sxv, neutron star kicks inferred from young pulsars
may be overestimated by $\sim 15$ % if pulsar beams are narrow and
preferentially aligned with the kick direction. The Mandel:2020qwb
prescription is applied equally to both core-collapse supernovae and electron-
capture supernovae. We do not allow for electron-capture supernovae to occur
in wide, non-interacting binaries, to avoid overproducing low-speed pulsars
which are not observed Willcox:2021kbg.

###  2.3 COMPAS Implementation of Pulsar Evolution

This study focuses on canonical pulsars, hence we only describe the evolution
of isolated pulsars. Here we follow magnetic braking for an isolated pulsar to
spin down, with braking index $n=3$. The pulsar evolution theory implemented
in COMPAS is described in Chattopadhyay2020. We assume an exponentially
decaying magnetic field of a neutron star after it spins down for a time
duration of $t$

$B_{f}=(B_{i}-B_{\textrm{min}})\exp(-t/\tau)+B_{\textrm{min}},$ (2)

where $B_{i}$ is the magnetic field strength of the neutron star before it is
evolved by time duration $t$ and $B_{f}$ the magnetic field strength
afterwards, $B_{\textrm{min}}$ is the minimum magnetic field of a pulsar which
is set to be $10^{8}$ G for the entire analysis, and $\tau$ is the time scale
of pulsar magnetic field decay.

During the time duration $t$, the change in the spin period of the neutron
star is described as

$P^{2}=\frac{16\pi^{2}R^{6}\sin^{2}{\alpha}}{3c^{3}I}[B_{\textrm{min}}^{2}t-\tau
B_{\textrm{min}}(B_{f}-B_{i})-\frac{\tau}{2}(B^{2}_{f}-B^{2}_{i})]+P_{0}^{2},$
(3)

where $R$ is the radius of the pulsar (we assume $R=10$ km for all pulsars at
all times in this work), $P_{0}$ is the spin period of the pulsar at the
beginning of the evolved time duration $t$, and $B_{i}$ and $B_{f}$ are the
pulsar surface magnetic field at the beginning and end of the time step as
described in Eq. 2, $c$ is speed of light, $\alpha$ is the alignment angle
between the magnetic and rotational axes, and

$I=(0.237\pm
0.008)mR^{2}[1+4.2\frac{\textrm{M}}{\textrm{M}_{\odot}}\frac{\textrm{km}}{\textrm{R}}+90(\frac{\textrm{M}}{\textrm{M}_{\odot}}\frac{\textrm{km}}{\textrm{R}})^{4}]$
(4)

being the moment of inertia of the pulsar lattimer2005. Currently in COMPAS,
$\alpha=\pi/2$ is set to be constant and not evolved (though see
Johnston2020). Evolution of alignment will be included in a future version of
COMPAS.

The evolution of period derivative, $\dot{P}$ , can be described after
updating magnetic field strength and spin period, as

$\dot{P}=\frac{8\pi^{2}R^{6}B^{2}\sin^{2}{\alpha}}{3c^{3}IP},$ (5)

###  2.4 Free Parameters of COMPAS Runs

In addition to the default binary parameters, we need to propose distributions
of pulsar properties at birth for subsequent modelling. We note that some of
the past population studies (Faucher-Giguere:2006ApJ; Popov:2012Ap&SS;
Szary2014; Johnston2020; Chattopadhyay2020; Dirson2022; Igoshev:2022MNRAS,
e.g.,) proposed and refined distributions of pulsar properties at birth, such
as surface magnetic field strength and spin period. In this work, we also
choose the distributions of these two quantities at birth, as well as magnetic
decay timescale ($\tau_{d}$) to be free parameters. By default, COMPAS sets
the birth distribution of NS magnetic field as a uniform distribution between
$10^{10}$ G to $10^{13}$ G, and the birth distribution of NS spin period as a
uniform distribution between 10 to 100 ms. These distributions are in general
agreement with those observed from young pulsar populations Manchester:2004bp.
We exclude extremely high magnetic field strengths with B $>10^{15}$ G in the
magnetar regime as we are uncertain if the magnetar formation and evolution
physics would follow those of a normal NS. In our default model, we assume
that the magnetic field decay time scale $\tau_{d}=500$ Myr Chattopadhyay2021.
We choose $\tau_{d}$ as a free parameter since there is no clear consensus on
whether the magnetic field of a pulsar decays at all, or how long it takes to
decay. For example, semi-analytical studies 1970ApJ...160..979G;
1981JApA....2..315V; 1985MNRAS.213..613L found that the magnetic field decays
on the timescale of the order of a million years, with Dirson2022 recently
finding $\tau_{d}=4.6\times 10^{5}$ yr. In comparison, other works
1977ApJ...215..885T; 1987A&A...178..143S; Kiel:2008MNRAS have shown the
opposite, finding that $\tau_{d}>100$ Myr is needed in order for theoretical
models to match the observed catalog. For this reason, we let $\tau_{d}$ vary
from 10 Myr to 1000 Myr in our models.

We list the details of all the models with different pulsar birth
distributions and parameters in Table. 1. The initial model (INIT) is as
described above. Following that, we vary our model parameters in such ways
that we cover different ranges, and in some cases different shapes of the
distribution. The D1, D2 and D3 models, combined with INIT, are specifically
targeted to examine the impact of the $\tau_{d}$ on the results of the
simulations, as only $\tau_{d}$ is varied in models D1, D2 and D3 compared to
INIT. For B1 and B2 models, only the range of the birth magnetic field
distributions is changed compared to INIT. Similarly for P1, P2 and P3 models,
only the ranges of the birth spin period distributions are changed and the
others are kept the same. For PN1, PN2 and PN3 models, the birth spin period
distribution is changed to a normal distribution (hence, PN models as in
period-normal models) with various means and standard deviations. In the BLN
series of models, the birth magnetic field strength distribution is changed to
a log normal distribution, with different means and standard deviations in
models BLN1–4. The remaining three models, BLN1-PN1, BLN4-PN1, and BLN1-PN5,
have a combination of the magnetic field distribution from the BL models and
the spin period distributions of the PN models, respectively. We specifically
note that the original form of the log-normal distribution of the birth
magnetic field strength is taken from Szary2014. The BLN1-PN5 model is exactly
as described in Szary2014. In each COMPAS model, 300,000 binary systems are
modelled.

Table 1: Description of models used for COMPAS simulations in this work.  Model |  $B_{\textrm{birth}}$ Range (G) |  $B_{\textrm{birth}}$ Distribution |  $P_{\textrm{birth}}$ Range (ms) |  $P_{\textrm{birth}}$ Distribution |  $\tau_{d}$ (Myrs)
---|---|---|---|---|---
Initial | ($10^{10}$ \- $10^{13}$) | Uniform | (10 - 100) | Uniform | 500
D1 | ($10^{10}$ \- $10^{13}$) | Uniform | (10 - 100) | Uniform | 100
D2 | ($10^{10}$ \- $10^{13}$) | Uniform | (10 - 100) | Uniform | 10
D3 | ($10^{10}$ \- $10^{13}$) | Uniform | (10 - 100) | Uniform | 1000
B1 | ($10^{11}$ \- $10^{13}$) | Uniform | (10 - 100) | Uniform | 500
B2 | ($10^{10}$ \- $10^{12}$) | Uniform | (10 - 100) | Uniform | 500
P1 | ($10^{10}$ \- $10^{13}$) | Uniform | (10 - 1000) | Uniform | 500
P2 | ($10^{10}$ \- $10^{13}$) | Uniform | (1 - 100) | Uniform | 500
P3 | ($10^{10}$ \- $10^{13}$) | Uniform | (1 - 1000) | Uniform | 500
PN1 | ($10^{10}$ \- $10^{13}$) | Uniform |  $\mu=75$,$\sigma=25$ | Normal | 500
PN2 | ($10^{10}$ \- $10^{13}$) | Uniform |  $\mu=50$,$\sigma=25$ | Normal | 500
PN3 | ($10^{10}$ \- $10^{13}$) | Uniform |  $\mu=50$,$\sigma=50$ | Normal | 500
BLN1 | $\mu=10^{12.6},\sigma=10^{0.55}$ | LogNormal | (10 - 100) | Uniform | 500
BLN2 | $\mu=10^{12.0},\sigma=10^{0.55}$ | LogNormal | (10 - 100) | Uniform | 500
BLN3 | $\mu=10^{12.6},\sigma=10^{0.25}$ | LogNormal | (10 - 100) | Uniform | 500
BLN4 | $\mu=10^{12.0},\sigma=10^{0.25}$ | LogNormal | (10 - 100) | Uniform | 500
PN1-B1 | ($10^{11}$ \- $10^{13}$) | Uniform |  $\mu=75$,$\sigma=25$ | Normal | 500
BLN1-PN1 | $\mu=10^{12.6},\sigma=10^{0.55}$ | LogNormal |  $\mu=75$,$\sigma=25$ | Normal | 500
BLN4-PN1 | $\mu=10^{12.0},\sigma=10^{0.25}$ | LogNormal |  $\mu=75$,$\sigma=25$ | Normal | 500
BLN1-PN5 | $\mu=10^{12.6},\sigma=10^{0.55}$ | LogNormal |  $\mu=300$,$\sigma=150$ | Normal | 500

###  2.5 Binary Birth, Pulsar Death, & Parameters

COMPAS begins evolving each binary system at $T=0$ Myr. The first stage of
post-processing is to assign a birth time of the binary system that
corresponds to the Milky Way star formation history. We consider a uniform
star formation history of the Milky Way gomez2018, enabling us to draw the
birth time of the binary systems from a uniform distribution. Without
recycling, pulsars rapidly spin down, losing their rotational energy to the
point that they can no longer sustain detectable radio or $\gamma$-ray
emission, at which point we define this as the death of pulsar. The typical
characteristic age of a canonical $\gamma$-ray pulsar is relatively small,
roughly $10^{5}$ yrs Fermi-LAT:2023zzt. For this reason, we assume that all
the neutron stars from the simulation that would potentially be $\gamma$-ray
emitting are formed in binary systems that are born in the recent history of
the Milky Way. The birth time of any given binary system is then drawn from a
uniform distribution with a minimum of 12.8 Gyr and a maximum of 13.0 Gyr. We
choose to neglect the population of pulsars produced through accretion-induced
collapse of white dwarfs as they represent only a small fraction of the
population ($\lesssim 10\%$), are formed on much longer timescales and may
potentially be born with different properties to those born in core-collapse
supernovae Bailyn:1990ApJ; Tauris:2013A&A.

At this point, the time stamp of the binary system in COMPAS output
corresponds to an actual point of time in the history of the Milky Way. We
make the observation of any neutron star at the present day, assuming the age
of the Milky Way is 13.0 Gyr. If the time stamp is not exactly at 13 Gyr, then
a linear interpolation of the two nearest data points is used as the recorded
parameter value. COMPAS provides $P$, $\dot{P}$ , $B$ and mass of the neutron
star at the observation, and one can estimate the spin-down luminosity via
Equation 1.

Before recording the parameters of the neutron star, it is essential to check
if the pulsar crosses the radio death-line rudakritter1994 on the
$P$-$\dot{P}$ diagram, which is defined by

$\log_{10}\dot{P}=3.29\times\log_{10}P-16.55,$ (6)

and

$\log_{10}\dot{P}=0.92\times\log_{10}P-18.65,$ (7)

If $\dot{P}$ of a neutron star from COMPAS is larger than both of these two
calculated values according to Eqs. 6, 7, then the pulsar has not yet crossed
the death-line and may be observable, and we record its parameters. Our
treatment of radio selection effects is given in Section 2.7.

###  2.6 Dynamical Evolution

Massive stars, which in later evolutionary stages form neutron stars, are
formed within an exponential disk around the Galactic plane. We assume that
these massive stars move along their Galactic orbits without any interruption
until the change in the evolutionary status of one of the stars in the binary
system. As described in § 2.1, we assume all the neutron stars are formed in
supernovae. Upon formation, the stellar remnant receives a natal kick as
described in § 2.2 which can disrupt the binary system, and in some cases, end
up with two individual stars.

Upon the supernova event that forms the neutron star, the binary system is
then assigned a random location in the Galaxy following a distribution of a
thin, exponential disk Paczynski1990:

$p_{z}(z)dz=\exp{-z/z_{\textrm{exp}}}\frac{dz}{z_{\textrm{exp}}},$ (8)

where $z$ is the vertical distance from the Galactic plane, and

$p_{R}(R)dR=a_{R}\exp{-R/R_{\textrm{exp}}}\frac{R}{R^{2}_{\textrm{exp}}}dR,$
(9)

where $a_{R}=1.0683$ and $R$ is the distance from the Galactic center on the
Galactic plane. We set $z_{\textrm{exp}}=75$ pc and $R_{\textrm{exp}}=4.5$ kpc
following Paczynski1990.

Knowing the location of the newly born neutron star, the time at which it is
born, and both the magnitude and direction of the natal kick velocity, we can
follow its subsequent motion in the Galactic potential. We use NIGO nigo2015,
a numerical galactic orbit integration tool, to evolve the motion of binary
systems after the supernova within the Galactic potential. NIGO describes the
Milky Way as a disc galaxy, with the Galactic potential comprised of three
components: a bulge, a disk and a dark matter halo. The Galactic bulge is
described as a Plummer sphere flynn1996; plummer1911:

$\Phi_{b}=-\frac{GM_{b}}{\sqrt{R^{2}+z^{2}+b_{b}^{2}}},$ (10)

where $G$ is the gravitation constant, $M_{b}=1.0\times 10^{10}$M⊙ is the mass
of the bulge and $b_{b}=0.4$kpc is the scale length and $R^{2}=x^{2}+y^{2}$.

The Galactic disc is modelled by an exponential disc formed by the
superposition of three Miyamoto-Nagai potentials MN1975:

$\Phi_{d}=-\sum_{n=1}^{3}\frac{GM_{d_{n}}}{\sqrt{R^{2}+[a_{d_{n}}+\sqrt{b_{d}^{2}+z^{2}}]^{2}}},$
(11)

where $M_{d_{n}}$ are the masses of each disk, $a_{d_{n}}$ are related to the
disk scale lengths of the three disk components (which are listed in Table. 2)
and $b_{d}=0.2376$ is related to the disk scale height.

Table 2: Parameters used for the three components of exponential disk | 1 | 2 | 3
---|---|---|---
$M_{d}$ (M⊙) | $9.0335\times 10^{10}$ | $-5.9112\times 10^{10}$ | $1.0590\times 10^{10}$
$a_{d}$ (kpc) | 4.5081 | 9.6483 | 1.5136

An NFW dark matter halo NFW1997 is used in the modelling of the Galactic
potential:

$\Phi_{h}=-\frac{GM_{h}}{r}\ln{(1+\frac{r}{a_{h}})},$ (12)

where $M_{h}=3.3\times 10^{12}$M⊙ is the mass of the halo component,
$a_{h}=45.02$ kpc is the length scale and $r^{2}=R^{2}+z^{2}$.

###  2.7 Radio Selection Effects

The observed radio pulsar population represents only a small fraction of the
total pulsar population in the Milky Way, and may be strongly biased due to
observational selection effects. In this section we describe our modelling of
the relevant radio selection effects.

####  2.7.1 Radio Luminosity

One of the primary factors determining whether a pulsar is observable is its
radio luminosity. Past population studies have shown that the radio luminosity
of a pulsar has no dependency on $P$ or $\dot{P}$ of the pulsar itself
Szary2014. Instead, in the unit of mJy kpc2 in 1400 MHz band, it follows a
log-normal distribution Szary2014

$\log L_{\textrm{1400}}\sim N(0.5,1.0),-3.0\leq\log L_{\textrm{1400}}\leq 4.0$
(13)

If a pulsar crosses the “death-line”, described in § 2.5, before the
observation is made, then the pulsar is not going to be observed.
Additionally, another similar criterion is to examine the radio emission
efficiency Szary2014 of a pulsar, defined as the fraction of rotational energy
transformed into radio emission

$\xi=\frac{L_{\textrm{r}}}{\dot{E}},$ (14)

where $L_{\textrm{r}}(\textrm{erg s}^{-1})\approx 7.4\times
10^{27}L_{\textrm{r}}$ (mJy). As the pulsar ages, its radio emission
efficiency will increase. For a pulsar to be detectable in radio, $\xi<0.01$
must be satisfied.

####  2.7.2 PSREvolve: Radio Sensitivity

To determine if a pulsar can be observed by a given radio survey, we use the
PSREvolve code Oslowski2021; Chattopadhyay2020 to calculate the radiometer
equation dewey1985; lm2004 and determine the minimum detectable flux,
$S_{\textrm{min}}$, at the sky location at a given signal-to-noise ratio
$(S/N_{\textrm{min}})$

$S_{\textrm{min}}=\beta\frac{(S/N_{\textrm{min}})(T_{\textrm{rec}}+T_{\textrm{sky}})}{G\sqrt{n_{p}t_{\textrm{int}}\Delta
f}}\sqrt{\frac{W_{e}}{P-W_{e}}},$ (15)

where $T_{\textrm{rec}}$ and $T_{\textrm{sky}}$ are the temperature of the
receiver and the direction in the sky, respectively; $\Delta f$ is the
bandwidth of the receiver, $G$ is the gain of the telescope, $n_{p}$ is the
number of polarisations in the detector, $W_{e}$ is pulse width and $P$ is the
period of the pulsar, and $\beta$ is a noise correction factor that increases
the noise with digitisation errors and bandpass distortion.

Pulsar detection has been carried out by almost every major radio facility in
the world, over a period of more than 50 years, using a variety of different
backends. This makes it complicated to individually model each pulsar survey
that has been conducted. However, one can focus on modelling some of the
largest, most successful surveys (Faucher-Giguere:2006ApJ; Graber:2023arXiv,
e.g.,). In this work, we choose pulsar survey parameters for Equation 15 to
represent the Parkes Multibeam Pulsar Survey manchester2001;
Lorimer:2006MNRAS; Sengar:2023MNRAS by setting $\beta=1$ and
$S/N_{\textrm{min}}\geq 10$. This survey discovered over half of the pulsars
detected. Additionally, when comparing the properties of radio pulsars
observed by Parkes multibeam surveys and the entire radio pulsar population,
the difference in their distributions is negligible, as shown in Fig. 1. For
these reasons, we assume that the radio selection is based on a “global”
Parkes survey that covers the entire sky, and we use the same specifications
as Parkes.

Figure 1: Cumulative distribution functions (CDFs) of $\dot{E}$ , radio flux,
period and $\dot{P}$ of all radio pulsars in the catalogs (red) and pulsars
observed in the Parkes Multibeam Surveys (blue).

The expression of $W_{e}$ below 2003Natur.426..531B describes how the ISM
impacts the observed pulses. Free electrons in the Galaxy broaden the
intrinsic pulses $W_{i}$ 2002astro.ph..7156C, and the interstellar medium
scatters the pulsar beam.

$W_{e}^{2}=W_{i}^{2}+\tau_{\textrm{samp}}^{2}+(\tau_{\textrm{samp}}\frac{\textrm{DM}}{\textrm{DM}_{0}})^{2}+\tau_{\textrm{scatt}}^{2},$
(16)

where $W_{i}$ is the intrinsic width of the pulse, $\tau_{\textrm{samp}}$ and
$\tau_{\textrm{scatt}}$ sampling and ISM scattering times, DM is the
dispersion measure in the direction of the pulsar and DM0 is the diagonal
dispersion measure of the survey. While the duty cycle of pulsars varies
widely, here we assume $\frac{W_{i}}{P}=0.05$ for all pulsars. If the radio
luminosity of a pulsar is larger than $S_{\textrm{min}}d^{2}$, then it can be
detected by the chosen radio survey.

####  2.7.3 Radio beaming

Pulsars are considered to be observable during the pulsations, which are
strongly beamed. A pulsar should only be observable when the beamed pulsation
crosses the line of sight from the Earth. To this end, we should consider the
beaming fraction of pulsars in radio, which depends on different emission
mechanisms and geometry of the pulsar.

The radio beaming fraction of a pulsar, $f_{r}$ Tauris1998, is considered to
be

$f_{r}=0.09\log{(P/10)}^{2}+0.03,$

where $P$ is the rotational period of the pulsar in seconds. We use a
probabilistic rejection scheme to use the beaming direction as an additional
selection criterion to determine the detectability of pulsars in radio and
$\gamma$-ray assuming a uniformly distributed viewing angle.

###  2.8 Gamma-ray Selection Effects

####  2.8.1 Gamma-ray Luminosity & Flux

An empirical description of the relationship between the $\gamma$-ray
luminosity of a pulsar and its physical properties as a fundamental plane has
been proposed Pulsar_FP2019; Pulsar_FP2022. In this work we estimate pulsar
$\gamma$-ray luminosity using a fundamental plane in four-dimensional (4D)
space that is fitted to 4FGL-DR3 pulsars Pulsar_FP2022,

$L_{\gamma}=10^{14.3\pm 1.3}E_{\textrm{cut}}^{1.39\pm 0.17}B^{0.12\pm
0.03}\dot{E}^{0.39\pm 0.05},$ (17)

where $E_{\textrm{cut}}$ is the cutoff energy of a pulsar’s $\gamma$-ray SED
when modelled by a power-law with the exponential cutoff model, $B$ is the
surface magnetic field strength, and $\dot{E}$ is the spin-down luminosity.
Pulsar $\gamma$-ray spectral parameters, cutoff energy $E_{\textrm{cut}}$ and
power-law index before cutoff $\Gamma$, are determined following Figs. 10 & 11
in Pulsar_FP2022. These two figures show correlations between
$E_{\textrm{cut}}$ and $\dot{E}$ and $\Gamma$ and $\dot{E}$ respectively when
the pulsar has $10^{33}$ erg/s $<$ $\dot{E}$ $<10^{38}$ erg/s. Within this
$\dot{E}$ range, we assume a power-law correlation of $E_{\rm cut}$ (in GeV)
with $\dot{E}$ as well as $\Gamma$ with $\dot{E}$ with these following
correlations

$\log_{10}E_{\rm cut}=\log_{10}0.5+0.26(\log_{10}\dot{E}-33)+U(0.05),$

and

$\Gamma=0.5+0.3(\log_{10}\dot{E}-33)+U(0.05),$

where $U(0.05)$ is a random number drawn from a uniform distribution between
-0.05 and 0.05. When $\dot{E}$ $<10^{33}$ erg/s, we set

$\log_{10}E_{\rm cut}=\log_{10}0.5+U(0.05),$

and

$\Gamma=0.5+U(0.05).$

When $\dot{E}$ $>10^{38}$ erg/s, $\log_{10}E_{\rm cut}=1.0+U(0.05)$,
$\Gamma=0.5+U(0.05)$.

We can then determine the pulsar’s energy flux by

$F_{\gamma}=\frac{L_{\gamma}}{4\pi d^{2}c_{g}},$

where $d$ is the distance of the pulsar to the Earth and $c_{g}$ is the flux
correction factor Watters2009, which depends on the emission mechanism of the
pulsar as well as geometry and viewing angle. Here we combine results from
Johnston2020 and Petri2011, and assume that $c_{g}$ follows a sigmoid function
versus $\dot{E}$ with a maximum of 0.8 Johnston2020 and a minimum of 0.22
Petri2011, while turning over at $10^{33}$erg s-1 with an added randomisation
drawn from a uniform distribution between -0.1 and 0.1.

####  2.8.2 Gamma-ray Sensitivity

The 12-year Fermi-LAT Source Catalog provides an all-sky map of the detection
threshold for the
catalog333https://fermi.gsfc.nasa.gov/ssc/data/access/lat/12yr_catalog/detthresh_P8R3_12years_PL22.fits.
By comparing the energy flux of a pulsar in the simulation to the detection
threshold at the corresponding sky location, we would be able to determine if
the simulated pulsar would be detectable or not.

In reality, pulsars that are already detected in radio surveys first are
easier to find in $\gamma$-rays when compared to $\gamma$-ray blind searches
Fermi-LAT:2013svs; Fermi-LAT:2023zzt. Fig. 17 of Fermi-LAT:2013svs clearly
indicates that for the detected pulsars in the 2nd Fermi-LAT pulsar catalog,
pulsars that are already detected in radio in general are easier to detect.
For this reason, we make the detection threshold five times lower at the
location of the pulsar if the simulated pulsar is detected in radio selection;
and 1.5 times higher at the location of the pulsar if the simulated pulsar is
not detected in radio.

####  2.8.3 Gamma-ray beaming

Gamma-ray pulsars have a beaming effect similar to that of radio pulsars. In
Johnston2020, a description of young and energetic (with $\dot{E}$ $>10^{35}$
erg/s) pulsars’ $\gamma$-ray beaming fraction ($f_{g}$) is in their Table. 3.
We adopt this description in determining the beaming fractions of our modelled
pulsars. As this work does not provide a beaming fraction for pulsars with
$\dot{E}$ $<10^{35}$ erg/s, here we assume that as $\dot{E}$ decreases, the
beaming fraction also decreases. The choices of beaming fraction can be
described as

$f_{g}=\begin{cases}0.92,&\text{if }\dot{E}>10^{38}\ {\rm erg\ s}^{-1}\\\
0.82,&\text{if }10^{37}<\dot{E}\leq 10^{38}\ {\rm erg\ s}^{-1}\\\
0.67,&\text{if }10^{36}<\dot{E}\leq 10^{37}\ {\rm erg\ s}^{-1}\\\
0.50,&\text{if }10^{35}<\dot{E}\leq 10^{36}\ {\rm erg\ s}^{-1}\\\
0.40,&\text{if }10^{33}<\dot{E}\leq 10^{35}\ {\rm erg\ s}^{-1}\\\
0.30,&\text{if }\dot{E}\leq 10^{33}\ {\rm erg\ s}^{-1}\end{cases}$

The assumptions above are made based on observations and outer gap model of
young, powerful pulsars Johnston2020. According to 2012A&A...545A..42P, this
choice can be dependent on the pulsar emission model and is consistent with
what we have chosen.

We also use a probabilistic rejection scheme same as that described in § 2.7.3
similarly for radio pulsars. If beaming is considered in the chosen emission
mechanism of the $\gamma$-ray pulsar in the simulation, then it is used as an
additional selection criterion to determine the detectability of this pulsar
in $\gamma$-ray , while also assuming a uniformly distributed viewing angle.

###  2.9 Reusing COMPAS Output

Depending on the model parameters, roughly 30% to 50% of binary systems would
produce at least one neutron star in each COMPAS simulation. However, due to
the selection effects introduced above, only an extremely small portion of
these simulated neutron stars are recorded as observable. Adhering to the
principle of effectively utilising the simulation, we develop a scheme to
reuse each model. Each time the simulation is re-used, the birth time and
location of a neutron star is re-assigned. We then make the star travel
through the Galactic potential described above, and make the observation again
at 13 Gyr, with re-assigned radio and $\gamma$-ray luminosity. This re-using
scheme is applied 9 times to each simulation, and along with the original
simulation, it is equivalent to 3,000,000 binary systems simulated in each
COMPAS model.

##  3 Results & Discussions

The results from the COMPAS simulations using the setups from § 2 are
presented in this section, as well as the discussion on the implications of
these results. The results are compared to the pulsar catalog to determine the
goodness of the different models. This will help us determine the model of
pulsar initial conditions that best fits both the radio and $\gamma$-ray
observations. We will first show the results of the simulation by comparing
the number of different types of pulsars produced in different simulation
models, then employ a statistical method to compare the distribution of pulsar
physical properties to those reported in the catalogs. We further discuss the
impact of the results from the selection effects of $\gamma$-ray pulsars. We
also aim to provide some insights into the $\gamma$-ray emission of low
$\dot{E}$ pulsars.

###  3.1 Star Formation Rate

In COMPAS simulations produced for this study, only stars within the initial
mass range of 5 M⊙ to 150 M⊙ are modelled. Because of this setup, we need to
account for the low-mass end of the initial mass function (IMF) when
calculating the star formation rate. Since all the models in this study use
the same parameters for binary population, they all have the same star
formation rate. Using the CosmicIntegration Python package in COMPAS, we are
able to account for the lower end of the IMF. By default, we use a
Kroupa:2000iv broken power-law IMF, with a slope of 0.3 between 0.01 M⊙ and
0.08 M⊙ ; slope of 1.3 between 0.08 M⊙ and 0.5 M⊙ ; slope of 2.3 between 0.5
M⊙ and 200 M⊙ . Each binary system generated in any given COMPAS run
corresponds to 234 M⊙ of star formation. Considering we set all the binary
systems to be born in the most recent 200 Myr, and each COMPAS output is
sampled 10 times, the star formation rate is $\sim 1.2$ M⊙ per year with the
current model selection. The observed star formation rate over the most recent
history of the Milky Way is around $1.65\pm 0.19$ 2015ApJ...806...96L or
$1.9\pm 0.4$ M⊙yr-1 2011AJ....142..197C. Other works, such as
2022ApJ...929L..18E, estimate the Milky Way star formation rate in the range
of 0.50 - 5.93 M⊙ yr-1. The star formation rate from our COMPAS simulations is
consistent with these observed and modelled values.

###  3.2 Results

The most straightforward and obvious way to determine if a COMPAS model is a
good fit for the observation is to compare the number of pulsars produced by
the model to different catalogs. We examine, for each model, the number of
radio pulsars, the number of $\gamma$-ray pulsars, the number of radio-quiet
$\gamma$-ray loud pulsars, and the number of radio-loud $\gamma$-ray quiet
pulsars, produced in the simulation. To fully eliminate the possible
randomness from the models, the observation part of post-processing is
repeated 10 times, and the average and standard deviation of these
realisations are recorded for each model in Table. 3.

Table 3: Results of the simulations: number of different types of pulsars. Columns left to right: COMPAS model; pulsars that are detected in radio ($N_{\rm R}$); pulsars detected in $\gamma$-ray ($N_{\rm G}$); radio-quiet $\gamma$-ray loud pulsars ($N_{\rm RQGL}$); radio-loud $\gamma$-ray quiet pulsars ($N_{\rm RLGQ}$); radio-loud $\gamma$-ray pulsars ($N_{\rm RLGL}$). The first row lists the total number of pulsars in the radio and $\gamma$-ray catalogs Manchester:2004bp; Fermi-LAT:2023zzt. MSPs are excluded as we focus on canonical pulsars. For the simulations, the values and uncertainties are obtained from the average and standard deviation of 10 different realisations of each model.  Model | NR | NG | NRQGL | NRLGQ | NRLGL
---|---|---|---|---|---
Catalogs | 2715 | 160 | 90 | 2645 | 70
INIT | $3414\pm 46$ | $195\pm 13$ | $113\pm 6$ | $3331\pm 47$ | $82\pm 10$
D1 | $3412\pm 55$ | $188\pm 7$ | $103\pm 8$ | $3327\pm 53$ | $85\pm 7$
D2 | $2533\pm 25$ | $208\pm 8$ | $127\pm 6$ | $2451\pm 29$ | $81\pm 9$
D3 | $3357\pm 40$ | $194\pm 6$ | $114\pm 6$ | $3277\pm 40$ | $80\pm 7$
B1 | $2803\pm 56$ | $185\pm 10$ | $108\pm 4$ | $2726\pm 57$ | $77\pm 8$
B2 | $20466\pm 89$ | $612\pm 18$ | $199\pm 7$ | $20053\pm 87$ | $413\pm 14$
P1 | $1953\pm 34$ | $41\pm 4$ | $22\pm 2$ | $1935\pm 35$ | $18\pm 2$
P2 | $3366\pm 45$ | $210\pm 10$ | $112\pm 5$ | $3269\pm 47$ | $98\pm 8$
P3 | $1993\pm 26$ | $35\pm 4$ | $18\pm 3$ | $1976\pm 27$ | $16\pm 3$
PN1 | $3263\pm 47$ | $138\pm 9$ | $75\pm 4$ | $3200\pm 45$ | $62\pm 6$
PN2 | $3436\pm 64$ | $190\pm 9$ | $109\pm 6$ | $3356\pm 60$ | $81\pm 7$
PN3 | $3337\pm 59$ | $188\pm 8$ | $102\pm 7$ | $3251\pm 53$ | $87\pm 8$
BLN1 | $3434\pm 58$ | $184\pm 10$ | $106\pm 7$ | $3356\pm 58$ | $78\pm 6$
BLN2 | $12020\pm 93$ | $451\pm 14$ | $179\pm 12$ | $11748\pm 93$ | $272\pm 13$
BLN3 | $1773\pm 28$ | $185\pm 8$ | $124\pm 7$ | $1711\pm 26$ | $62\pm 7$
BLN4 | $7918\pm 74$ | $373\pm 13$ | $177\pm 7$ | $7723\pm 70$ | $196\pm 10$
PN1-B1 | $2736\pm 27$ | $146\pm 7$ | $79\pm 5$ | $2669\pm 26$ | $67\pm 6$
BLN1-PN1 | $3446\pm 59$ | $152\pm 5$ | $75\pm 4$ | $3369\pm 59$ | $77\pm 5$
BLN1-PN5 | $2586\pm 34$ | $52\pm 5$ | $25\pm 3$ | $2559\pm 33$ | $27\pm 5$
BLN4-PN1 | $7727\pm 79$ | $262\pm 13$ | $96\pm 4$ | $7562\pm 80$ | $166\pm 13$

We then tabulate the difference of numbers of pulsars in Table. 4. For many of
the models, the predicted number of pulsars is largely different from the
observed pulsars. For example, models INIT and D1 have 26% differences in
radio pulsars, and 21% and 18% differences in $\gamma$-ray pulsars,
respectively. An even more extreme case is B2 model, with 7.5 times of radio
pulsars and 3.8 times of $\gamma$-ray pulsars predicted compared to the
observed populations. Due to these large differences, we scale the predicted
numbers according to the radio pulsar population. For each model, the numbers
of all types of pulsars are scaled by the factor ${\rm N_{\rm R}}/2715$, with
2715 is the number of observed radio pulsars from the catalogs. These scaled
numbers of $\gamma$-ray pulsar, radio-quiet and radio-loud $\gamma$-ray
pulsars, and radio-loud $\gamma$-ray-quiet pulsars are then compared to those
from the catalogs to calculate the percent difference. The B1 and PN1-B1
models are not scaled as the predicted numbers of radio pulsars are already
close to the observed value. Scaling the number of pulsars in a model
effectively changes the modelled star formation rate. We scale the default
star formation rate of 1.2M⊙ /yr, as described in § 3.1, by the same scaling
factor for each model, to obtain the appropriate number as listed in the last
column in Table. 4. Another indicator of the goodness of the models is the
ratio between different types of pulsars. Here we check the ratio between
radio and $\gamma$-ray pulsars, and ratio between radio-quiet and radio-loud
$\gamma$-ray pulsars for each model, and compare them to those from the
catalogs. The percent differences between these predicted and observed ratios
are also listed in Table. 4. Both the INIT and PN1-B1 models are able to
reproduce the observed numbers and ratios of different types of pulsars to
within $\approx 10\%$.

Table 4:  Percentage difference (PD) between the numbers of different types of pulsars generated in COMPAS models and the radio and $\gamma$-ray catalogs. For all the models, except for B1 and PN1-B1, the comparisons are carried out after scaling the number of radio pulsars to be the same as the catalog. First five columns: percent differences of the numbers of radio, $\gamma$-ray , radio-quiet $\gamma$-ray-loud, radio-loud $\gamma$-ray-quiet, and radio-loud $\gamma$-ray-loud pulsars between each model and the catalogs. Sixth column: percent difference of the ratio between radio and $\gamma$-ray pulsars between each model and the catalogs. Seventh column: percent difference of the ratio between radio and $\gamma$-ray pulsars between each model and the catalogs. Last column: star formation rate of each model after scaling by the number of radio pulsars.  Model | PDR | PDG | PDRQGL | PDRLGQ | PDRLGL | PDQ1 | PDQ2 | SFR (M⊙/yr)
---|---|---|---|---|---|---|---|---
INIT | - | 3.1% | 0.0% | 0.11% | 7.14% | 7.01% | 1.2% | 0.95
D1 | - | 6.25% | 8.9% | 3.97% | 2.8% | 6.95% | 5.5% | 0.95
D2 | - | 39.4% | 51.1% | 0.68% | 24.3% | 28.2% | 22.6% | 1.3
D3 | - | 1.9% | 2.2% | 0.2% | 7.1% | 1.9% | 11.3% | 0.97
B1∗ | 3.24% | 15.6% | 20% | 3.1% | 10% | 10.72% | 9.37% | 1.2
B2 | - | 49.3% | 71.1% | 0.6% | 21.4% | 97% | 62.3% | 0.16
P1 | - | 64.4% | 65.6% | 1.7% | 64.3% | 180.5% | 4.7% | 1.67
P2 | - | 5.6% | 0.0% | 0.3% | 12.8 | 5.6% | 11% | 0.97
P3 | - | 70% | 72.2% | 1.7% | 68.6% | 235.3% | 12.5% | 1.6
PN1 | - | 28% | 31% | 0.64% | 27% | 39.30% | 5.47% | 1.0
PN2 | - | 6.9% | 5% | 0.26% | 9.3% | 6.54% | 4.69% | 0.95
PN3 | - | 5% | 7.8% | 0.0% | 1.4% | 4.6% | 8.6% | 0.98
BLN1 | - | 9.4% | 6.7% | 0.3% | 11.4% | 9.9% | 6.2% | 0.95
BLN2 | - | 36.9% | 55.5% | 0.3% | 12.8% | 57% | 56.2% | 0.27
BLN3 | - | 76.9% | 111.1% | 0.9% | 35.7% | 43.5% | 56.2% | 1.84
BLN4 | - | 20% | 32.2% | 0.1% | 3.4% | 24.9% | 29.7% | 0.41
PN1-B1∗ | 0.7% | 8.75% | 12.2% | 0.9% | 4.3% | 10.4% | 7.8% | 1.2
BLN1-PN1 | - | 25% | 34.4% | 0.3% | 12.8% | 33.6% | 24.2% | 0.94
BLN1-PN5 | - | 65.6% | 71.1% | 1.5% | 60% | 193% | 28.1% | 1.26
BLN4-PN1 | - | 42.5% | 62.2% | 0.03% | 17.1% | 73.8% | 54.7% | 0.42
∗ models not scaled.

Beyond simply comparing the total number of pulsars between observations and
our simulations, we can also compare the predicted and observed distributions
of various physical quantities of the pulsar populations between the catalog
and the simulations. We first plot the cumulative distribution functions
(CDFs) for each pulsar parameter for each of the models listed in Table 1 in
Figure 2, and compare them to the observed distributions. In general, some of
the COMPAS model results produce distributions in most physical quantities
that match the observations. But there are a few quantities that cannot be
reproduced by any of the COMPAS models (e.g., $\gamma$-ray flux).

Figure 2: CDFs of physical quantities produced by the COMPAS models described
in Table. 1 (coloured), compared to the data from the catalogs (black). CDFs
from a given model are drawn from a randomised sampling from 10 realisations.
From top left to bottom right, the panels show the distributions of
$\gamma$-ray pulsar $\dot{E}$, radio pulsar $\dot{E}$, $\gamma$-ray flux,
radio flux, $\gamma$-ray pulsar period, radio pulsar period, $\gamma$-ray
pulsar $\dot{P}$, radio pulsar $\dot{P}$, $\gamma$-ray pulsar galactic
latitude, radio pulsar galactic latitude, $\gamma$-ray pulsar galactic
longitude and radio pulsar galactic longitude.

Next, we use statistical tests to examine how well the predicted parameter
distributions match the observed distributions for each physical quantity.
Given that the simulation and catalog have different-sized samples on each
quantity, and the quantity of each pulsar should be independent of each other,
we use Mann-Whitney U test (MWUtest, MWU test;) to make the comparison. The
MWU test is chosen because it is a nonparametric test independent of
underlying distributions. The null hypothesis is that both the predicted and
observed samples are drawn from the same underlying distribution. Given the
model and simulation have different size, and there is no certain theory on
the shape of the distribution of each parameter, it is important to use a
statistical test that is independent of the distribution itself. We record the
$p$-value of the test for a given physical quantity between the comparison of
simulation and catalog. If $p<0.05$, it means we reject the null hypothesis
that the two compared samples are drawn from the same distribution. When
$p>0.05$, we cannot reject the null hypothesis, which means we consider the
two distributions to be the same.

Table 5: P-value of Mann-Whitney U Test, comparing the physical quantities obtained from simulations and catalog. Quantities listed in the second to the third last columns are, respectively, radio flux, $\gamma$-ray flux, period of radio pulsars, period of $\gamma$-ray pulsars, $\dot{P}$ of radio pulsars, $\dot{P}$ of $\gamma$-ray pulsars, $\dot{E}$ of radio pulsars, $\dot{E}$ of $\gamma$-ray pulsars, Galactic latitude of radio pulsars, Galactic latitude of $\gamma$-ray pulsars, Galactic longitude of radio pulsars, and Galactic longitude of $\gamma$-ray pulsars. The last column is the number of the quantities that have p-values $>0.05$ for each model. For a given model, the p-values for each parameter is determined using the random realisation of the model that is shown in Fig. 2 Model |  $\dot{E}$(G) |  $\dot{E}$(R) | Flux (G) | Flux(R) | Period (G) | Period(R) |  $\dot{P}$(G) |  $\dot{P}$(R) | b(G) | b(R) | l(G) | l(R) | NO
---|---|---|---|---|---|---|---|---|---|---|---|---|---
INIT | 0.059 | 0.09 | 0.0 | 0.0 | 0.001 | 0.0 | 0.999 | 0.0 | 0.767 | 0.487 | 0.144 | 0.052 | 7
D1 | 0.58 | 0.526 | 0.0 | 0.0 | 0.0 | 0.0 | 0.443 | 0.0 | 0.085 | 0.541 | 0.048 | 0.011 | 5
D2 | 0.0 | 0.0 | 0.0 | 0.0 | 0.001 | 0.0 | 0.05 | 0.0 | 0.548 | 0.728 | 0.426 | 0.017 | 4
D3 | 0.113 | 0.528 | 0.0 | 0.0 | 0.003 | 0.0 | 0.644 | 0.0 | 0.461 | 0.085 | 0.112 | 0.047 | 6
B1 | 0.019 | 0.094 | 0.0 | 0.0 | 0.003 | 0.4 | 0.266 | 0.0 | 0.117 | 0.579 | 0.041 | 0.006 | 5
B2 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.845 | 0.572 | 0.058 | 0.064 | 4
P1 | 0.524 | 0.0 | 0.0 | 0.0 | 0.161 | 0.0 | 0.497 | 0.0 | 0.819 | 0.515 | 0.988 | 0.003 | 6
P2 | 0.758 | 0.363 | 0.0 | 0.0 | 0.0 | 0.0 | 0.241 | 0.0 | 0.242 | 0.893 | 0.388 | 0.073 | 7
P3 | 0.097 | 0.0 | 0.0 | 0.0 | 0.029 | 0.0 | 0.646 | 0.0 | 0.04 | 0.405 | 0.003 | 0.124 | 4
PN1 | 0.379 | 0.003 | 0.0 | 0.0 | 0.66 | 0.0 | 0.181 | 0.0 | 0.405 | 0.074 | 0.267 | 0.013 | 6
PN2 | 0.312 | 0.09 | 0.0 | 0.0 | 0.006 | 0.0 | 0.295 | 0.0 | 0.315 | 0.589 | 0.019 | 0.008 | 5
PN3 | 0.193 | 0.858 | 0.0 | 0.0 | 0.008 | 0.0 | 0.031 | 0.0 | 0.159 | 0.265 | 0.624 | 0.031 | 5
BLN1 | 0.191 | 0.022 | 0.0 | 0.0 | 0.007 | 0.005 | 0.048 | 0.0 | 0.868 | 0.047 | 0.53 | 0.084 | 4
BLN2 | 0.038 | 0.042 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.268 | 0.035 | 0.014 | 0.005 | 1
BLN3 | 0.006 | 0.0 | 0.0 | 0.0 | 0.908 | 0.0 | 0.0 | 0.0 | 0.292 | 0.876 | 0.121 | 0.107 | 5
BLN4 | 0.84 | 0.029 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.535 | 0.159 | 0.102 | 0.019 | 4
PN1-B1 | 0.618 | 0.004 | 0.0 | 0.0 | 0.687 | 0.266 | 0.065 | 0.0 | 0.296 | 0.836 | 0.189 | 0.007 | 7
BLN1-PN1 | 0.626 | 0.002 | 0.0 | 0.0 | 0.54 | 0.001 | 0.017 | 0.0 | 0.476 | 0.109 | 0.151 | 0.034 | 5
BLN1-PN5 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.759 | 0.0 | 0.654 | 0.178 | 0.283 | 0.001 | 4
BLN4-PN1 | 0.0 | 0.002 | 0.0 | 0.0 | 0.102 | 0.0 | 0.0 | 0.0 | 0.224 | 0.086 | 0.031 | 0.02 | 3

In Table. 5, the p-values of the MWU test for a selection of physical
quantities of the pulsars compared between each model and the catalog are
shown. For each model, the p-values shown in Table. 5 are determined from the
random realisation of the model shown in Fig. 2. We note a few of the
prominent features in this table below:

  1. 1.

None of the given models can reproduce the distribution of $\dot{P}$ of radio
pulsars. While the CDFs of radio pulsar $\dot{P}$ look similar between the
catalog and some of the models as shown in Fig. 2, the MWU p-values are always
zero.

  2. 2.

None of the given models can reproduce the distributions of radio flux. This
is surprising, as the distribution of radio luminosities that we assume is
based on radio observations Szary2014. This is likely due to the fact that the
models we explored in this study do not cover a larger or fine enough grid
within the parameter space, which is unlike some of the commonly used Bayesian
analysis methods or machine learning Graber:2023arXiv; Ronchi:2021ApJ that can
explore a larger range of the parameter space with the ability to comb through
different parts finely. It is likely that a slightly different choice of radio
luminosity distribution would yield a better fit to the flux distribution
since this is not something that we varied in our models. Additionally, the
catalog is made up of pulsars from various pulsar surveys, which each
individually would pose different selection effects on the radio pulsar
population. In our work, we only chose to implement selection effects based on
the Parkes multi-beam pulsar survey (cf. Section 2.7.2), and this should also
have an impact on the resulting radio flux distribution.

  3. 3.

Similarly for $\gamma$-ray flux, all of the models report a MWU test p-value
of 0. When examining the CDFs of the $\gamma$-ray flux (e.g. in Fig. 3), we
find that the simulated $\gamma$-ray flux of the pulsars peaks at a lower
value compared to those of the catalog. This effect is to some level expected
as the changes from the sensitivity level mean that our simulation might be
sensitive to $\gamma$-ray pulsars with lower flux. Given that the difference
in the peak value is roughly a factor of 2, we think this is an acceptable
result as it is important to keep the selection of $\gamma$-ray sensitivity as
described above to reflect the observations. Additionally, a recent conference
abstract 2024AAS...24343803T suggested that the 3D fundamental plane might not
work well for $\gamma$-ray pulsars. As described in § 2.8.1, we use 3D
fundamental plane in this work to model the $\gamma$-ray luminosity of
pulsars. This might also cause problem in matching model to the catalogs.

  4. 4.

Only one model, which is the PN1-B1 model, is able to reproduce the period
distribution of radio pulsars that resemble that of the catalog.

  5. 5.

In general, our models reproduce $\gamma$-ray pulsar populations that have a
better resemblance to the catalog, compared to the modelled radio pulsar
populations. This means that more physical quantities for $\gamma$-ray pulsars
have p-values $>0.05$ compared to radio pulsars.

  6. 6.

The Galactic coordinates of the pulsar populations seem to be better fitted
compared to the other physical quantities. This seems to confirm that the
Galactic model, kick distribution and other related models used in this work
is good enough. However, the large distance uncertainties of the pulsars in
the catalog might be a restraint on this modelling.

Since none of the COMPAS models produce results that completely agree with
observations, we will not attempt to choose a best-fit model. Here we examine
and discuss a few models that perform better than the rest by checking the
number of quantities that have p-value $>0.05$. The PN1-B1, P2, and INIT
models have 7 quantities with large p-values; and the P1, PN1 and D3 models
have 6 quantities with large p-values. We examine these models in detail
below, combining the results from Table. 3 in the discussion. In addition, we
also pay attention to the B1 and BLN1-PN4 models, as discussed below. From
Table 4, it is obvious that some models are badly matched to the observations,
including the D2, B2, P1, P3, PN1, BLN2, BLN3, BLN4, BLN1-PN1, and BLN4-PN4.
The bad match for model D2 means that short magnetic decay timescale is not
favoured. The bad match for models B2, BLN2, BLN3, and BLN4, means that low
birth magnetic field is not favored. The bad match for models P1, P3 and PN1
means that large birth rotational periods are not favoured. On the other hand,
simply by examining this table, it appears that the model’s INIT, D1, D3, B1,
P2, BLN1, and PN1-B1 match the catalogs. Combining the results in Tables. 4
and 5, it appears that models INIT, D3, PN1-B1 can best describe the observed
properties of canonical pulsar populations. Examining the CDFs on the physical
parameters of all 10 realisations of the INIT, D3, and PN1-B1 models, it
appears that visually, the realisations of the PN1-B1 model best match to the
distributions of the physical parameters out of the three models.

Figure 3: CDFs of various physical quantities predicted/produced by the PN1-B1
model (red), compared to the data from the catalogs (blue). The red lines in
these plots are from all ten realisations of the PN1B1 model that are sampled
randomly as described in § 2.9. From top left to bottom right, the panels show
the distributions of $\gamma$-ray pulsar $\dot{E}$, radio pulsar $\dot{E}$,
$\gamma$-ray flux, radio flux, $\gamma$-ray pulsar period, radio pulsar
period, $\gamma$-ray pulsar $\dot{P}$, radio pulsar $\dot{P}$, $\gamma$-ray
pulsar galactic latitude, radio pulsar galactic latitude, $\gamma$-ray pulsar
galactic longitude and radio pulsar galactic longitude.

Overall, as discussed above, we determine that the PN1-B1 model best describes
the pulsar initial conditions which produce pulsar populations in agreement
with the catalogs. We first plot the CDFs of the quantities listed in Table. 5
for PN1-B1 model in Fig. 3, showing the results of all 10 realisations.
Overall, the CDFs of a given physical quantity between the model and the
catalog match well with each other. We then visualise the results with the
P-$\dot{P}$ diagrams of the $\gamma$-ray pulsars and radio pulsars generated
from one randomly chosen realisation of the PN1-B1 model, and overlay them
with the canonical pulsars from the catalog, as shown in Fig. 4. These
diagrams, for pulsar researchers, would serve as a better illustration for the
comparison between the model and the catalog. It is obvious that for the
$\gamma$-ray pulsars, the simulation and catalog pulsar populations occupy the
same part of the parameter space. For the radio pulsar populations, however,
the shapes are slightly different between the catalog and simulation. There
are fewer radio pulsars from the simulation that occupy the top left corner of
the P-$\dot{P}$ diagram, indicating that the simulation produces fewer young,
more energetic pulsars with higher $\dot{E}$ . There are more radio pulsars
from the simulation that occupy the parameter space in the bottom left part of
the plot, indicating there are more pulsars with lower magnetic field strength
compared to the catalog.

Figure 4: Top: P-$\dot{P}$ diagram of Fermi detected canonical pulsars (blue)
and detected $\gamma$-ray pulsars from the characteristic PN1-B1 model (red);
bottom: P-$\dot{P}$ diagram of radio detected canonical pulsars as recorded in
the ATNF catalog (blue) and detected radio pulsars from the characteristic
PN1-B1 model (red). The two black dashed lines represent the death-lines
described in Eqs. 6 & 7. The gray dashed lines labeled $10^{3}$ yr, $10^{6}$
yr, and $10^{9}$ yr represent constant characteristic age each respectively;
and those labeled with $10^{31}$ erg s-1 , $10^{33}$ erg s-1 , and $10^{37}$
erg s-1 represent constant $\dot{E}$ each respectively.

In Smith:2019ApJ, a discussion of distance biases on the pulsar catalog was
presented. Lower $\dot{E}$ pulsars, the majority of which are off the Galactic
plane, cannot be detected easily due to their low efficiency, small sample,
and possible large distance. However, Smith:2019ApJ suggest that if low
$\dot{E}$ pulsars emit in $\gamma$-rays, it would be useful to use the
population synthesis approach to help understand if a large population of such
pulsars would contribute to a part of the diffuse Galactic $\gamma$-ray
emission. This work does not provide a full population modelling of pulsars,
especially old pulsars, as only the recent star formation history is
considered. Therefore we are unable to determine if the diffuse background can
be attributed to old pulsars. This will be discussed in a follow-up study.

Here we show the heuristic flux of detected and stacked pulsars in Fig. 5, as
for the detected pulsars, both those from the simulation and the Fermi-LAT
catalog occupy similar parts of the parameter space. The heuristic flux is
defined as $\sqrt{\dot{E}}/4\pi d^{2}$, and is often used in observation to
estimate the potential $\gamma$-ray flux of a pulsar, especially if it is not
detectable. In Fig. 5, excluding those with estimated distance 25 kpc (which
is the flagged value based on the dispersion measure YMW17), the stacked
pulsars from the simulation, in general, have larger distances compared to
those from the catalog. This might also stem from the fact that the distance
measurements of the stacked pulsars from the catalog are largely uncertain.

Figure 5: Expected $\gamma$-ray flux of pulsars (the heuristic flux) as a
function of the spindown luminosity $\dot{E}$ (similar to Fig. 7 in
Smith:2019ApJ). We plot the expected $\gamma$-ray flux of pulsars (the
heuristic flux) when assuming $L_{\gamma}=\sqrt{10^{33}\dot{E}}$ erg s-1. Blue
squares are Fermi-LAT detected pulsars, red triangles are the detected pulsars
from the simulation using the PN1-B1 model, green dots are the pulsars from
the sample list of Song2023, and the orange crosses are the pulsars that
satisfy the criteria of being stacked from the simulation using the PN1-B1
model. The gray dotted lines labelled as 0.3 kpc, 1 kpc, 5 kpc and 15 kpc
represent those of constant distance respectively. Most of the pulsars stacked
in Song2023 with $\dot{E}$ $<10^{33}$ erg s-1 remain undetected in spite of
having similar heuristic fluxes to those detected with higher $\dot{E}$. Most
of the pulsars being stacked from the simulation have larger distances
compared to those from Song2023.

We show the spatial coordinates of detected pulsars in Figs. 6 & 7. These
plots show a similar story as the P-$\dot{P}$ diagram, where the distributions
of $\gamma$-ray pulsars across the Galaxy seem to agree better with the
catalogs, and the radio pulsars less so. We note the main difference in the
spatial distribution of radio pulsars between the catalog and simulation is a
large population of radio pulsars around the Galactic Centre is present in the
simulation. There is however no such population present in the catalog. There
are two possible causes for this issue: 1). The underlying stellar population
that produced the pulsars is extremely concentrated around the Galactic
Centre, making it the dominant effect. It is reasonable for this to happen,
even with the selection effect accounting for high sky temperature around this
direction of observation. 2). It is also likely that the distance measurements
of the detected radio pulsars are largely uncertain, especially towards the
direction of the Galactic Centre, due to high DM/scatter in this direction,
causing the catalog pulsars to have inaccurate coordinates.

Figure 6: Galactocentric coordinates of radio pulsars from the catalog (blue)
and from PN1-B1 model (red). The yellow star in all panels except for the Y-Z
planes indicates the location of the Sun.
Figure 7: Galactocentric coordinates of $\gamma$-ray pulsars from the catalog
(blue) and from PN1-B1 model (red). The yellow star in all panels except for
the Y-Z planes indicates the location of the Sun.

###  3.3 Impact of $\gamma$-ray Luminosity Model and Sensitivity

We understand that there are impacts on this work from various model choices
implemented in this study. Here we discuss the impact on the simulation
results from two model choices: $\gamma$-ray luminosity of pulsars, and the
$\gamma$-ray sensitivity limit.

Luminosity determination can impact the simulation results from the
uncertainty in the fitting of the fundamental plane itself. As described above
in § 2.8, we use the fundamental plane relations with the best-fit parameters
without any scatter. Introducing the uncertainties back into the parameters
would end up with slightly different fundamental plane relations for each
pulsar. We attempt to set the fundamental plane with the uncertainties
included in Eq. 17, and re-analyse the PN1-B1 model.

In Fig. 8, we plot the $\gamma$-ray luminosity of the observed pulsars in
PN1-B1 model and catalog vs. their $\dot{E}$ for both the original analysis
and the $\gamma$-ray luminosity with the scattered fundamental plane. Overall,
when taking the fundamental plane without any uncertainties, the Fermi
detected pulsars and pulsars from the simulation follow the fundamental plane
relationship. The obvious effect of using a scattered fundamental plane is
that the luminosity in the model does not follow the apparent fundamental
plane relations where $L\sim\dot{E}^{0.39}$. As shown in the right panel of
the figure, the correlation is closer to the $L\sim\dot{E}$ line. This is
simply due to the fact that the scatter from the power law indices on various
quantities have uncertainties that are too big, causing this up scatter.
Additionally, some pulsars have luminosity larger than their $\dot{E}$ both in
simulation and catalog. When this occurs in the simulation, it is due to the
observational effect of the flux correction factor which is caused by the
beaming of the pulsar. There are other reasons for this to happen in observed
pulsars, most likely due to the large errors in the estimation of pulsar
distances and beaming geometry of the pulsars.

Figure 8: Left: $\gamma$-ray luminosity of all detected pulsars vs. $\dot{E}$
, assuming for fixed fundamental plane. Right: $\gamma$-ray luminosity of all
detected pulsars vs. their $\dot{E}$ for the scattered fundamental plane. In
both panels, the blue data points with error bars are the pulsars detected by
Fermi-LAT and the red data points are from the simulation using PN1-B1 model.
The magenta dashed line represents the heuristic relation
$L_{\gamma}=\sqrt{10^{33}\dot{E}}$ and the red solid line represents the limit
of energy conservation as $L_{\gamma}=\dot{E}$.

We also discovered that the choice of the sensitivity limit for the radio-
detected pulsars, as described in § 2.8.2, can impact the final results.
Specifically, when we change our assumption about the sensitivity floor, the
detected population in the model will change too. As stated in § 2.8.2, we
deviate from the standard 4FGL sensitivity based on the radio detectability of
the pulsar. We define $\alpha$ as the factor to raise the sensitivity floor
from the default 4FGL-DR4 sensitivity when the pulsar is not detected in the
radio; and $\beta$ as the factor to lower the sensitivity floor when the
pulsar is detected in the radio. Here we set a few different options for the
combination of $\alpha$ and $\beta$ as described in Table. 6, then re-run the
analysis pipelines. We compare the results by checking the number of
$\gamma$-ray pulsars produced in each instance, as well as the ratio between
radio-quiet and radio-loud $\gamma$-ray pulsars ($q$). This analysis is
performed on the PN1-B1 model with 6 different combinations of $\alpha$ and
$\beta$. Table. 6 shows the results of the number of $\gamma$-ray pulsars
detected in the model (NG) as well as the ratio $q$ between the $\gamma$-ray
loud and $\gamma$-ray quiet radio pulsars. When $\alpha$ increases, more
$\gamma$-ray pulsars are detected and the ratio between radio-quiet and radio-
loud $\gamma$-ray pulsars increases. In this case, most of the increase of NG
can be attributed to more radio-quiet pulsars being detected in $\gamma$-rays.
When $\beta$ increases, more $\gamma$-ray pulsars are detected and the ratio
between radio-quiet and radio-loud $\gamma$-ray pulsars decreases. In this
case, most of the increase of NG can be attributed to more radio pulsars being
detected in $\gamma$-rays.

Table 6: Results of different setup for the $\gamma$-ray sensitivity. $\alpha$ is the factor to raise the sensitivity floor from the 4FGL value at the location if the pulsar is not detected in radio, and $\beta$ is the factor to lower the sensitivity from the 4FGL value at the location if the pulsar is detected in radio. Model-0 is the same that adopted for the results above. Models 1-6 are additional models. NG is the number of $\gamma$-ray pulsars produced by the model given the sensitivity setup, and the q is the ratio between the radio-quiet and radio-loud $\gamma$-ray pulsars. Cat denotes the observed values for NG and q.  Model | $\alpha$ | $\beta$ | NG | q
---|---|---|---|---
Cat | – | – | 160 | 90:70
0 | 1.5 | 5.0 | 152 | 87:65
1 | 1.0 | 3.0 | 171 | 122:49
3 | 1.0 | 5.0 | 171 | 114:57
3 | 1.25 | 4.0 | 160 | 99:61
4 | 2.0 | 1.0 | 81 | 62:19
5 | 2.0 | 2.5 | 98 | 57:41
6 | 2.0 | 4.0 | 128 | 68:60

Overall, this further confirms that there are multiple factors that impact the
results of this study. Future development in better physical understanding of
pulsar emission mechanism and stronger correlation between the radio and
$\gamma$-ray detection limit from the current population can further improve
these results.

###  3.4 Low $\dot{E}$ Pulsars in $\gamma$-rays

Theoretical studies of $\gamma$-ray emission from pulsars, such as the
equatorial current sheet model current_sheet0; current_sheet1; current_sheet2
or the outer gap model outer_cheng1986a; outer_cheng1986b; outer_romani1996;
Takata2011, seem to indicate that pulsars with $\dot{E}$ $<10^{33}$ erg/s do
not emit in $\gamma$-ray s. Indeed, if we investigate Eq. 17, it shows that at
$\dot{E}$ $=10^{33}$ erg/s, the $\gamma$-ray luminosity of a pulsar is equal
to $\dot{E}$ , reaching maximum efficiency. With $\dot{E}$ $<10^{33}$, if we
still follow the fundamental plane relation, the $\gamma$-ray luminosity of
the pulsar becomes larger than $\dot{E}$ , breaking the conservation of
energy. It is also worth noting that according to Pulsar_FP2019;
Pulsar_FP2022, the fundamental plane relations are rooted in the equatorial
current sheets model with theoretical foundations.

However, the most up-to-date 4FGL/3PC catalogs include a handful of
$\gamma$-ray pulsars that are detected with $\dot{E}$ $<10^{33}$ erg/s. Given
the low number of detections of $\gamma$-ray pulsars with $\dot{E}$ $<10^{33}$
erg/s, there is no clear understanding of their emission mechanism. In fact,
according to Song2023, pulsars might have a universal, weak, and isotropic
$\gamma$-ray emission. We try to replicate the stacking analysis performed in
Song2023 to show further evidence of the possible emission mechanism for low
$\dot{E}$ pulsars.

We choose pulsars from the PN1-B1 model that are radio loud, $\gamma$-ray
quiet that are at least $20^{\circ}$ away from the Galactic plane for this
study, consistent with Song2023. There are 590 such pulsars in the simulation.
We then assign each pulsar with two different $\gamma$-ray luminosities, one
that follows the fundamental plane if $\dot{E}>10^{33}$ erg s-1and
$L=0.8\dot{E}$ when $\dot{E}<10^{33}$ erg s-1; and the other one is assigned a
luminosity drawn from a log normal distribution $\log_{10}\sim N(32,0.5)$
(ergs/s). The former luminosity assignment corresponds to beamed $\gamma$-ray
emission correlated to the $\dot{E}$ of each pulsar, and the latter
corresponds to the universal, weak and isotropic $\gamma$-ray emission. For
the first emission mechanism, we still account for the $\gamma$-ray beaming
effect discussed in § 2.8.3 when being stacked, while the isotropic emission
does not need to be beamed.

We then calculate the average luminosity of the pulsars in five different
$\dot{E}$ bins, with each bin having the same number of pulsars. We check
these results in a similar way to Fig. 6 in Song2023. Here in this work, we
choose the canonical pulsars from Song2023 and plot them in Fig. 9, and use
this as the observational results to be compared with the model. We then plot
in Fig. 9 the average luminosity vs. $\dot{E}$ from the model with three
different scenarios: only the fundamental plane is considered, only isotropic
emission is considered, and combining both. Unfortunately, because the
observational results are largely uncertain at this point, we are unable to
determine which would be a better fit to describe the emission mechanism of
lower $\dot{E}$ pulsars. We either need a better understanding of pulsar
emission in this regime, or more data to improve the observational results.

Figure 9: Left: $\gamma$-ray luminosity vs. $\dot{E}$ for stacked pulsars. The
green data points with error bars in both directions are from Song2023 and the
orange data points with only horizontal error bars are from the simulation
using PN1-B1 model. The magenta dashed line represents the heuristic relation
$L_{\gamma}=\sqrt{10^{33}\dot{E}}$ and the red solid line represents the limit
of energy conservation as $L_{\gamma}=\dot{E}$. The left panel calculates the
stacked luminosity from the simulation by assuming the stacked pulsars have a
luminosity that follows the fundamental plane while $\dot{E}>10^{33}$ erg s-1
and $L_{\gamma}=0.8\dot{E}$ when $\dot{E}<10^{33}$ erg s-1Ṫhe middle panel is
calculating the luminosity assuming all these pulsars have a weak, isotropic
$\gamma$-ray emission follows a log-normal distribution where
$\log{L_{\gamma}}~{}N(32,0.5)$. The right panel combines both.

Another way to examine the results and potential emission mechanism is to
compare the spectral energy distributions (SEDs) for various scenarios. We
calculate the stacked SED of the 590 pulsars for the fundamental plane
emission and isotropic emission as described above. The $\gamma$-ray spectral
parameters of the pulsars, $E_{\textrm{cut}}$ and PL index before cutoff, are
drawn from the described text above. When considering the fundamental plane
emission, a pulsar needs to be beamed towards the observer. When considering
the isotropic emission, it does not need to be beamed. Averaging the SEDs, we
can produce the stacked SEDs for both scenarios as shown in Fig. 10. We
compare this result to Song2023 as well as the averaged pulsar SED from
2015ApJ...804...86M. As expected, if considering the FP emission, this stacked
pulsar population is relatively older compared to the young, energetic
pulsars, hence it is somewhere in between the averaged young pulsar and MSP
SED. That does mean that this is much more luminous than the stacked pulsars
from Song2023. When considering the isotropic emission, however, the SED is
much less luminous compared to the stacked pulsars from the data. Neither of
our models here can completely explain the observation. It might be a mixture
of both scenarios that can account for the observation, averaging out to the
observed data points.

Figure 10: SED of stacked pulsars from the simulation using the PN1-B1 model.
The orange solid curve is the averaged SED from the stacked pulsars in the
simulation assuming they all follow the empirical fundamental plane (FP)
emission mechanism in $\gamma$-rays, and the red dashed curve is for the same
population of pulsars assuming they all emit weakly isotropic $\gamma$-ray
emission (ISO). The blue data points with error bars are the averaged SED of
the sampled canonical pulsars from Song2023, The green dotted-dashed curve
represents the averaged SED of young pulsars, and the green dotted curve SED
of MSPs 2015ApJ...804...86M.

Many low $\dot{E}$ pulsars are located outside of the Galactic plane Fermi-
LAT:2023zzt. It is commonly believed that the high latitude pulsars are a
result of pulsars receiving larger kicks when they went through supernova
explosions. The results from this work can be used to confirm if that is
indeed the case, as the kick velocities are one of the modelled components in
our population synthesis models. In Fig. 11, we show the distributions of kick
velocities received by pulsars from the model that are either detected or
stacked, as selected above. Consistent with our hypothesis, we find that the
high latitude, low $\dot{E}$ pulsars in the stacked sample received larger
natal kicks than the typical detected population.

Figure 11: Kick distributions of simulated pulsars. The blue histogram shows
the kick distribution of the pulsars detected in the simulation, and the red
histogram shows the kick distribution of the stacked pulsars from the
simulation. The selection of pulsars off of the Galactic plane preferentially
selects pulsars born with large kick velocities.

##  4 Future Development

As the first of a series of studies using COMPAS to study $\gamma$-ray pulsar
populations, there are some caveats as discussed above, and these are
summarised below:

  1. 1.

This work currently does not include MSP populations. Pulsar recycling physics
will be implemented in a direct follow-up of this work so that a full
population study of pulsars including MSPs can be carried out. MSPs are
products of binary evolution, and so comparing to binary evolution models will
be important to fully understand the population.

  2. 2.

Some NSs are formed through accretion-induced collapses, and these NSs are not
accounted for in this work. This is mainly because it takes much longer than
200 Myrs for a NS to form through this channel, and they are naturally left
out from this study as we only choose to include binary systems formed within
the most recent 200 Myrs. Additionally, neutron stars formed through
accretion-induced collapse may directly form MSPs (Hurley:2010MNRAS;
Tauris:2013A&A; Ruiter:2019MNRAS, e.g.,), which are not the focus of this
work.

  3. 3.

We currently lack a physical model for the birth spin periods and magnetic
fields for pulsars as a function of the progenitor properties (Mandel:2020qwb,
cf.). This work, and any follow-up in the foreseeable future will only be able
to use some empirically assumed distribution of these quantities.

  4. 4.

We have neglected the population of magnetars with long spin periods and high
magnetic fields (Popov:2010MNRAS; Gullon:2015MNRAS; Makarenko:2021MNRAS,
e.g.,). As with canonical pulsars, most magnetars are isolated, whilst their
massive star progenitors were presumably born in binaries. Binary interactions
such as stellar mergers may be responsible for producing highly magnetic
neutron stars Schneider:2019Nature.

After addressing the above-mentioned caveats, we would be able to simulate
full populations of NSs and broaden the scope of our research with several
different follow-up studies:

  1. 1.

A full scope of pulsar population synthesis can enable the study of NS
contribution to the diffuse $\gamma$-ray background 2017arXiv170500009F

  2. 2.

The inclusion of MSPs in the population synthesis models would provide further
insight into the possible MSP origin of the Galactic Centre GeV Excess
hooperm2016; ploeg2017; eckner2018.

  3. 3.

A comparison of pulsar populations between the Galactic field and globular
clusters with results from dynamical simulations, such as those from NBODY and
CMC (Ye:2019ApJ; Ye:2024ApJ, e.g.,).

  4. 4.

After adding in the formation of MSPs and NSs formed through accretion-induced
collapse, the model will be able to perform on the Hubble time scale. This
will enable the modelling of merger events for gravitational wave detections,
such as black hole-NS binaries Chattopadhyay2021, double NSs Chattopadhyay2020
and NS-white dwarf binaries.

  5. 5.

Include magnetar physics in future versions of COMPAS and synthesize magnetar
populations.

  6. 6.

Investigate more sophisticated model comparison methods to fully explore the
parameter space and identify a best-fit model that produces the population of
pulsars observed to date. Several approaches have been taken in the
literature, including hierarchical Bayesian parameter estimation
(Ploeg:2020JCAP, e.g.,) and machine learning approaches (Ronchi:2021ApJ;
Graber:2023arXiv, e.g.,). At present, such a large parameter space exploration
is computationally infeasible using COMPAS. In a follow-up study, we will
optimise, parameterise and automate our COMPAS simulations and post-processing
pipelines so that the use of more sophisticated model exploration and
comparison techniques becomes computationally feasible.

In addition to varying the birth distributions of the pulsar properties (birth
magnetic fields and spin periods), in future work, we aim to study how
different binary properties, including initial orbital properties, different
mass transfer theories, neutron star kick velocities and common envelope
physics can impact the pulsar population, and how this population can be used
to constrain the underlying physics of massive binary evolution.

##  5 Conclusions

We performed a population synthesis of the canonical Galactic pulsar
population using the rapid binary population synthesis code COMPAS COMPAS on
both radio and $\gamma$-ray canonical pulsars. We generated 20 different
models for pulsar birth properties for COMPAS runs incorporating a magnetic
braking model with a braking index of 3 for pulsar spindown. We discover that
the PN1-B1 model, where the birth period of pulsar follows a normal
distribution with $\mu=75$ ms and $\sigma=25$ ms, the birth magnetic field of
pulsars follows a uniform distribution between $10^{11}$ and $10^{13}$ G, and
magnetic decay time scale $\tau_{d}=500$ Myr best describes the observed
pulsar population, based on the facts that:

  1. 1.

PN1-B1 model produces the number of radio and $\gamma$-ray canonical pulsars
that, given the star formation rate of the model, best fits the catalog.

  2. 2.

The distribution of physical properties of radio and $\gamma$-ray pulsars
produced from this model matches the catalog the best.

We then compare in detail the PN1-B1 model results to the catalog by examining
the P-$\dot{P}$ plots, the Galactocentric coordinates, and the heuristic flux
between the simulation and the catalog. We also discussed the impact of the
selection effect on the population synthesis results, demonstrating that the
slightest change in pulsar physics and observation sensitivity can bring
significant differences to the final simulation.

Using the results from PN1-B1 model, we also attempted to explore the possible
emission mechanism from low $\dot{E}$ pulsars ($\dot{E}$ $<10^{33}$ erg s-1,
and failed to make a distinction between a beamed emission mechanism or an
isotropic emission mechanism.

We will follow up this study to include the Galactic MSP population and
utilise the results to answer other physical questions, including the
$\gamma$-ray diffuse background, Galactic Centre GeV excess, and gravitational
wave events from a NS binary.

The authors appreciate coding assistance from R. Wilcox and J. Riley, and
helpful discussions from T. Paglione, J. Tan and J. Hurley.

##### Funding Statement

This research was conducted by the Australian Research Council Centre of
Excellence for Gravitational Wave Discovery (OzGrav), through project numbers
CE170100004 and CE230100016. S. Stevenson is a recipient of an ARC Discovery
Early Career Research Award (DE220100241). D. Chattopadhyay is supported by
the UK’s Science and Technology Facilities Council grant ST/V005618/1. This
work used the OzSTAR high-performance computer at Swinburne University of
Technology. OzSTAR is funded by Swinburne University of Technology and the
National Collaborative Research Infrastructure Strategy (NCRIS).

##### Competing Interests

None.

##### Data Availability Statement

COMPAS is open source and hosted on Github. Due to the size of COMPAS output
files, the authors can share configuration files used to produce the pop-synth
results upon request.

Post-processing scripts will be hosted on Github as well. Before being made
available on Github, they will be available upon request provided that this
work is properly cited.