1. Introduction
Pulsars are rotating neutron stars most commonly observed by their periodic radio emission. Around 3 500 pulsars are known to date (Manchester et al. Reference Manchester, Hobbs, Teoh and Hobbs2005). Modern radio facilities such as MeerKAT, ASKAP and FAST are uncovering new pulsars (e.g. Han et al. Reference Han2021; Padmanabh et al. Reference Padmanabh2023), whilst modern search techniques are still discovering pulsars in decades-old archival surveys (e.g. Sengar et al. Reference Sengar2023, Reference Sengar2025). 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 (e.g. Kramer et al. Reference Kramer2021).
Whilst most pulsars are discovered and studied through radio surveys, some are also observed at other wavelengths, including optical (Cocke, Disney, & Taylor Reference Cocke, Disney and Taylor1969; Nather, Warner, & Macfarlane Reference Nather, Warner and Macfarlane1969), X-rays (Bradt, Rappaport, & Mayer Reference Bradt, Rappaport and Mayer1969; Fritz et al. Reference Fritz, Henry, Meekins, Chubb and Friedman1969; Oegelman, Finley, & Zimmerman Reference Oegelman, Finley and Zimmerman1993) and
$\gamma$
-rays (Kniffen et al. Reference Kniffen, Hartman, Thompson, Bignami and Fichtel1974; Tumer et al. Reference Tumer1984; Thompson et al. Reference Thompson1992). A growing number of pulsars are being detected in
$\gamma$
-rays by the Large Area Telescope (LAT) onboard the Fermi space telescope (Smith et al. Reference Smith2023), 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 & Holt Reference Halpern and Holt1992; Bertsch et al. Reference Bertsch1992; Malofeev & Malov Reference Malofeev and Malov1997).
The spindown luminosity,
$\dot{E}$
, is a useful parameter describing a pulsar, and is defined as
\begin{equation} \dot{E} = 4 \pi^{2} I \dot{P} P^{-3}, \end{equation}
where P is the spin period of the pulsar,
$\dot{P}$
is its spin derivative, and I is its moment of inertia. Most of the pulsars observed in
$\gamma$
-rays (e.g. Abdo et al. Reference Abdo2009a; Abdo et al. Reference Abdo2013; Pletsch et al. Reference Pletsch2013; Smith et al. Reference Smith2023) have high spin down luminosities (
$\dot{E} \gtrsim 10^{33}$
erg s
$^{-1}$
). This includes both young pulsars (mostly found at low Galactic latitudes) and millisecond pulsars (MSPs). 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 (with this fraction
$\eta \gtrsim 10\%$
; Abdo et al. Reference Abdo2013). On the other hand, comparing the radio luminosity of pulsars to their spindown luminosities suggests that only a small fraction of the energy is radiated away in radio wavelengths. The radio emission efficiency
$\xi$
, defined as the fraction of rotational energy transformed into radio emission
\begin{equation} \xi = \frac{L_{\textrm{r}}}{\dot{E}}, \end{equation}
where
$L_{\textrm{r}} (\textrm{erg s}^{-1}) \approx 7.4 \times 10^{27} L_{\textrm{r}}$
(mJy), will always be smaller than 1% (Szary et al. Reference Szary, Zhang, Melikidze, Gil and Xu2014). 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 (Abdo et al. Reference Abdo2008, Reference Abdo2009a, Reference Abdo2013; Smith et al. Reference Smith2023). The distinction between these two populations likely depends on whether their radio emission is beamed towards the Earth (Abdo et al. Reference Abdo2009a).
The current sample of
$\gamma$
-ray detected pulsars is around 300 (Abdo et al. Reference Abdo2013; Smith et al. Reference Smith2019, Reference Smith2023)Footnote
a
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 (Abdo et al. Reference Abdo2013). An extensive search of
$\gamma$
-ray pulsations from more than one thousand pulsars (Smith et al. Reference Smith2019) 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, Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) performed a stacking analysis of 12 yr 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| \gt 20^\circ$
). This probed a population of relatively old, low
$\dot{E}$
pulsars, different to the usual population of
$\gamma$
-ray loud pulsars. Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) report the significant detection of
$\gamma$
-rays from this set of pulsars compared to the expected background. Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) 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. For example, Gonthier et al. (Reference Gonthier, Ouellette, Berrier, O’Brien and Harding2002) used Monte Carlo simulations to predict the radio and
$\gamma$
-ray pulsar populations that would be observed by Fermi in the context of the polar cap model for
$\gamma$
-ray emission. In a seminal work, Faucher-Giguère & Kaspi (Reference Faucher-Giguère and Kaspi2006) 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 (e.g. Takata, Wang, & Cheng Reference Takata, Wang and Cheng2011; Dirson, Pétri, & Mitra Reference Dirson, Pétri and Mitra2022) that include radio and
$\gamma$
-ray selection effects and pulsar evolution theories. Modern pulsar population synthesis studies (e.g. Gullón et al. Reference Gullón, Miralles, Viganò and Pons2014; Shi & Ng Reference Shi and Ng2024; Graber et al. Reference Graber, Ronchi, Pardo-Araujo and Rea2024; Pardo-Araujo et al. Reference Pardo-Araujo, Ronch, Graber and Rea2025) make use of updated assumptions regarding the birth distributions, magnetic field decay, and pulsar emission geometry, and some of them use sophisticated machine learning inference techniques to compare model predictions to observations.
Typically, existing pulsar population synthesis works, such as those done by using PsrPopPy (Bates et al. Reference Bates, Lorimer, Rane and Swiggum2014), do not model the details of the full evolution of the binary stellar systems. Neutron stars – and hence pulsars – originate from massive stars, and the majority of massive stars are born in binaries (e.g. Sana et al. Reference Sana2012; Moe & Di Stefano Reference Moe and Di Stefano2017). 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 (e.g. Brandt & Podsiadlowski Reference Brandt and Podsiadlowski1995; Kalogera Reference Kalogera1996; Renzo et al. Reference Renzo2019). 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 (Riley et al. Reference Riley2022), 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.
Some previous studies have combined pulsar population synthesis with binary stellar evolution. The study of Kiel et al. (Reference Kiel, Hurley, Bailes and Murray2008) is closest to ours in terms of methodology. They use binary population synthesis to study the Galactic pulsar population, but they focus on radio pulsars and do not explore the
$\gamma$
-ray pulsar population. More recently, Titus et al. (Reference Titus2020) used binary population synthesis to study the radio pulsar population in the Small Magellanic Cloud. In this work, we attempt to explore the canonical radio and
$\gamma$
-ray pulsar populations produced through massive binary stellar evolution and compare them to the pulsar populations observed in the Milky Way. In a follow-up work, we will attempt to understand the full effects on pulsar populations from stellar evolution physics.
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 (Wijnands & van der Klis Reference Wijnands and van der Klis1998), spider pulsars (Fruchter, Stinebring, & Taylor Reference Fruchter, Stinebring and Taylor1988) and double neutron star binaries (e.g. Hulse & Taylor Reference Hulse and Taylor1975; Abbott et al. Reference Abbott2017), 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 et al. Reference Stevenson2017; Vigna-Gómez et al. Reference Vigna-Gómez2018; Riley et al. Reference Riley2022; Mandel et al. Reference Mandel2025). COMPAS is based on the Single Star Evolution (SSE, Hurley, Pols, & Tout Reference Hurley, Pols and Tout2000) code and the Binary Star Evolution (BSE, Hurley, Tout, & Pols Reference Hurley, Tout and Pols2002) code with updates and modifications made as described in Riley et al. (Reference Riley2022). Binary interactions such as mass transfer and common envelope evolution are included, but we neglect the effects of tides (Riley et al. Reference Riley2022).
We use default binary property parameters for our COMPAS simulations (Riley et al. Reference Riley2022). 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 (IMF; Kroupa Reference Kroupa2001) within the range of 5–150 M
$_{\odot}$
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 \lt 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
$_\odot$
(Stevenson et al. Reference Stevenson2019). We assume initially circular orbits (Hurley et al. Reference Hurley, Tout and Pols2002). We draw the orbital separation from a uniform distribution in log space with a minimum of 0.01 and maximum of 1 000 AU (Sana et al. Reference Sana2012). 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 et al. Reference Asplund, Grevesse, Sauval and Scott2009). A list of default COMPAS option values, including mass transfer rate and common envelope evolution, can be found in the link in footnote.Footnote
b
We assume a binary fraction of 100% for massive stars and neglect single stars. This is appropriate for massive O stars (
$M \gt 16$
M
$_\odot$
). For example, Sana et al. (Reference Sana2012) find that 70% of massive stars are in close binaries with
$\log(P_\mathrm{orb}/\mathrm{d}) \lt 3.5$
, and the binary fraction reaches 100% when including wider binaries (Sana et al. Reference Sana2014) which we include in our COMPAS models (see also discussion in de Mink & Belczynski Reference de Mink and Belczynski2015). Many of these stars are effectively single, in the sense that they do not interact with their wide binary companion. Assuming a binary fraction of 100% is likely an overestimate for B-type stars, which are the progenitors of most neutron stars and have more typical binary fractions of 70% (Moe & Di Stefano Reference Moe and Di Stefano2017), see also discussion in Zapartas et al. (Reference Zapartas2017). However, this effect is likely subdominant to many of the other approximations in this work. We leave a thorough exploration of the binary fraction as a function of mass to future work.
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 electron-capture supernovae that are not in wide, non-interacting binaries (Willcox et al. Reference Willcox2021), or through core collapse supernovae are included in the following analysis. The reasoning for this decision is discussed in Section 2.2.
2.2. Neutron star kicks
Neutron stars receive natal kicks at birth (Lyne & Lorimer Reference Lyne and Lorimer1994; Hobbs et al. Reference Hobbs, Lorimer, Lyne and Kramer2005; Verbunt, Igoshev, & Cator Reference Verbunt, Igoshev and Cator2017; Igoshev Reference Igoshev2020; Igoshev et al. Reference Igoshev, Chruslinska, Dorozsmai and Toonen2021; Disberg & Mandel Reference Disberg and Mandel2025). It is common in the pulsar population synthesis literature to use a distribution inferred from the velocities of isolated pulsars as the kick distribution (e.g. Faucher-Giguère & Kaspi Reference Faucher-Giguère and Kaspi2006; Graber et al. Reference Graber, Ronchi, Pardo-Araujo and Rea2024). A common choice is the distribution from Hobbs et al. (Reference Hobbs, Lorimer, Lyne and Kramer2005).Footnote c However, as we are simulating binary progenitors of pulsars, we need to simulate the underlying kick distribution, not the observed kick distribution. Several classes of supernovae are expected to give rise to reduced kicks, including electron-capture supernovae (Gessner & Janka Reference Gessner and Janka2018), ultra-stripped supernovae (e.g. Tauris, Langer, & Podsiadlowski Reference Tauris, Langer and Podsiadlowski2015; Suwa et al. Reference Suwa, Yoshida, Shibata, Umeda and Takahashi2015; Müller et al. Reference Müller2019) and low-mass iron core-collapse supernovae (Stockinger et al. Reference Stockinger2020; Müller et al. Reference Müller2019).
We use the physical remnant prescription from Mandel & Müller (Reference Mandel and Müller2020) 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 et al. Reference Kapil, Mandel, Berti and Müller2023); that is, the velocity distribution of young isolated pulsars matches those inferred from very long baseline interferometry observations (Deller et al. Reference Deller2019). 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 et al. (Reference Kapil, Mandel, Berti and Müller2023) 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 & Igoshev (Reference Mandel and Igoshev2023), 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 & Müller (Reference Mandel and Müller2020) 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 et al. Reference Willcox2021). This choice is also motivated theoretically as binary interactions are expected to widen the parameter space for the production of electron-capture supernovae compared to single stars (e.g. Podsiadlowski et al. Reference Podsiadlowski2004).
2.3. COMPAS implementation of pulsar evolution
This study focuses on canonical pulsars, hence we only describe the evolution of spin and magnetic field of pulsars that are not experiencing any mass transfer. 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 Chattopadhyay et al. (Reference Chattopadhyay, Stevenson, Hurley, Rossi and Flynn2020), and is briefly introduced below.
The surface magnetic field strength of isolated pulsars may decay over time. The exact timescale on which this decay occurs (if at all) is uncertain (see further discussion in Section 2.4). The field decay is thought to be due to the Hall effect leading to a redistribution of the magnetic energy, leading to enhanced Ohmic dissipation (Romani Reference Romani1990; Konar & Bhattacharya Reference Konar and Bhattacharya1997; Konar & Bhattacharya Reference Konar and Bhattacharya1999a; Konar & Bhattacharya Reference Konar and Bhattacharya1999b; Pons et al. Reference Pons, Link, Miralles and Geppert2007; Aguilera, Pons, & Miralles Reference Aguilera, Pons and Miralles2008). In COMPAS, for simplicity we assume an exponentially decaying magnetic field of a neutron star similar to many recent studies (Kiel et al. Reference Kiel, Hurley, Bailes and Murray2008; Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Rossi and Flynn2020; Dirson et al. Reference Dirson, Pétri and Mitra2022; Shi & Ng Reference Shi and Ng2024). The decay of the magnetic field after it spins down for a time duration of t
\begin{equation} B_f = (B_i - B_{\textrm{min}}) \exp({-}t / \tau) + B_{\textrm{min}}, \end{equation}
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 (Zhang & Kojima Reference Zhang and Kojima2006) for the entire analysis, and
$\tau$
is the time scale of pulsar magnetic field decay.
The rate of change of the spin period of a neutron star over time is known as the spin period derivative,
$\dot{P}$
. Assuming spindown due to a rotating magnetic dipole in a vacuum,
$\dot{P}$
is given by
\begin{equation} \dot{P} = \frac{8 \pi ^ 2 R ^ 6 B^2 \sin^2{\alpha}}{3 c^3 I P} , \end{equation}
where R is the radius of the pulsar (we assume
$R = 10$
km for all pulsars at all times in this work), c is speed of light,
$\alpha$
is the alignment angle between the magnetic and rotational axes, and
\begin{equation} I = (0.237 \pm 0.008) M R^2 \left[1 + 4.2\frac{\textrm{M}}{\textrm{M}_{\odot}} \frac{\textrm{km}}{\textrm{R}} + 90\left(\frac{\textrm{M}}{\textrm{M}_{\odot}}\frac{\textrm{km}}{\textrm{R}}\right)^4\right] , \end{equation}
being the moment of inertia of the pulsar (Lattimer & Schutz Reference Lattimer and Schutz2005). Substituting B from Eq. (3) into Eq. (4) and integrating over time t gives
\begin{align} P^2 & = \frac{16 \pi ^ 2 R ^ 6 \sin^2{\alpha}}{3 c^3 I} \Big[B_{\textrm{min}}^2 t - \tau B_{\textrm{min}} (B_f - B_i)\nonumber\\[5pt]& \quad \left.- \frac{\tau}{2}\left(B^2_f - B^2_i\right)\right] + P_0^2 . \end{align}
$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. (3). The evolution of P is dependent on
$\alpha$
(Eq. 4). As pointed out in other pulsar population synthesis studies (e.g. Dirson et al. Reference Dirson, Pétri and Mitra2022), the evolution of
$\alpha$
is also dependent on B, as the field decay impacts the evolution of
$\alpha$
. A simplification in COMPAS is made so that
$\alpha$
is set constant and not evolved. This decouples the evolution between P/B and
$\alpha$
. Some works have suggested that
$\alpha$
tends to decrease as the pulsar ages, and the old pulsar populations have near-alignment between the magnetic field and rotational axes (Young et al. Reference Young, Chan, Burman and Blair2010). However, results from Johnston et al. (Reference Johnston, Smith, Karastergiou and Kramer2020) suggest that large
$\alpha$
is needed to model both
$\gamma$
-ray and radio pulsars well, especially for younger populations. Hence in this work we set
$\alpha = \pi/2$
. Evolution of alignment will be included in a future version of COMPAS.
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. Previous pulsar population studies (e.g. Faucher-Giguère & Kaspi Reference Faucher-Giguère and Kaspi2006; Popov & Turolla Reference Popov and Turolla2012; Szary et al. Reference Szary, Zhang, Melikidze, Gil and Xu2014; Johnston et al. Reference Johnston, Smith, Karastergiou and Kramer2020; Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Rossi and Flynn2020; Dirson et al. Reference Dirson, Pétri and Mitra2022; Igoshev et al. Reference Igoshev2022; Du et al. Reference Du2024) have proposed and refined distributions of pulsar properties at birth, such as the 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 field 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 et al. Reference Manchester, Hobbs, Teoh and Hobbs2005). We exclude extremely high magnetic field strengths with B
$\gt 10^{15}$
G in the magnetar regime as magnetars are not powered by rotation like typical pulsars, and may have also have a different formation history.
In our default model, we assume that the magnetic field decay time scale
$\tau_d = 500$
Myr (Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Bailes and Broekgaarden2021). 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 (Gunn & Ostriker Reference Gunn and Ostriker1970; Vivekanand & Narayan Reference Vivekanand and Narayan1981; Lyne, Manchester, & Taylor Reference Lyne, Manchester and Taylor1985) found that the magnetic field decays on the timescale of the order of a million years, with Dirson et al. (Reference Dirson, Pétri and Mitra2022) recently finding
$\tau_d = 4.6\times10^{5}$
yr. In comparison, other works (Taylor & Manchester Reference Taylor and Manchester1977; Stollman Reference Stollman1987; Kiel et al. Reference Kiel, Hurley, Bailes and Murray2008) have shown the opposite, finding that
$\tau_d \gt 100$
Myr is needed in order for theoretical models to match the observed catalogue. Recent work incorporating more sophisticated models of magnetic field decay, pulsar beaming and magnetic field alignment find
$\tau_d$
in the range 1–10 Myr (e.g. Graber et al. Reference Graber, Ronchi, Pardo-Araujo and Rea2024; Shi & Ng Reference Shi and Ng2024). Given the uncertainty, we let
$\tau_d$
vary from 1 Myr to 10 000 Myr in some of 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. Models D1–5, 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–5 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. Model PLN1 gives the birth spin period distribution in a log uniform distribution. 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. Models P2-D1–5, PN1-B1-D1–5 and PN5-B1-D1–5 are set to explore the impact of
$\tau_d$
on different initial distributions of spin period and magnetic field, whilst also varying
$\tau_d$
in the same way as in models D1–5. 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. The original form of the log-normal distribution of the birth magnetic field strength is taken from Szary et al. (Reference Szary, Zhang, Melikidze, Gil and Xu2014). The BLN1-PN5 model is intended to reproduce the model described in Szary et al. (Reference Szary, Zhang, Melikidze, Gil and Xu2014). The birth magnetic field and spin period distributions in BLN5-PLN2-D4 are found in Graber et al. (Reference Graber, Ronchi, Pardo-Araujo and Rea2024). The birth distributions for BLN6-PLN3-D4 are found in Pardo-Araujo et al. (Reference Pardo-Araujo, Ronch, Graber and Rea2025). In Graber et al. (Reference Graber, Ronchi, Pardo-Araujo and Rea2024), Pardo-Araujo et al. (Reference Pardo-Araujo, Ronch, Graber and Rea2025) the authors model the magnetic field decay based on numerical simulations of magnetothermal evolution in the neutron star crust (Viganò et al. Reference Viganò, Garcia-Garcia, Pons, Dehman and Graber2021) up to 1 Myr, and with a power law decay model beyond 1 Myr. Their results are equivalent to a short decay timescale. Hence for these two models, the magnetic field decay time scale is set to 1 Myr. In each COMPAS model, 300 000 binary systems are modelled.
Table 1. Description of models used for COMPAS simulations in this work.
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 (Vigna-Gómez et al. Reference Vigna-Gómez2018), 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}$
yr (Smith et al. Reference Smith2023). The median characteristic age of radio pulsars selected from the catalogue in this study (as described later in Section 3.2) is
$3.1\times10^{6}$
yr. While this age is higher than typical
$\gamma$
-ray pulsars, it still suggests that a typical radio pulsar should be formed in a recent star formation episode. 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. Readers should note that in general, the characteristic age is a poor primer for the true age of a pulsar when the neutron star is young or when field decay is included and the star is older than 1 Myr.
The birth time of any given binary system is then drawn from a uniform distribution with a minimum of 12.0 Gyr and a maximum of 13.0 Gyr. We choose to neglect the population of pulsars produced through accretion-induced collapse (AIC) of white dwarfs as they represent only a small fraction of the population (
$\lesssim 10\%$
, as indicated by our COMPAS models), are formed on much longer timescales and may potentially be born with different properties to those born in core-collapse supernovae (Bailyn & Grindlay Reference Bailyn and Grindlay1990; Tauris et al. Reference Tauris, Sanyal, Yoon and Langer2013).
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 Eq. (1).
The emission of a pulsar is turned off when it can no longer produce electron-positron pairs in its magnetosphere. More realistic modelling of pulsar polar cap have suggested a broad ‘valley’ like region in which pulsars cease their emission (Chen & Ruderman Reference Chen and Ruderman1993; Beskin & Istomin Reference Beskin and Istomin2022), which is starting to be implemented in some recent pulsar population works (e.g. Sautron et al. Reference Sautron, Pétri, Mitra and Dirson2024). In this work, we choose a more simplistic description with a death-line model. Before recording the parameters of the neutron star, it is essential to check if the pulsar crosses the radio death-line (Rudak & Ritter Reference Rudak and Ritter1994) on the P-
$\dot{P}$
diagram, which is defined by
\begin{equation} \log_{10} \dot{P} = 3.29 \times \log_{10} P - 16.55, \end{equation}
and
\begin{equation} \log_{10} \dot{P} = 0.92 \times \log_{10} P - 18.65, \end{equation}
where Eq. (7) is commonly referred to as the death-line for canonical pulsars, and Eq. (8) for MSPs. In the simulation, we decide if a pulsar crosses the death-line (and hence becomes unobservable) by comparing its
$\dot{P}$
to the
$\dot{P}$
values of Eqs. (7) and (8) given the current spin period of the neutron star. If
$\dot{P}$
of a neutron star from COMPAS is larger than both
$\dot{P}$
values calculated by Eqs. (7) and (8), then the pulsar has not yet crossed the death-line and may be observable, and its parameters are recorded. Otherwise, it is discarded from the simulation. There are other pulsar population models that do not need to explicitly include a death-line, e.g. when assuming both radio and
$\gamma$
-ray luminosity are dependent on
$\dot{E}$
(Shi & Ng Reference Shi and Ng2024; Graber et al. Reference Graber, Ronchi, Pardo-Araujo and Rea2024; Pardo-Araujo et al. Reference Pardo-Araujo, Ronch, Graber and Rea2025). Our treatment of radio selection effects is given in Section 2.7, and
$\gamma$
-ray selection effects in Section 2.8.
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 Section 2.1, we assume all the neutron stars are formed in supernovae. Upon formation, the stellar remnant receives a natal kick as described in Section 2.2 which may or may not disrupt the binary system, in some case leading to 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 (Paczynski Reference Paczynski1990):
\begin{equation} p_z(|z|) dz = \exp{({-}|z|/z_{\textrm{exp}})}\frac{dz}{z_{\textrm{exp}}},\end{equation}
where
$|z|$
is the absolute value of the vertical distance from the Galactic plane, and
\begin{equation} p_R(R) dR = a_R \exp{({-}R/R_{\textrm{exp}})} \frac{R}{R^2_{\textrm{exp}}} dR,\end{equation}
where
$a_R = 1.0683$
and R is the distance from the Galactic centre on the Galactic plane. We set
$z_{\textrm{exp}} = 75$
pc and
$R_{\textrm{exp}} = 4.5$
kpc following Paczynski (Reference Paczynski1990) and Dirson et al. (Reference Dirson, Pétri and Mitra2022). Quintana et al. (Reference Quintana, Wright and Martnez Garca2025) recently performed a census of massive OB stars (neutron star progenitors) within 1 kpc of the Sun and found a scale height close to 75 pc. It is worth noting that a recent study of OB stars in Gaia survey (Li et al. Reference Li2019) has suggested a varying scale height with respect to the Galactic radius. Previous pulsar population synthesis studies have used both larger and smaller scale heights (Gullón et al. Reference Gullón, Miralles, Viganò and Pons2014; Shi & Ng Reference Shi and Ng2024).
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 (Rossi Reference Rossi2015), 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 (Flynn, Sommer-Larsen, & Christensen Reference Flynn, Sommer-Larsen and Christensen1996; Plummer Reference Plummer1911):
\begin{equation} \Phi_b = - \frac{GM_b}{\sqrt{R^2 + z^2 + b_b^2}}, \end{equation}
where G is the gravitation constant,
$M_b = 1.0\times10^{10}$
M
$_{\odot}$
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 (Miyamoto & Nagai Reference Miyamoto and Nagai1975):
\begin{equation} \Phi_d = -\sum_{n=1}^{3} \frac{G M_{d_n}}{\sqrt{R^2 + \left[a_{d_n}+\sqrt{b_d^2 + z^2}\right]^2}}, \end{equation}
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.
An NFW dark matter halo (Navarro, Frenk, & White Reference Navarro, Frenk and White1997) is used in the modelling of the Galactic potential:
\begin{equation} \Phi_h = -\frac{GM_h}{r}\ln{\left(1+\frac{r}{a_h}\right)}, \end{equation}
where
$M_h = 3.3\times10^{12}$
M
$_{\odot}$
is the mass of the halo component,
$a_h = 45.02$
kpc is the length scale and
$r^2 = R^2 + z^2$
.
Table 2. Parameters used for the three components of exponential disk.
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. Theoretical work on modelling pulsar magnetosphere has suggested a dependence of radio luminosity on spin-down luminosity
$\dot{E}$
(Lorimer & Kramer Reference Lorimer and Kramer2004). Recent measurements made by the MeerKAT Thousand Pulsar Array have revealed evidence for such correlation, although it is a weak correlation (Posselt et al. Reference Posselt2023). Some population studies have also found no significant dependence of radio luminosity on P,
$\dot{P}$
or
$\dot{E}$
of the pulsar itself (Szary et al. Reference Szary, Zhang, Melikidze, Gil and Xu2014). For simplicity, we opt to assume no correlation between radio luminosity and spin-down luminosity following Szary et al. (Reference Szary, Zhang, Melikidze, Gil and Xu2014). We assume that the distribution of pulsar radio luminosities of a pulsar (in units of mJy kpc
$^2$
in 1 400 MHz band) follows a log-normal distribution (Szary et al. Reference Szary, Zhang, Melikidze, Gil and Xu2014)
\begin{equation} \log L_{\textrm{1 400}} \sim N(0.5, 1.0), -3.0 \leq \log L_{\textrm{1 400}} \leq 4.0\end{equation}
If a pulsar crosses the ‘death-line’, described in Section 2.5, before the observation is made, then the pulsar is not going to be observed. As the pulsar ages, its radio emission efficiency
$\xi$
as defined in Eq. (2) will increase. As per Szary et al. (Reference Szary, Zhang, Melikidze, Gil and Xu2014), for a pulsar to be detectable in radio,
$\xi \lt 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 (Osłowski et al. 2011; Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Rossi and Flynn2020) to calculate the radiometer equation (Dewey et al. Reference Dewey, Taylor, Weisberg and Stokes1985; Lorimer & Kramer Reference Lorimer and Kramer2004) and determine the minimum detectable flux,
$S_{\textrm{min}}$
, at the sky location at a given signal-to-noise ratio
$(S/N_{\textrm{min}})$
\begin{equation} 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}}, \end{equation}
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,
$t_\textrm{int}$
is the integration time,
$W_e$
is the 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 yr, 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 (e.g. Faucher-Giguère & Kaspi Reference Faucher-Giguère and Kaspi2006; Graber et al. Reference Graber, Ronchi, Pardo-Araujo and Rea2024). In this work, we choose pulsar survey parameters for Eq. (15) to represent the Parkes Multibeam Pulsar Survey (Manchester et al. Reference Manchester2001) and previous binary population synthesis work (Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Bailes and Broekgaarden2021) by setting
$\beta = 1$
,
$S/N_{\textrm{min}} \geq 10$
,
$T_{\textrm{rec}} = 24$
K,
$\Delta f =$
288 MHz,
$G = 0.65$
K/Jy,
$n_p = 2$
, and
$t_\textrm{int} = 2\,100$
s. The sky temperature
$T_\textrm{sky}$
is determined by the location of the pulsar calculated by NIGO as described in Section 2.6, and is calculated by scaling the values of sky temperature at 408 MHz from the all-sky survey by Haslam et al. (Reference Haslam, Salter, Stoffel and Wilson1982) to the frequency of Parkes Multibeam survey of 1 374 MHz, assuming the sky temperature spectral index of −2.5 (Reich & Reich Reference Reich and Reich1988; Hobbs et al. Reference Hobbs2004; Guzmán et al. Reference Guzmán, May, Alvarez and Maeda2011). The pulsar period P is determined by the COMPAS simulation. The Parkes Multibeam Survey discovered around half of the pulsars detected (Manchester et al. Reference Manchester2001; Lorimer et al. Reference Lorimer2006; Sengar et al. Reference Sengar2023).
The expression of
$W_e$
below (Burgay et al. Reference Burgay2003) describes how the interstellar medium (ISM) impacts the observed pulses. Free electrons in the Galaxy broaden the intrinsic pulses
$W_i$
(Cordes & Lazio Reference Cordes and Lazio2002), and the ISM scatters the pulsar beam.
\begin{equation} W_e^2 = W_i^2 + \tau_{\textrm{samp}}^2 + \left(\tau_{\textrm{samp}}\frac{\textrm{DM}}{\textrm{DM}_0}\right)^2 + \tau_{\textrm{scatt}}^2,\end{equation}
where
$W_i$
is the intrinsic width of the pulse,
$\tau_{\textrm{samp}} = 250$
$\unicode{x03BC}$
s the sampling time of Parkes Multibeam Survey, and
$\tau_{\textrm{scatt}}$
ISM scattering times based on a fit with respect to DM from Bhat et al. (Reference Bhat, Cordes, Camilo, Nice and Lorimer2004), DM is the dispersion measure in the direction of the pulsar and DM
$_0$
is the diagonal dispersion measure of the survey. While the duty cycle of pulsars varies significantly (Lyne et al. Reference Lyne, Manchester and Taylor1985), here we assume
$W_i/{P} = 0.05$
for all pulsars for simplicity. 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.
Additionally, only pulsars that are located in the regions observable by the Parkes Multibeam survey are selected. Their Galactic latitude (b) and longitude (l) would satisfy
\begin{equation*} (|b| \leq 5), \& (0 \leq l \leq 50 \ \textrm{or}\ 260 \leq l \lt 360)\end{equation*}
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$
(Tauris & Manchester Reference Tauris and Manchester1998), is considered to be
\begin{equation*} f_r = 0.09 \log{(P/10)}^2 + 0.03,\end{equation*}
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 (Kalapotharakos et al. Reference Kalapotharakos, Harding, Kazanas and Wadiasingh2019, Reference Kalapotharakos, Wadiasingh, Harding and Kazanas2022). In this work we estimate pulsar
$\gamma$
-ray luminosity using a fundamental plane in four-dimensional (4D) space that is fitted to pulsars in the Fourth Fermi-LAT Source Catalogue (4FGL; Abdollahi et al. Reference Abdollahi2020) using 12-yr survey data (4FGL-DR3; Abdollahi et al. Reference Abdollahi2022), which contains 190
$\gamma$
-ray detected pulsars with well-defined distances in the ATNF catalogue. (Kalapotharakos et al. Reference Kalapotharakos, Wadiasingh, Harding and Kazanas2022) find that the
$\gamma$
-ray luminosity is given by,
\begin{equation} L_{\gamma} = 10^{14.3\pm1.3} E_{\textrm{cut}}^{1.39\pm0.17} B^{0.12\pm0.03}\dot{E}^{0.39\pm0.05}, \end{equation}
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 Figures 10 & 11 in Kalapotharakos et al. (Reference Kalapotharakos, Wadiasingh, Harding and Kazanas2022). 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
$\lt$
$\dot{E}$
$\lt 10^{38}$
erg/s. Within this
$\dot{E}$
range, we assume a power-law correlation of
$E_\textrm{cut}$
(in GeV) with
$\dot{E}$
as well as
$\Gamma$
with
$\dot{E}$
simply by following the trends in Figures 10 & 11 in Kalapotharakos et al. (Reference Kalapotharakos, Wadiasingh, Harding and Kazanas2022) with these following correlations
\begin{equation*} \log_{10}E_\textrm{cut} = \log_{10}0.5 + 0.26 (\log_{10}\dot{E}-33) + U(0.05),\end{equation*}
and
\begin{equation*} \Gamma = 0.5 + 0.3 \left(\log_{10} \dot{E} - 33\right) + U(0.05),\end{equation*}
where
$U(0.05)$
is a random number drawn from a uniform distribution between −0.05 and 0.05. These relations were estimated visually to fit the approximate trends. Outside this
$\dot{E}$
range, we assume constants for
$E_\textrm{cut}$
and
$\Gamma$
. When
$\dot{E}$
$\lt 10^{33}$
erg/s, we set
\begin{equation*} \log_{10}E_\textrm{cut} = \log_{10}0.5 + U(0.05),\end{equation*}
and
\begin{equation*} \Gamma = 0.5 + U(0.05).\end{equation*}
When
$\dot{E}$
$\gt 10^{38}$
erg/s,
$\log_{10}E_\textrm{cut} = 1.0 + U(0.05)$
,
$\Gamma = 0.5 +U(0.05)$
.
We can then determine the pulsar’s energy flux by
\begin{equation*} F_{\gamma} = \frac{L_{\gamma} }{4\pi d^2 c_g},\end{equation*}
where d is the distance of the pulsar to the Earth and
$c_g$
is the flux correction factor (Watters et al. Reference Watters, Romani, Weltevrede and Johnston2009), which depends on the emission mechanism of the pulsar as well as geometry and viewing angle. Here we combine results from Johnston et al. (Reference Johnston, Smith, Karastergiou and Kramer2020) and Pétri (Reference Pétri2011), and assume that
$c_g$
follows a sigmoid function versus
$\dot{E}$
with a maximum of 0.8 (Johnston et al. Reference Johnston, Smith, Karastergiou and Kramer2020) and a minimum of 0.22 (Pétri Reference Pétri2011), 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 14-year Fermi-LAT Source Catalogue provides an all-sky map of the detection threshold for the catalogue.Footnote d By comparing the energy flux of a pulsar in the simulation with 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 (Abdo et al. Reference Abdo2013; Smith et al. Reference Smith2023). Figure 17 of Abdo et al. (Reference Abdo2013) clearly indicates that for the detected pulsars in the 2nd Fermi-LAT pulsar catalogue, pulsars that are already detected in radio in general are easier to detect. In this work, 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. Evaluating Figure 17 in Abdo et al. (Reference Abdo2013), we choose
$\alpha = 3.0$
and
$\beta = 1.5$
.
2.8.3. Gamma-ray beaming
Gamma-ray pulsars have a beaming effect similar to that of radio pulsars. In Johnston et al. (Reference Johnston, Smith, Karastergiou and Kramer2020), a description of young and energetic (with
$\dot{E}$
$\gt 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 Johnston et al. (Reference Johnston, Smith, Karastergiou and Kramer2020) does not provide a beaming fraction for pulsars with
$\dot{E}$
$\lt 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
\begin{align*} f_g =\begin{cases} 0.92,& \text{if }\ \dot{E} \gt 10^{38}\,\textrm{erg s}^{-1}\\[3pt] 0.82,& \text{if }\ 10^{37} \lt \dot{E} \leq 10^{38} \textrm{erg s}^{-1}\\[3pt] 0.67,& \text{if }\ 10^{36} \lt \dot{E} \leq 10^{37} \textrm{erg s}^{-1}\\[3pt] 0.50,& \text{if }\ 10^{35} \lt \dot{E} \leq 10^{36} \textrm{erg s}^{-1}\\[3pt] 0.40,& \text{if }\ 10^{33} \lt \dot{E} \leq 10^{35} \textrm{erg s}^{-1}\\[3pt] 0.30,& \text{if }\ \dot{E} \leq 10^{33} \textrm{erg s}^{-1}\end{cases}\end{align*}
The assumptions above are made based on observations and outer gap model of young, powerful pulsars (Johnston et al. Reference Johnston, Smith, Karastergiou and Kramer2020). According to Pierbattista et al. (Reference Pierbattista, Grenier, Harding and Gonthier2012), 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 Section 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.
Table 3. Columns 2–7: p-values of Mann-Whitney U Test, comparing the distribution of physical quantities obtained from simulations and catalogue for radio pulsars. Quantities listed in the second to the second last columns are, respectively, radio flux, period,
$\dot{P}$
,
$\dot{E}$
, Galactic latitude, and Galactic longitude. The last column: the averaged number of the physical quantities with p-values
$\gt 0.05$
out of all 10 realisations of each model. For a given model, the p-value for each parameter is determined using the random realisation of the model that is shown in Figure 1.
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 and discussions
The results of the COMPAS simulations using the setups from Section 2 are presented in this section, as well as the discussion on the implications of these results. The results are compared to the pulsar catalogue 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 employ a statistical method to compare the distribution of pulsar physical properties to those reported in the catalogues, then compare the number of different types of pulsars produced in different simulation models. These results are combined to choose a model that best describe the birth distributions of pulsar properties. 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. Pulsar birth rate
In this section, we calculate the pulsar birth rate (PBR) in our COMPAS simulations. Only binaries where the primary (initially more massive) star is within the initial mass range of 5–150 M
$_{\odot}$
are modelled using COMPAS, since lower mass stars are not expected to produce pulsars. Because of this setup, we need to account for the low-mass end of the initial mass function (IMF) when calculating the pulsar birth rate. By default, we use a Kroupa (Reference Kroupa2001) broken power-law IMF, with a slope of 0.3 between 0.01 and 0.08 M
$_{\odot}$
; slope of 1.3 between 0.08 and 0.5 M
$_{\odot}$
; slope of 2.3 between 0.5 and 200 M
$_{\odot}$
. To correct for the lower mass stars (0.08 M
$_{\odot}$
$\lt M \lt$
5.0 M
$_{\odot}$
), we use the CosmicIntegration Python package in COMPAS to calculate a correction factor, which is the ratio of the true average mass produced per binary and the averaged mass produced per binary in COMPAS. Using the implementation of the Kroupa IMF in CosmicIntegration, we draw a sample of
$10^{9}$
binaries with the mass range between 0.08 to and calculate the average mass per binary to be
$\bar{M}_\textrm{true} = 1.33$
M
$_{\odot}$
. Next we calculate the average mass per binary in COMPAS simulations. Since all the COMPAS models in this study use the same parameters for initial mass and orbit for the binaries and the same evolution physics, they all have the same pulsar birth rate. In each COMPAS simulation, the total mass simulated in binaries within the 5–150 M
$_{\odot}$
range is
$M_\textrm{tot, COMPAS} = 8.05 \times 10^6$
M
$_{\odot}$
. The average mass formed per binary in one COMPAS simulation is
$\bar{M}_\textrm{COMPAS} = 20.3$
M
$_{\odot}$
. The correction factor
$\kappa$
is calculated as
$\kappa = \bar{M}_\textrm{COMPAS} / \bar{M}_\textrm{true} = 15.14$
.
With
$N_\mathrm{psr} = 1.2 \times 10^5$
neutron stars produced in each simulation, the pulsar yield (the number of pulsars born per Solar mass of star formation) is calculated as,
\begin{equation} \zeta_\mathrm{PSR} = N_\mathrm{psr} / (M_\textrm{tot, COMPAS} \times \kappa) = 1.3 \times 10^{-3}\,\mathrm{M}_\odot^{-1} . \end{equation}
To translate the pulsar yield into a PBR, we have to assume a star formation history. The observed star formation rate (SFR) over the most recent history of the Milky Way is around
$1.65 \pm 0.19$
M
$_{\odot}$
yr
$^{-1}$
(Licquia & Newman Reference Licquia and Newman2015),
$1.9 \pm 0.4$
M
$_{\odot}$
yr
$^{-1}$
(Chomiuk & Povich Reference Chomiuk and Povich2011) or
$3.3^{+0.7}_{-0.6}$
M
$_{\odot}$
yr
$^{-1}$
(Zari, Frankel, & Rix Reference Zari, Frankel and Rix2023). Other works, such as Evans et al. (Reference Evans, Kim and Ostriker2022), estimate the Milky Way star formation rate in the range of 0.50–5.93 M
$_{\odot}$
yr
$^{-1}$
. The pulsar birth rate scales linearly with the assumed SFR. Here we assume the star formation rate in the recent Milky Way history is
$\mathrm{SFR} = 1.65$
M
$_{\odot}$
yr
$^{-1}$
, and the pulsar birth rate is then,
\begin{equation} \mathrm{PBR} = 2.1 \times 10^{-3}\, \mathrm{yr}^{-1} \left( \frac{\mathrm{SFR}}{1.65\, \mathrm{M}_\odot \, \mathrm{yr}^{-1}} \right) \left( \frac{\zeta_\mathrm{PSR}}{1.3 \times 10^{-3}\, \mathrm{M}_\odot^{-1}} \right) \end{equation}
or 1 pulsar born every 466 yr. This value is lower than most of the literature, e.g. 1/50 yr (Keane & Kramer Reference Keane and Kramer2008), 1/46–1/58 yr (Graber et al. Reference Graber, Ronchi, Pardo-Araujo and Rea2024), 1/70 yr (Dirson et al. Reference Dirson, Pétri and Mitra2022), and 1/100 yr (Johnston et al. Reference Johnston, Smith, Karastergiou and Kramer2020). The main reason for this difference could be due to some of the formation channels not included in this work, such as AIC, stellar mergers and double white dwarf mergers. Previous COMPAS models (e.g. Stevenson et al. Reference Stevenson2019) have also found the predicted rate of core collapse supernova (CCSN) to be somewhat low compared to the observed rate. Increasing the SFR would also increase the pulsar birth rate in our work, with a detailed exploration on this effect described in Section 3.5.
Figure 1. CDFs of physical quantities for radio pulsars produced by the COMPAS models described in Table 1 (coloured), compared to the data from the catalogues (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
$\dot{E}$
, radio flux, period,
$\dot{P}$
, Galactic latitude, and Galactic longitude.
3.2. Results: Radio pulsars
We start the discussion of the results of the simulations by comparing the predicted and observed distributions of various physical quantities of the pulsar populations between the catalogue and the simulations. We first discuss radio pulsars, and then discuss
$\gamma$
-ray pulsars in Section 3.3.
For model comparison purposes, we select Galactic radio pulsars that are not associated with any globular clusters, have positive radio flux at 1 400 MHz and have an association with Parkes. Only pulsars with
$P \gt 0.03$
s and
$\dot{P} \gt 10^{-17}$
s/s are chosen, to ensure that the sample does not include pulsars that are or have been going through binary interactions/recycling. All pulsars are chosen from the ATNF Pulsar catalogue v2.5.1.
The cumulative distribution functions (CDFs) for each parameter of radio pulsars for each of the models listed in Table 1 are plotted in Figure 1.Footnote
e
Overall, several COMPAS models yield distributions of most physical quantities that align well with the observations. The largest variation is seen in the distributions of P and
$\dot{P}$
.
We use statistical tests to examine how well the predicted parameter distributions match the observed distributions for each physical parameter of the pulsar populations. These physical parameters include: radio flux,
$\gamma$
-ray flux, spin period P, spin-down rate
$\dot{P}$
, spin-down luminosity
$\dot{E}$
, Galactic longitude (l) and latitude (b). Given that the observed data and simulations have different-sized samples, and the physical parameters for each pulsar are independent measurements that is not influenced by measuring a different pulsar, we use Mann-Whitney U test (MWU test; Mann & Whitney Reference Mann and Whitney1947) 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 catalogue. If
$p \lt 0.05$
, it means we reject the null hypothesis that the two compared samples are drawn from the same distribution. When
$p \gt 0.05$
, we cannot reject the null hypothesis, which means we consider the two distributions to be the same.
To account for the stochastic variation due to the limited sample size from the models, the observation part of post-processing is repeated 10 times. In Table 3, the p-values of the MWU test are shown for a selection of the physical quantities of the pulsars compared between each model and the catalogue. The p-values listed in Table 3 are based on the random realisation of the model shown in Figure 1. To further quantify model performance, we introduce the metric NO, defined as the average number of physical quantities (out of six considered) for which the MWU test p-value exceeds 0.05 across the 10 realisations. A higher NO value indicates better agreement with the observed distributions, with a maximum possible score of 6. Notable trends and features from Table 3 are discussed below.
The spatial distributions of the simulated radio pulsars closely match those in the catalogue. This is evident as the p-values for b across all models, and for l in all but the D4 and P2-D3 models, are larger than 0.05, meaning the simulation results align well with the catalogue distributions shown in Figure 1. This is also an indication that the Galactic models are well-suited for representing radio pulsars.
Only six models successfully reproduce the distributions of radio flux: D4, P2-D4, PN1-B1-D4, PN5-B1-D4, BLN5-PLN2-D4, and BLN6-PLN3-D4. Seven models capture the spin period distribution, nine for the
$\dot{P}$
distribution, and four for the
$\dot{E}$
distribution. Further examination of the CDF plots for these properties indicates that models B1, BLN1, BLN1-PN1, PN1, D1, D5, P2, PN1-B1, PN1-B1-D2, PN3, PN5-B1-D4, BLN5-PLN2-D4, and BLN6-PLN3-D4 visually exhibit a good fit.
Based solely on the MWU test p-values, model PN5-B1-D4 best matches the observed radio population. It achieves an average NO score of 3.9, meaning that on average, 3.9 out of 6 physical quantities have p-values exceeding 0.05 across the 10 realisations. The next best-performing models are BLN6-PLN3-D4 with NO score of 3.6, BLN5-PLN2-D4 with NO score of 3.4, then BLN4-PN1 and PN1-B1-D4 each with NO score of 3.0. Close behind are models PN5-B1-D5, PN1-B1-D2, P2-D4, P3, P1, and BLN4, each with an average of 2.9 quantities showing p-values above 0.05. Interestingly, despite good visual comparisons, models PN1-B1-D4, BLN4-PN1, P2-D4, P3, P1 and BLN4 do not exhibit a high number of physical quantities p-values exceeding 0.05. All the models mentioned earlier in this paragraph present a mixture of long (PN5 series models and P3) and short (PN1 series models) initial spin periods; large (INIT,
$\tau_d = 500$
Myr series models) and small (D2/D4 series models) magnetic decay timescale. All the models mentioned earlier in this paragraph favour larger initial magnetic field strengths.
3.3. Results:
$\gamma$
-ray pulsars
The results of modelling
$\gamma$
-ray pulsars, and the comparisons to the catalogue are presented in this section. The catalogue
$\gamma$
-ray pulsars are chosen from 4FGL-DR4, the 14-Year Version of the 4th Fermi-LAT Source catalogueFootnote
f
(Ballet et al. Reference Ballet, Bruel, Burnett and Lott2023). Only Galactic
$\gamma$
-ray pulsars not associated with any globular clusters that have
$P \gt 0.03$
s and
$\dot{P} \gt 10^{17}$
s/s are chosen for model comparison. To ensure consistency, both the catalogue and simulated
$\gamma$
-ray pulsars are selected from the regions covered by the Parkes Multibeam survey. The simulated
$\gamma$
-ray pulsars are compared to the catalogue
$\gamma$
-ray pulsars using the same approach applied to radio pulsars, as described in Section 3.2. For the
$\gamma$
-ray pulsars, the same random realisations of each model, as used in Section 3.2, are selected to plot the CDFs in Figure 2 and to list the MWU test p-values in Table 4.
Figure 2. CDFs of physical quantities for
$\gamma$
-ray pulsars produced by the COMPAS models described in Table 1 (coloured), compared to the data from the catalogues (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
$\dot{E}$
, radio flux, period,
$\dot{P}$
, Galactic latitude (b), and Galactic longitude (l).
Table 4. Columns 2–7: p-values of Mann-Whitney U Test, comparing the physical quantities obtained from simulations and catalogue for
$\gamma$
-ray pulsars. Quantities listed in the second to the second last columns are, respectively,
$\gamma$
-ray flux, period,
$\dot{P}$
,
$\dot{E}$
, Galactic latitude, and Galactic longitude. The last column: the averaged number of the physical quantities with p-values
$\gt 0.05$
out of all 10 realisations of each model. For a given model, the p-value for each parameter is determined using the random realisation of the model that is shown in Figure 2.
The first notable observation is that the simulated
$\gamma$
-ray pulsar populations appear to more closely resemble the catalogue, compared to the modelled radio pulsar populations. This is reflected in the greater number of physical quantities for
$\gamma$
-ray pulsars with p-values larger than 0.05. Specifically, 31 models achieve a fit to
$\gamma$
-ray flux with p-values
$\gt 0.05$
, while 23 models do so for spin period, 14 for
$\dot{P}$
, and 18 for
$\dot{E}$
. These values are all higher than those reported for the results from the same models on radio pulsars in Table 3. While it is possible that our
$\gamma$
-ray models are better than the radio ones, another possibility is that the smaller observational sample of
$\gamma$
-ray pulsars and the lower number of
$\gamma$
-ray pulsars produced in these models (see Section 3.5, especially Table 5) make it easier to fit the simulation to observation, resulting in higher p-values due to small number statistics. As with the radio pulsars, most models produce distributions of physical parameters that are close to those of the observation. However, greater scatter is present in Figure 2, often attributable to low number statistics. Examining Table 4 indicates that many of these models provide a good fit to the catalogue.
Based solely on the MWU test p-values, models BLN1-PN5, P3 and PN5-B1-D2 exhibit average of more than 5.5 quantities with p-values larger than 0.05. Models P1, PN5-B1-D5, PN5-B1-D3 and PN5-B1-D1 have an average of over 5.0, while models D5, PN5-B1-D4, PN1-B1-D1, P2 and PN1 have averages larger than 4.5, and models INIT, D1, BLN1PN1, PN1B1D4, PN1B1, BLN5-PLN2-D4 and BLN6-PLN3-D4 have averages above 4.0. Once again, similar to the case of radio pulsars, these models cover both long and short initial spin periods, as well as large and small
$\tau_d$
.
3.4. Comparing predicted and observed numbers of pulsars
Another straightforward way to determine if a COMPAS model is a good fit for the observation is to compare the number of pulsars predicted by the model to those listed in the catalogues. We evaluate, for each model, the following quantities: the number of radio pulsars (N
$_\textrm{R}$
), the number of
$\gamma$
-ray pulsars in the Parkes Multibeam Survey covered regions (N
$_\textrm{G}$
), the number of
$\gamma$
-ray pulsars in the entire sky (N
$_\textrm{G, all-sky}$
), the number of radio-loud
$\gamma$
-ray pulsars in the Parkes Multibeam Survey covered regions (N
$_\textrm{RLGL}$
). These quantities are listed in columns 2–5 in Table 5. As discussed before, we repeat the observation part of post-processing 10 times in order to account for stochastic variation. The average and standard deviation of these realisations are recorded for each model in Table 5.
For many of the models, the predicted number of pulsars differs significantly from the observations. For example, model B1 produces 18.3% fewer radio pulsars, and 49.4% fewer
$\gamma$
-ray pulsars, respectively. A more extreme example is model B2, predicting 4.2 times as many radio pulsars as observed. Due to these large differences, we do not calculate the percent difference with the numbers reported in columns 2–5 in Table 5. Instead, we calculate a scaled percent difference (SPD) using scaled numbers of all types of pulsars. The scaling factor is defined as
$F = 1\,060/{\textrm{N}_\textrm{R}}$
, with 1 060 is the number of observed radio pulsars from the Parkes Multi Beam survey, as listed in the sixth column in Table 5 for each model. These scaled numbers of
$\gamma$
-ray pulsar in Parkes covered regions,
$\gamma$
-ray pulsars observed in the entire sky, and radio-loud
$\gamma$
-ray pulsars are then compared to those from the catalogues to calculate the SPD. Here we give the equations on how these SPDs are calculated:
\begin{equation*} {\textrm{SPD}_{\textrm{R}}} = 100 \times \frac{\textrm{N}_{\textrm{R}} \times F - \textrm{N}_{\textrm{R,0}}}{\textrm{N}_{\textrm{R,0}}},\end{equation*}
where
${\textrm{N}_{\textrm{R,0}}}= 1\,060$
is the number of radio pulsars in Parkes covered regions in the catalogue;
\begin{equation*} {\textrm{SPD}_{\textrm{G}}} = 100 \times \frac{\textrm{N}_{\textrm{G}} \times F - \textrm{N}_{\textrm{G,0}}}{\textrm{N}_{\textrm{G,0}}},\end{equation*}
where
${\textrm{N}_{\textrm{G,0}}} = 85$
is the number of
$\gamma$
-ray pulsars in Parkes covered regions in the catalogue;
\begin{equation*} {\textrm{SPD}_{\textrm{G,all-sky}}} = 100\times \frac{\textrm{N}_{\textrm{G,all-sky}} \times F - \textrm{N}_{\textrm{G,all-sky,0}}}{\textrm{N}_{\textrm{G,all-sky,0}}},\end{equation*}
where
${\textrm{N}_{\textrm{G,0}}} = 156$
is the number of
$\gamma$
-ray pulsars in the entire sky in the catalogue;
\begin{equation*} {\textrm{SPD}_{\textrm{RLGL}}} = 100\times \frac{\textrm{N}_{\textrm{RLGL}} \times F - \textrm{N}_{\textrm{RLGL,0}}} {\textrm{N}_{\textrm{RLGL,0}}},\end{equation*}
where
${\textrm{N}_{\textrm{RLGL,0}}} = 38$
is the number of radio-loud
$\gamma$
-ray pulsars in the Parkes covered region in the catalogue. The SPD values are tabulated in columns 7–10 in Table 5. By design, all models would have a SPD
$_\textrm{R}$
= 0 and their SPD
$_\textrm{R}$
’s are hence not listed in Table 5. This scaling practice is equivalent of setting the SFR as a free parameter. As a result, the SFR for each model will be changed to
$F \times 1.65$
M
$_{\odot}$
yr
$^{-1}$
and the PBR is changed to
$F \times 2.15\times10^{-3}$
yr
$^{-1}$
, or 1 in every
$446/F$
yr, which are listed in the second last and last columns in Table 5 respectively.
Another indicator of the goodness of the models is the ratio between different types of pulsars. We check the ratio between radio and
$\gamma$
-ray pulsars predicted by the models in Parkes covered region,
${\mathrm{Q_1} = {\textrm{N}_{\textrm{R}}}}/{\textrm{N}_{\textrm{G}}}$
; and the ratio between all
$\gamma$
-ray pulsars and radio-loud
$\gamma$
-ray pulsars predicted by the models in the Parkes covered region,
${\mathrm{Q_2} = {\textrm{N}_{\textrm{G}}}}/{\textrm{N}_{\textrm{RLGL}}}$
. Percent differences for
${\mathrm{Q_1}}$
and
${\mathrm{Q_2}}$
when compared to the observed values respectively are calculated and present in columns 11 and 12 in Table 5.
A striking observation is that the PN5 series models, including PN5-B1-D1–5 and BLN1-PN5 significantly underproduce
$\gamma$
-ray pulsars. This indicates that the initial spin periods of
$\gamma$
-ray pulsars should not be as long as described in PN5 models.
Another prominent feature is that all models, except D4 and P2-D4, underproduce
$\gamma$
-ray pulsars to a certain degree. It could be because of the selection criteria for
$\gamma$
-ray pulsars is too stringent, and this is discussed in Section 3.7. A strange phenomenon is that the degree of underproduction of
$\gamma$
-ray pulsars decreases when comparing all-sky
$\gamma$
-ray pulsars, instead of Parkes covered
$\gamma$
-ray pulsars.
Based on the percent differences across the various pulsar types, we consider models with all percent differences below 60% to be satisfactory. With this criterion we identify the following models to be reasonable: D2, D4, PN2, PLN1, BLN1, BLN3, P2-D2, P2-D3, P2-D4, P2-D5, PN1-B1-D2, PN1-B1-D4, BLN5-PLN2-D4 and BLN6-PLN1-D4. Among these, PN1-B1-D4 and BLN5-PLN2-D4 exhibit the lowest SPDs across all categories of pulsar types, followed by D2, P2-D4, PN1-B1-D2, and BLN6-PLN3-D4. These results indicate that models with lower
$\tau_d$
produces the right number of pulsars. However, some of the other models with larger
$\tau_d$
also perform reasonably well in this regard. The other models listed above with reasonable numbers of different types of pulsars suggest that models with high birth magnetic field and short birth spin periods are more likely to produce the right numbers of different types of pulsars.
Table 5. Simulation results showing numbers of different types of pulsars, the corresponding scaled percent differences, the star formation rate and the pulsar birth rate. Columns from left to right: COMPAS model; pulsars that are detected in radio (
$N_\textrm{R}$
); pulsars detected in
$\gamma$
-ray in Parkes Multibeam Survey covered region (
$N_\textrm{G}$
); pulsars detected in
$\gamma$
-ray in the entire sky (
$N_\textrm{G, all-sky}$
); pulsars detected in both radio and
$\gamma$
-ray in Parkes Multibeam Survey covered region (
$N_\textrm{RLGL}$
); scaling factor for the model as described in text (F); scaled percent difference between catalogue and simulated radio pulsars (SPD
$_\textrm{R}$
); scaled percent difference between catalogue and simulated
$\gamma$
-ray pulsars (SPD
$_\textrm{G}$
); scaled percent difference between catalogue and simulated all-sky
$\gamma$
-ray pulsars (SPD
$_\textrm{G, all-sky}$
); scaled percent difference between catalogue and simulated radio-loud
$\gamma$
-ray pulsars (SPD
$_\textrm{RLGL}$
); percent difference between catalogue and simulated ratio of radio over
$\gamma$
-ray pulsars (PD
$_\textrm{Q1}$
); percent difference between catalogue and simulated ratio of all
$\gamma$
-ray pulsars over radio-loud
$\gamma$
-ray pulsars (PD
$_\textrm{Q2}$
); star formation rate of the model; inverse of pulsar birth rate (yr). The first row lists the total number of pulsars in the radio and
$\gamma$
-ray catalogues (Manchester et al. Reference Manchester, Hobbs, Teoh and Hobbs2005; Smith et al. Reference Smith2023). MSPs are excluded as we focus on canonical pulsars. Details on the calculations of percent differences are given in Section 3.4. For the simulations, the values and uncertainties are obtained from the average and standard deviation of 10 different realisations of each model.
Figure 3. CDFs of various physical quantities predicted/produced by the PN1-B1-D4 model (red), compared to the data from the catalogues (blue). The red lines in these plots are from all ten realisations of the PN1-B1-D4 model that are sampled randomly as described in Section 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. The p-values for each individual parameter as listed in Tables 3 and 4 for one random realisation are plotted on each respective panel.
3.5. Best fit model
Combining the results presented in Tables 3, 4 and 5, we identify models that best describe both the radio and
$\gamma$
-ray pulsar populations. Ideal models need to explain the observed distributions of the physical properties and reproduce the numbers of different pulsar types. However, no single model satisfies all three criteria perfectly. Among the models, PN1-B1-D4 is the best fit, as it is listed as one of the top models for both the radio and
$\gamma$
-ray populations. This model produces the numbers of different pulsar types and their ratios that are comparable to the observations. The star formation rate is 3.43 M
$_{\odot}$
yr s
$^{-1}$
, agreeing with the values listed in Section 3.1. The pulsar birth rate of 1 per 224 yr, while low compared to the literature values, is comparable to the rate of CCSN in COMPAS models, as described in Section 3.1. A degeneracy between the number of pulsars formed in the model and the fraction of observable ones is expected. With higher pulsar birth rate, a smaller fraction is needed to produce the same amount of observable pulsars. This effect will be explored in the future once the issue of the low CCSN rate has been addressed in COMPAS. Model BLN5-PLN2-D4 is a strong contender for the best fit model, describing the distributions of physical parameters of both radio and
$\gamma$
-ray pulsars well, and having the lowest percent differences in almost all types of pulsars and the two ratios. However, the SFR at 7.49 M
$_{\odot}$
yr s
$^{-1}$
is considerably high compared to the literature values, even though the resulting pulsar birth rate at 1 per 103 yr is much closer to the literature values as listed in Section 3.1. This rules it out for consideration of the best fit model.
Figure 4. Top: P-
$\dot{P}$
diagram of Fermi detected canonical pulsars (blue) and detected
$\gamma$
-ray pulsars from the characteristic PN1-B1-D4 model (red); bottom: P-
$\dot{P}$
diagram of radio detected canonical pulsars as recorded in the ATNF catalogue (blue) and detected radio pulsars from the characteristic PN1-B1-D4 model (red). The two black dashed lines represent the death-lines described in Eqs. (7) and (8). The gray dashed lines labelled
$10^3$
,
$10^6$
, and
$10^9$
yr represent constant characteristic age each respectively; and those labelled with
$10^{31}$
,
$10^{33}$
, and
$10^{37}$
erg s
$^{-1}$
represent constant
$\dot{E}$
each respectively.
To further examine model PN1-B1-D4, we plot the CDFs of the quantities listed in Tables 3 and 4 in Figure 3, showing the results from all 10 realisations. The CDFs of some of the parameters, such as spin period of
$\gamma$
-ray pulsars, are well matched to the observations, while some others such as spin period of radio pulsars, are not. This is expected when examining Tables 3 and 4.
Additionally, we visualise the results with the P-
$\dot{P}$
diagrams for
$\gamma$
-ray and radio pulsars generated from one randomly chosen realisation of the PN1-B1-D4 model, overlaid with the canonical pulsars from the catalogue, as shown in Figure 4. These diagrams provide a clearer comparison between the modelled and observed populations. For
$\gamma$
-ray pulsars, the simulation and catalogue pulsar populations occupy the same part of the parameter space. The simulated radio pulsars, compared to the observation, are shifted towards the bottom left of the P-
$\dot{P}$
diagram. This suggests that the simulated radio pulsars have lower magnetic fields compared to observation. This is anticipated given the PN1-B1-D4 model has a small
$\tau_d$
, letting the magnetic field to decay faster, resulting in pulsars with lower field and shorter spin period.
We show the spatial coordinates of detected pulsars in Figures 5 and 6. In both Figures 5 and 6, the spatial distributions of pulsars in the
$X-Y$
,
$X-Z$
and
$Y-X$
planes have good overlap between the observation and the simulation, with a slight over density around the Galactic Centre region. This over density is the most obvious in the R-Z plot at the bottom right panel in both Figures 5 and 6. We attribute this to the initial spatial distributions of massive stars, but overall the Galactic model used in this work is sufficiently good.
3.6. Magnetic decay timescale
Magnetic decay timescale (
$\tau_d$
) is a key parameter for pulsar evolution models. As mentioned in Section 2.4, there is no consensus on the value of
$\tau_d$
. Some studies, mentioned in Section 2.4, suggest short
$\tau_d$
, indicating rapid magnetic field decay; while others prefer long
$\tau_d$
, indicating a slow-changing or even constant magnetic field throughout the life time of a NS. In this section, we use population synthesis methods to examine both perspectives. As shown in Table 1, four sets of models are examined: INIT, P2, PN1-B1 and PN5-B5, each with specific birth magnetic field and spin period distributions and
$\tau_d$
ranging from 1 to 10 000 Myr.
For radio pulsars, shorter
$\tau_d$
values (1–10 Myr) generally produce better matches to observations across most model sets, while longer
$\tau_d$
values (100–10 000 Myr) tend to yield worse fits. The caveat of having shorter
$\tau_d$
for radio pulsars is that the models tend to generate more low field pulsars, resulting with a population with shorter spin periods. For
$\gamma$
-ray pulsars, moderate to long
$\tau_d$
values (e.g.
$\gt100$
Myr) better match the observed population, as slow-decaying or effectively non-decaying magnetic fields are more consistent with the high
$\dot{E}$
s and low characteristic ages of these pulsars. Notably, model set PN1-B1 stands out for its strong performance in describing radio pulsar populations, with
$\tau_d = 1$
Myr emerging as the best fit.
Overall, the results favour shorter
$\tau_d$
values (1–10 Myr) for describing the population of neutron stars, particularly radio pulsars, as these models consistently align with observed distributions of pulsar properties, as well as the numbers and ratios of different pulsar types. However,
$\gamma$
-ray pulsars require longer
$\tau_d$
values for optimal fits.
One caveat is that this result is based on the assumption that the magnetic field always decays exponentially with time, as described in Section 2.3. The exponentially decaying field in our model suggests that different decay timescales are needed to explain older radio and younger
$\gamma$
-ray pulsars. However, accounting for the Hall effect at earlier times, as mentioned in Section 2.3, would bridge this issue. This would require modelling the magnetic field decay with different prescriptions (e.g. Colpi, Geppert, & Page Reference Colpi, Geppert and Page2000; Dall’Osso, Granot, & Piran Reference Dall’Osso, Granot and Piran2012) and it is planned for future work (see Section 4).
Figure 5. Galactocentric coordinates of radio pulsars from the catalogue (blue) and from PN1-B1-D4 model (red). The yellow star in all panels except for the
$Y-Z$
planes indicates the location of the Sun.
Figure 6. Galactocentric coordinates of
$\gamma$
-ray pulsars from the catalogue (blue) and from PN1-B1-D4 model (red). The yellow star in all panels except for the Y-Z planes indicates the location of the Sun.
3.7. Impact of
$\gamma$
-ray beaming and sensitivity
Two model choices can impact the simulation results:
$\gamma$
-ray beaming and the
$\gamma$
-ray sensitivity limit.
Our initial suspicion for the underproduction of
$\gamma$
-ray pulsars was that it was due to the beaming factor (as described in Section 2.8.3). Since
$\gamma$
-ray beaming is applied as a probabilistic rejection scheme, variations in the beaming will impact the numbers of
$\gamma$
-ray pulsars. We made modifications to the beaming prescription outlined as follows: Beaming 1 option is the one described in Section 2.8.3. Beaming 2 assumes the beaming fraction remains high to a floor of f
$_g$
= 0.7 for all pulsars with
$\dot{E} \lt 10^{35}$
erg s
$^{-1}$
. For pulsars with
$10^{35}$
erg s
$^{-1} \lt \dot{E} \lt 10^{36}$
erg s
$^{-1}$
, f
$_g$
= 0.75. For pulsars with
$10^{36}$
erg s
$^{-1} \lt \dot{E} \lt 10^{37}$
erg s
$^{-1}$
, f
$_g$
= 0.8. Beaming 3 decreases the beaming fraction of pulsars by steps in the same way as Beaming 1 until a floor at f
$_g$
= 0.4 at
$\dot{E}$
= 10
$^{35}$
erg s
$^{-1}$
. These beaming models are simply ad hoc generalisation made to compare different prescriptions. All three beaming options are plotted in Figure 7.
The sensitivity limit for radio-detected
$\gamma$
-ray pulsars, as described in Section 2.8.2, can also affect the final results. Specifically, when we change our assumption about the sensitivity limit, the detected population in the model will also change. As stated in Section 2.8.2, we deviate from the standard 4FGL sensitivity based on the radio detectability of the
$\gamma$
-ray pulsar. Here we explore the effects of varying
$1.0 \leq \alpha \leq 4.0$
and
$1.0 \leq \beta \leq 4.0$
, then re-run the analysis pipelines with all three beaming options on the PN1-B1-D4 model. We compare the resulting number of
$\gamma$
-ray pulsars produced in each instance (N
$_{G}$
) and the ratio of
$\gamma$
-ray pulsars and radio-loud
$\gamma$
-ray pulsars (q) to the catalogue values in Figure 8.
The results of this analysis are plotted in Figure 8 with contours of N
$_G$
and q for different combinations of
$\alpha$
and
$\beta$
for all three beaming options. Assuming the catalogue contains the same number of radio pulsars as predicted by the PN1-B1-D4 model (509; see Table 5), the number of catalogue
$\gamma$
-ray pulsars is scaled down proportionally from 85 to 41, reflecting the same reduction ratio as applied to the radio pulsars (from 1 060 to 509), to ensure a fair comparison. N
$_G$
= 41, and the catalogue value of
$q = 2.24$
are both plotted with dashed lines in the contours. Firstly, the impact of beaming is examined. Beaming 2 increases number of
$\gamma$
-ray pulsars for all combinations of
$\alpha$
and
$\beta$
, whilst showing a slight decrease in q in most cases. Beaming 3 shows marginal difference from Beaming 1.
For a given beaming option, when
$\alpha$
increases, fewer
$\gamma$
-ray pulsars are detected, and q decreases. In this case, the decrease of N
$_G$
can be attributed to radio-quiet pulsars being harder to detect in
$\gamma$
-rays without existing timing solutions. When
$\beta$
increases, more
$\gamma$
-ray pulsars are detected, but q decreases. This can be attributed to more radio pulsars being detected in
$\gamma$
-rays. In all three different beaming options, high
$\alpha \gtrsim 2.8$
is required to achieve the catalogue values of q and
$N_G$
. This agrees with the initial assumption that according to the 2PC/3PC catalogues (Abdo et al. Reference Abdo2013; Smith et al. Reference Smith2023), Fermi-LAT is less sensitive to pulsars discovered without known radio counterparts (e.g. through blind searches). In beaming options 1 and 3, high
$\beta$
is required to achieve the catalogue values, which suggests that pulsars with prior radio detection is much easier to be detected in Fermi-LAT. In beaming option 2, while only a low
$\beta$
is required, the beaming fractions for pulsars with different
$\dot{E}$
are all raised significantly. This is another indication that radio pulsars should be easier to detect in
$\gamma$
-ray to agree with the catalogues. In comparison, if choosing the unaltered Fermi-LAT sensitivity with
$\alpha = \beta = 1$
, the models overproduce
$\gamma$
-ray pulsars, and also produce more radio quiet pulsars indicated by higher q.
Future development in better physical understanding of
$\gamma$
-ray beaming is required. It is also important to establish a correlation between the radio and
$\gamma$
-ray detection limit from the current population, which can further improve results for future population study.
Figure 7. Three models for the beaming of
$\gamma$
-ray pulsars as a function of
$\dot{E}$
as described in Section 3.7.
3.8. Low
$\dot{\boldsymbol{E}}$
pulsars in
$\gamma$
-rays
Theoretical studies of
$\gamma$
-ray emission from pulsars, such as the equatorial current sheet model (Pétri Reference Pétri2012; Kalapotharakos et al. Reference Kalapotharakos, Brambilla, Timokhin, Harding and Kazanas2018; Hakobyan, Philippov, & Spitkovsky Reference Hakobyan, Philippov and Spitkovsky2023 or the outer gap model (Cheng, Ho, & Ruderman Reference Cheng, Ho and Ruderman1986a, b; Romani Reference Romani1996; Takata et al. Reference Takata, Wang and Cheng2011), seem to indicate that pulsars with
$\dot{E}$
$\lt 10^{33}$
erg/s do not emit in
$\gamma$
-rays. 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}$
$\lt 10^{33}$
erg s
$^{-1}$
, if we still follow the fundamental plane relation, the
$\gamma$
-ray luminosity of the pulsar becomes larger than
$\dot{E}$
, breaking conservation of energy. It is also worth noting that according to Kalapotharakos et al. (Reference Kalapotharakos, Harding, Kazanas and Wadiasingh2019, Reference Kalapotharakos, Wadiasingh, Harding and Kazanas2022), the fundamental plane relations are rooted in the equatorial current sheets model with theoretical foundations.
Figure 8. Results showing
$N_{G}$
(yellow-brown contours) and q (blue-green contours) with different combination of
$\alpha$
and
$\beta$
, under Beaming 1 (left), Beaming 2 (middle) and Beaming 3 (right) options, respectively. The catalogue values of
$N_G$
= 41 and
$q = 2.24$
are labelled as dashed lines in the contours. The fuchsia cross in each panel represents the default choices of
$\alpha = 3.0$
and
$\beta = 1.5$
used in the analysis as stated in Section 2.8.2. The red star in each panel denotes the intersection where the simulation results match both
$N_G$
and q.
Figure 9.
$\gamma$
-ray luminosity vs.
$\dot{E}$
for stacked pulsars. The green data points with error bars in both directions are from Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) 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} \gt 10^{33}$
erg s
$^{-1}$
and
$L_{\gamma} = 0.8 \dot{E}$
when
$\dot{E} \lt 10^{33}$
erg s
$^{-1}$
${\dot{\rm{T}}}$
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.
However, the most up-to-date 4FGL/3PC catalogues include a handful of pulsars that are detected with
$\dot{E}$
$\lt10^{33}$
erg/s. Given the low number of detections of
$\gamma$
-ray pulsars with
$\dot{E}$
$\lt10^{33}$
erg/s, there is no clear understanding of their emission mechanism. In Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023), a stacked signal is detected from low
$\dot{E}$
pulsars at 4.8
$\sigma$
. Applying the updated stacking methods from Henry et al. (Reference Henry2024) boosts the detection significance to 7
$\sigma$
. However, an emission mechanism was not determined due to large uncertainty in the results. We replicate the stacking analysis performed in Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) for the simulated pulsars, in an attempt to determine the possible emission mechanism for low
$\dot{E}$
pulsars. The
$\gamma$
-ray stacking analysis in Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) is performed on the undetected Galactic pulsars outside of globular clusters. The target pulsars are more than 20 away from the Galactic plane. No selections were made based on pulsar type. For this part of the discussion, we also need the population synthesis models to produce pulsars over the entire sky. For this reason, we assume that the radio pulsars are discovered with a ‘global’ Parkes survey that has the same sensitivity but covers both the northern and southern hemispheres.
$\gamma$
-ray pulsars are selected as described in Section 2.8.2.
We choose pulsars from the PN1-B1-D4 model that are radio-loud,
$\gamma$
-ray quiet, that are at least
$20^{\circ}$
away from the Galactic plane for this study, consistent with Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023). There are 46 such pulsars in the randomly chosen realisation of the simulation. We then assign each pulsar with two different
$\gamma$
-ray luminosities, one that follows the fundamental plane if
$\dot{E} \gt 10^{33}$
erg s
$^{-1}$
and
$L = 0.8 \dot{E}$
when
$\dot{E} \lt 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 Section 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 four different
$\dot{E}$
bins, with the lowest two
$\dot{E}$
bins have the same number of pulsars and highest two
$\dot{E}$
bins having the same number of pulsars. We check these results in a similar way to Figure 6 in Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023). Here in this work, we choose the canonical pulsars from Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) and plot them in Figure 9, and use this as the observational results to be compared with the model. We then plot in Figure 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. While the observational results are largely uncertain, they are consistent with all three scenarios proposed in this work.
We further 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 46 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 Figure 10. We compare this result to Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) as well as the averaged pulsar SED from McCann (Reference McCann2015). If considering the FP emission, this stacked pulsar population is relatively older compared to the young, energetic pulsars, hence the SED is much lower. Interestingly, when considering the isotropic emission, the SED is more luminous compared to the FP emission, and roughly agrees with the stacked pulsars from the data. By examining the stacked SEDs of modelled pulsars, it seems to indicate that the isotropic emission is favoured over the fundamental plane emission. However the caveat here is that the model produces a lot fewer pulsars that satisfy the stacking criteria.
Figure 10. SED of stacked pulsars from the simulation using the PN1-B1-D4 model. The orange solid curve is the averaged SED from the stacked pulsars in the simulation assuming they all follow the empirical fundamental plane emission description 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 Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023), The green dotted-dashed curve represents the averaged SED of young pulsars, and the green dotted curve SED of MSPs (McCann Reference McCann2015).
Many low
$\dot{E}$
pulsars are located outside of the Galactic plane (Smith et al. Reference Smith2023). 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 Figure 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.
Figure 12. Expected
$\gamma$
-ray flux of pulsars (the heuristic flux) as a function of the spindown luminosity
$\dot{E}$
(similar to Figure 7 in Smith et al. Reference Smith2019). 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 Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023), 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, 1, 5 kpc and 15 kpc represent those of constant distance respectively. Most of the pulsars stacked in Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023) with
$\dot{E}$
$\lt 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 Song et al. (Reference Song, Paglione, Tan, Lee-Georgescu and Herrera2023).
In Smith et al. (Reference Smith2019), a discussion of distance biases on the pulsar catalogue 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 et al. (Reference Smith2019) 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. Diffuse Galactic emission has been observed in GeV energies by Fermi-LAT (Abdo et al. Reference Abdo2009b) and TeV energies (e.g. by HESS; Abramowski et al Reference Abramowski2014), with its origin still under research (Di Mauro Reference Di Mauro2016) and unresolved pulsars being potential contributors (Smith et al. Reference Smith2019). 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
$\gamma$
-ray flux of detected and stacked pulsars in Figure 12, as for the detected pulsars, both those from the simulation and the Fermi-LAT catalogue occupy similar parts of the parameter space. The heuristic
$\gamma$
-ray 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. We exclude pulsars that are estimated to have distances equal to 25 kpc, where the distances are estimated based on the dispersion measure of radio observations, and rely on an underlying electron density model (Yao, Manchester, & Wang Reference Yao, Manchester and Wang2017). In Figure 12, the stacked pulsars from the simulation, in general, occupy similar ranges of distances compared to those from the catalogue. However the distribution of the distances from the stacked pulsars in the simulation peaks at larger distance compared to those from the catalogue.
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. 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. Some NSs are formed through AICs, and these NSs are not accounted for in this work. We exclude them as neutron stars formed through AIC may directly form MSPs (e.g. Hurley et al. Reference Hurley, Tout, Wickramasinghe, Ferrario and Kiel2010; Tauris et al. Reference Tauris, Sanyal, Yoon and Langer2013; Ruiter et al. Reference Ruiter2019), which are not the focus of this work. Our models also do not include neutron stars formed from stellar merger products. These populations, whilst not expected to produce the majority of neutron stars, may go some way to explaining the low pulsar birth rates found in this study.
-
3. Our model currently assumes a binary fraction of 100% for massive stars. Whilst this is appropriate for massive O-type stars (Sana et al. Reference Sana2012, Reference Sana2014), B-type stars likely have somewhat lower binary fractions around 70% (Moe & Di Stefano Reference Moe and Di Stefano2017; Zapartas et al. Reference Zapartas2017). Future work should incorporate a mass dependent binary fraction and the contribution of single stars. This will be particularly important hen considering millisecond pulsars and AICs, as binary interactions play a more important role in the evolution of these systems.
-
4. Our model of an exponentially decaying magnetic field (Section 2.3) is meant to provide a simple starting point for exploring the possibility of pulsar magnetic fields decaying on a broad range of timescales, as well as allowing for comparison to previous work (Kiel et al. Reference Kiel, Hurley, Bailes and Murray2008; Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Rossi and Flynn2020; Dirson et al. Reference Dirson, Pétri and Mitra2022). Recent work has introduced more sophisticated models for magnetic field decay motivated by detailed theoretical work (e.g. Viganò et al. Reference Viganò, Garcia-Garcia, Pons, Dehman and Graber2021; Sarin, Brandenburg, & Haskell Reference Sarin, Brandenburg and Haskell2023; Graber et al. Reference Graber, Ronchi, Pardo-Araujo and Rea2024). We leave the implementation of these models, and an exploration of their impact to future work.
-
5. We currently lack a physical model for the birth spin periods and magnetic fields for pulsars as a function of the progenitor properties (cf. Mandel & Müller Reference Mandel and Müller2020, for neutron star masses and kicks). This work, and any follow-up in the foreseeable future will only be able to use some empirically assumed distribution of these quantities.
-
6. We have neglected the population of magnetars with long spin periods and high magnetic fields (e.g. Ferrario & Wickramasinghe Reference Ferrario and Wickramasinghe2006, Reference Ferrario and Wickramasinghe2008; Popov et al. Reference Popov, Pons, Miralles, Boldin and Posselt2010; Gullón et al. Reference Gullón2015; Makarenko, Igoshev, & Kholtygin Reference Makarenko, Igoshev and Kholtygin2021). 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 (Ferrario et al. Reference Ferrario, Pringle, Tout and Wickramasinghe2009; Schneider et al. Reference Schneider2019; Frost et al. Reference Frost2024).
-
7. Our model assumes that the angle between the magnetic and spin axes of pulsars is constant for all pulsars and does not evolve with time. This is a simplification, and we will explore the implications of a varying angle in future work (e.g. Johnston et al. Reference Johnston, Smith, Karastergiou and Kramer2020).
-
8. Updated treatment of both radio and
$\gamma$
-ray selection effects, including the dependence of pulsar radio luminosity and beam geometry on pulsar spin properties (cf. Posselt et al. Reference Posselt2023) will improve the fidelity of our models.
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. A full scope of pulsar population synthesis can enable the study of NS contribution to the diffuse
$\gamma$
-ray background (Fermi-LAT Collaboration 2017). -
2. The inclusion of MSPs in the population synthesis models would provide further insight into the possible MSP origin of excess of GeV
$\gamma$
-ray emission from the Galactic Centre (Hooper & Mohlabeng Reference Hooper and Mohlabeng2016; Ploeg et al. Reference Ploeg, Gordon, Crocker and Macias2017; Eckner et al. Reference Eckner2018; Gautam et al. Reference Gautam2022). -
3. A comparison of pulsar populations between the Galactic field and globular clusters with results from dynamical simulations, such as those from NBODY simulations (e.g., Ye et al., Reference Ye, Kremer, Chatterjee, Rodriguez and Rasio2019, Reference Ye, Kremer, Ransom and Rasio2024).
-
4. After adding in the formation of MSPs and NSs formed through AIC, 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 (Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Bailes and Broekgaarden2021), double NSs (Chattopadhyay et al. Reference Chattopadhyay, Stevenson, Hurley, Rossi and Flynn2020) and NS-white dwarf binaries.
-
5. Including stellar mergers and the physics of magnetar formation and evolution in future versions of COMPAS will allow us to synthesise magnetar populations and compare them to observations. This will open up avenues for comparing COMPAS simulations to observations of fast radio bursts.
-
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 (e.g. Ploeg et al. Reference Ploeg, Gordon, Crocker and Macias2020) and machine learning approaches (e.g. Ronchi et al. Reference Ronchi, Graber, Garcia-Garcia, Rea and Pons2021; Graber et al. Reference Graber, Ronchi, Pardo-Araujo and Rea2024; Pardo-Araujo et al. Reference Pardo-Araujo, Ronch, Graber and Rea2025). 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 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 (Riley et al. Reference Riley2022) on both radio and
$\gamma$
-ray canonical pulsars. We generated 40 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-D4 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 = 1$
Myr best describes the observed pulsar population, based on the facts that:
-
(a) PN1-B1-D4 model produces the number of radio and
$\gamma$
-ray canonical pulsars that, given the star formation rate of the model, best fits the catalogue. -
(b) The distribution of physical properties of radio and
$\gamma$
-ray pulsars produced from this model matches the catalogue the best.
We then compare in detail the PN1-B1-D4 model results to the catalogue by examining the P-
$\dot{P}$
plots, the Galactocentric coordinates, and the heuristic flux between the simulation and the catalogue. 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-D4 model, we also attempted to explore the possible emission mechanism from low
$\dot{E}$
pulsars (
$\dot{E}$
$\lt10^{33}$
erg s
$^{-1}$
. Our results seem to favour the isotropic emission mechanism over the beamed 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.
Acknowledgement
The authors appreciate coding assistance from R. Wilcox and J. Riley, and helpful discussions from J. Hurley.
Data availability statement
COMPAS is open source and hosted on Github (https://github.com/TeamCOMPAS/COMPAS). Due to the size of COMPAS output files, the authors can share configuration files used to produce the pop-synth results upon request.
Our post-processing scripts are available on Github at https://github.com/yuzhesong/compas_pulsar_paper.git.
Funding statement
This research was conducted by the Australian Research Council (ARC) 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 and Ngarrgu Tindebeek (NT) high-performance computers at Swinburne University of Technology. OzSTAR and NT are funded by Swinburne University of Technology and the National Collaborative Research Infrastructure Strategy (NCRIS), and the Victorian Higher Education State Investment Fund (VHESIF).
Competing interests
None.
Open access
